Numerical Analysis of a Contact Problem with Wear

TL;DR

This paper develops and analyzes a numerical scheme for a quasistatic elastic contact problem with wear, providing optimal error bounds and confirming convergence through simulations.

Contribution

It introduces a more general fully discrete numerical scheme for contact problems with wear and establishes optimal error bounds with validation via simulations.

Findings

Numerical convergence orders match theoretical predictions.

Optimal error bounds are derived for the scheme.

Simulations confirm the effectiveness of the numerical method.

Abstract

This paper represents a sequel to the previous one, where numerical solution of a quasistatic contact problem is considered for an elastic body in frictional contact with a moving foundation. The model takes into account wear of the contact surface of the body caused by the friction. Some preliminary error analysis for a fully discrete approximation of the contact problem was provided in the previous paper. In this paper, we consider a more general fully discrete numerical scheme for the contact problem, derive optimal order error bounds and present computer simulation results showing that the numerical convergence orders match the theoretical predictions.

Tables

Table 1. Table 1: Numerical errors

$h+k$	$1$	$1/2$	$1/4$	$1/8$	$1/16$
${\Vert 𝒖-{𝒖}^{hk}\Vert }_V/{\Vert 𝒖\Vert }_V$	$4.1698e^{-1}$	$2.6840e^{-1}$	$1.4360e^{-1}$	$7.3979e^{-2}$	$3.4882e^{-2}$
Convergence order		0.6355	0.9022	0.9569	1.0846
${\Vert w-w^{hk}\Vert }_W/{\Vert w\Vert }_W$	$2.9009e^{-1}$	$1.0328e^{-1}$	$3.8385e^{-2}$	$1.4694e^{-2}$	$5.0891e^{-3}$
Convergence order		1.4898	1.4280	1.3853	1.5297

Equations

$$\mbox{\boldmath{$\varepsilon$}}(\mbox{\boldmath{$u$}})=(\varepsilon_{ij}(\mbox{\boldmath{$u$}})),\quad\varepsilon_{ij}(\mbox{\boldmath{$u$}})=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right).$$

$$\mbox{\boldmath{$n$}}^{*}(t)=-\mbox{\boldmath{$v$}}^{*}(t)/\|\mbox{\boldmath{$v$}}^{*}(t)\|,\quad\alpha(t)=\kappa\,\|\mbox{\boldmath{$v$}}^{*}(t)\|,$$

$$\displaystyle\mbox{\boldmath{$\sigma$}}(t)={\mathcal{F}}\mbox{\boldmath{$\varepsilon$}}(\mbox{\boldmath{$u$}}(t))\quad$$

$$\displaystyle{\rm Div}\,\mbox{\boldmath{$\sigma$}}(t)+\mbox{\boldmath{$f$}}_{0}(t)=\mbox{\boldmath{$0$}}\quad$$

$$\displaystyle\mbox{\boldmath{$u$}}(t)=\mbox{\boldmath{$0$}}\quad$$

$$\displaystyle\mbox{\boldmath{$\sigma$}}(t)\mbox{\boldmath{$\nu$}}=\mbox{\boldmath{$f$}}_{N}(t)\quad$$

$$\displaystyle\left.\begin{array}[]{ll}u_{\nu}(t)\leq g,\ \sigma_{\nu}(t)+p(u_{\nu}(t)-w(t))\leq 0,\\ (u_{\nu}(t)-g)(\sigma_{\nu}(t)+p(u_{\nu}(t)-w(t)))=0\end{array}\right\}\quad$$

$$\displaystyle-\mbox{\boldmath{$\sigma$}}_{\tau}(t)=\mu\,p(u_{\nu}(t)-w(t))\,\mbox{\boldmath{$n$}}^{*}(t)\quad$$

$$w^{\prime }(t)=\alpha (t)p(u_{\nu }(t)-w(t))$$

$$w(0)=0$$

$${\cal H}=L^{2}(\Omega;\mathbb{S}^{d}),\quad V=\{\mbox{\boldmath{$v$}}\in H^{1}(\Omega)^{d}\mid\mbox{\boldmath{$v$}}=\mbox{\boldmath{$0$}}\ {\rm on}\ \Gamma_{D}\}$$

$$(\mbox{\boldmath{$\sigma$}},\mbox{\boldmath{$\tau$}})_{\cal H}=\int_{\Omega}\sigma_{ij}\tau_{ij}dx,\quad(\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}})_{V}=(\mbox{\boldmath{$\varepsilon$}}(\mbox{\boldmath{$u$}}),\mbox{\boldmath{$\varepsilon$}}(\mbox{\boldmath{$v$}}))_{\cal H}$$

$$U=\{\mbox{\boldmath{$v$}}\in V\mid v_{\nu}\leq g\ {\rm on}\ \Gamma_{C}\}.$$

$$\|\mbox{\boldmath{$v$}}\|_{L^{2}(\Gamma_{C})^{d}}\leq c_{0}\|\mbox{\boldmath{$v$}}\|_{V}\quad\forall\,\mbox{\boldmath{$v$}}\in V.$$

$$|p(\mbox{\boldmath{$x$}},r)|\leq L_{p}|r|\quad\forall\,r\in\mathbb{R},\ {\rm a.e.}\ \mbox{\boldmath{$x$}}\in\Gamma_{C}.$$

$$\mbox{\boldmath{$f$}}_{0}\in C([0,T];L^{2}(\Omega)^{d}),\quad\mbox{\boldmath{$f$}}_{N}\in C([0,T];L^{2}(\Gamma_{N})^{d}).$$

$$\mbox{\boldmath{$n$}}^{*}\in C([0,T];\mathbb{R}^{d}),\quad\alpha\in C([0,T];L^{\infty}(\Gamma_{C})),$$

$$\displaystyle\langle F\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}}\rangle_{V^{*}\times V}=({\cal F}(\mbox{\boldmath{$\varepsilon$}}(\mbox{\boldmath{$u$}})),\mbox{\boldmath{$\varepsilon$}}(\mbox{\boldmath{$v$}}))_{\cal H},$$

$$\displaystyle\langle\mbox{\boldmath{$f$}}(t),\mbox{\boldmath{$v$}}\rangle_{V^{*}\times V}=\int_{\Omega}\mbox{\boldmath{$f$}}_{0}(t)\cdot\mbox{\boldmath{$v$}}\,dx+\int_{\Gamma_{N}}\mbox{\boldmath{$f$}}_{N}(t)\cdot\mbox{\boldmath{$v$}}\,da,$$

$$\displaystyle\varphi(t,w,\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}})=\int_{\Gamma_{C}}p(u_{\nu}-w)\left[v_{\nu}+\mu\,\mbox{\boldmath{$n$}}^{*}(t)\cdot\mbox{\boldmath{$v$}}_{\tau}\right]da.$$

$$\displaystyle\langle F\mbox{\boldmath{$u$}}(t),\mbox{\boldmath{$v$}}-\mbox{\boldmath{$u$}}(t)\rangle_{V^{*}\times V}+\varphi(t,w(t),\mbox{\boldmath{$u$}}(t),\mbox{\boldmath{$v$}})-\varphi(t,w(t),\mbox{\boldmath{$u$}}(t),\mbox{\boldmath{$u$}}(t))$$

$$\displaystyle\qquad\qquad\qquad{}\geq\langle\mbox{\boldmath{$f$}}(t),\mbox{\boldmath{$v$}}-\mbox{\boldmath{$u$}}(t)\rangle_{V^{*}\times V}\quad\forall\,\mbox{\boldmath{$v$}}\in U,$$

$$w(t)={\int }_0^t\alpha (s)p(u_{\nu }(s)-w(s))ds.$$

$$\mbox{\boldmath{$u$}}\in C([0,T];V),\quad\mbox{\boldmath{$\sigma$}}\in C([0,T];{\cal H}),\quad w\in C^{1}([0,T];W).$$

$$\alpha \in W^{1,\infty }(0,T;L^{\infty }({\Gamma }_C)).$$

$$\displaystyle\langle F\mbox{\boldmath{$u$}}^{hk}_{n},\mbox{\boldmath{$v$}}^{h}-\mbox{\boldmath{$u$}}^{hk}_{n}\rangle_{V^{*}\times V}+\varphi(t_{n},w^{hk}_{n},\mbox{\boldmath{$u$}}^{hk}_{n},\mbox{\boldmath{$v$}}^{h})-\varphi(t_{n},w^{hk}_{n},\mbox{\boldmath{$u$}}^{hk}_{n},\mbox{\boldmath{$u$}}^{hk}_{n})$$

$$\displaystyle\qquad\qquad\qquad\geq\langle\mbox{\boldmath{$f$}}_{n},\mbox{\boldmath{$v$}}^{h}-\mbox{\boldmath{$u$}}^{hk}_{n}\rangle_{V^{*}\times V}\quad\forall\,\mbox{\boldmath{$v$}}^{h}\in U^{h},$$

$$w_n^{hk}=\sum _{j=0}^{n-1}{k}_j{\alpha }_{j}p(u_{j,\nu }^{hk}-w_j^{hk}).$$

$$e_n\le c{g}_n+c\sum _{j=1}^{n-1}{k}_j{e}_j,n=1,\dots ,N.$$

$$e_{n}\leq c\,\Big{(}g_{n}+\sum_{j=1}^{n-1}k_{j}g_{j}\Big{)},\quad n=1,\dots,N.$$

$$\underset{1\le n\le N}{\max }{e}_n\le c\underset{1\le n\le N}{\max }{g}_n.$$

$$\displaystyle\max_{1\leq n\leq N}\left(\|\mbox{\boldmath{$u$}}_{n}-\mbox{\boldmath{$u$}}^{hk}_{n}\|_{V}^{2}+\|w_{n}-w^{hk}_{n}\|_{W}^{2}\right)$$

$$\displaystyle\quad{}+c\,\max_{1\leq n\leq N}\left(\|\mbox{\boldmath{$u$}}_{n}-\mbox{\boldmath{$v$}}^{h}_{n}\|_{V}^{2}+|R_{n}(w_{n},\mbox{\boldmath{$u$}}_{n},\mbox{\boldmath{$v$}}^{h}_{n})|\right)$$

$$\displaystyle R_{n}(w_{n},\mbox{\boldmath{$u$}}_{n},\mbox{\boldmath{$v$}}^{h}_{n})$$

$$\displaystyle\quad{}-\langle\mbox{\boldmath{$f$}}_{n},\mbox{\boldmath{$v$}}^{h}_{n}-\mbox{\boldmath{$u$}}_{n}\rangle_{V^{*}\times V}.$$

$$\displaystyle\|\mbox{\boldmath{$u$}}_{n}-\mbox{\boldmath{$u$}}^{hk}_{n}\|_{V}^{2}+\|w_{n}-w^{hk}_{n}\|_{W}^{2}$$

$$\displaystyle\quad{}+c\sum_{j=1}^{n-1}k_{j}\left(\|\mbox{\boldmath{$u$}}_{j}-\mbox{\boldmath{$u$}}^{hk}_{j}\|_{V}^{2}+\|w_{j}-w^{hk}_{j}\|_{W}^{2}\right).$$

$$\displaystyle\varphi(t,w_{1},\mbox{\boldmath{$u$}}_{1},\mbox{\boldmath{$v$}}_{2})+\varphi(t,w_{2},\mbox{\boldmath{$u$}}_{2},\mbox{\boldmath{$v$}}_{1})-\varphi(t,w_{1},\mbox{\boldmath{$u$}}_{1},\mbox{\boldmath{$v$}}_{1})-\varphi(t,w_{2},\mbox{\boldmath{$u$}}_{2},\mbox{\boldmath{$v$}}_{2})$$

Full text

Numerical Analysis of a Contact Problem with Wear

Danfu Han111Department of Mathematics, Hangzhou Normal University, Hangzhou, China. Email: mhdf@hznu.edu.cn, Weimin Han222Program in Applied Mathematical and Computational Sciences (AMCS) & Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA. Email: weimin-han@uiowa.edu, Michal Jureczka333Jagiellonian University in Krakow, Faculty of Mathematics and Computer Science, Lojasiewicza 6, 30-348 Krakow, Poland. Email: michal.jureczka@uj.edu.pl and  Anna Ochal444Jagiellonian University in Krakow, Faculty of Mathematics and Computer Science, Lojasiewicza 6, 30-348 Krakow, Poland. Email: ochal@ii.uj.edu.pl

Abstract. This paper represents a sequel to [13] where numerical solution of a quasistatic contact problem is considered for an elastic body in frictional contact with a moving foundation. The model takes into account wear of the contact surface of the body caused by the friction. Some preliminary error analysis for a fully discrete approximation of the contact problem was provided in [13]. In this paper, we consider a more general fully discrete numerical scheme for the contact problem, derive optimal order error bounds and present computer simulation results showing that the numerical convergence orders match the theoretical predictions.

Keywords. Quasistatic contact problem, variational inequality, numerical methods, optimal order error estimate.

AMS Classification. 65N30

1 Introduction

Contact phenomenon is common in engineering applications. Mathematical studies and numerical analysis of contact problems are most suitably carried out within the frameworks of variational inequalities or hemivariational inequalities, which have attracted the attention of many researchers. The related mathematical literature grows rapidly. Some representative comprehensive references in this area are [5, 14, 12, 10, 21, 6, 22, 8] in the context of variational inequalities and [19, 17, 11, 16, 23] in the context of hemivariational inequalities.

For a contact problem, the mathematical model is constructed based on considerations of various aspects of the contact process. Factors to be taken into account include the type of the contact process (static, quasi-static, or dynamic), constitutive relations of the deformable bodies, contact conditions of various application-specific forms. In certain applications, it is important to consider heating or thermo effects ([18]), or piezoelectricity effects ([25]). Since the contact process inevitably causes material wear or even damage, it is not surprising that the wear effect has been built into mathematical models for a variety of contact processes, cf. [3, 9, 20, 15, 2, 26, 7]. In a recent paper [24], a mathematical model is proposed and studied for contact with wear described by Archard’s law of surface wear. In this model, the friction between a deformable body and the foundation leads to wear of the contact surface of the body over time. Solution existence and uniqueness for the model are provided in [24]. Numerical approximation of the contact problem is the subject of [13] where some error bounds are derived for a fully discrete scheme. In this paper, we take a further step by considering a more general fully discrete numerical scheme for the contact problem that allows an arbitrary partition of the time interval, providing optimal order error estimates of the fully discrete scheme to solve the contact problem. Moreover, we present numerical results showing deformation of the contact body and numerical convergence orders of the fully discrete solutions that confirm the theoretical error bounds.

The remainder of this paper is organized as follows. In Section 2 we introduce the contact problem and its variational formulation. In Section 3, we study a fully discrete numerical scheme and derive optimal order error bounds. In Section 4, we present computational simulation results for numerical convergence orders that match the theoretical predictions.

2 The contact problem and its variational formulation

First, we describe the physical setting of the contact problem. Consider a deformable body that occupies a domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ in application. The body is subject to the action of volume forces with a total density \mbox{\boldmath{f}}_{0}. The boundary $\Gamma$ of the domain $\Omega$ is assumed to be Lipschitz continuous and is divided into three disjoint measurable parts $\Gamma_{D}$, $\Gamma_{N}$ and $\Gamma_{C}$, with ${\rm meas}\,(\Gamma_{D})>0$. Denote by $\nu$ the unit outward normal vector on $\Gamma$ that is defined a.e. on $\Gamma$. The body is clamped on $\Gamma_{D}$, i.e., the displacement is equal to [math] on $\Gamma_{D}$. Surface transactions of a total density \mbox{\boldmath{f}}_{N} act on the boundary $\Gamma_{N}$. The contact boundary is $\Gamma_{C}$ where the contact is modeled by a normal compliance condition with a unilateral constraint and Coulomb’s law of dry friction. Following [24], we assume that the body is elastic, in contact with a moving obstacle (foundation) made of a hard perfectly rigid material, and assume that the contact surface of the body $\Gamma_{C}$ is covered by a layer of soft material. This layer is deformable and the foundation may penetrate it. Frictional contact with the foundation may cause this layer to wear over time.

We assume that the acceleration of the body is negligible and so the problem is quasistatic. In our model, the framework of the small strain theory is employed. We are interested in the body displacement and foundation wear in a time interval $[0,T]$, with $T>0$. We denote by “$\cdot$” and $\|\cdot\|$ the scalar product and the Euclidean norm in $\mathbb{R}^{d}$ or $\mathbb{S}^{d}$, respectively, where $\mathbb{S}^{d}$ is the space of symmetric matrices of order $d$. The indices $i$ and $j$ run from $1$ to $d$ and the index after a comma represents the partial derivative with respect to the corresponding component of the independent variable. Summation convention over repeated indices is adopted. We denote the divergence operator by {\rm Div}\mbox{\boldmath{\sigma}}=(\sigma_{ij,j}) for an $\mathbb{S}^{d}$-valued field $\sigma$. Standard Lebesgue and Sobolev spaces will be used, such as $L^{2}(\Omega)^{d}=L^{2}(\Omega;\mathbb{R}^{d})$ and $H^{1}(\Omega)^{d}=H^{1}(\Omega;\mathbb{R}^{d})$. Recall that the linearized strain tensor of a displacement field \mbox{\boldmath{u}}\in H^{1}(\Omega)^{d} is

[TABLE]

Let u_{\nu}=\mbox{\boldmath{u}}\cdot\mbox{\boldmath{\nu}} and \sigma_{\nu}=\mbox{\boldmath{\sigma}}\mbox{\boldmath{\nu}}\cdot\mbox{\boldmath{\nu}} be the normal components of $u$ and $\sigma$, respectively, and let \mbox{\boldmath{u}}_{\tau}=\mbox{\boldmath{u}}-u_{\nu}\mbox{\boldmath{\nu}} and \mbox{\boldmath{\sigma}}_{\tau}=\mbox{\boldmath{\sigma}}\mbox{\boldmath{\nu}}-\sigma_{\nu}\mbox{\boldmath{\nu}} be their tangential components, respectively. To simplify the notation, we will usually not indicate explicitly the dependence of various functions on the spatial variable $x$.

Denote by \mbox{\boldmath{v}}^{*}(t)\not=\mbox{\boldmath{0}} the velocity of the foundation. Let

[TABLE]

where $\kappa$ represents the wear coefficient, and let $\mu$ be the friction coefficient. The classical formulation of the contact problem with wear is as follows.

Problem 1

Find a displacement field \mbox{\boldmath{u}}\colon\Omega\times[0,T]\to\mathbb{R}^{d}, a stress field \mbox{\boldmath{\sigma}}\colon\Omega\times[0,T]\to\mathbb{S}^{d}, and a wear function $w\colon\Gamma_{C}\times[0,T]\to\mathbb{R}_{+}=[0,\infty)$ such that for all $t\in[0,T]$,

[TABLE]

In Problem 1, equation (2.2) represents an elastic constitutive law with an elasticity operator ${\cal F}$. Equation (2.3) is the equilibrium equation. The equality (2.4) describes the fact that body is clamped on $\Gamma_{D}$ and (2.5) represents external forces acting on $\Gamma_{N}$. The relations in (2.8) describe the damping response of the foundation, $g>0$ being the thickness of a soft layer covering $\Gamma_{C}$. The friction is modeled by equation (2.9). Here, the size of \mbox{\boldmath{v}}^{*} is assumed to be significantly larger than that of the tangential body velocity \mbox{\boldmath{u}}^{\prime}_{\tau}. Equations (2.10) and (2.11) govern the evolution of the wear function. Detailed derivation of this model is presented in [24].

The contact problem will be studied in its variational formulation. For this purpose, we introduce function spaces and hypotheses on the problem data. We recall that for a normed space $X$, $C([0,T];X)$ is the space of continuous functions from $[0,T]$ to $X$. We will use the following Hilbert spaces:

[TABLE]

endowed with the inner scalar products

[TABLE]

with the corresponding norms. Denote by $\langle\cdot,\cdot\rangle_{V^{*}\times V}$ the duality pairing between a dual space $V^{*}$ and $V$. The set of admissible displacements is

[TABLE]

For a function \mbox{\boldmath{v}}\in V, we use the same symbol $v$ for its trace on the boundary $\Gamma$. By the Sobolev trace theorem, there exists a constant $c_{0}>0$ depending only on $\Omega$, $\Gamma_{D}$ and $\Gamma_{C}$ such that

[TABLE]

Now we introduce the hypotheses on the data needed in the study of Problem 1.

$H({\cal F})$: For the elasticity operator ${\cal F}\colon\Omega\times\mathbb{S}^{d}\to\mathbb{S}^{d}$,

a (a) {\cal F}(\cdot,\mbox{\boldmath{\varepsilon}}) is measurable on $\Omega$ for all \mbox{\boldmath{\varepsilon}}\in\mathbb{S}^{d}, {\cal F}(\cdot,\mbox{\boldmath{0}})\in{\cal H};

a (b) $\exists$ $L_{\cal F}>0$ s.t. \|{\cal F}(\mbox{\boldmath{x}},\mbox{\boldmath{\varepsilon}}_{1})-{\cal F}(\mbox{\boldmath{x}},\mbox{\boldmath{\varepsilon}}_{2})\|\leq L_{\cal F}\|\mbox{\boldmath{\varepsilon}}_{1}-\mbox{\boldmath{\varepsilon}}_{2}\| \forall\,\mbox{\boldmath{\varepsilon}}_{1},\mbox{\boldmath{\varepsilon}}_{2}\in\mathbb{S}^{d}, a.e. \mbox{\boldmath{x}}\in\Omega;

a (c) $\exists$ $m_{\cal F}>0$ s.t. ({\cal F}(\mbox{\boldmath{x}},\mbox{\boldmath{\varepsilon}}_{1})-{\cal F}(\mbox{\boldmath{x}},\mbox{\boldmath{\varepsilon}}_{2}))\cdot(\mbox{\boldmath{\varepsilon}}_{1}-\mbox{\boldmath{\varepsilon}}_{2})\geq m_{\cal F}\|\mbox{\boldmath{\varepsilon}}_{1}-\mbox{\boldmath{\varepsilon}}_{2}\|^{2} \forall\,\mbox{\boldmath{\varepsilon}}_{1},\mbox{\boldmath{\varepsilon}}_{2}\in\mathbb{S}^{d}, a.e. \mbox{\boldmath{x}}\in\Omega.

$H(p)$: For the normal compliance function $p\colon\Gamma_{C}\times\mathbb{R}\to\mathbb{R}_{+}$,

a (a) $p(\cdot,r)$ is measurable on $\Gamma_{C}$ $\forall\,r\in\mathbb{R}$;

a (b) $\exists$ $L_{p}>0$ s.t. |p(\mbox{\boldmath{x}},r_{1})-p(\mbox{\boldmath{x}},r_{2})|\leq L_{p}|r_{1}-r_{2}| $\forall\,r_{1},r_{2}\in\mathbb{R}$, a.e. \mbox{\boldmath{x}}\in\Gamma_{C};

a (c) (p(\mbox{\boldmath{x}},r_{1})-p(\mbox{\boldmath{x}},r_{2}))(r_{1}-r_{2})\geq 0 $\forall\,r_{1},r_{2}\in\mathbb{R}$, a.e. \mbox{\boldmath{x}}\in\Gamma_{C};

a (d) p(\mbox{\boldmath{x}},r)=0 $\forall\,r\leq 0$, a.e. \mbox{\boldmath{x}}\in\Gamma_{C}.

Note that $H(p)$ (b) and (d) imply

[TABLE]

H(\mbox{\boldmath{f}}): For the densities of body and traction forces,

[TABLE]

$H_{0}$: For the friction and wear coefficients, and the foundation velocity,

a (a) $\mu\in L^{\infty}(\Gamma_{C})$, \mu(\mbox{\boldmath{x}})\geq 0 a.e. \mbox{\boldmath{x}}\in\Gamma_{C};

a (b) $\kappa\in L^{\infty}(\Gamma_{C})$, \kappa(\mbox{\boldmath{x}})\geq 0 a.e. \mbox{\boldmath{x}}\in\Gamma_{C};

a (c) \mbox{\boldmath{v}}^{*}\in C([0,T];\mathbb{R}^{d}), \|\mbox{\boldmath{v}}^{*}(t)\|\geq v_{0}>0 $\forall\,t\in[0,T]$.

We notice that hypotheses $H_{0}$ implies the following regularities:

[TABLE]

where \mbox{\boldmath{n}}^{*} and $\alpha$ are defined in (2.1).

Finally, we will need a smallness assumption on the combined effect of the Lipschitz constant of the normal compliance function $p$ and the friction coefficient $\mu$. Recall that $c_{0}$ is the constant in the inequality (2.12).

$H_{s}$: $c^{2}_{0}L_{p}\|\mu\|_{L^{\infty}(\Gamma_{C})}<m_{\cal F}$.

Now we define some operators and functions needed in the variational formulation of Problem 1. Let $F\colon V\to V^{*}$, \mbox{\boldmath{f}}\colon[0,T]\to V^{*} and $\varphi\colon[0,T]\times L^{2}(\Gamma_{C})\times V\times V\to\mathbb{R}$ be defined for all \mbox{\boldmath{u}},\mbox{\boldmath{v}}\in V, $w\in L^{2}(\Gamma_{C})$, $t\in[0,T]$ as follows:

[TABLE]

Let $W=L^{2}(\Gamma_{C})$ be the space for the wear variable $w$. Using the standard procedures in the mathematical theory of contact mechanics, we obtain the week formulation of Problem 1.

Problem 2

Find \mbox{\boldmath{u}}\colon[0,T]\to U and $w\colon[0,T]\to W$ such that for all $t\in[0,T]$,

[TABLE]

We recall the following existence and uniqueness result for Problem 2 from [24].

Theorem 3

Assume $H({\cal F})$, $H(p)$, H(\mbox{\boldmath{f}}), $H_{0}$ and $H_{s}$. Then Problem 2 has a unique solution with the regularity

[TABLE]

In addition, $w(t)\geq 0$ for all $t\in[0,T]$, a.e. on $\Gamma_{C}$.

3 Numerical analysis

We turn to the numerical solution of Problem 2. Let $V^{h}\subset V$ and $W^{h}\subset W$ be two families of finite dimensional subspaces with a discretization parameter $h>0$. Then define $U^{h}=U\cap V^{h}$. Let $0=t_{0}<t_{1}<\cdots<t_{N}=T$ be a partition of the time interval $[0,T]$. Denote $k_{n}=t_{n+1}-t_{n}$, $0\leq n\leq N-1$, and $k=\max_{0\leq n\leq N-1}k_{n}$ for the time step size. For a function $z$ continuous in $t$, we write $z_{n}=z(t_{n})$.

We make the following additional assumptions on the solution $u$ to Problem 2 and the velocity of the foundation \mbox{\boldmath{v}}^{*}.

$H_{1}$: \mbox{\boldmath{u}}\in H^{1}(0,T;V), \mbox{\boldmath{v}}^{*}\in W^{1,\infty}(0,T;\mathbb{R}^{d}).

Note that assumptions $H_{1}$ and $H_{0}$ (b) imply that

[TABLE]

Consider the following fully discrete scheme for solving Problem 2.

Problem 4

Find \mbox{\boldmath{u}}^{hk}=\{\mbox{\boldmath{u}}^{hk}_{n}\}_{n=0}^{N}\subset U^{h} and $w^{hk}=\{w^{hk}_{n}\}_{n=0}^{N}\subset W^{h}$, $w^{hk}_{0}=0$, such that for $0\leq n\leq N$,

[TABLE]

and for $1\leq n\leq N$,

[TABLE]

We remark that existence of a unique solution to Problem 4 follows from an application of discrete version of Theorem 3. We also remark that the numerical scheme considered in [13] is a special case of Problem 4 where a uniform partition of the time interval $[0,T]$ is used. For a uniform partition of $[0,T]$ into $N$ equal size sub-intervals, we let $k=T/N$ be the time step and $t_{n}=n\,k$, $0\leq n\leq N$, the node points.

We will make use of the following discrete Gronwall inequality ([10, Lemma 7.25]).

Lemma 5

Assume $\{g_{n}\}_{n=1}^{N}$ and $\{e_{n}\}_{n=1}^{N}$ are two sequences of non-negative numbers satisfying

[TABLE]

Then

[TABLE]

Therefore,

[TABLE]

We have Ceá’s inequality useful for error estimation.

Theorem 6

Under the assumptions stated in Theorem 3 and the additional hypothesis $H_{1}$, there exists a constant $c>0$ such that for any \mbox{\boldmath{v}}^{h}_{n}\in U^{h}, $1\leq n\leq N$,

[TABLE]

where

[TABLE]

Proof. By modifying the proof of Theorem 4 in [13], we can establish the inequality

[TABLE]

Applying Lemma 5 on (3.6), we get the inequality (3.4).

Note that from $H(p)$ and $H_{0}$, we have (cf. [13, (27)]), for $t\in[0,T]$,

[TABLE]

The inequality (3.4) is the starting point for further error estimation. For simplicity, we assume $\Omega$ is a polygonal/polyhedral domain. Then $\Gamma_{D}$, $\Gamma_{N}$ and $\Gamma_{C}$ can be expressed as unions of flat components (line segments for $d=2$ and polygons for $d=3$) that have pairwise disjoint interiors. In particular, we write $\overline{\Gamma_{C}}=\cup_{i=1}^{i_{0}}\Gamma_{C,i}$, where each component $\Gamma_{C,i}$ is a line segment if $d=2$ or a polygon if $d=3$. Consider a regular family of finite element partitions $\{{\cal T}^{h}\}$ of the domain $\overline{\Omega}$ into triangular or tetrahedral elements such that if the intersection of one side/face of an element with one flat component of the boundary has a positive relative measure, then the side/face lies entirely in that flat component. Corresponding to ${\cal T}^{h}$, we define the linear element space

[TABLE]

Then we define the discrete admissible finite element set

[TABLE]

We assume $g$ is a concave function. Then, $U^{h}=V^{h}\cap U\subset U$. We proceed to derive an optimal order error estimate for the finite element solution defined by Problem 4.

Theorem 7

Keep the assumptions stated in Theorem 6. Assume further the solution regularities

[TABLE]

Then we have the optimal order error estimate

[TABLE]

Proof. By following the arguments presented in [10, Section 8.1], it can be shown that under the stated regularity assumptions, the solution of Problem 2 satisfies, for $t\in[0,T]$,

[TABLE]

where

[TABLE]

Using these relations we find that

[TABLE]

Thus,

[TABLE]

This provides an upper bound for the term \left|R_{n}(w_{n},\mbox{\boldmath{u}}_{n},\mbox{\boldmath{v}}^{h}_{n})\right| on the right hand side of (3.4).

Now we bound the error \|\mbox{\boldmath{u}}_{0}-\mbox{\boldmath{u}}^{hk}_{0}\|_{V}. For simplicity, we denote

[TABLE]

Write

[TABLE]

From (2.15) with $t=0$,

[TABLE]

From (3.2) with $n=0$,

[TABLE]

Take \mbox{\boldmath{v}}=\mbox{\boldmath{u}}^{hk}_{0} in (3.15), and use the resulting inequality and the inequality (3.16) in (3.14) to obtain

[TABLE]

By $H({\cal F})$ (c),

[TABLE]

By $H({\cal F})$ (b),

[TABLE]

Then, for an arbitrarily small $\epsilon>0$, there is a constant $c$ depending on $\epsilon$ such that

[TABLE]

By (3.13),

[TABLE]

By (3.7),

[TABLE]

Since

[TABLE]

for the arbitrarily small $\epsilon>0$, there is a constant $c$ depending on $\epsilon$ such that

[TABLE]

Using these relations in (3.17), we obtain

[TABLE]

Recall the condition $H_{s}$; choosing $\epsilon=\left(m_{\cal F}-c^{2}_{0}L_{p}\|\mu\|_{L^{\infty}(\Gamma_{C})}\right)/4$ we obtain from the above inequality that

[TABLE]

Using (3.18) and (3.13) in (3.4), we have

[TABLE]

for any \mbox{\boldmath{v}}^{h}_{n}\in U^{h}.

Thus, by applying the finite element interpolation theory (e.g., [1, 4]), we have the optimal order error bound (3.12) from (3.19), under the solution regularities (3.10) and (3.11).

We comment that if {\cal F}(\mbox{\boldmath{x}},\mbox{\boldmath{\varepsilon}}) is a smooth function of $x$, in particular if {\cal F}(\mbox{\boldmath{x}},\mbox{\boldmath{\varepsilon}}) does not depend on $x$, then (3.11) follows from (3.10) and thus there is no need to assume (3.11).

4 Numerical results

In this section, we report computer simulation results on a numerical example. Let $d=2$ and consider a square-shaped set $\Omega=(0,1)\times(0,1)$ with the following partition of the boundary

[TABLE]

The linear elasticity operator $\mathcal{F}$ is defined by

[TABLE]

Here $I$ denotes the identity matrix, tr denotes the trace of the matrix, $\lambda>0$ and $\eta>0$ are the Lame coefficients. In our simulations, we choose $\lambda=\eta=4$, $T=1$ and take the following data

[TABLE]

We use the linear finite element space $V^{h}$ defined in (3.8) and its subset $U^{h}$ defined in (3.9), based on uniform triangular partitions of $\overline{\Omega}$. We use the uniform partition of the time interval $[0,1]$ with the time step size $k=1/N$ for a positive integer $N$.

We first demonstrate the effect of some input data on the deformation of the body. In all cases, we show the shape of the body at final time $t=1$, and the numerical solutions correspond to the time step size $1/16$ and where the boundary $\Gamma_{C}$ of the body is divided into $16$ equal parts.

In Figure 2 we show the deformed configuration for $\mu(\bm{x})=0.3$, $\kappa(\bm{x})=0.04$ and $\bm{v^{*}}(\bm{x},t)=(1,0)$. We push the body down towards the moving foundation with a force $\bm{f}_{0}$, and as a result of friction, the soft layer of material covering $\Gamma_{C}$ wears out allowing the body to move downward. We observe that in this case coefficient $\kappa$, governing the rate of wear, is not big enough to cause the body to touch the foundation. Because of the friction, the body moves in the same direction as the foundation, i.e. to the right.

We then increase the wear coefficient $\kappa$ to $\kappa(\bm{x})=0.08$. The deformed configuration is shown in Figure 2. We observe that the layer of soft material on part of the boundary $\Gamma_{C}$ completely wears out, allowing the body to rest on the rigid foundation as it cannot penetrate it further.

In Figure 4, we show the deformed configuration for $\mu(\bm{x})=1$, $\kappa(\bm{x})=0.04$ and $\bm{v^{*}}(\bm{x},t)=(1,0)$. We observe that the body moves further to the right, which is a result of increased friction between soft layer of material covering $\Gamma_{C}$ and the rigid foundation.

The result in Figure 4 corresponds to $\mu(\bm{x})=0.3$, $\kappa(\bm{x})=0.02$ and $\bm{v^{*}}(\bm{x},t)=(-1,0)$. Note that the direction of the motion of the foundation is reversed. As a result, the lower part of the body squeezes to the left and we observe that the boundary $\Gamma_{C}$ is slightly curled. We conclude that all those modifications lead to results that can be expected.

Finally, we explore the numerical convergence orders of the numerical method on the model problem with $\bm{f}_{N}(\bm{x},t)=(-0.5,-0.5)$, $\bm{f}_{0}(\bm{x},t)=(-0.5,-2)$, $\mu(\bm{x})=1.0$, $\kappa(\bm{x})=0.05$, and $\bm{v^{*}}(\bm{x},t)=(1,0)$. We present a comparison of numerical errors $\|\bm{u}-\bm{u}^{hk}\|_{V}$ and $\|w-w^{hk}\|_{W}$ computed for a sequence of solutions to discretized problems. We use a uniform discretization of the problem domain and time interval according to the spatial discretization parameter $h$ and time step size $k$, respectively. The boundary $\Gamma_{C}$ of $\Omega$ is divided into $1/h$ equal parts. We start with $h=1/2$ and $k=1/2$, which are successively halved. The numerical solution corresponding to $h=1/64$ and $k=1/64$ is taken as the “exact” solution $\bm{u}$ and $w$ with $\|\bm{u}\|_{V}\doteq 0.19131$ and $\|w\|_{W}\doteq 0.08192$. The results are presented in Table 1 and Figures 6 and 6, where the dependence of the relative error estimates $\|\bm{u}-\bm{u}^{hk}\|_{V}/\|\bm{u}\|_{V}$ and $\|w-w^{hk}\|_{W}/\|w\|_{W}$ with respect to $h+k$ are plotted on a log-log scale. A first order convergence is clearly observed for the numerical solutions of the displacement. The numerical convergence orders for the numerical solutions of the wear function are somewhat higher than 1.

**Acknowledgments

**The project leading to this application has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 823731 CONMECH.