# Contact constraints within coupled thermomechanical analysis-A finite element model

Peter Wriggers, Christian Miehe

► To cite this version:

Peter Wriggers, Christian Miehe. Contact constraints within coupled thermomechanical analysis-A finite element model. Computer Methods in Applied Mechanics and Engineering, 1994, 113 (3-4), pp.301 - 319. 10.1016/0045-7825(94)90051-5. hal-01574175

HAL Id: hal-01574175

https://hal.science/hal-01574175

Submitted on 11 Aug 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Contact constraints within coupled thermomechanical analysis – A finite element model

P. Wriggers*

Institut für Mechanik, Technische Hochschule Darmstadt, Hochschulstrasse 1, D-6100 Darmstadt, Germany

C. Miehe

Institut für Baumechanik und Numerische Mechanik, Universität Hannover, Appelstrasse 9a, D-3000 Hannover 1, Germany

Based on a sophisticated interface model a thermomechanical contact formulation is derived. The development addresses theoretical and numerical aspects. It starts from a continuous formulation of the local contact geometry and includes constitutive equations for the contact stresses, the contact heat flux and the frictional dissipation. Based on these considerations, a finite element model for large deformation processes is developed. An operator split technique concerning the solution of the global coupled equations is applied for the algorithmic treatment of the thermomechanical coupling. By means of three numerical examples, the theoretical and numerical formulations are validated.

## 1. Introduction

In engineering problems such as shrink fitting, deep drawing or hot rolling, the effects of heat generation and heat transfer have to be considered during the process. Since many of these simulations involve contact problems with dry friction, it is of interest to develop contact laws that are also able to predict, besides normal and tangential stresses, the amount of heat generated due to friction and the heat transfer across the contact surface.

During recent years, a series of publication was devoted to the derivation of thermomechanical interface models that take into account the microgeometrical shape of macroscopically smooth surfaces (see e.g. [1]). The associated theoretical studies consider a statistical characterization of the surface profile and of the related thermomechanical phenomena and usually predefine the shape and the mechanical behaviour of the asperities. The derived constitutive relations for the thermomechanical contact interface also include effects due to the finishing treatments of the surface layer which involve different characteristics from the bulk material. Hence, the effective yield pressure (determined by hardness tests) depends on the depth of the plastic zone [2].

Connected to these sophisticated interface models for the heat transfer within the contact area are constitutive relations for the mechanical contact which are also based on statistical models (see e.g. [3] for an overview). Here we have to distinguish between constitutive equations for the normal contact and the tangential contact. For the normal contact, statistical models have been stated in [4]. For the
tangential contact, which splits into the stick and slip part, several constitutive relations have been
developed (see e.g. [4]).

Numerical formulations of contact problems within the finite element method have been given in many papers. The standard approach is to model the contact surface as an idealized smooth interface which then leads to a condition describing the non-penetration by a pure geometrical constraint. Within this assumption, several approaches are possible for the variational formulation of contact problems in case of finite deformations. The so-called Lagrangian multiplier technique enforces the contact constraints exactly, (see e.g. [5, 6]). A consistent formulation for implicit computations using penalty methods is given in e.g. [7] based on a discrete model discussed in [8]. Recently also constitutive equations within the contact interface have been considered in finite element approximations (see e.g. [9]).

Finite element formulations for the frictional contact have been discussed e.g. in [3, 9–12]. Finite element models which include thermomechanical coupling for contact problems have been presented in e.g. [13] for lubricated contact, in [14] for large deformation frictionless contact and for small deformations in [15] including frictional heating.

In this paper we focus on a finite element formulation of a sophisticated interface model that includes frictional contact, heat exchange in the contact area and frictional heating. For this purpose, we give (i) a continuous formulation of the local contact geometry based on a parametrization of a moving master surface, where we restrict our consideration to two-dimensional problems. This formulation extends the formulation proposed in [12] for the rigid obstacle problem. Then we postulate (ii) constitutive equations for the contact stress, the contact heat flux and the frictional dissipation on the contact surface. These equations have model character but exhibit essential features of known experimental and microgeometrical investigations of contact interfaces. After having formulated (iii) the structure of the initial boundary value problem of finite strain thermoelasticity combined with thermomechanical contact, we discuss (iv) two main algorithmic aspects. Firstly we discuss the solution algorithm for the coupled global equations. Global solution algorithms for coupled thermomechanical analysis have been formulated in [16–18]. Here we use an algorithm that is based on an operator split of the global thermomechanical evolution operator. This strategy was recently introduced in the context of finite strain thermoplasticity without contact constraints in [19, 20] and leads to an algorithmic decoupling of the coupled problem within the time step by preserving the symmetry properties of the subproblems. Secondly we consider an extension of the return map algorithm proposed in [21] for the update of the contact stress and frictional dissipation. Finally we discuss (v) the spatial discretization of the contact surfaces and show as an interesting aspect the transition of the continuous contact geometry to the discrete contact geometry of the three node contact element used e.g. in [9] for frictional contact. Several numerical examples show the performance and the accuracy of the derived finite numerical formulation.

## 2. Contact geometry

Let $\mathcal{B}^\alpha$, $\alpha = 1, 2$, denote two bodies of interest and $\varphi_i^\alpha : \mathcal{B}^\alpha \to \mathbb{R}^3$ the associated deformation maps at time $t \in \mathbb{R}_+$. $\varphi_i^\alpha$ maps points $X^\alpha \in \mathcal{B}^\alpha$ of the reference configuration onto points $x^\alpha = \varphi_i^\alpha(X^\alpha)$ of the current configuration.

Assume that both bodies come into contact. Motivated by micromechanical investigations of contact problems, we view mechanical contact as a (microscopical) penetration of the current mathematical boundaries $\varphi_i^\alpha(\Gamma_c^\alpha)$ where $\Gamma_c^\alpha \subset \partial\mathcal{B}^\alpha$ are possible contact surfaces of the bodies $\mathcal{B}^\alpha$; see Fig. 1 for an illustration of this idea. In what follows we denote $\varphi_i^1(\Gamma_c^1)$ as the current slave surface which penetrates in the case of contact into the current master surface $\varphi_i^2(\Gamma_c^2)$. Within our formulation of the contact geometry, the latter plays the role of a (moving) reference surface. We parameterize the master surface in its reference and current configuration by the natural parameters $\xi^1, \xi^2$ (for 2-D problems we only need $\xi^1$), i.e. we consider material surfaces $(\xi^1, \xi^2) \mapsto X^2 = \hat{X}^2(\xi^1, \xi^2) \subset \Gamma_c^2$ and $(\xi^1, \xi^2) \mapsto x^2 = \hat{x}^2(\xi^1, \xi^2) \subset \varphi_i^2(\Gamma_c^2)$ (see Fig. 1).
<figure><img src="image_2.png" /><figcaption>Fig. 1. Basic contact geometry.</figcaption></figure>

Then the local deformation gradient of the master surface is given by $F_t^2 := a_α^2 \otimes A^{2α}$ based on the tangent vectors of the contact surface $a_α^2 := \hat{x}_{t,α}^2(\xi^1, \xi^2)$ and $A_α^2 := \tilde{X}_{,α}^2(\xi^1, \xi^2)$ with the standard relations $a_α^2 \cdot a^{2β} = δ_α^β$ and $A_α^2 \cdot A^{2β} = δ_α^β$.

## 2.1. Penetration

As the first relevant function for the contact geometry, we define a penetration function on the current slave surface $\varphi_t^1(\Gamma_c^1)$ by setting

$$g_N = \begin{cases} \|\mathbf{x}^1 - \hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)\|, & \text{for } [\mathbf{x}^1 - \hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \bar{\mathbf{n}}^2 < 0, \\ 0, & \text{otherwise.} \end{cases} \quad (2.1)$$

Here $(\bar{\xi}^1, \bar{\xi}^2)$ is the minimizer of the distance function (see Fig. 2),

$$\hat{d}^1(\xi^1, \xi^2) = \|\mathbf{x}^1 - \hat{\mathbf{x}}_t^2(\xi^1, \xi^2)\| \rightarrow \text{MIN} \quad (2.2)$$

for a given slave point $\mathbf{x}^1$. $(\bar{\xi}^1, \bar{\xi}^2)$ are obtained by writing the necessary condition

$$\frac{d}{d\xi^\alpha} \hat{d}^1(\xi^1, \xi^2) = - \frac{\mathbf{x}^1 - \hat{\mathbf{x}}_t^2(\xi^1, \xi^2)}{\|\mathbf{x}^1 - \hat{\mathbf{x}}_t^2(\xi^1, \xi^2)\|} \cdot \hat{\mathbf{x}}_{t,\alpha}^2(\xi^1, \xi^2) = 0 \quad (2.3)$$

for the minimum of the distance function (2.2). Clearly the first term in Eq. (2.3) denotes the normal $\bar{\mathbf{n}}^2$ and $\hat{\mathbf{x}}_{t,α}^2(\xi^1, \xi^2)$ is the tangent vector $\bar{\mathbf{a}}_α^2$. Thus, we have the condition $-\bar{\mathbf{n}}^2 \cdot \bar{\mathbf{a}}_α^2 = 0$ which means that the current master point $\hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)$ is the orthogonal projection of a given slave point $\mathbf{x}^1$ onto the current master surface $\varphi_t^2(\Gamma_c^2)$.

$\bar{\mathbf{n}}^2 := (\bar{\mathbf{a}}_1^2 \times \bar{\mathbf{a}}_2^2) / \|\bar{\mathbf{a}}_1^2 \times \bar{\mathbf{a}}_2^2\|$ is the outward unit normal on the current master surface at the master point where $\bar{\mathbf{a}}_α^2$ are tangent vectors at $\hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)$. Here and in the following we denote by a bar over a quantity its evaluation at the minimal distance point $(\bar{\xi}^1, \bar{\xi}^2)$.

The penetration function gives two pieces of informations: Firstly it serves as a local contact check, i.e. we set contact $\Leftrightarrow g_N > 0$. Secondly it enters for $g_N > 0$ as a local kinematical variable the constitutive function for the contact pressure (see Section 3.1).

By taking the time derivative of Eq. (2.2) at the minimal distance point $(\bar{\xi}^1, \bar{\xi}^2)$, one obtains, in the case of contact, the rate of penetration

$$-\dot{g}_N = [\mathbf{v}_t^1 - \hat{\mathbf{v}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \bar{\mathbf{n}}^2 \quad (2.4)$$

for given spatial velocities $\mathbf{v}_t^1$ and $\hat{\mathbf{v}}_t^2(\bar{\xi}^1, \bar{\xi}^2)$ at the slave and master points.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Penetration: solution of minimum distance problem.</figcaption></figure>

## 2.2. Tangential relative velocity and tangential relative slip

The tangential relative slip between two bodies is related to the change of the solution point $(\bar{\xi}^1, \bar{\xi}^2)$ of the minimal distance problem. Thus, we can compute the time derivative of $\xi^\alpha$ from Eq. (2.3). This yields the following result:

$$ \frac{d}{dt} [\mathbf{x}_t^1 - \hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \bar{\mathbf{a}}_\alpha^2 = [\mathbf{v}_t^1 - \hat{\mathbf{v}}_t^2(\bar{\xi}^1, \bar{\xi}^2) - \bar{\mathbf{a}}_\beta \dot{\bar{\xi}}^\beta] \cdot \bar{\mathbf{a}}_\alpha^2 + [\mathbf{x}_t^1 - \hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \hat{\mathbf{a}}_\alpha^2 = 0, \quad (2.5) $$

with $\hat{\mathbf{a}}_\alpha^2 = \hat{\mathbf{v}}_{t,\alpha}^2(\bar{\xi}^1, \bar{\xi}^2) + \hat{\mathbf{x}}_{t,\alpha\beta}^2(\bar{\xi}^1, \bar{\xi}^2)\dot{\bar{\xi}}^\beta$, we obtain $\dot{\bar{\xi}}^\beta$ from the following system of equations:

$$ \bar{H}_{\alpha\beta} \dot{\bar{\xi}}^{\beta} = \bar{R}_{\alpha}, \quad (2.6) $$

with

$$ \begin{aligned} \bar{H}_{\alpha\beta} &= [\bar{a}_{\alpha\beta} - g_N \hat{\mathbf{x}}_{t,\alpha\beta}^2(\bar{\xi}^1, \bar{\xi}^2) \cdot \bar{\mathbf{n}}^2], \\ \bar{R}_{\alpha} &= [\mathbf{v}_t^1 - \hat{\mathbf{v}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \bar{\mathbf{a}}_{\alpha}^2 + [\mathbf{x}_t^1 - \hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \hat{\mathbf{v}}_{t,\alpha}^2(\bar{\xi}^1, \bar{\xi}^2) \\ &= [\mathbf{v}_t^1 - \hat{\mathbf{v}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \bar{\mathbf{a}}_{\alpha}^2 - g_N \bar{\mathbf{n}}^2 \cdot \hat{\mathbf{v}}_{t,\alpha}^2(\bar{\xi}^1, \bar{\xi}^2). \end{aligned} \quad (2.7) $$

Note that $[\mathbf{x}_t^1 - \hat{\mathbf{x}}_t^2(\bar{\xi}^1, \bar{\xi}^2)]$ has been written with Eqs. (2.1) and (2.2) as $-g_N \bar{\mathbf{n}}^2$. A well known result from differential geometry introduces for $\hat{\mathbf{x}}_{t,\alpha\beta}^2(\bar{\xi}^1, \bar{\xi}^2) \cdot \bar{\mathbf{n}}^2$ the curvature tensor $\bar{b}_{\alpha\beta}$. Thus, we can rewrite $\bar{H}_{\alpha\beta} = [\bar{a}_{\alpha\beta} - g_N \bar{b}_{\alpha\beta}]$.

We define as the second important kinematical function (with a view to problems with friction) a tangential relative velocity function on the current slave surface $\varphi_t^1(\Gamma_c^1)$ by setting

$$ \mathcal{L}_v g_T := \dot{\bar{\xi}}^\alpha \bar{a}_\alpha. \quad (2.8) $$

Eq. (2.8) determines per definition the evolution of the tangential slip $g_T$ which enters as a local kinematical variable the constitutive function for the contact tangential stress (see Section 3.2). The rate $\dot{\bar{\xi}}^\alpha$ in Eq. (2.8) at the solution point $(\bar{\xi}^1, \bar{\xi}^2)$ has already been computed in Eq. (2.6).
REMARK 1. Note that the second term on the right-hand side of Eq. (2.7) depends on the penetration $g_N$. Thus, in the case of a strong enforcement of the non-penetration condition ($g_N = 0$) with Lagrangian multipliers, this term vanishes. Then the evolution $\mathcal{L}_v$ in Eq. (2.8) is given by the projection of the spatial velocities $v_i^1$ and $\hat{v}_i^2(\xi)$ evaluated at the slave and master points onto the tangential direction of the master surface at the master point.

REMARK 2. If we have a flat contact surface the curvature tensor $\bar{b}_{\alpha\beta}$ is zero.

For a penetration $g_N > 0$, we have to take into account the second term in Eq. (2.8) and the scaling factors $\bar{H}_{\alpha\beta}$, both consequences of the time dependence of $\bar{a}_\alpha^2$.

### 2.3. Frame invariance

Consider the vectors $g_N := -g_N \bar{n}^2$ and $g_T := \dot{\xi}^\alpha \bar{a}_\alpha$ normal and tangential to the master surface at the master point. Clearly these vectors transform objectively under rigid body motions superposed onto the spatial configuration $(x_i^2)^+ = Q(t)x_i^2 + c(t) \forall Q(t) \in SO(3)$ where SO(3) is the proper orthogonal group. Example: $(g_N)^+ = Qg_N$ is a consequence of $(\bar{n}^2)^+ = Q\bar{n}^2$. Observe furthermore that the (a priori objective) Lie derivative of the tangential vector $g_T$ has the representation

$$ \mathcal{L}_v g_T = F_i^2 \left\{ \frac{d}{dt} [F_i^{2-1}(g_T)] \right\} = \dot{\xi}^\alpha \bar{a}_\alpha $$

based on the deformation gradient $F_i^2$ of the master surface defined above. Thus, Eq. (2.3) represents an evolution equation for the objective rate $\mathcal{L}_v g_T$ of the tangential vector introduced above.

## 3. Constitutive modelling on the contact surface

When the micromechanical behaviour of the contact area is studied, different phenomena have to be considered for the mechanical as well as for the thermal interface description. In this paper, we use constitutive models which have been derived based on micromechanical observations of physical contact surfaces. These models can be related to formulations relative to mathematical contact surfaces by an averaging process as symbolically indicated in Fig. 3. The goal of this section is to formulate local constitutive equations for the pressure, tangential stress and normal (outward) heat flux on the slave surface at the slave point $x^1$ relative to the bases $\{\bar{a}_1^2, \bar{a}_2^2, \bar{n}^2\}$ acting on body $\mathcal{B}^1$; see general structure in Eq. (4.7).

<figure><img src="image_4.png" /><figcaption>Fig. 3. Averaging of micromechanical contact relations.</figcaption></figure>
### 3.1. Contact normal stress

It is well known that the contact pressure is related to the approach of the physical surfaces which come into contact, i.e. the penetration of the mathematical surfaces results from the deformation of the micro-asperities (see Fig. 3). Based on experimental investigations Kragelsky et al. [4] formulated the following non-linear elastic constitutive equation for the contact pressure:

$$p_N = c_N(g_N)^m \qquad (3.1)$$

in terms of the penetration $g_N$ defined in Eq. (2.1). Here $c_N$ and $m$ are material parameters which have to be determined by experiment.

### 3.2. Contact tangential stress and plastic tangential slip

In this paper we restrict our consideration to classical Coulomb-type friction problems. Coulomb's constitutive equation can formally be formulated in the framework of elastoplasticity. This has been considered by several authors (e.g. [22, 23]) and within a finite element formulation by Wriggers [21]. The key idea of this approach is a split of the tangential slip component $g_T$ into an elastic part $g_T^e$ (microdisplacement describing the stick behaviour) and a plastic part $g_T^p$ (frictional slip), (see Eq. (3.2b)). We assume a linear elastic constitutive equation for the tangential contact stress component

$$t_T - c_T g_T^e \quad \text{with} \quad g_T^e := g_T - g_T^p, \qquad (3.2a,b)$$

where $c_T$ is a material parameter. The plastic tangential slip $g_T^p$ is governed by a constitutive evolution equation which can be formally derived by using standard concepts of elastoplasticity theory. Let

$$\mathcal{D}^p := t_T \cdot \mathcal{L}_v g_T^p \ge 0 \qquad (3.3)$$

be the dissipation due to the plastic slip. Now consider an elastic domain $E_t := \{t_T \in \mathbb{R}^2 | \hat{f}_p(t_T) \le 0\}$ in the space of the contact tangential stress. Here

$$\hat{f}_p(t_T) = \|t_T\| - \mu p_N \le 0 \qquad (3.4)$$

is the plastic slip criterion function for given contact pressure $p_N$ with material parameter $\mu$. Holding $p_N$ fixed, one obtains from the so-called maximum dissipation principle, $(t_T - t_T^*) \cdot \mathcal{L}_v g_T^p \ge 0 \quad \forall t_T^* \in E_t$, the constitutive evolution equation for the plastic slip

$$\mathcal{L}_v g_T^p = \lambda \frac{\partial \hat{f}_p(t_T)}{\partial t_T} = \lambda n_T \quad \text{with} \quad n_T = \frac{t_T}{\|t_T\|} \qquad (3.5)$$

(normality rule for fixed contact pressure) along with the loading-unloading conditions in Kuhn-Tucker form

$$\lambda \ge 0, \quad \hat{f}_p(t_T) \le 0, \quad \lambda \hat{f}_p(t_T) = 0, \qquad (3.6)$$

which determine the plastic parameter $\lambda$.

### 3.3. Contact heat flux

Different techniques for the computation of the thermal contact resistance, taking into account the dependence on various parameters have been proposed (see e.g. [1]). We assume the following structure for the constitutive equation for the heat flux $q_N$ (outflux through slave surface):

$$q_N = \hat{q}_N(\Theta^1, \Theta^2, p_N) = \hat{h}(\Theta^1, \Theta^2, p_N)(\Theta^1 - \Theta^2). \qquad (3.7)$$

The contact heat transfer coefficient $h = \hat{h}(\Theta^1, \Theta^2, p_N)$ depends on the current temperatures $\Theta^1$ and $\Theta^2$ at the slave and master points, $x^1$ and $\hat{x}_i^2(\xi)$, respectively, as well as on the contact pressure $p_N$ defined in Eq. (3.1). The resistance against heat transfer is mainly due to the low percentage of physical surface area which is really in contact. The presence of a reduced set of spots surrounded by microcavities
characterizes the contact zone, hence heat exchange is possible by heat conduction through the spots (S), heat conduction through the gas (G) contained in the cavities and radiation (R) between microcavity surfaces. These mechanisms act in parallel which induces the additive decomposition

$$ \hat{h} = \hat{h}_S + \hat{h}_G + \hat{h}_R \quad (3.8) $$

of the contact heat transfer function. Based on statistical micromechanical models, different relations for the gas and spot conductance have been developed (see e.g. [2]).

To obtain the relation representing the heat transfer associated with the spots, the behaviour of one spot is analyzed in detail [24]. Here, a heat flux tube having a narrowing in the contact zone is considered around each spot. Then the mean value of the heat transfer through the spots is obtained from a statistical model. The additional effect of taking into account the hardness variation with the mean planes approach has been suggested in [2] leading to the following interface conductivity for the spots:

$$ h_s = \frac{1.25\hat{k}\bar{m}}{\sigma} \left[ \frac{p_N}{c_1} \left( 1.62 \frac{10^6 \sigma}{\bar{m}} \right)^{-c_2} \right]^{0.95/(1+0.071c_2)} \quad (3.9) $$

Here $c_1$, $c_2$ describe the harness variation and $\hat{k}$ is the mean thermal conductivity depending on the conductivities of the two bodies contact. $\bar{m}$ is the mean absolute asperity slope and $\sigma$ denotes the surface roughness. All values have to be deduced from experimental data.

The heat transfer due to gas or liquid within the microcavities is mainly governed by conduction. This fact results from the small height of the cavities which do not allow for convective flow. Based on this assumption, Yovanovich [25] derived a relationship, which also takes into account the change of cavity height by the contact pressure $p_N$. With the current gas temperature $\theta_G$, a constitutive constant $c_{PC}$ for the gas and the Vickers hardness $H_e$, the gas conductivity $k_g$ and the surface roughness $\sigma$, the heat transfer through the gas yields

$$ h_G = \frac{k_g}{1.36\sigma\sqrt{-\ln\left(5.59\frac{p_N}{H_e}\right) + c_{PC}\theta_G}} \quad (3.10) $$

In this paper, we use a simplified relation which neglects radiation effects and gas conductance but takes into account pressure dependency of the contact conductance

$$ h_c(p_N) = h_{c0} \left[ \frac{p_N}{H_e} \right]^\epsilon \quad \text{and} \quad h_G = h_R = 0, \quad (3.11) $$

based on the material parameters $h_{c0}$, $H_e$, $\epsilon$. This model is valid for high pressures. A more sophisticated constitutive model is implemented in [14] for the frictionless case.

## 4. Initial boundary value problem

In the previous sections we have discussed the contact geometry and the constitutive equations associated with contact interface. Let us now formulate the initial boundary value problem for non-linear thermoelasticity combined with thermomechanical (frictional) contact. In order to obtain a compact structure, we introduce locally at the material point $X^\alpha \in \mathcal{B}^\alpha$ and time $t \in \mathbb{R}_+$ a vector of primary variables $\mathcal{S}^\alpha(X^\alpha, t) := \{\varphi^\alpha(X^\alpha, t), V^\alpha(X^\alpha, t), \Theta^\alpha(X^\alpha, t)\}$ (which contains the configuration $\varphi^\alpha$, the material velocity $V^\alpha$ and the material temperature $\Theta^\alpha$) and locally at points $X^\alpha \in \Gamma_c^\alpha$ of the contact surfaces a vector of history variables

$$ \mathcal{S}(X^\alpha, t) := \{g_T^P(X^\alpha, t)\} $$

(which contains the plastic tangential slip $g_T^P$). Internal variables in the domains $\mathcal{B}^\alpha$ do not appear since we restrict our presentation to thermoelastic constitutive response. The thermomechanical initial boundary value problem is governed by the local field equations in the domains $\mathcal{B}^\alpha$,
$$ \frac{\partial}{\partial t} \boldsymbol{\varphi}^{\alpha} = \boldsymbol{V}^{\alpha} + \mathbf{0}, \qquad (4.1a) $$

$$ \frac{\partial}{\partial t} \boldsymbol{V}^{\alpha} = \frac{1}{\rho_{0}^{\alpha}} \mathrm{DIV} \boldsymbol{P}^{\alpha} + \mathbf{0}, \qquad (4.1b) $$

$$ \frac{\partial}{\partial t} \boldsymbol{\Theta}^{\alpha} = 0 - \frac{1}{c^{\alpha}} \mathrm{DIV} \boldsymbol{Q}^{\alpha} + \frac{1}{c^{\alpha}} \boldsymbol{S}^{\alpha}, \qquad (4.1c) $$

which represent the definition of the material velocity, the balance of linear momentum and the balance of internal energy in form of the temperature evolution equation. $\boldsymbol{P}^\alpha$ is the first Piola-Kirchhoff stress tensor, $\boldsymbol{Q}^\alpha$ denotes the material heat flux vector and $\boldsymbol{S}^\alpha$ is a heat source which describes the Gough-Joule coupling effect in the framework of thermoelasticity. The field equations (4.1) form a coupled first order evolution system for the primary variable vector $\mathcal{S}^\alpha$ introduced above, i.e.

$$ \frac{\partial}{\partial t} \mathcal{S}^{\alpha} = \hat{\boldsymbol{A}}^{\alpha}(\mathcal{S}^{\alpha}), \qquad (4.2) $$

where the non-linear evolution operator $\hat{\boldsymbol{A}}$ represents the right-hand side of Eq. (4.1) by taking into account the dependency of the thermoelastic constitutive equations in the domains $\mathcal{B}^\alpha$

$$ \boldsymbol{P}^{\alpha} = \hat{\boldsymbol{P}}^{\alpha}(\mathcal{S}^{\alpha}), \quad \boldsymbol{Q}^{\alpha} = \hat{\boldsymbol{Q}}^{\alpha}(\mathcal{S}^{\alpha}), \quad \boldsymbol{S}^{\alpha} = \hat{\boldsymbol{S}}^{\alpha}(\mathcal{S}^{\alpha}), \qquad (4.3) $$

on $\mathcal{S}^\alpha$. Next the initial conditions for the primary and history variables for a time interval $[t_n, t_{n+1}] \in \mathbb{R}_+$ of interest are given by

$$ \mathcal{S}^{\alpha}(\boldsymbol{X}^{\alpha}, t_n) = \mathcal{S}_{n}^{\alpha}(\boldsymbol{X}^{\alpha}) \quad \text{and} \quad \mathcal{I}^{\alpha}(\boldsymbol{X}^{\alpha}, t_n) = \mathcal{I}_{n}^{\alpha}(\boldsymbol{X}^{\alpha}). \qquad (4.4) $$

We consider boundary conditions for the deformation and the temperature on $\Gamma_{\varphi}^{\alpha} \subset \partial\mathcal{B}^{\alpha}$ and $\Gamma_{\theta}^{\alpha} \subset \partial\mathcal{B}^{\alpha}$

$$ \begin{aligned} \boldsymbol{\varphi}^{\alpha}(\boldsymbol{X}^{\alpha}, t) &= \bar{\boldsymbol{\varphi}}^{\alpha}(\boldsymbol{X}^{\alpha}, t) && \text{on } \Gamma_{\varphi}^{\alpha}, \\ \boldsymbol{\Theta}^{\alpha}(\boldsymbol{X}^{\alpha}, t) &= \bar{\boldsymbol{\Theta}}^{\alpha}(\boldsymbol{X}^{\alpha}, t) && \text{on } \Gamma_{\theta}^{\alpha}. \end{aligned} \qquad (4.5) $$

as well as for the traction vector and the heat flux on $\Gamma_{t}^{\alpha} \subset \partial\mathcal{B}^{\alpha}$ and $\Gamma_{q}^{\alpha} \subset \partial\mathcal{B}^{\alpha}$

$$ \begin{aligned} t^{\alpha}(\boldsymbol{\varphi}_{t}^{\alpha}(\boldsymbol{X}^{\alpha}), t) &= t^{\alpha}(\boldsymbol{\varphi}_{t}^{\alpha}(\boldsymbol{X}^{\alpha}), t) && \text{on } \boldsymbol{\varphi}_{t}^{\alpha}(\Gamma_{t}^{\alpha}), \\ q^{\alpha}(\boldsymbol{\varphi}_{t}^{\alpha}(\boldsymbol{X}^{\alpha}), \boldsymbol{\Theta}_{t}^{\alpha}(\boldsymbol{X}^{\alpha}), t) &= q^{\alpha}(\boldsymbol{\varphi}_{t}^{\alpha}(\boldsymbol{X}^{\alpha}), \boldsymbol{\Theta}_{t}^{\alpha}(\boldsymbol{X}^{\alpha}), t) && \text{on } \boldsymbol{\varphi}_{t}^{\alpha}(\Gamma_{q}^{\alpha}). \end{aligned} \qquad (4.6) $$

The initial boundary value problem is completed by the thermomechanical constitutive equations on the current slave contact surface $\boldsymbol{\varphi}_t^1(\Gamma_c^1)$,

$$ \begin{aligned} p_N &= \hat{p}_N[\mathcal{S}^{\alpha, \beta}, \mathcal{I}, \lambda], & t_T &= \hat{t}_T[\mathcal{S}^{\alpha, \beta}, \mathcal{I}, \lambda], \\ q_N &= \hat{q}_N[\mathcal{S}^{\alpha, \beta}, \mathcal{I}, \lambda], & \mathcal{D}^P &= \hat{\mathcal{D}}^P[\mathcal{S}^{\alpha, \beta}, \mathcal{I}, \lambda], \end{aligned} \qquad (4.7) $$

for the contact pressure, the contact tangential stress, the contact heat flux, the frictional dissipation and the evolution equation for the plastic slip with structure,

$$ \frac{\partial}{\partial t} \mathcal{I} = \lambda \hat{\mathcal{E}}[\mathcal{S}^{\alpha, \beta}, \mathcal{I}], \qquad (4.8) $$

which is constrained by the Kuhn-Tucker loading-unloading conditions

$$ \lambda \ge 0, \quad \hat{f}_p[\mathcal{S}^{\alpha, \beta}, \mathcal{I}] \le 0, \quad \lambda \hat{f}_p[\mathcal{S}^{\alpha, \beta}, \mathcal{I}] = 0. \qquad (4.9) $$

The superscript $(\alpha, \beta)$ indicates the dependency of the constitutive functions on both contacting surfaces. Eqs. (4.7) and (4.8) are of central interest within this paper and have already been specified for model problems in Section 3; cf. Eqs. (3.1)-(3.7).
## 5. Solution algorithms

The main solution strategies for coupled problems are monolithic schemes where the differential equations for the different variables are all solved together, or staggered schemes where the different variables are solved separately (see [26] for an overview). Global solution algorithms for coupled thermomechanical analysis have been formulated by Argyris and Doltsinis [16], Miehe [17], Doltsinis [18], Simo and Miehe [19], among others. In what follows, we adopt the staggered scheme proposed in [19] based on an operator split of the global thermomechanical evolution operator discussed in Section 4. This strategy yields an algorithmic decoupling of the thermomechanical equations within a time step on the basis of (in the frictionless case) symmetric subproblems.

A further objective of this section is to formulate the local update algorithm for the contact tangential stress and the frictional dissipation. Here we extend the return mapping algorithm proposed e.g. in [21] for isothermal frictional contact to the non-isothermal case.

### 5.1. Global thermomechanical operator split algorithm

The central idea is a split of the evolution operator $\hat{A}^\alpha$ in Eq. (4.2) into its natural mechanical and thermal parts as indicated in Eq. (4.1): $\hat{A}^\alpha(\mathcal{S}^\alpha) = \hat{A}_M^\alpha(\mathcal{S}^\alpha) + \hat{A}_T^\alpha(\mathcal{S}^\alpha)$. This split defines within the time step $\Delta t_{n+1} := t_{n+1} - t_n$ two subproblems

$$ (M): \quad \frac{\partial}{\partial t} \bar{\mathcal{F}}^{\alpha} = \hat{A}_{M}^{\alpha}(\bar{\mathcal{F}}^{\alpha}) \quad \text{and} \quad (T): \quad \frac{\partial}{\partial t} \mathcal{F}^{\alpha} = \hat{A}_{T}^{\alpha}(\mathcal{F}^{\alpha}), \qquad (5.1) $$

a purely mechanical subproblem (M) at frozen temperature along with the mechanical initial and boundary conditions in Eqs. (4.4)-(4.7) followed by a purely thermal subproblem (T) at frozen configuration along with the thermal initial and boundary conditions in Eqs. (4.4)-(4.7). Both subproblems are constrained by the evolution of the plastic variables Eqs. (4.8)-(4.9). We consider the quasistatic problem and integrate both phases of Eq. (5.1) with the backward Euler algorithm. Then the algorithmic counterpart $\text{ALGO}_M$ for the mechanical subproblem (M) takes the form (integration of Eq. (4.1a,b))

$$ \text{DIV } \hat{P}^{\alpha}(\varphi_{n+1}^{\alpha}; \Theta_{n}^{\alpha}) = 0 \quad \text{and} \quad V_{n+1}^{\alpha} = \frac{1}{\Delta t_{n+1}} (\varphi_{n+1}^{\alpha} - \varphi_{n}^{\alpha}) \qquad (5.2a,b) $$

in $\mathcal{B}^\alpha$ at frozen thermal primary variable $\Theta_n^\alpha$. The thermal subproblem (T) obtains the algorithmic form $\text{ALGO}_T$ (integration of Eq. (4.1c))

$$ \frac{c^{\alpha}}{\Delta t_{n+1}} (\Theta_{n+1}^{\alpha} - \Theta_{n}^{\alpha}) = -\text{DIV } \hat{Q}^{\alpha}(\Theta_{n+1}^{\alpha}; \varphi_{n+1}^{\alpha}) + \hat{S}^{\alpha}(\Theta_{n+1}^{\alpha}; \varphi_{n+1}^{\alpha}, V_{n+1}^{\alpha}) \qquad (5.3) $$

in $\mathcal{B}^\alpha$ at frozen mechanical primary variables $\varphi_{n+1}^\alpha, V_{n+1}^\alpha$. Thus we solve within a time step $\Delta t_{n+1}$ first the mechanical problem (5.2) for the actual configuration field $\varphi_{n+1}^\alpha$. Next we compute the actual temperature field $\Theta_{n+1}^\alpha$ by solving the thermal problem (5.3). The overall global thermomechanical solution algorithm within a typical time step can be regarded as a composition $\text{ALGO}_{\text{TM}} = \text{ALGO}_T \circ \text{ALGO}_M$ of the two subalgorithms $\text{ALGO}_T$ and $\text{ALGO}_M$. Thus the operator split algorithm results in an algorithmic decoupling of the coupled thermomechanical equations within the time interval. The proposed split (5.1) yields an identical coupling algorithm to that suggested in [16], characterized by an isothermal deformation predictor followed by a heat conduction corrector. Within their work, this algorithm has been interpreted as a one-pass Gauss-Seidel scheme in terms of the variables $\varphi$ and $\Theta$ (see also [18]).

In view of the finite element treatment, we construct the weak forms of the time discrete field equations (5.2a) and (5.3) by standard Galerkin procedures. The weak form of the mechanical subproblem $\text{ALGO}_M$ takes the form
$$
\begin{equation}
\begin{split}
\hat{G}_M(\varphi_{n+1}^\alpha, \dot{\boldsymbol{u}}^\alpha; \Theta_n^\alpha) := & - \int_{\varphi_{n+1}^1(\Gamma_c^1)} [\delta g_N(-p_N) + \delta g_T \cdot t_T] \, da^1 \\
& + \sum_{\alpha=1}^2 \left\{ \int_{\mathcal{B}^\alpha} \nabla \boldsymbol{u}^\alpha \cdot \boldsymbol{\tau}^\alpha \, dV - \int_{\varphi_{n+1}^\alpha(\Gamma_t^\alpha)} \dot{\boldsymbol{u}}^\alpha \cdot \bar{\boldsymbol{t}}^\alpha \, da \right\} = 0,
\end{split}
\tag{5.4}
\end{equation}
$$

where we have dropped for convenience the subscript $n+1$ on the right hand side. $\dot{\boldsymbol{u}}^\alpha$ denotes a
mechanical test function field defined on the current configuration (virtual displacements). The first
term in Eq. (5.4) describes the contributions due to the contact with the relative (Lie-type) variations
$-\delta g_N = [\dot{\boldsymbol{u}}_t^1 - \dot{\boldsymbol{u}}_t^2(\bar{\xi}^1, \bar{\xi}^2)] \cdot \bar{\boldsymbol{n}}^2$ analogous to Eq. (2.4) and $\mathcal{L}_v g_T := \bar{\xi}^\alpha \bar{a}_\alpha$ analogous to Eqs. (2.6)-(2.8).
The second term contains the standard domain contributions with the Kirchhoff stress tensor $\boldsymbol{\tau}^\alpha$. The
weak form of the thermal subproblem ALGO$_T$ takes the form

$$
\begin{equation}
\begin{split}
\hat{G}_T(\Theta_{n+1}^\alpha, \vartheta^\alpha; \varphi_{n+1}^\alpha) := & - \int_{\varphi_{n+1}^1(\Gamma_c^1)} [\delta\vartheta_q(-q_N) + \delta\vartheta_{\mathcal{D}}^\alpha \mathcal{D}^P] da^1 \\
& + \sum_{\alpha=1}^2 \left\{ \int_{\mathcal{B}^\alpha} \left[ \nabla\vartheta^\alpha \cdot (-q^\alpha) - \vartheta^\alpha \left( S^\alpha - \frac{c^\alpha}{\Delta t} (\Theta^\alpha - \Theta_n^\alpha) \right) \right] dV + \int_{\varphi_{n+1}^\alpha(\Gamma_q^\alpha)} \vartheta^\alpha (-\bar{q}^\alpha) da \right\} = 0.
\end{split}
\tag{5.5}
\end{equation}
$$

Here $\vartheta^\alpha$ denotes a thermal test function field defined on the current configuration (virtual temperatures). The first term in Eq. (5.5) describes the contributions due to the thermal contact conductance and frictional heating with the variations $\delta\vartheta_q = \vartheta^1 - \vartheta^2$ and $\delta\vartheta_{\mathcal{D}} = \frac{1}{2}(\vartheta^1 + \vartheta^2)$ on $\varphi_t^1(\Gamma_c^1)$. The second term characterizes the standard domain contributions with Kirchhoff heat flux vector $q^\alpha$. The non-linear problems (5.4) and (5.5) have to be solved by an iterative solver, e.g. Newton's algorithm.

5.2. Local update algorithm for contact tangential stress and dissipation

The algorithmic update of the tangential stress $t_{Tn+1}$ and dissipation $\mathcal{D}_{n+1}^p$ in Eqs. (5.4) and (5.5) is
performed by the return algorithm based on an objective (backward Euler) integration of the evolution
equation (3.5) for the plastic slip (see e.g. [9, 21, 27, 28]). The results can be summarized as follows.
Integration of Eq. (2.3) gives the increment of the total slip within the time step $\Delta t_{n+1}$

$$
\Delta g_{Tn+1} = (\bar{\xi}_{n+1} - \bar{\xi}_n)^\beta \bar{a}_{\beta n+1}, \quad (5.6)
$$

which has to be decomposed into an elastic and a plastic part. Thus in the case of contact, i.e. for
$g_{Nn+1} > 0$, we compute the elastic trial state

$$
\begin{align}
t_{Tn+1}^{\text{trial}} &:= t_{Tn} + c_T \Delta g_{Tn+1} \nonumber \\
f_{pn+1}^{\text{trial}} &:= \|t_{Tn+1}^{\text{trial}}\| - \mu p_{Nn+1} \tag{5.7}
\end{align}
$$

and perform the return algorithm

$$
t_{Tn+1} = \begin{cases} t_{Tn+1}^{\text{trial}}, & \text{for } f_{pn+1}^{\text{trial}} \le 0, \\ t_{Tn+1}^{\text{trial}} (\mu p_{Nn+1} / \|t_{Tn+1}^{\text{trial}}\|), & \text{otherwise,} \end{cases} \quad (5.8)
$$

where we have chosen $t_{Tn}$ as history variable. The dissipation due to the plastic slip is given by

$$
\mathcal{D}_{n+1}^{p} = \begin{cases} 0, & \text{for } f_{pn+1}^{\text{trial}} \le 0, \\ \frac{t_{Tn+1}^{\text{trial}} (t_{Tn+1}^{\text{trial}} - t_{Tn+1})}{c_T \Delta t_{n+1}}, & \text{otherwise.} \end{cases} \quad (5.9)
$$
## 6. Finite element discretization

The discretization of the domain contributions in Eqs. (5.4) and (5.5) is not the objective of this work. Within this context, we refer to the finite element implementations of coupled finite thermoelastic initial boundary value problems documented in [17] and coupled finite strains thermoplastic problems in [19, 29] (see also references therein).

In this paper, we focus on the contribution due to the presence of the contact constraints, i.e. the discretization of the integrals in Eqs. (5.4) and (5.5) over the moving master surface. The basis for our subsequent finite element treatment provides a two-dimensional contact element formulation which has been proposed in [7] for the frictionless case and extended in [9] to finite deformations including frictional slip. In what follows, we extend this element formulation to the thermomechanically coupled case under consideration involving contact heat exchange and frictional heating. Within the context of formulating the contact element geometry, it is shown that both strategies, (i) starting directly from the discrete formulation and (ii) starting from the continuous formulation outlined in Section 2, result in the same discretized formulation of the local contact geometry.

Consider the contact element in Fig. 4 and assume the case of contact, i.e. the discrete slave point (s) penetrates into the master segment (1)–(2). $l := \|\mathbf{x}_2^2 - \mathbf{x}_1^2\|$ denotes the current length of the master segment and $\{a_1^2 := (\mathbf{x}_2^2 - \mathbf{x}_1^2)/l; n^2 := E_3 \times a_1^2\}$ a local orthonormal base attached to the master segment. Here we have defined the base vector $\mathbf{a}^1$ as a unit vector. Its connection to the tangent vector of the continuum case in (2.3) is given by $\bar{\mathbf{a}}_1^2 = [(1 - \bar{\xi}^1)\mathbf{x}_1^2 + \bar{\xi}^1 \mathbf{x}_2^2]_1 = l\mathbf{a}_1^2$.

$\bar{\xi}^1$ and $g_N$ are given by the solution of the minimal distance problem, i.e. by the projection of the slave node (s) onto the master segment (1)–(2)

$$ \bar{\xi}^1 = \frac{1}{l} (\mathbf{x}_s^1 - \mathbf{x}_1^2) \cdot \mathbf{a}_1^2 \quad \text{and} \quad g_N = \|\mathbf{x}_s^1 - (1 - \bar{\xi}^1)\mathbf{x}_1^2 - \bar{\xi}^1 \mathbf{x}_2^2\|. \quad (6.1) $$

From these equations and the local formulation (2.4), we compute directly by assuming the linear interpolation $\hat{\mathbf{x}}_i^2(\xi) = \mathbf{x}_i^2 + \xi^1(\mathbf{x}_2^2 - \mathbf{x}_1^2)$ on the straight master segment (1)–(2).

$$ -\dot{g}_N = [\hat{\mathbf{x}}_s^1 - (1 - \bar{\xi}^1)\hat{\mathbf{x}}_1^2 - \bar{\xi}^1 \hat{\mathbf{x}}_2^2] \cdot \mathbf{n}^2. \quad (6.2a) $$

The local equation (2.7) yields the expression for $\bar{\xi}^1$. With the interpolation $\hat{v}_i^2(\xi^1) = \hat{\mathbf{x}}_i^2 + \xi^1(\hat{\mathbf{x}}_2^2 - \hat{\mathbf{x}}_1^2)$ on the straight master segment (1)–(2), we specialize

$$ \begin{aligned} \bar{H}_{\alpha\beta} &= (a_{\alpha\beta} - g_N b_{\alpha\beta}) \Rightarrow \bar{H}_{11} = a_{11} = l^2, \\ \bar{R}_1 &= [\mathbf{v}_i^1 - \hat{v}_i^2(\bar{\xi}, \bar{\xi}^2)] \cdot \bar{\mathbf{a}}_1^2 - g_N \bar{\mathbf{n}}^2 \cdot \hat{v}_{i,1}^2(\bar{\xi}^1), \end{aligned} $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Approach in the contact area.</figcaption></figure>
which leads to

$$ \dot{\xi}^1 l = [\dot{\mathbf{x}}_s^1 - (1 - \bar{\xi}^1) \dot{\mathbf{x}}_1^2 - \bar{\xi}^1 \dot{\mathbf{x}}_2^2] \cdot \mathbf{a}_1^2 - \frac{g_N}{l} [\dot{\mathbf{x}}_2^2 - \dot{\mathbf{x}}_1^2] \cdot \mathbf{n}^2. \qquad (6.2b) $$

Eqs. (6.1) and (6.2) characterize the main kinematical relations of the contact element in Fig. 4.

In what follows, we formulate the constitutive relationships for the normal force $P_N$, tangential force $T_T$, the discrete heat flux $Q_N$ and the discrete dissipation $\mathcal{D}^p$ at the discrete contact point(s) of the contact element under consideration. Assume that $\bar{\xi}_{n+1}$ and $g_{Nn+1}$ are known from the exploitation of Eq. (6.1). Then analogous to Eq. (3.1), we compute the normal force

$$ P_{Nn+1} = C_N (g_{Nn+1})^m. \qquad (6.3) $$

The algorithmic update of the tangential force is performed by the return algorithm as described in Section 5.2. Note that we have only one component in the two-dimensional case. Analogous to Eqs. (5.6)-(5.9), we compute the elastic trial state

$$ T_{Tn+1}^{\text{trial}} := T_{Tn} + C_T \Delta g_{Tn+1}, \qquad F_{Pn+1}^{\text{trial}} := |T_{Tn+1}^{\text{trial}}| - \mu P_{Nn+1}, \qquad (6.4) $$

perform the return algorithm

$$ T_{Tn+1} = \begin{cases} T_{Tn+1}^{\text{trial}}, & \text{for } F_{Pn+1}^{\text{trial}} \le 0, \\ T_{Tn+1}^{\text{trial}} (\mu P_{Nn+1} / |T_{Tn+1}^{\text{trial}}|), & \text{otherwise,} \end{cases} \qquad (6.5) $$

and compute the frictional dissipation

$$ \mathcal{D}_{n+1}^{p} = \begin{cases} 0, & \text{for } F_{Pn+1}^{\text{trial}} \le 0, \\ T_{Tn+1}(T_{Tn+1}^{\text{trial}} - T_{Tn+1}) / C_T \Delta t_{n+1}, & \text{otherwise.} \end{cases} \qquad (6.6) $$

Finally the discrete heat flux takes the form analogous to Eq. (3.7)

$$ Q_{Nn+1} = H_{S0} \left[ \frac{P_{Nn+1}}{H} \right]^\epsilon \left[ \Theta_{Sn+1}^1 - (1 - \bar{\xi}_{n+1}) \Theta_{1n+1}^2 - \bar{\xi}_{n+1} \Theta_{2n+1}^2 \right]. \qquad (6.7) $$

The contributions of the thermomechanical contact in the mechanical and thermal weak forms (5.4) and (5.5) take the contact element of the form

$$ G_M^c := \delta g_{Nn+1} P_{Nn+1} - \delta g_{1n+1} T_{1n+1}, \qquad (6.8a) $$

$$ G_T^c := \delta g_{Qn+1} Q_{Nn+1} - \delta g_{\mathcal{D}n+1} \mathcal{D}_{n+1}^p, \qquad (6.8b) $$

for the discrete contact point (s) with the mechanical relative (Lie-type) variations analogous to Eq. (6.2),

$$ \delta g_{Nn+1} = [\mathbf{u}_s^* - (1 - \bar{\xi}_{n+1}) \mathbf{u}_1^* - \bar{\xi}_{n+1} \mathbf{u}_2^*] \cdot \mathbf{e}_{Nn+1}^2 \qquad (6.9a) $$

$$ \delta g_{Tn+1} = [\mathbf{u}_s^* - (1 - \bar{\xi}_{n+1}) \mathbf{u}_1^* - \bar{\xi}_{n+1} \mathbf{u}_2^*] \cdot \mathbf{e}_{Tn+1}^2 - \frac{g_{Nn+1}}{l_{n+1}} [\mathbf{u}_2^* - \mathbf{u}_1^*] \cdot \mathbf{e}_{Nn+1}^2 \qquad (6.9b) $$

and thermal variations

$$ \begin{aligned} \delta g_{Qn+1} &= [\vartheta_s^1 - (1 - \bar{\xi}_{n+1}) \vartheta_1^2 - \bar{\xi}_{n+1} \vartheta_2^2], \\ \delta g_{\mathcal{D}n+1} &= [\vartheta_s^1 + (1 - \bar{\xi}_{n+1}) \vartheta_1^2 + \bar{\xi}_{n+1} \vartheta_2^2]/2. \end{aligned} \qquad (6.10) $$

These equations can now be cast into a matrix formulation. For the mechanical part Eq. (6.8a), we set for the variation (6.9a) of the penetration
$$ \delta g_{Nn+1} = \dot{u}^* \cdot N_{n+1}. \qquad (6.11) $$

With the same notation, we can express the variation (6.9b) of the tangential gap

$$ \delta g_{Tn+1} = \dot{u}^* \cdot \left( T_{n+1} - \frac{g_{Nn+1}}{l_{n+1}} N_{0n+1} \right). \qquad (6.12) $$

In Eqs. (6.11) and (6.12), the following vectors have been used:

$$ \dot{u}^* := [u_s^1 \quad u_s^2 \quad u_1^2 \quad u_1^2 \quad u_2^2 \quad u_2^2]^t, \qquad (6.13) $$

$$ N_{n+1} := \begin{Bmatrix} e_N^2 \\ -(1 - \xi^1)e_N^2 \\ -\bar{\xi}^1 e_N^2 \end{Bmatrix}_{n+1}, \quad N_{0n+1} := \begin{Bmatrix} \mathbf{0} \\ -e_N^2 \\ e_N^2 \end{Bmatrix}_{n+1}, \qquad (6.14) $$

and

$$ T_{n+1} := \begin{Bmatrix} e_T^2 \\ -(1 - \bar{\xi}^1)e_T^2 \\ -\bar{\xi}^1 e_T^2 \end{Bmatrix}_{n+1}, \quad T_{0n+1} := \begin{Bmatrix} \mathbf{0} \\ -e_T^2 \\ e_T^2 \end{Bmatrix}_{n+1}. \qquad (6.15) $$

Thus the virtual mechanical work, Eq. (6.8a), of the contact element can be written in the matrix formulation $G_M^c = \dot{u}^* \cdot R_{Mn+1}^c$ with the mechanical contact element residual,

$$ R_{Mn+1}^c := P_{Nn+1} N_{n+1} - T_{Tn+1} \left( T_{n+1} - \frac{g_{Nn+1}}{l_{n+1}} N_{0n+1} \right). \qquad (6.16) $$

In Eq. (6.16), the contact normal force $P_{Nn+1}$ follows from Eq. (6.3) whereas the tangential force $T_{Tn+1}$ is given by the return algorithm (6.5). Due to this approach, a pure displacement formulation of the contact problem is possible which is in contrast to the Lagrangian multiplier technique often used to enforce Eq. (2.1). For a global algorithmic treatment using Newton's method, we have to linearize Eq. (6.16). The associated formulation for this discretization can be found in [9].

The matrix formulation of the thermal part, Eq. (6.8b), is similar to the mechanical part. As a consequence of the global operator split algorithm discussed in Section 5 and the assumed simplified constitutive equation, Eq. (6.7) is the thermal part linear.

A formulation completely analogous to Eqs. (6.11)–(6.15) yields the matrix representation $G_T^c = \dot{\mathbf{\theta}}^* \cdot \mathbf{R}_{Tn+1}^c$ for the virtual thermal work Eq. (6.8b) with the thermal contact element residual

$$ \mathbf{R}_{Tn+1}^c := Q_{Nn+1} \mathbf{H}_{n+1} - \mathcal{D}_{n+1}^p \mathbf{H}_{0n+1}. \qquad (6.17) $$

Here we have introduced the matrices

$$ \dot{\mathbf{\theta}}^* := [\dot{\vartheta}_s^1 \quad \dot{\vartheta}_1^2 \quad \dot{\vartheta}_2^2]^t, \qquad (6.18) $$

analogous to Eq. (6.13) and

$$ \mathbf{H}_{n+1} := \begin{Bmatrix} 1 \\ -(1 - \bar{\xi}^1) \\ -\bar{\xi}^1 \end{Bmatrix}_{n+1}, \quad \mathbf{H}_{0n+1} := \frac{1}{2} \begin{Bmatrix} 1 \\ (1 - \bar{\xi}^1) \\ \bar{\xi}^1 \end{Bmatrix}_{n+1}, \qquad (6.19) $$

for the matrix formulation of the thermal variations (6.10). In Eq. (6.19) the discrete contact heat flux $Q_{Nn+1}$ follows from Eq. (6.7) whereas the frictional dissipation $\mathcal{D}_{n+1}^p$ is given by the return algorithm (6.6). The nodal contact contribution Eqs. (6.16) and (6.19) have to be added to the global system of equations.
## 7. Numerical examples

The proposed contact-model has been implemented in the finite element code FEAP and was tested by means of several examples. For the following examples, we use finite element formulations for large strain coupled thermo-elastic and thermo-plastic response, which are documented in detail in [17] and [19], respectively.

### 7.1. Heat transfer through the contact surface

The first example is a test for the developed and implemented pressure dependent thermomechanical constitutive equation. The system and its finite element discretization is depicted in Fig. 5. A temperature of 100 K and a pressure $p$ is applied to the upper surface of $\mathcal{B}^2$ whereas the bottom of $\mathcal{B}^1$ is fixed and has a temperature of 0 K, the vertical boundaries are thermally isolated. Body $\mathcal{B}^1$ is assumed to be rigid but can conduct heat. Body $\mathcal{B}^2$ is thermoelastic. The thermoelastic constitutive parameters for aluminium are given in Table 1.

Fig. 5 shows the comparison with the analytical solution which is given for this one-dimensional steady state example by

<figure><img src="image_6.png" /><figcaption>Fig. 5. System and discretization. Contact pressure versus contact temperature.</figcaption></figure>
<table><caption>Table 1
Material parameters for examples in Sections 7.1 and 7.2</caption><tr><td>Bulk modulus</td><td>κ</td><td>58333 N/mm<sup>2</sup></td></tr><tr><td>Shear modulus</td><td>G</td><td>26926 N/mm<sup>2</sup></td></tr><tr><td>Density</td><td>ρ</td><td>2.7 × 10<sup>-9</sup> N s<sup>2</sup>/mm<sup>4</sup></td></tr><tr><td>Expansion coefficient</td><td>α</td><td>23.86 × 10<sup>-6</sup> K<sup>-1</sup></td></tr><tr><td>Conductivity</td><td>k</td><td>150 N/s K</td></tr><tr><td>Capacity</td><td>c</td><td>0.9 × 10<sup>9</sup> mm<sup>2</sup>/s<sup>2</sup> K</td></tr><tr><td>Resistance coefficient</td><td>h<sub>so</sub></td><td>150 N/s K</td></tr><tr><td>Vickers hardness</td><td>H</td><td>932 N/mm<sup>2</sup></td></tr><tr><td>Exponent</td><td>ε</td><td>0.95</td></tr></table>

$$ \vartheta_c^1 = \frac{(1 + \eta)\vartheta_A + \eta\vartheta_B}{(1 + 2\eta)}, \quad \vartheta_c^2 = \frac{(1 + \eta)\vartheta_B + \eta\vartheta_A}{(1 + 2\eta)} \quad \text{with } \eta = h_c \frac{1}{k}. $$

These results clearly show the pressure dependency of the constitutive law for the heat transfer in the contact area. They depict also that the finite element approximation is in very good accordance with the analytical solution.

## 7.2. Frictional heating of a block on a rigid surface

The mechanism of heating due to friction is studied in this example. For this purpose, we consider an elastic block $B^2$ which is slid over a rigid block $B^1$. Both bodies are heat conductors. The upper body moves within $3.75 \times 10^{-3}$ s from the left to the right end of the lower block, the displacements are prescribed at the top of $B^2$. During this process a pressure of $p = 10 \text{ N/mm}^2$ is applied on block $B^2$. The system and its finite element discretization are shown in Fig. 6. In order to check the numerical result, we consider idealized thermally isolated surfaces of both bodies. Again we assume that both bodies consist of aluminium with the material parameters given in Table 1. The frictional coefficient is assumed to be $\mu = 0.2$.

A total of 100 time steps have been used to compute the evolution of the temperature during the process. The temperature distributions at time $t = 0.75 \times 10^{-3}$ s and $t = 2.25 \times 10^{-3}$ s are shown in Fig. 7. After that time, the movement of the upper body was stopped.

Another 50 time steps with increasing step length were applied to compute the temperature evolution up to the final homogeneous steady state at the process time of 10 s. After that time, the temperature is balanced and thus equal at all elements of the bodies $B^a$. Since in this example, we wanted to evaluate

<figure><img src="image_7.png" /><figcaption>Fig. 6. System of frictional heating problem.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Temperature distribution at different time steps.</figcaption></figure>

only the frictional heating mechanism of idealized thermally insulated surfaces of both bodies were considered within each stage of the process. This made it possible to compare the finite element result with a global energy check. We compute the dissipative mechanical work: $W_{\text{diss}} = \max uF_{\text{T}} = \max$
<table><caption>Table 2</caption><thead><tr><th></th><th>Finite element result<br>(K)</th><th>Global energy check<br>(K)</th></tr></thead><tbody><tr><td>ϑ<sup>1</sup></td><td>0.317</td><td>0.309</td></tr><tr><td>ϑ<sup>2</sup></td><td>1.267</td><td>1.235</td></tr></tbody></table>

<table><caption>Table 3
Material parameters for the example in Section 7.3</caption><tr><td>Bulk modulus</td><td>κ</td><td>58333 N/mm<sup>2</sup></td></tr><tr><td>Shear modulus</td><td>μ</td><td>26926 N/mm<sup>2</sup></td></tr><tr><td>Flow stress</td><td>y<sub>0</sub></td><td>70 N/mm<sup>2</sup></td></tr><tr><td>Hardening modulus</td><td>h</td><td>210 N/mm<sup>2</sup></td></tr><tr><td>Density</td><td>ρ</td><td>2.7 × 10<sup>-9</sup> N s<sup>2</sup>/mm<sup>4</sup></td></tr><tr><td>Expansion coefficient</td><td>α</td><td>23.8 × 10<sup>-6</sup> K<sup>-1</sup></td></tr><tr><td>Conductivity</td><td>k</td><td>150 N/s K</td></tr><tr><td>Capacity</td><td>c</td><td>0.9 × 10<sup>9</sup> mm<sup>2</sup>/s<sup>2</sup> K</td></tr><tr><td>Reference temperature</td><td>Θ<sub>0</sub></td><td>293 K</td></tr><tr><td>Dissipation factor</td><td>χ</td><td>0.9</td></tr><tr><td>Flow stress softening</td><td>ω<sub>0</sub></td><td>3 × 10<sup>-4</sup> K<sup>-1</sup></td></tr><tr><td>Hardening softening</td><td>ω<sub>h</sub></td><td>0 K<sup>-1</sup></td></tr></table>

$u[μAp] = 9.375$ N mm. From this quantity we obtain the average temperature of the isolated system $\mathcal{B}^\alpha$ for an infinite process time $\vartheta_m^\alpha = W_{\text{diss}} / (2\rho^\alpha c^\alpha V^\alpha)$ which yields the results shown in Table 2. Again we observe good agreement of the global energy check with the finite element solution based on the proposed model.

### 7.3. Upsetting of a billet

A thermoplastic upsetting process of an aluminium block is considered under plane strain conditions in this example. The block is pressed between two rigid plates which are able to conduct heat. Within this process, heating of the block occurs due to plastic dissipation within the body and frictional dissipation within the contact surface. The material parameters are listed in Table 3. We have used a temperature-dependent yield stress, given by the function $\hat{y}(e^P, \Theta) = [1 - \omega_0(\Theta - \Theta_0)]y_0 + [1 - \omega_h(\Theta - \Theta_0)]he^P$. The finite element mesh is shown in Fig. 8. Due to symmetry, only one quarter of the system was discretized by 340 finite elements with a total number of 591 unknowns.

The billet is deformed within $T = 0.0035$ s. A total of 100 time steps were needed to arrive at the final reduction of 42%. Fig. 9 shows the temperature distribution and the deformed system at the final

<figure><img src="image_9.png" /><figcaption>Fig. 8. Finite element mesh of the system.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Temperature distribution at reduction of 42%.</figcaption></figure>

reduction. The maximum temperature increase in the middle of the specimen is 28.2 K whereas the average temperature increase is 16 K within the contact interface.

References

[1] C.V. Madhusudana and L.S. Fletcher, Contact heat transfer – the last decade, AIAA J. 24 (1986) 510–523.

[2] S. Song and M.M. Yovanovich, Explicit relative contact pressure expression: dependence upon surface roughness parameters and Vickers microhardness coefficients, AIAA Paper 87-0152 (1987).

[3] J.T. Oden and J.A.C. Martins, Models and computational methods for dynamic friction phenomena, Comput. Methods Appl. Mech. Engrg. 52 (1985) 527–634.

[4] I.V. Kragelsky, M.N. Dobychin and V.S. Kombalov, Friction and Wear – Calculation Methods (translated from Russian by N. Standen) (Pergamon Press, Oxford, 1982).

[5] T.J.R. Hughes, R.L. Taylor, J.L. Sackman, A. Curnier and W. Kanoknukulchai, A finite element method for a class of contact-impact problems, Comput. Methods Appl. Mech. Engrg. 8 (1976) 249–276.

[6] B. Nour-Omid and P. Wriggers, A two-level iteration method for solution of contact problems, Comput. Methods Appl. Mech. Engrg. 54 (1986) 131–144.

[7] P. Wriggers and J.C. Simo, A note on tangent stiffness for fully nonlinear contact problems, Comm. Appl. Numer. Methods 1 (1985) 199–203.

[8] J. Hallquist, NIKE2D: An implicit, finite deformation, finite-element code for analysing the static and dynamic response of two-dimensional solids, Rept. UCRL-52678, UC-Lawrence Livermore National Laboratory, 1979.

[9] P. Wriggers, T. Vu Van and E. Stein, Finite element formulation of large deformation impact-contact problems with friction, Comput. & Structures 37 (1990) 319–331.

[10] A.B. Chaudhary and K.J. Bathe, A solution method for static and dynamic analysis of three-dimensional contact problems with friction, Comput. & Structures 24 (1985) 855–873.

[11] I.St. Doltsinis, J. Luginsland and S. Nölting, Some developments in the numerical simulation of metal forming processes, Engrg. Comput. 4 (1987) 266–280.

[12] T.A. Laursen and J.C. Simo, On the formulation and numerical treatment of finite deformation frictional contact problems, in P. Wriggers and W. Wagner, eds., Computational Methods in Nonlinear Mechanics (Springer, Berlin, 1991).

[13] A. Curnier and R.L. Taylor, A thermomechanical formulation and solution of lubricated contacts between deformable bodies, J. Lubrication Technol. 104 (1982).
[14] G. Zavarise, P. Wriggers, E. Stein and B.A. Schrefler, A numerical model for thermomechanical contact based on microscopic interface laws, Mech. Res. Comm. 19 (3) (1992) 173-182.

[15] L. Johannson and A. Klarbring, Thermoelastic frictional contact problems: modelling, finite element approximation and numerical realization (1992) preprint.

[16] J.H. Argyris and J.S. Doltsinis, On the natural formulation and analysis of large deformation coupled thermomechanical problems, Comput. Methods Appl. Mech. Engrg. 25 (1981) 195-253.

[17] C. Miehe, On the numerical treatment of thermomechanical processes (in German), Forschungs- und Seminarberichte aus dem Bereich der Mechanik der Universität Hannover, Report-No. F 88/6, 1988.

[18] I.St. Doltsinis, Aspects of modelling and computation in the analysis of metal forming, Comput. 7 (1990) 2-20.

[19] J.C. Simo and C. Miehe, Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementation, Comput. Methods Appl. Mech. Engrg. 98 (1992) 41-104.

[20] J.C. Simo and Armero, A priori stability estimates and unconditionally stable product formula algorithms for nonlinear coupled thermoplasticity, Internat. J. Plasticity 9 (1993) 749-782.

[21] P. Wriggers, On consistent tangent matrices for frictional contact problems, in J. Middleton and G.N. Pande, eds., Proc. NUMETA 87 Conf. (1987).

[22] R. Michalowski and Z. Mroz, Associated and non-associated sliding rules in contact friction problems, Arch. Mech. 30 (1978) 259-276.

[23] A. Curnier, A theory of friction, Internat. J. Solids Structures 20(7) (1984) 637-647.

[24] I. Shai and M. Santo, Heat transfer with contact resistance, Internat. J. Heat Mass Transfer 24 (1982) 465-470.

[25] M.M. Yovanovich, Thermal contact correlations, AIAA Paper 81-1164 (1981).

[26] W.L. Wood, Practical Time-Stepping Schemes, Oxford Applied Mathematics and Computing Science Series (Clarendon Press, Oxford, 1990).

[27] W. Ju and R.L. Taylor, A perturbed Lagrangian formulation for the finite element solution of nonlinear frictional contact problems, J. Theoret. Appl. Mech. 7 (1988) 1-14.

[28] A.E. Giannokopoulos, The return mapping method for the integration of friction constitutive relations, Comput. & Structures 32 (1989) 157-168.

[29] C. Miehe, A model of finite strain thermo-plasticity: formulation and aspects of the numerical implementation, IBNM Report 92/1, 1992.

[30] N. Back, M. Burdekin and A. Cowley, Review of the research on fixed and sliding joints, in S.A. Tobias and F. Koenigsberger, eds., Proc. 13th Internat. MTDR Conf. (Macmillan, London, 1973) 87-97.

[31] J.R. Barber, Indentation of the semi-infinite elastic solid by a hot sphere, Internat. J. Mech. Sci. 15 (1973) 813-819.

[32] M.G. Cooper, B.B. Mikic and M.M. Yovanovich, Thermal contact conductance, Internat. J. Heat Mass Transfer 12 (1969) 279-300.

[33] C. Miehe and E. Stein, A canonical model of multiplicative elasto-plasticity. Formulation and aspects of the numerical implementation, European J. Mech. (A/Solids) (1992) in press.

[34] C. Panek, A thermomechanical example of auto-oscillations, J. Appl. Mech. 47 (1980) 875-878.

[35] J.C. Simo and R.L. Taylor, Consistent tangent operators for rate-independent elastoplasticity, Comput. Methods Appl. Mech. Engrg. 48 (1985) 101-118.

[36] P. Wriggers, C. Miehe, M. Kleiber and J.C. Simo, On the coupled thermomechanical treatment of necking problems via finite element methods, Internat. J. Numer. Methods Engrg. 33 (1992) 869-883.