Numerical analysis of a non-clamped dynamic thermoviscoelastic contact problem

TL;DR

This paper presents a comprehensive numerical analysis of a complex thermoviscoelastic contact problem involving thermal effects and non-monotone friction laws, providing optimal error estimates and numerical validation.

Contribution

It introduces a fully discrete approximation for the problem and derives optimal error estimates without smallness restrictions on the data.

Findings

Optimal error estimates achieved for the numerical scheme

Numerical simulations confirm theoretical results

Handles non-monotone friction relations effectively

Abstract

In this work, we analyze a non-clamped dynamic viscoelastic contact problem involving thermal effect. The friction law is described by a non-monotone relation between the tangential stress and the tangential velocity. This leads to a system of second-order inclusion for displacement and a parabolic equation for temperature. We provide a fully discrete approximation of the problem and find optimal error estimates without any smallness assumption on the data. The theoretical result is illustrated by numerical simulations.

Equations

$${\phi }^0(x;v)=\underset{y\to x,t\downarrow 0}{\mathrm{limsup}}\frac{\phi (y+tv)-\phi (y)}{t}.$$

$$\partial \phi (x)=\{\zeta \in X^{\ast }\mid {\phi }^0(x;v)\ge {⟨\zeta ,v⟩}_{X^{\ast }\times X}\text{mbox}\mathrm{forall}v\in X\}.$$

$$u\cdot v$$

$$\sigma :\tau$$

$$\mathrm{div}u=u_{i,i},\epsilon (u)=({\epsilon }_{ij}(u)),{\epsilon }_{ij}(u)=\frac{1}{2}(u_{i,j}+u_{j,i}).$$

$$\mathrm{Div}\sigma =({\sigma }_{ij,j}).$$

$$H=L^2(\Omega ;{\mathbb{R}}^d),Q=\{\sigma =({\sigma }_{ij})|{\sigma }_{ij}={\sigma }_{ji}\in L^2(\Omega )\}.$$

$$(u,v{)}_H={\int }_{\Omega }{u}_i{v}_{i}dx,(\sigma ,\tau {)}_Q={\int }_{\Omega }{\sigma }_{ij}{\tau }_{ij}dx.$$

$${\int }_{\Omega }(u\mathrm{div}v+\nabla u\cdot v)dx={\int }_{\Gamma }u(v\cdot \nu )d\Gamma \text{mbox}\mathrm{forall}u\in W^{1,p}(\Omega ),v\in W^{1,q}(\Omega ;{\mathbb{R}}^d),$$

$$(\sigma,\varepsilon(v))_{Q}+(\mathop{\rm Div}\nolimits\sigma,\,v)_{H}=\int_{\Gamma}\sigma\nu\cdot v\,d\,\Gamma\ \ \mbox{\rm for all}\ v\in H^{1}(\Omega;\mathbb{R}^{d}),\,\,\sigma\in C^{1}(\Omega;\mbox{{\boldmath{${\mathbb{S}}$}}${}^{d}$}).$$

$$\Vert x{\Vert }_Y\le \epsilon \Vert x{\Vert }_X+C(\epsilon )\Vert x{\Vert }_Z\text{mbox}\mathrm{forall}x\in X.$$

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

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

$$\sigma =\mathcal{A}(\epsilon (\dot{u}))+\mathcal{B}\epsilon (u)+C\theta \text{mbox}in\Omega \times (0,T).$$

$$\rho \ddot{u}=\mathrm{Div}\sigma +f_0\text{mbox}in\Omega \times (0,T),$$

$$\sigma \nu =f_2\text{mbox}on{\Gamma }_N\times (0,T).$$

$$-{\sigma }_{\nu }=p(u_{\nu })\text{mbox}on{\Gamma }_C\times (0,T),$$

$$-{\sigma }_{\tau }\in \partial j({\dot{u}}_{\tau })\text{mbox}on{\Gamma }_C\times (0,T),$$

$$\rho c_p\dot{\theta }-\mathrm{div}(K\nabla \theta )=c_{ij}{\dot{u}}_{i,j}+g\text{mbox}on\Omega \times (0,T),$$

$$-K\nabla \theta \cdot \nu =-r(\theta )-h(\Vert {\dot{u}}_{\tau }{\Vert }_{{\mathbb{R}}^d})\text{mbox}on{\Gamma }_C\times (0,T).$$

$$u(0)=u_0,\dot{u}(0)=v_0,\theta (0)={\theta }_0\text{mbox}in\Omega .$$

$$f_0\in C(0,T;L^2(\Omega ;{\mathbb{R}}^d)),f_2\in C(0,T;L^2({\Gamma }_N;{\mathbb{R}}^d)).$$

$$g\in C(0,T;L^2(\Omega )).$$

$$E=H^1(\Omega ),V=H^1(\Omega ;{\mathbb{R}}^d),$$

$$(u,v{)}_V=(u,v{)}_H+(\epsilon (u),\epsilon (v){)}_Q\text{mbox}forallu,v\in V$$

$$\mathbb{V}=\{v\in L^2(0,T;V)\mid \dot{v}\in L^2(0,T;V^{\ast })\},\mathbb{E}=\{v\in L^2(0,T;E)\mid \dot{v}\in L^2(0,T;E^{\ast })\},$$

$$⟨A(t,u),v{⟩}_V=(\mathcal{A}(t,\epsilon (u)),\epsilon (v){)}_Q\mathrm{for}\mathrm{all}u,v\in V,t\in (0,T),$$

$$⟨Bu,v{⟩}_V=(\mathcal{B}\epsilon (u),\epsilon (v){)}_Q\mathrm{for}\mathrm{all}u,v\in V,$$

$$⟨C_1\theta ,v{⟩}_V={\int }_{\Omega }{c}_{ij}{v}_{i,j}\theta dx\mathrm{for}\mathrm{all}v\in V,\theta \in E,$$

$$⟨C_2(t,\theta ),\eta {⟩}_E=(K(t,\nabla \theta ),\nabla \eta {)}_{L^2(\Omega )}\mathrm{for}\mathrm{all}\theta ,\eta \in E,$$

$$⟨C_{3}v,\eta {⟩}_E=-{\int }_{\Omega }{c}_{ij}{v}_{i,j}\eta dx\mathrm{for}\mathrm{all}v\in V,\eta \in E.$$

$$⟨f(t),v{⟩}_{V^{\ast }\times V}=⟨f_0(t),v{⟩}_H+⟨f_2(t),v{⟩}_{L^2({\Gamma }_N;{\mathbb{R}}^d)},$$

Full text

Numerical analysis of a non-clamped dynamic thermoviscoelastic contact problem

Piotr Bartman* 1*, Krzysztof Bartosz* 1*,

Michał Jureczka* 1*, Paweł Szafraniec* 1*,  

 

1 Jagiellonian University, Faculty of Mathematics and Computer Science

ul. Lojasiewicza 6, 30348 Krakow, Poland

The paper is dedicated to Professor Weimin Han

on the occasion of his 60th birthday

Abstract. In this work, we analyze a non-clamped dynamic viscoelastic contact problem involving thermal effect. The friction law is described by a non-monotone relation between the tangential stress and the tangential velocity. This leads to a system of second-order inclusion for displacement and a parabolic equation for temperature. We provide a fully discrete approximation of the problem and find optimal error estimates without any smallness assumption on the data. The theoretical result is illustrated by numerical simulations.

Keywords: dynamic contact, non-clamped contact, thermoviscoelastic material, non-monotone friction law, finite element method, error estimate, numerical simulations.

2020 Mathematics Subject Classification: 65J08, 65M15, 65M50, 74B05, 74F05, 74H15, 74S05, 80M10.

1 Introduction

Physical contact processes appear in many fields of industry and everyday life. Hence, it is reasonable to study mathematical models concerning various properties and parameters of physical bodies in contact with each other. In particular, thermal effects play a special role in such models. A physical body may change its properties and shape when heated. On the other hand, contact between two physical bodies is directly related to a heat exchange. Additionally, frictional contact usually leads to heat generation. It is clear that a well-formulated mathematical model should take into account all of the above phenomena.

There are numerous publications concerning mathematical analysis of contact processes with thermal effect, which essentially deal with the existence and uniqueness of a solution for the stated problems. We recommend [1, 2, 3, 4, 5, 6, 7] as representative references in this field.

Apart from the theoretical aspect of mathematical modeling, numerical analysis also plays a significant role. Typically, it includes approximation schemes based on temporal and spatial discretization, an estimate of the error between the approximate and exact solutions, and numerical simulations. This kind of approach carried out for various mechanical contact problems with or without thermal effect is presented, for instance, in [8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

In this paper we deal with the dynamic contact of a viscoelastic body with a foundation involving a thermal effect. The contact is frictional, and the friction law is modeled by a non-monotone relation between the tangential velocity and the tangential stress. More precisely, the relation has the form of an inclusion involving Clarke subdifferential of a locally Lipschitz non-convex potential. It is worth mentioning that such inclusions are closely related to the concept of hemivariational inequalities (HVIs), which are used as a variational formulation of various non-monotone contact problems of mechanics. The theory and applications of HVIs have been broadly developed in the last few decades. For example, we refer to [18] for a comprehensive review of recent achievements in this field. The numerical approach to stationary HVIs based on the Finite Element Method (FEM) was first investigated in [19]. In recent years, similar techniques supported by numerical simulations have been applied to the numerical analysis of various contact processes in mechanics modeled by Clarke subdifferential inclusions or HVIs. We refer to [20, 21, 22, 23, 24, 25, 26] for further details.

The non-monotone character of the friction law introduces one of the main difficulties in our model. Despite this issue, we provide a numerical analysis of the problem based on the fully discrete scheme. The main result consists of the estimate of the error between the discrete and the exact solution. We show that the error estimate is linear with respect to the temporal and spatial discretization parameters in the case where the spatial discretization is based on the first-order FEM. We also present results from numerical simulations that illustrate and validate the linear rate of convergence with respect to the discretization parameters.

This paper provides a generalization of the result obtained in [23]. In contrast to the model studied there, we deal with the non-clamped body, which introduces additional difficulty to the problem. Furthermore, our results are obtained without smallness assumptions, which are required in [23]. These advances are possible due to the arguments contained in Corollary 3.

The remainder of the paper is organized as follows. In Section 2 we introduce some preliminary material. In Section 3 we formulate the mathematical model of the dynamic thermoviscoelastic contact problem and present its variational formulation. The main results on the error estimates for the fully discrete numerical scheme are presented in Section 4. Finally, in Section 5, we present numerical simulations of a two- and three-dimensional contact problem and provide numerical evidence of optimal order convergence for the linear finite elements.

2 Mathematical preliminaries

In this section, we briefly present the definitions and notation used in the paper. In what follows, we denote by $\langle\cdot,\cdot\rangle_{X}$ the duality pairing between a Banach space $X$ and its dual $X^{*}$. Throughout the paper, we denote by $C$ a generic constant whose value may change from line to line. The symbol $C(0,T;X)$ represents the space of continuous functions on $[0,T]$ with values in $X$. For a function $u\colon[0,T]\to X$ we denote by $\dot{u}$ and $\ddot{u}$ its first and second time derivatives, respectively.

We recall the definitions of the generalized directional derivative and the generalized Clarke subgradient for a locally Lipschitz function $\varphi\colon X\to\mathbb{R}$, where $X$ is a Banach space (see [27]). The generalized directional derivative of $\varphi$ at the point $x\in X$ in direction $v\in X$, denoted by $\varphi^{0}(x;v)$, is defined by

[TABLE]

The generalized subgradient of $\varphi$ at $x$, denoted by $\partial\varphi(x)$, is a subset of a dual space $X^{*}$ given by

[TABLE]

Let $\mathbb{R}^{d}$ denote an Euclidean space of dimension $d$, with the scalar product and the norm given by

[TABLE]

Hereafter, all indices $i,j,k,l$ run between $1$ and $d$. Moreover, we apply the so-called summation convention over repeated indices, e.g. the notation $a_{i}b_{i}$ means $\sum_{i=1}^{d}a_{i}b_{i}$ and $a_{ij}b_{ij}$ means $\sum_{i,j=1}^{d}a_{ij}b_{ij}$. We denote by ${\mathbb{S}}$d the space of second-order symmetric tensors in $\mathbb{R}^{d}$, or equivalently the space of symmetric matrices $d\times d$. The scalar product and the norm in ${\mathbb{S}}$d are given by

[TABLE]

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded domain with a Lipschitz boundary $\Gamma=\partial\Omega$. For any function $u\colon\Omega\to\mathbb{R}^{d}$ such that $u=(u_{i})$, $u_{i}\colon\Omega\to\mathbb{R}$, we use the symbol $u_{i,j}$ for its partial derivative $\frac{\partial u_{i}}{\partial x_{j}}$. Furthermore, we use symbols $\mathop{\rm div}\nolimits u$ and $\varepsilon(u)$ to denote its divergence operator and the deformation operator, respectively, defined by

[TABLE]

For a function $\sigma\colon\Omega\to\mathbb{S}^{d}$, such that $\sigma=(\sigma_{ij})$, $\sigma_{ij}\colon\Omega\to\mathbb{R}$, we use the symbol $\mathop{\rm Div}\nolimits\sigma$ to denote its divergence operator given by

[TABLE]

In the sequel, we will use the following spaces of functions defined in $\Omega$.

[TABLE]

The spaces $H$ and $Q$ are real Hilbert spaces with inner products given by

[TABLE]

The norms induced by the scalar products mentioned above are denoted by $\|\cdot\|_{H}$ and $\|\cdot\|_{Q}$, respectively.

Let ${\nu}=(\nu_{i})$ be the outward unit normal vector at the boundary $\partial\Omega$. We denote by $v_{\nu}$ and $v_{\tau}$ the normal and tangential parts of the vector field $v\colon\partial\Omega\to\mathbb{R}^{d}$, defined by $v_{\nu}=v\cdot\nu$ and $v_{\tau}=v-v_{\nu}\nu$, respectively. Similarly, $\sigma_{\nu}$ and $\sigma_{\tau}$ represent the normal and tangential parts of the tensor field \sigma\colon\partial\Omega\to\mbox{{\boldmath{{\mathbb{S}}}}{}^{d}} defined by $\sigma_{\nu}=(\sigma\nu)\cdot\nu$ and $\sigma_{\tau}=\sigma-\sigma_{\nu}\nu$, respectively.

We now recall two Green formulas (cf. [18], (2.6) and (2.7), respectively) needed to obtain a variational formulation of the problem of contact mechanics.

[TABLE]

where $1\leq p<\infty$, $\frac{1}{p}+\frac{1}{q}=1$ and

[TABLE]

We complete this section with the following two lemmas, which will be used throughout the paper.

Lemma 1** ([28], Lemma 8.4.12)**

Let $X$, $Y$ and $Z$ be Banach spaces such that $X$ is compactly embedded in $Y$, and $Y$ is continuously embedded in $Z$. Then, for every $\varepsilon>0$, there exists a constant $C(\varepsilon)>0$ such that

[TABLE]

Lemma 2** ([29], Lemma 7.25)**

Let $T>0$ be fixed, $N>1$ and $k=T/N$. Assume $\{e_{n}\}_{n=1}^{N}$ and $\{g_{n}\}_{n=1}^{N}$ be two sequences of non-negative numbers satisfying

[TABLE]

with a positive constant $c$ independent of $N$. Then there exists a positive constant $\hat{c}$, independent of $N$, such that

[TABLE]

3 Mechanical problem and variational formulation

In this section, we describe the mathematical model of the thermoviscoelastic dynamic contact problem in its classical and variational formulations.

We consider a body that, in its reference configuration, occupies a bounded domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ with a Lipschitz boundary $\Gamma$, consisting of two parts $\Gamma_{N}$ and $\Gamma_{C}$, such that $\Gamma_{N}\cap\Gamma_{C}=\emptyset$. The part of the boundary $\Gamma_{C}$ is in contact with a foundation. We are interested in a mathematical model that describes the evolution of body displacement, stress, and temperature. Let $T>0$ be given, and let $[0,T]$ represent the time interval considered. In what follows, we denote by $x\in\Omega$ and $t\in[0,T]$ the spatial and temporal variables, respectively. We also denote by $u=(u_{i}(x,t))$, $\sigma=(\sigma_{ij}(x,t))$ and $\theta=\theta(x,t)$ the displacement field, the stress field, and the temperature, respectively. These quantities are assumed to satisfy the following constitutive law (for simplicity, we omit the symbols $x$ and $t$ below)

[TABLE]

The symbols $\mathcal{A}$, $\mathcal{B}$ and $C$ in (3) represent the viscosity operator, the elasticity operator, and the heat expansion tensor, respectively. For specific examples of these operators, we refer to Section 5. The process is assumed to be dynamic, and hence the equation of motion takes the form of

[TABLE]

where $\rho$ is the density of the material and $f_{0}$ is the density of the volume forces. For simplicity of notation, we assume that $\rho=1$. We assume that the traction force is applied on $\Gamma_{N}$, that is,

[TABLE]

The normal compliance condition is described by relation

[TABLE]

where $p_{\nu}$ is a given non-negative function that vanishes when its argument is negative. Moreover, we assume that the frictional law takes the following multivalued, possibly non-monotone form

[TABLE]

where $j$ is a locally Lipschitz function, and $\partial j$ stands for its Clarke subgradient; see Section 2. In our model, we assume a simplified and linearized version of the energy law that takes the form

[TABLE]

where $c_{p}$ is the heat capacity, $K$ is the thermal conductivity and $g$ represents the heat sources. We assume that the heat transfer between the body and the foundation is described by the law

[TABLE]

Finally, we impose the initial conditions

[TABLE]

The mechanical problem reads as follows

Problem ${\cal P_{M}}$. *Find the displacement $u\colon\Omega\times[0,T]\rightarrow\mathbb{R}^{d}$, the stress $\sigma\colon\Omega\times[0,T]\rightarrow\mathbb{S}^{d}$ and the temperature $\theta\colon\Omega\times[0,T]\to\mathbb{R}$ that satisfy conditions $(\ref{classical_sigma})$-$(\ref{classical_intial})$.

In the study of Problem ${\cal P_{M}}$ we need the following assumptions on its data.

\underline{H(\mathcal{A})}\: The viscosity operator $\mathcal{A}:\Omega\times[0,T]\times\mathbb{S}^{d}\to\mathbb{S}^{d}$ satisfies

• (a)

$\mathcal{A}(\cdot,\cdot,\varepsilon)$ is continuous on $\Omega\times[0,T]$ for all $\varepsilon\in\mathbb{S}^{d}$;

• (b)

\big{(}\mathcal{A}(x,t,\varepsilon_{1})-\mathcal{A}(x,t,\varepsilon_{2})\big{)}:(\varepsilon_{1}-\varepsilon_{2})\geq m_{\mathcal{A}}\|\varepsilon_{1}-\varepsilon_{2}\|_{\mathbb{S}^{d}}^{2} for all $\varepsilon_{1},\varepsilon_{2}\in\mathbb{S}^{d}$, a.e. $(x,t)\in\Omega\times(0,T)$ with $m_{\mathcal{A}}>0$;

• (c)

$\|\mathcal{A}(x,t,\varepsilon_{1})-\mathcal{A}(x,t,\varepsilon_{2})\|_{\mathbb{S}^{d}}\leq L_{\mathcal{A}}\|\varepsilon_{1}-\varepsilon_{2}\|_{\mathbb{S}^{d}}$ for all $\varepsilon_{1},\varepsilon_{2}\in\mathbb{S}^{d}$, a.e. $(x,t)\in\Omega\times(0,T)$ with $L_{\mathcal{A}}>0$.

\underline{H(\mathcal{B})}\: The elasticity operator $\mathcal{B}:\Omega\times\mathbb{S}^{d}\to\mathbb{S}^{d}$ is bounded, symmetric, positive fourth-order tensor, i.e.

• (a)

$\mathcal{B}_{ijkl}\in L^{\infty}(\Omega)$, $1\leq i,j,k,l\leq d$;

• (b)

$\mathcal{B}\sigma:\tau=\sigma:\mathcal{B}\tau$ for all $\sigma,\tau\ \in\mathbb{S}^{d}$, a.e. in $\Omega$;

• (c)

$\mathcal{B}\tau:\tau\geq 0$ for all $\tau\in\mathbb{S}^{d}$, a.e. in $\Omega$.

\underline{H(C)}\: The thermal expansion tensor $C=(c_{ij})$ satisfies

• (a)

$c_{ij}\in L^{\infty}(\Omega)$, $1\leq i,j\leq d$;

• (b)

$c_{ij}=c_{ji}$, $1\leq i,j\leq d$.

$\underline{H(K)}:$ The thermal conductivity operator $\displaystyle K\colon\Omega\times[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ satisfies

• (a)

$K(\cdot,\cdot,\xi)$ is continuous on $\Omega\times[0,T]$ for all $\xi\in\mathbb{R}^{d}$;

• (b)

\displaystyle\big{(}K(x,t,\xi_{1})-K(x,t,\xi_{2})\big{)}\cdot(\xi_{1}-\xi_{2})\geq m_{K}\,\|\xi_{1}-\xi_{2}\|_{\mathbb{R}^{d}}^{2} for all $\xi_{1}$, $\xi_{2}\in\mathbb{R}^{d}$, a.e. $(x,t)\in\Omega\times(0,T)$ with $m_{K}>0$;

• (c)

$\|K(x,t,\xi_{1})-K(x,t,\xi_{2})\|_{\mathbb{R}^{d}}\leq L_{K}\|\xi_{1}-\xi_{2}\|_{\mathbb{R}^{d}}$ for all $\xi_{1},\xi_{2}\in\mathbb{R}^{d}$, a.e. $(x,t)\in\Omega\times(0,T)$ with $L_{K}>0$.

$\underline{H(p)}:$ The function $\displaystyle p\colon\Gamma_{C}\times[0,T]\times\mathbb{R}\to\mathbb{R}$ satisfies

• (a)

$p(\cdot,\cdot,\xi)$ is measurable on $\Gamma_{C}\times(0,T)$ for all $\xi\in\mathbb{R}^{d}$;

• (b)

$\displaystyle|p(x,t,\xi_{1})-p(x,t,\xi_{2})|\leq L_{p}\,|\xi_{1}-\xi_{2}|$ for all $\xi_{1}$, $\xi_{2}\in\mathbb{R}$, a.e. $(x,t)\in\Gamma_{C}\times(0,T)$ with $L_{p}>0$.

\underline{H(j)}\: The potential $j\colon\Gamma_{C}\times[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ satisfies

• (a)

$j(\cdot,\cdot,\xi)$ is measurable on $\Gamma_{C}\times(0,T)$ for all $\xi\in\mathbb{R}^{d}$;

• (b)

$j(x,t,\cdot)$ is locally Lipschitz on $\mathbb{R}^{d}$ for a.e. $(x,t)\in\Gamma_{C}\times(0,T)$;

• (c)

$\|\partial j(x,t,\xi)\|_{\mathbb{R}^{d}}\leq c_{j}$ for all $\xi\in\mathbb{R}^{d}$, a.e. $(x,t)\in\Gamma_{C}\times(0,T)$ with $c_{j}\geq 0$;

• (d)

$(\zeta_{1}-\zeta_{2})\cdot(\xi_{1}-\xi_{2})\geq-m_{j}\|\xi_{1}-\xi_{2}\|^{2}_{\mathbb{R}^{d}}$ for all $\zeta_{i}\in\partial j(x,t,\xi_{i})$, $\xi_{i}\in\mathbb{R}^{d}$, $i=1$, $2$, a.e. $(x,t)\in\Gamma_{C}\times(0,T)$ with $m_{j}\geq 0$.

$\underline{H(r)}:$ The function $\displaystyle r\colon\Gamma_{C}\times\mathbb{R}\to\mathbb{R}$ satisfies

• (a)

$r(\cdot,\zeta)\in L^{2}(\Gamma_{C})$ for all $\zeta\in\mathbb{R}$;

• (b)

$\displaystyle|r(x,\zeta_{1})-r(x,\zeta_{2})|\leq L_{r}\,|\zeta_{1}-\zeta_{2}|$ for all $\zeta_{1}$, $\zeta_{2}\in\mathbb{R}$, a.e. $x\in\Gamma_{C}$ with $L_{r}>0$.

$\underline{H(h)}:$ The function $\displaystyle h\colon\Gamma_{C}\times[0,\infty)\to[0,\infty)$ satisfies

• (a)

$h(\cdot,\zeta)\in L^{2}(\Gamma_{C})$ for all $\zeta\in[0,\infty)$;

• (b)

$\displaystyle|h(x,\zeta_{1})-h(x,\zeta_{2})|\leq L_{h}\,|\zeta_{1}-\zeta_{2}|$ for all $\zeta_{1}$, $\zeta_{2}\in[0,\infty)$, a.e. $x\in\Gamma_{C}$ with $L_{h}>0$.

\underline{H(f)}\: The volume forces and traction densities satisfy

[TABLE]

\underline{H(g)}\: The heat source satisfies

[TABLE]

We remark that the result presented in the sequel can be obtained under a weaker assumption than $H(j)(c)$. Namely it is enough to assume a sublinear growth of the subgradient $\partial j$ with respect to the last variable. However, for the sake of simplicity we assume that it is bounded.

In order to provide a variational formulation of the problem, we define spaces

[TABLE]

where $E$ is equipped with the classical $H^{1}$ norm, and in $V$ we define the inner product

[TABLE]

and the corresponding norm $\|u\|^{2}:=(u,u)_{V}$ for all $u\in V$. It follows that the norms $\|\cdot\|_{V}$ and $\|\cdot\|_{H^{1}(\Omega;\mathbb{R}^{d})}$ are equivalent. We know that the embedding $i\colon E\to H^{1-\delta}(\Omega)$ is compact for $\delta\in(0,\frac{1}{2})$, and the trace operator $\overline{\gamma}\colon H^{1-\delta}(\Omega)\to L^{2}(\Gamma)$ is continuous. Hence, we see that $\gamma=\overline{\gamma}i\colon E\to L^{2}(\Gamma)$ is compact. We omit the notation of the embedding and the trace operator, and write simply $v$ instead of $\gamma v$ for $v\in E$. We proceed similarly for the space $V$. We introduce the following spaces of vector-valued functions

[TABLE]

where $\dot{v}$ denotes time derivative of $v$ in the sense of distributions.

Next, we define operators $A\colon(0,T)\times V\to V^{*}$, $B\colon V\to V^{*}$, $C_{1}\colon E\to V^{*}$, $C_{2}\colon(0,T)\times E\to E^{*}$, $C_{3}\colon V\to E^{*}$ by

[TABLE]

We define the functions $f\colon(0,T)\to V^{*}$ and $\tilde{g}\colon(0,T)\to E^{*}$ for all $v\in V$, $\eta\in E$, a.e. $t\in(0,T)$ by

[TABLE]

We also introduce an additional hypothesis on the data

$\underline{H_{0}}\ :$ $u_{0}\in V$, $v_{0}\in H$, $\theta_{0}\in E$.

Now we are in a position to present a weak formulation of Problem ${\cal P_{M}}$.

Problem ${\cal P_{E}}$. Find a displacement $u\in L^{2}(0,T;V)$ with $\dot{u}\in\mathbb{V}$, and a temperature $\theta\in\mathbb{E}$ such that

[TABLE]

Problem ${\cal P_{E}}$ can be obtained from equations (4) and (8) by multiplying them by test functions $w\in V$ and $\eta\in E$, respectively, performing the Green formulas (1) and (2), applying conditions (3)–(10), using definitions of operators above and denoting $\xi=\sigma_{\tau}$. Every solution of Problem ${\cal P_{E}}$ is called a weak solution of Problem ${\cal P_{M}}$.

Existence and uniqueness of the solution to Problem ${\cal P_{E}}$ can be proved following the proof of Theorem 9 in [30] and using Corollary 3. Since this paper is devoted to numerical analysis, we omit this proof. For more details on non-clamped problems see e.g. [31] and for problems with no smallness assumption see e.g. [32].

Corollary 3

Lemma 1 holds for spaces $X=H^{1}(\Omega)$, $Y=H^{1-\delta}(\Omega)$, $Z=L^{2}(\Omega)$ for $\delta\in(0,1/2)$. Moreover, from the trace theorem, we have $\|v\|_{L^{2}(\Gamma)}\leq C\|v\|_{H^{1-\delta}(\Omega)}$, hence for all $\varepsilon>0$ there exists $C(\varepsilon)>0$ such that

[TABLE]

Similarly, we obtain

[TABLE]

4 Discretization and error estimates

In this section we consider a fully discrete approximation of Problem ${\cal P_{E}}$ and examine error estimates for the approximate solution and the exact one.

Let $E^{h}$ and $V^{h}$ be finite-dimensional subspaces of $E$ and $V$, respectively, with $h>0$ denoting the spatial discretization parameter. Let $u_{0}^{h},v_{0}^{h}\in V^{h}$ and $\theta_{0}^{h}\in E^{h}$ be suitable approximations of $u_{0},v_{0},\theta_{0}$. We define a uniform partition of $[0,T]$ denoted by $0=t_{0}<t_{1}<\ldots<t_{N}=T$, where $N\in\mathbb{N}$, $t_{j}=jk$ for $j=1,...,N$ and $k=\frac{T}{N}$ is the time step. For a continuous function $\varphi$ we use the notation $\varphi_{n}=\varphi(t_{n})$. For a sequence $\{z_{n}\}_{n=0}^{N}$ we denote by $\delta z_{n}=(z_{n}-z_{n-1})/k$ for $n=1,\ldots,N$ its backward divided difference.

Let $(u,\theta)$ be a solution of Problem ${\cal P_{E}}$ and $v=\dot{u}$. The fully discrete approximation of Problem ${\cal P_{E}}$ is the following.

Problem ${\cal P}_{\cal E}^{hk}$. Find sequences of velocity $\{v_{n}^{hk}\}_{n=1}^{N}\subset V^{h}$ and temperature $\{\theta_{n}^{hk}\}_{n=1}^{N}\subset E^{h}$ such that there exists $\{\xi_{n}^{hk}\}_{n=1}^{N}\subset L^{2}(\Gamma_{C};\mathbb{R}^{d})$ satisfying for $j=1,\ldots,N$

[TABLE]

where $v_{0}^{hk}=v_{0}^{h},\theta_{0}^{hk}=\theta_{0}^{h}$ and $\{u_{n}^{hk}\}_{n=1}^{N}\subset V^{h}$ is given by

[TABLE]

*with $u_{0}^{hk}=u_{0}^{h}$. *

Now we deal with the error analysis for the discrete scheme, Problem ${\cal P}_{\cal E}^{hk}$.

Let $\{u_{n}^{hk}\}_{n=1}^{N}$, $\{v_{n}^{hk}\}_{n=1}^{N}$ and $\{\theta_{n}^{hk}\}_{n=1}^{N}$ denote the discrete approximation of displacement $u$, velocity $v$ and temperature $\theta$, respectively, derived by the numerical scheme in Problem ${\cal P}_{\cal E}^{hk}$. We prove the following estimate.

Theorem 4

Let the assumptions $H({\cal A})$, $H({\cal B})$, $H(C)$, $H(K)$, $H(j)$, $H(r)$, $H(h)$, $H(f)$, $H(g)$, and $H_{0}$ hold. Moreover, we impose the following additional regularity conditions $u\in C^{2}(0,T;H)\cap C^{1}(0,T;V)$, $\dot{u}_{\tau}\in C(0,T;L^{2}(\Gamma_{C};\mathbb{R}^{d}))$, $\theta\in C(0,T;L^{2}(\Gamma_{C}))$. Then, there exists $C>0$ such that for all $\{w_{n}^{h}\}_{n=1}^{N}\subset V^{h}$ and $\{\eta_{n}^{h}\}_{n=1}^{N}\subset E^{h}$ we have

[TABLE]

Proof. Take $w^{h}\in V^{h}$ and $\eta^{h}\in E^{h}$ as test functions in both Problem ${\cal P_{E}}$ and Problem ${\cal P}_{\cal E}^{hk}$, respectively, and subtract relevant equations to obtain

[TABLE]

for $j=1,\ldots,N.$ Let $\{w_{j}^{h}\}_{j=1}^{N}\subset V$ and $\{\eta_{j}^{h}\}_{j=1}^{N}\subset E$ be arbitrary. Applying (19) with $w^{h}=v_{j}^{hk}$, $\eta^{h}=\theta_{j}^{hk}$ and then with $w^{h}=w_{j}^{h}$, $\eta^{h}=\eta_{j}^{h}$, we find

[TABLE]

Using the fact that $\langle C_{1}\eta,w\rangle_{V}+\langle C_{3}w,\eta\rangle_{E}=0$ for all $w\in V$, $\eta\in E$, we get

[TABLE]

From the formula $2\left\langle a-b,a\right\rangle_{H}=\|a-b\|_{H}^{2}+\|a\|_{H}^{2}-\|b\|_{H}^{2}$ for $a,b\in H$, we get

[TABLE]

By $H({\cal A})(b)$, $H(j)(d)$ and Corollary 3, we get

[TABLE]

Similarly, we obtain

[TABLE]

Next, from the Lipschitz continuity of ${\cal A}$ and $K$, we have

[TABLE]

[TABLE]

By the linearity of $B$, $C_{1}$ and $C_{3}$, we obtain

[TABLE]

From $H(p)$, and the continuity of the trace operator $\gamma\colon V\to L^{2}(\Gamma;\mathbb{R}^{d})$, we deduce

[TABLE]

By $H(j)(c)$, $H(r)(b)$, $H(h)(b)$, and the continuity of the trace operator $\gamma\colon E\to L^{2}(\Gamma)$, we get

[TABLE]

The symbol ${\rm meas}(\Gamma_{C})$ in (32) represents the Lebesgue surface measure of $\Gamma_{C}$. Finally, we estimate

[TABLE]

Applying (21)-(36) in (20), we deduce

[TABLE]

Summing up the last inequality from $j=1$ to $j=n$, and taking into account that $v_{0}=\dot{u}(0)=v_{0}$, we obtain the following

[TABLE]

for all $\{w_{j}^{h}\}_{j=1}^{n}\subset V^{h}$ and $\{\eta_{j}^{h}\}_{j=1}^{n}\subset E^{h}$. We have

[TABLE]

and

[TABLE]

Next, we estimate

[TABLE]

where $I_{j}$ is the integration error

[TABLE]

Since $I_{j}\leq k\|u\|_{H^{2}(0,T;V)}$, we get

[TABLE]

Using the fact that $Nk=T$, we have the estimate

[TABLE]

We denote

[TABLE]

and

[TABLE]

Then, we find that

[TABLE]

with $C>0$. Using Lemma 2, we obtain the desired estimate (4). $\hfill{\Box}$

Now we consider a special situation in which the spaces $E^{h}$ and $V^{h}$ are based on the piecewise affine finite element approximation. Assume that $\Omega$ is a polygonal domain in $\mathbb{R}^{2}$ or a polyhedral domain in $\mathbb{R}^{3}$. Then each of parts of the boundary $\Gamma_{N}$ and $\Gamma_{C}$ is a sum of segments (in the case $\Omega\subset\mathbb{R}^{2}$) or polygons (if $\Omega\subset\mathbb{R}^{3}$), having mutually disjoint interiors i.e.

[TABLE]

Let $\{\mathcal{T}^{h}\}$ be a family of regular triangulations of the set $\overline{\Omega}$ of diameter $h$ (cf. [33], Chapter 3, § 3.1) which introduces a division of $\overline{\Omega}$ into triangles/tetrahedrons compatible with the division of the boundary $\partial\Omega$ into parts $\Gamma_{j,i}$, $1\leq i\leq N_{j}$, $j\in\{N,C\}$, in that sense that if any edge of a triangle/any face of tetrahedron that belongs to $\mathcal{T}^{h}$ has an intersection of positive boundary measure with any of sets $\Gamma_{j,i}$, then the whole edge/face is contained in $\Gamma_{j,i}$. For an element $T\in\mathcal{T}^{h}$, let $P_{1}(T)$ denote the space of first-order polynomials on $T$. Then $E^{h}$ and $V^{h}$ are defined as the spaces of continuous piecewise affine functions on elements $T\in T^{h}$, that is,

[TABLE]

Corollary 5

Suppose that the solution $(u,\theta)$ of Problem ${\cal P_{E}}$ satisfies the following regularity assumptions

[TABLE]

Moreover, we assume that the initial data satisfy

[TABLE]

Then, for $h<1$, we have the optimal order error estimate

[TABLE]

Proof. To prove (5), we will use Theorem 4 and estimate each of the terms on the right-hand side of (4). We start with the initial conditions. Let us define elements $u_{0}^{h}=\Pi^{h}u_{0}$, $v_{0}^{h}=\Pi^{h}v_{0}$, $\theta_{0}^{h}=\Pi^{h}\theta_{0}$ as finite element interpolants of $u_{0}$, $v_{0}$ and $\theta_{0}$, respectively (see [33], (2.3.29)) By the standard finite element interpolation error estimates ([34, 35, 33]) we have the following

[TABLE]

Let $w^{h}(t)=\Pi^{h}\dot{u}(t)\in V^{h}$ be the finite element interpolant of $\dot{u}(t)$, for $t\in(0,T]$. Note that $w_{\tau}^{h}(t)=(\Pi^{h}\dot{u}(t))_{\tau}$ is the continuous piecewise-linear interpolant of $\dot{u}_{\tau}(t)$ on $\Gamma_{C}$. Similarly, let $\eta^{h}(t)=\Pi^{h}\theta(t)\in E^{h}$ be the finite element interpolant of $\theta(t)$, a.e. $t\in[0,T]$. Then we have

[TABLE]

Applying (4) with $w_{n}^{h}:=w^{h}(t_{n})$ and $\eta_{n}^{h}:=\eta^{h}(t_{n})$, $n=1,\dots,N$, we get the following estimates

[TABLE]

Moreover, similarly to Section 5 of [14], we get

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Then the error bound (5) follows from (4) combined with (37). The proof of the corollary is complete. $\hfill{\Box}$

5 Numerical simulations

In this section, we present the results of our numerical simulations. Such simulations are conducted to compare the empirical and theoretical convergence results of the proposed fully discrete scheme and visualize the role of various factors in the considered model. In Section 5.1, an academic two-dimensional sample problem is solved to obtain empirical error estimates and compare them with the theoretical estimates presented in Corollary 5. Next, in Section 5.2 we consider the three-dimensional problem with a more complex shape of the body.

We employ FEM and use the definition of spaces $E^{h}$ and $W^{h}$ presented in (40) and (41), respectively. For the sake of clarity, we assume dimensionless models with mass density and heat capacity $\rho=c_{p}=1$. In our simulations, we use the original open source package conmech (https://github.com/KOS-UJ/conmech). The code is written in Python and uses just-in-time compilation mechanisms provided by packages numba [36] and JAX [37].

For two-dimensional simulations, we use the direct optimization method with the Schur complement method to decrease the dimension of the system [38]. We use Powell’s conjugate direction method from the SciPy [39] package since the minimized functional is not necessarily differentiable.

In the case of three-dimensional simulations, to speed up the calculations, we perform computations on GPU using the L-BFGS solver implemented using JAX. We use TetWild [40] to generate the volumetric mesh from the surface mesh obtained from the Stanford 3D Scanning Repository. The results are visualized using Blender ([41]).

5.1 Two-dimensional numerical example

We consider an isotropic homogeneous domain $\Omega=[0,1.5]\times[0,1]$, set $T=1$, and assume the following boundary conditions

[TABLE]

The domain $\Omega$ represents the cross-section of a three-dimensional viscoelastic body. We use the Kelvin-Voigt type short-memory viscoelastic law. The viscosity operator $\mathcal{A}$, elasticity operator $\mathcal{B}$, heat expansion tensor $C$ and the thermal conductivity operator $K$ are defined by

[TABLE]

Here, $I\in\mathbb{R}^{2\times 2}$ is the identity matrix, $\text{tr}(\cdot)$ is the trace of a square matrix, $\lambda$ and $\mu$ are the Lame coefficients, while $\xi$ and $\phi$ represent the viscosity coefficients, $\lambda,\mu,\xi,\phi,c,k>0$. In the initial state, the body rests on the foundation, that is, the initial displacement $u_{0}$ and the initial velocity $v_{0}$ are equal to zero at any point in $\Omega$, and the gap between the body and the foundation $g_{0}=0$ at any point on the contact boundary $\Gamma_{C}$. The viscoelastic properties of the body are characterized by $\phi=4.5,\xi=10.5,\nu=45,\lambda=105$. The thermal effect is parameterized by the heat expansion coefficient $c=0.5$ and the thermal conductivity coefficient $\kappa=0.033$. The internal heat source $g$ and the heat exchange between the body and the foundation $r$ are modeled by the functions

[TABLE]

where $r_{F}=0.1$ and $\theta_{F}=0.27$ denote the heat exchange coefficient and the constant foundation temperature, respectively. We assume the initial temperature $\theta_{0}=0.25$, constant density of the volume forces $f_{0}=(0,\,-1)^{T}$ inside the body, and non-zero traction force applying squeezing from both sides of the body, namely

[TABLE]

On the contact boundary $\Gamma_{C}$ we set

[TABLE]

as conditions governing the response in the normal direction, friction, and heat generated by friction. Note that $j(\cdot)$ is in the form of a non-differentiable and non-convex function.

Figure 1 demonstrates the solution of the sample problem in which we squeeze the body that lies on the foundation. The traction force from the left side is greater than from the right side, so we can observe slipping of the body on the foundation and therefore friction and heat generated by friction. On the right side, we can see the corresponding relative temperature error computed in comparison with the reference solution with $h=2^{-7}$, $k=2^{-9}$ and with 148419 degrees of freedom.

Figure 2 depicts the empirical error estimation for the sample problem considered plotted on a log-log scale. The color of the plot corresponds to the time step used to compute the solution. As we can see, for any fixed spatial step $h$ we can observe that the numerical error goes to zero as the time step goes to zero. The gray line shows the linear order of convergence for comparison with the theoretical results.

5.2 Three-dimensional numerical example

To visualize the influence of various modeled factors, we set $d=3$ and consider a volumetric mesh consisting of 42684 nodes and 230213 tetrahedrons. We set $\Gamma_{N}=\emptyset$ and $\Gamma_{C}=\Gamma$. The final time $T=2.5$, the time step $k=0.01$, and the results are presented at times $t\in\{0.5,\,1,\,1.5,\,2,\,2.5\}$. The position and shape of the foundation can be deduced from the attached illustrations. Body temperature is once again represented by its color, where blue is the coldest region and yellow is the hottest. We apply analogous constitutive laws as in a two-dimensional example and set

[TABLE]

The initial displacement $u_{0}$ and the initial velocity $v_{0}$ are equal to zero at any point in $\Omega$. The internal force, the internal heat source and the heat exchange between the body and the foundation are modeled by

[TABLE]

On the contact boundary $\Gamma_{C}$ we take the following data

[TABLE]

The values of the hardness of the foundation $\hat{p}$, the friction coefficient $\hat{j}$ and the coefficient of heat generated by friction $\hat{h}$ vary between examples.

We now provide a brief overview of the results. The first simulation shows the data $\hat{p}=200$, $\hat{j}=0.1$ and $\hat{h}=0.2$ and is presented in figure 6. This is the base configuration that will be modified in the following examples. We push the body down and to the right with force ${f}_{0}$. The body falls to the obstacle, rebounds, and slides to the right. We can also observe the dissipation of heat generated by friction.

The second simulation shows the data $\hat{p}=200$, $\hat{j}=1.0$ and $\hat{h}=0.2$ and is presented in Figure 6. The only change compared to the base configuration is the increase in friction coefficient. This causes the body to cover a shorter distance and, as a result of the lower speed of the nodes at the contact boundary, the heat generated by friction is decreased.

The third simulation shows the data $\hat{p}=200$, $\hat{j}=0.1$ and $\hat{h}=0.6$ and is presented in Figure 6. In this case, the coefficient of heat generated by friction is increased compared to the base configuration. This increases the heat generated at the contact boundary, which is then diffused inside the body.

The fourth simulation shows the data $\hat{p}=50$, $\hat{j}=0.1$ and $\hat{h}=0.2$ and is presented in Figure 6. In this case, the boundary hardness coefficient is decreased. The softer foundation causes the body to initially dive further, and the resulting rebound produces a less stable trajectory and body rotation. We remark that the infinitesimal strain tensor considered in this paper results in a constitutive law that is not rotationally invariant. As a result, its application in the case of large body rotation would cause the volume of the body to increase unnaturally.

6 Acknowledgements

The research was supported by the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Sklodowska-Curie grant agreement No. 823731 CONMECH, the Ministry of Science and Higher Education of Republic of Poland under Grant Nos. 4004/GGPJII/H2020/2018/0 and 440328/PnH2/2019, and the National Science Centre of Poland under Project No. 2021/41/B/ST1/01636.