# Isogeometric Analysis and thermomechanical Mortar contact problems

M Dittmann, M Franke, I Temizer, C Hesch

► To cite this version:

M Dittmann, M Franke, I Temizer, C Hesch. Isogeometric Analysis and thermomechanical Mortar contact problems. Computer Methods in Applied Mechanics and Engineering, 2014, 274, pp.192 - 212. 10.1016/j.cma.2014.02.012. hal-01403835

HAL Id: hal-01403835

https://hal.science/hal-01403835

Submitted on 27 Nov 2016

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
# Isogeometric Analysis and thermomechanical Mortar contact problems

M. Dittmann<sup>a</sup>, M. Franke<sup>b</sup>, I. Temizer<sup>c</sup>, C. Hesch<sup>b</sup>

<sup>a</sup>Chair of Computational Mechanics, University of Siegen, Paul-Bonatz-Str. 9-11, Siegen 57068, Germany

<sup>b</sup>Institute of Mechanics, Karlsruhe Institute of Technology (KIT), Otto-Ammann-Platz 9, Karlsruhe 76131, Germany

<sup>c</sup>Department of Mechanical Engineering, Bilkent University, Ankara, Turkey

Thermomechanical Mortar contact algorithms and their application to NURBS based Iso-geometric Analysis are investigated in the context of nonlinear elasticity. Mortar methods are applied to both the mechanical field and the thermal field in order to model frictional contact, the energy transfer between the surfaces as well as the frictional heating. A series of simplifications are considered so that a wide range of established numerical techniques for Mortar methods such as segmentation can be transferred to IGA without modification. The performance of the proposed framework is illustrated with representative numerical examples.

## 1. Introduction

In this contribution, we address transient isogeometric thermomechanical contact and impact problems in the context of nonlinear elasticity. Thermoelastic material models have been investigated in detail over the past two decades, see Reese and Govindjee [32], Miehe [25] and Holzapfel and Simo [15] among many others. Corresponding formulations for thermomechanical frictional contact and impact interfaces have also been investigated and analysed in, e.g., Zavarise et al. [42], Strömberg et al. [35], Saracibar [6] and Laursen [18]. Based on the first and second laws, the thermodynamic foundations of frictional interfaces govern the formulation of appropriate constitutive laws and additionally play a major role in the construction of the accompanying numerical schemes. The emphasis of the present contribution will be on the numerical aspects.

The spatial discretisation of the bodies in contact will be carried out in the context of NURBS based Isogeometric Analysis (IGA), see Cottrell et al. [4] for a comprehensive review. Contact problems for IGA have been addressed in a series of papers throughout the past years, see De Lorenzis et al. [20,21] and Temizer et al. [38,36,39]. In these works, a Knot-to-Surface (KTS) method has been developed and extended to Mortar based contact formulations. Recently, Matzen et al. [24] have proposed a collocation based approach, analogous to the well-known Node-to-Surface (NTS) method. See also Kim and Youn [17] for a Mortar approach and Lu [22] for an alternative contact treatment.

Mortar formulations in the context of IGA domain decomposition problems have been presented in Hesch and Betsch [14]. In the present contribution, we extend the ideas developed therein to thermomechanical contact in order to achieve a variationally consistent Mortar formulation for the discrete contact interface. In particular, the Mortar projections will be calculated via a newly developed segmentation procedure of the surface intersections, see Puso et al. [31] for a discussion
about different spatial integration schemes. For the thermal contributions, the Mortar concept will be applied by introducing triple Mortar integrals to accurately capture the frictional dissipation contribution to the contact heat flux and to establish a correct thermal interaction among the contacting surfaces. See also Hüeber and Wohlmuth [16] for Mortar methods applied to thermoelasticity.

For transient impact problems, the bodies will be discretised in time using structure preserving integration schemes, following the approach in Hesch and Betsch [10]. This allows us to investigate the effect of different integration schemes for the contact contributions, following Franke et al. [8] in the context of purely mechanical NTS methods. Specifically, we will devote particular attention to the conservation of angular momentum and its possible violation.

An outline of the paper is as follows. The underlying thermomechanical framework for the bodies and the frictional contact interface is presented in Section 2. The spatial discretisation using IGA as well as the Mortar formulation for the semi-discrete system is developed in Section 3. The temporal discretisation is outlined in Section 4, followed by representative numerical examples in Section 5 and conclusive remarks in Section 6.

## 2. Governing equations

In this section we summarize the variational form of the thermomechanical theory along with a most general description of frictional contact contributions, embedded within the thermomechanical framework. In particular, we consider Lipschitz bounded domains $B_0^{(i)} \subset \mathbb{R}^n$, $n \in [2, 3]$ in their reference configuration, where the upper index (i) will denote the respective body in the remainder of this article.<fnref n="1" /> Furthermore, we introduce the mapping

$$ \varphi(\mathbf{X}, t) : B_0 \times \mathcal{I} \to \mathbb{R}^n, \qquad (1) $$

to characterise the time dependent deformation along with the absolute temperature

$$ \theta(\mathbf{X}, t) : B_0 \times \mathcal{I} \to \mathbb{R} \qquad (2) $$

for the time interval $t \in \mathcal{I} = [0, T]$ elapsed during the motion. Here, $\mathbf{X} \in B_0$ labels material points in the reference configuration and both fields are assumed to be sufficiently smooth.

### 2.1. Finite strain thermoelastodynamics

In a first step we postulate that the material behaviour is governed by a Helmholtz energy density function $\Psi(\mathbf{C}, \theta)$, where $\mathbf{C} = \mathbf{F}^T \mathbf{F}$ denotes the right Cauchy-Green tensor and $\mathbf{F} : B_0 \times \mathcal{I} \to \mathbb{R}^{n \times n}$, $\mathbf{F} = D\varphi$ the deformation gradient. Accordingly, we define the local constitutive relations

$$ \Sigma = 2 \frac{\partial \Psi}{\partial \mathbf{C}}, \qquad (3) $$

$$ \eta = - \frac{\partial \Psi}{\partial \theta}. \qquad (4) $$

Therein, $\Sigma$ represents the second Piola-Kirchhoff stress tensor and $\eta$ the local entropy density. A third constitutive relation is required to account for the heat transfer

$$ \mathbf{Q} = -\hat{\mathbf{K}}(\mathbf{C}, \theta)\nabla_{\mathbf{X}}(\theta), \qquad (5) $$

known as Duhamel's law of heat conduction, where $\hat{\mathbf{K}}(\mathbf{C}, \theta)$ is the material thermal conductivity tensor. Introducing the space of virtual or admissible test functions for the deformation as well as for the absolute temperature

$$ \mathcal{V}^{\varphi} = \{ \delta\varphi \in \mathcal{H}^1(B_0) | \delta\varphi = \mathbf{0} \text{ on } \partial\mathcal{B}_0^{\varphi} \}, \qquad (6) $$

$$ \mathcal{V}^{\theta} = \{ \delta\theta \in \mathcal{H}^1(B_0) | \delta\theta = 0 \text{ on } \partial\mathcal{B}_0^{\theta} \}, \qquad (7) $$

where $\mathcal{H}^1$ denotes the Sobolev functional space of square integrable functions and derivatives, the weak form of the balance of linear momentum and the energy balance equation reads

$$ G_{\varphi} := \int_{B_0} \rho_0 \delta\varphi \cdot \ddot{\varphi} + \mathbf{\Sigma} : \mathbf{F}^T \nabla_{\mathbf{X}}(\delta\varphi) \, dV - \int_{B_0} \delta\varphi \cdot \bar{\mathbf{B}} \, dV - \int_{\partial\mathcal{B}_0^g} \delta\varphi \cdot \bar{\mathbf{T}} \, dA = 0, \qquad (8) $$

$$ G_{\theta} := \int_{B_0} \delta\theta\dot{\eta} - \mathbf{Q} \cdot \mathrm{Grad}(\delta\theta) \, dV - \int_{B_0} \delta\theta \bar{\mathbf{R}} \, dV - \int_{\partial\mathcal{B}_0^Q} \delta\theta \bar{\mathbf{Q}} \, dA = 0. \qquad (9) $$

Note that $\dot{\eta}(\mathbf{C}, \theta)$ can be decomposed into components denoted as heat capacity and latent heat, where the latter one is responsible for the Gough-Joule effect. The external contributions at the boundary are specified by Dirichlet and Neumann boundary conditions on the mechanical and thermal field, respectively

<footnote n="1">¹ If convenient and unique the superscript index will be omitted for the ease of exposition. Moreover, we make use of the summation convention for repeated indices.</footnote>
$$
\begin{align}
\boldsymbol{\varphi} &= \boldsymbol{\varphi}' \text{ on } \partial\mathcal{B}_0^{\varphi} \times [0, T], \quad \theta = \tilde{\theta} \text{ on } \partial\mathcal{B}_0^{\theta} \times [0, T], \tag{10} \\
\boldsymbol{PN} &= \boldsymbol{T} \text{ on } \partial\mathcal{B}_0^{\sigma} \times [0, T], \quad \boldsymbol{Q} \cdot \boldsymbol{N} = -\tilde{Q} \text{ on } \partial\mathcal{B}_0^{Q} \times [0, T],
\end{align}
$$

whereas $\bar{\mathbf{B}}$ and $\bar{\mathbf{R}}$ represents the material descriptions of prescribed body force and heat supply per unit volume. To complete the set of equations for the initial boundary value problem (IBVP), initial conditions in $\mathcal{B}_0$ at time $t=0$ are specified by

$$
\varphi(\mathbf{X}, 0) = \varphi_0, \quad \dot{\varphi}(\mathbf{X}, 0) = \mathbf{v}_0, \quad \theta(\mathbf{X}, 0) = \theta_0 \quad \text{in } \mathcal{B}_0. \tag{11}
$$

The IBVP defined above has to be thermodynamical consistent, i.e. the IBVP has to obey the fundamental first and second law of thermodynamics. The second law is of constitutive nature and is consistently represented in the definitions of (3)–(5), i.e. we will not enforce the second law explicitly throughout this paper. The first law postulates the validity of the global energy balance. To prove this for the problem at hand, we introduce suitable substitutions of the variations $\delta\varphi = \dot{\varphi}$ in (8) and obtain

$$
\int_{\mathcal{B}_0} \frac{d}{dt} \left( \frac{1}{2} \rho_0 \dot{\varphi} \cdot \dot{\varphi} \right) + \dot{\Psi}(\mathbf{C}, \eta) - \frac{\partial \Psi}{\partial \theta} \dot{\theta} dV = \int_{\mathcal{B}_0} \dot{\varphi} \cdot \bar{\mathbf{B}} dV + \int_{\partial \mathcal{B}_0^Q} \dot{\varphi} \cdot \bar{\mathbf{T}} dA. \quad (12)
$$

A Legendre transformation is applied next to recast the system in terms of the inner energy density $e(\mathbf{C}, \theta) = \Psi(\mathbf{C}, \theta) + \eta\theta$, such that

$$
\frac{d}{dt} T + \int_{\mathcal{B}_0} \dot{e} - \dot{\eta}\theta \, dV = P^{\text{ext}}, \quad (13)
$$

where $T$ denotes the kinetic energy and $P^{\text{ext}}$ the power of the external mechanical contributions. Substituting $\delta\theta = \mu$ in (9), where $\mu \in \mathbb{R}$ is arbitrary and constant, gives

$$
\int_{\mathcal{B}_0} \theta \dot{\eta} dV = \int_{\mathcal{B}_0} \tilde{R} dV - \int_{\partial \mathcal{B}_0^Q} \mathbf{Q} \cdot \mathbf{N} dA. \quad (14)
$$

Insertion in (13) yields the global energy balance equation

$$
\frac{d}{dt}[T + E] = P^{\text{ext}} + Q, \tag{15}
$$

where $Q$ is the total net heating of the continuum body and $E = \int_{\mathcal{B}_0} e dV$. The last equation represents the fundamental first law of thermodynamics for the coupled thermomechanical system.

2.2. Contact formulation

Assuming two bodies ($i \in [1, 2]$) in a contact situation, the boundaries are further subdivided and satisfy

$$
\partial \mathcal{B}_0^{(i,c)} \cup \partial \mathcal{B}_0^{(i),\varphi} \cup \partial \mathcal{B}_0^{(i),\sigma} = \partial \mathcal{B}_0^{(i),c} \cup \partial \mathcal{B}_0^{(i),\theta} \cup \partial \mathcal{B}_0^{(i),Q} = \partial \mathcal{B}_0^{(i)} \quad (16)
$$

and

$$
\begin{align}
& \partial \mathcal{B}_0^{(i),c} \cap \partial \mathcal{B}_0^{(i),\varphi} = \partial \mathcal{B}_0^{(i),\varphi} \cap \partial \mathcal{B}_0^{(i),\sigma} = \partial \mathcal{B}_0^{(i),\sigma} \cap \partial \mathcal{B}_0^{(i),c} = \emptyset, \tag{17} \\
& \partial \mathcal{B}_0^{(i),c} \cap \partial \mathcal{B}_0^{(i),\theta} = \partial \mathcal{B}_0^{(i),\theta} \cap \partial \mathcal{B}_0^{(i),Q} = \partial \mathcal{B}_0^{(i),Q} \cap \partial \mathcal{B}_0^{(i),c} = \emptyset. \nonumber
\end{align}
$$

In this connection, $\partial\mathcal{B}_0^{(i),c}$ represents the potential contact area on surface (i). Assuming that the local linear momentum bal-
ance equation across the contact interface

$$
\mathbf{t}^{(1)} \, dA = - \mathbf{t}^{(2)} \, dA, \tag{18}
$$

where $t^{(i)}$ denote the Piola tractions$^2$ associated with surface $\partial \mathcal{B}_0^{(i),c}$, is valid, the global variational statements (8) and (9) can be
augmented as follows

$$
G_{\varphi} + G_{\varphi}^{c} = 0, \tag{19}
$$

$$
G_{\theta} + G_{\theta}^{c} = 0, \tag{20}
$$

where the contact contributions are given by

$$
G_{\varphi}^{c} = \int_{\partial \mathcal{B}_{0}^{(1),c}} \mathbf{t}^{(1)} \cdot (\delta\boldsymbol{\varphi}^{(1)} - \delta\boldsymbol{\varphi}^{(2)}) dA, \quad (21)
$$

$$
G_{\theta}^{c} = - \int_{\partial \mathcal{B}_{0}^{(1),c}} \delta\theta^{(1)} Q_{c}^{(1)} + \delta\theta^{(2)} Q_{c}^{(2)} dA. \quad (22)
$$

<footnote>$^2$ Note that the traction $\mathbf{t}^{(i)}$ as well as the heat flux $Q_c^{(i)}$ can be considered as Cauchy traction and flux associated with the actual surfaces $\partial \mathcal{B}^{(2),c} = \partial \mathcal{B}^{(1),c}$ and pulled back onto $\partial \mathcal{B}_0^{(1),c}$.</footnote>
Here, $Q_c^{(i)}$ represent the heat supply at the respective contact interface (i). To account for different contact reactions in normal and tangential direction, we may decompose the Piola tractions

$$ \mathbf{t}^{(1)} = t_N \mathbf{n} + \mathbf{t}_T, \quad \mathbf{t}_T \cdot \mathbf{n} = 0, \qquad (23) $$

where $\mathbf{n} := \mathbf{n}^{(1)}$ denotes the current outward normal vector on $\partial B_0^{(1),c}$, and define the local normal and tangential gap functions

$$ g_N = \mathbf{n} \cdot (\boldsymbol{\varphi}^{(1)} - \boldsymbol{\varphi}^{(2)}), \quad g_T = (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) \cdot (\boldsymbol{\varphi}^{(1)} - \boldsymbol{\varphi}^{(2)}), \qquad (24) $$

see Puso & Laursen [30]. Note that the components of the tangential tractions and the tangential gap function can be written in terms of co-/contravariant basis vectors

$$ \mathbf{t}_T = t_{T,\alpha} \mathbf{a}^{\alpha}, \quad \mathbf{g}_T = g_T^{\alpha} \mathbf{a}_{\alpha}, \qquad (25) $$

which is often useful for the definition of the constitutive frictional law.

## 2.3. Constitutive equations for contact

The contact reactions are in general characterised by constitutive laws for the mechanical as well as for the thermal field. In particular, Coulomb's friction law is imposed on the mechanical side, whereas a local dissipation potential is introduced to account for the contact heat flux.

### Normal contact

Standard Karush-Kuhn-Tucker conditions are employed in normal directions

$$ g_N \le 0, \quad t_N \ge 0, \quad t_N g_N = 0, \qquad (26) $$

to enforce the non-penetration condition. Further constitutive relations to account for, e.g., adhesion effects, can be embedded without changing the considered framework.

### Friction law

Given the definitions above and following standard arguments as outlined, for example, in Yang & Laursen [41] and in Puso & Laursen [30], Coulomb's friction law can be written as

$$ \hat{\phi}_c := \|\mathbf{t}_T\| - \mu|t_N| \le 0, \qquad (27) $$

$$ \dot{\zeta} \ge 0, \qquad (28) $$

$$ \hat{\phi}_c \dot{\zeta} = 0. \qquad (29) $$

Here, $\mu := \mu(\theta^1, \theta^2)$ denote the friction coefficient, depending on the temperatures at the contact interface and $\dot{\zeta}$ a consistency parameter, where $\dot{\zeta} = 0$ represents stick and $\dot{\zeta} > 0$ slip, in analogy to the plastic multiplier in plasticity. Now, the Lie derivative $\mathcal{L}\mathbf{t}_T = t_{T,\alpha}\mathbf{a}^{\alpha}$ of the frictional tractions reads

$$ \mathcal{L}\mathbf{t}_T = \epsilon_T \left( \dot{\mathbf{g}}_T^s - \dot{\zeta} \frac{\mathbf{t}_T}{\|\mathbf{t}_T\|} \right), \qquad (30) $$

where $\epsilon_T$ is a tangential penalty parameter. Note that the subsequently used additive split of the tangential gap into a reversible (elastic) part $g_T^e$ and an irreversible (inelastic) part $g_T^s$ follows from the penalisation of the stick condition.

**Remark 1.** In contrast to the normal contributions, where we enforce a geometrical non-penetration condition via Lagrange multipliers, the constitutive relations in tangential direction are enforced via a penalty method. Other choices are possible and can be applied without changing the proposed framework.

### Thermal contact interface

The thermodynamic analysis [35,6] of the deformed contact interface suggests an additive decomposition

$$ q_c^{(i)} = q_f^{(i)} + q_t^{(i)}, \qquad (31) $$

of the contact heat flux $q_c^{(i)}$ on $\partial B^{(1),c}$ into frictional dissipation ($q_f^{(i)}$) and thermal interaction ($q_t^{(i)}$) contributions. In particular,

$$ q_f^{(1)} + q_f^{(2)} = \mathcal{F} \qquad (32) $$

and

$$ q_t^{(1)} + q_t^{(2)} = 0 \qquad (33) $$
are admitted. Here, $\mathcal{F}$ is the frictional dissipation such that $\mathcal{F} da = \mathcal{F}_o dA$ with $\mathcal{F}_o = \mathbf{t}_T \cdot \dot{\mathbf{g}}_T^s = \mu t_N \|\dot{\mathbf{g}}_T^s\|$. The constitutive modeling of $q_f^{(i)}$ is based on a phenomenological partitioning

$$q_f^{(i)} = \gamma^{(i)} \mathcal{F}, \qquad (34)$$

of $\mathcal{F}$, which implies $Q_f^{(i)} = \gamma^{(i)} \mathcal{F}_o$, where the relative effusivities $\gamma^{(i)} \ge 0$ satisfy $\gamma^{(1)} + \gamma^{(2)} = 1$. The relative effusivities are interface constitutive parameters that could depend on a number of variables, including the temperature and the pressure. Presently, they will be assumed to be constants.

The constitutive modeling of $q_t^{(1)} = -q_t^{(2)}$ must be based on the consideration of various heat exchange mechanisms across surfaces in vicinity [23,2,27]. In a continuum setting, primarily three mechanisms are of concern: conduction, convection and radiation. Among these, conduction is a contact interaction due to the temperature jump $\vartheta_c = \theta_c^{(1)} - \theta_c^{(2)}$ across the contact interface whereas the latter two are of non-contact type. The radiation contribution is often negligible unless the temperatures involved are too high. Macroscopically contacting surfaces, on the other hand, may simultaneously involve the remaining two mechanisms due to inherent roughness since the microscopic interface topography not only has contacting asperities but also gaps containing an interstitial material which serves as a medium for convection. Moreover, in view of the multiscale nature of roughness, the thermal interaction of the contacting asperities is also governed by both mechanisms. Due to the small length scales associated with the gaps, the convective mechanism degenerates into conduction through the interstitial medium. Consequently, the overall heat exchange across the contact interface can be modeled as a function of $\vartheta_c$. Finally, an increasing contact pressure leads to higher microscopic conformity among the contacting surfaces such that the macroscopic interface resistance decreases. In order to reflect this effect, constitutive models for the overall heat exchange across the contact interface typically admit a nonlinear dependence on the contact pressure in the small deformation regime. In the present work where large deformations are of concern [9,3], the constitutive choice is based on referential quantities in the simplified form

$$Q_t^{(1)} = -Q_t^{(2)} = -k_\theta |t_N| \vartheta_c, \qquad (35)$$

where $k_\theta$ is assumed to be a constant. Such a formulation may alternatively be written in terms of a dissipation potential $\mathcal{H} = -k_\theta |t_N| \vartheta_c^2/2$ such that $Q_t^{(1)} = \partial \mathcal{H} / \partial \vartheta_c$, which can assist in the derivation of alternative algorithms that are suitable for the modeling and simulation of microscale contact. Combining all contributions, the final expressions for the contact heat fluxes $q_c^{(i)}$ pulled back onto $\partial \mathcal{B}_o^{(1),c}$ are

$$Q_c^{(1)} = \gamma^{(1)} \mathcal{F}_o - k_\theta |t_N| \vartheta_c, \quad Q_c^{(2)} = \gamma^{(2)} \mathcal{F}_o + k_\theta |t_N| \vartheta_c. \qquad (36)$$

### Thermodynamical consistency

Assuming that the first and second law of thermodynamics holds for the contact interface as well, we introduce a local boundary Helmholtz energy density function $\Psi_c := \Psi_c(\varphi, \theta)$ and state the global balance equations across the interface

$$\int_{\partial \mathcal{B}_o^{(1),c}} \delta_{\varphi} \Psi_c \, dA = - \int_{\partial \mathcal{B}_o^{(1),c}} \mathbf{t}^{(1)} \cdot (\delta \varphi^{(1)} - \delta \varphi^{(2)}) \, dA, \qquad (37)$$

$$\int_{\partial \mathcal{B}_o^{(1),c}} \delta \theta_c \dot{\eta}_c \, dA = \int_{\partial \mathcal{B}_o^{(1),c}} \delta \theta^{(1)} Q_c^{(1)} + \delta \theta^{(2)} Q_c^{(2)} \, dA. \qquad (38)$$

Note that we do not account for the heat capacity of the contact interface, although this can be included in the left hand side of (38). Substituting $\delta\varphi = \dot{\varphi}$, we obtain for the contact interface

$$\int_{\partial \mathcal{B}_o^{(1),c}} \dot{\Psi}_c + \eta_c \dot{\theta}_c \, dA = - \int_{\partial \mathcal{B}_o^{(1),c}} t_N \dot{g}_N + \mathbf{t}_T \cdot (\dot{\mathbf{g}}_T^e + \dot{\mathbf{g}}_T^s) \, dA. \qquad (39)$$

Applying again a Legendre transformation $\dot{\Psi}_c + \eta_c \dot{\theta}_c = \dot{e}_c - \dot{\eta}_c \dot{\theta}_c$, substituting $\delta\theta = \mu, \mu \in \mathbb{R}$ in (38) and subsequent insertion yields the local energy balance contributions of the contact boundary

$$\int_{\partial \mathcal{B}_o^{(0),c}} \dot{e}_c \, dA = \int_{\partial \mathcal{B}_o^{(0),c}} -t_N \dot{g}_N - \mathbf{t}_T \cdot (\dot{\mathbf{g}}_T^e + \dot{\mathbf{g}}_T^s) + Q_c^{(1)} + Q_c^{(2)} \, dA. \qquad (40)$$

Here, $\dot{e}_c$ is the stored energy per unit area on the contact interface. An important issue has been shown here: Since $\gamma_1 + \gamma_2 = 1$ the dissipative energy $\mathbf{t}_T \cdot \dot{\mathbf{g}}_T^s$ drops out of the last equation and the mechanical energy is directly dissipated into the thermal field without changing the inner energy. By enforcing (26) exactly, we can state that the time derivative $\dot{g}$ vanishes as well. Following additionally the arguments in Oancea & Laursen [26], stating that the elastic part of the tangential gap is small enough to be neglected, we obtain

$$\int_{\partial \mathcal{B}_o^{(0),c}} \dot{e}_c \, dA = 0, \qquad (41)$$

which is the expected result if the heat capacity of the interface is omitted. With regard to (38), the local entropy production rate at the interface can be written as follows
<table><caption>**Table 1**

Weak form of the coupled thermomechanical contact problem.</caption>
  <tbody>
    <tr>
      <td>(1) Mechanical field</td>
      <td>(43)</td>
    </tr>
    <tr>
      <td>(2) Thermal field</td>
      <td>(44)</td>
    </tr>
    <tr>
      <td>(3) Interface conditions</td>
      <td></td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&bull; Normal contact</td>
      <td>(45)</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&bull; Tangential contact</td>
      <td>(46)</td>
    </tr>
    <tr>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&bull; Thermal contact</td>
      <td>(47)</td>
    </tr>
  </tbody>
</table>

$$
\dot{\eta}_c = \frac{Q_c^{(1)}}{\theta_c^{(1)}} + \frac{Q_c^{(2)}}{\theta_c^{(2)}} \ge 0. \tag{42}
$$

Assuming that the tangential penalty parameter $\epsilon_T$ is high enough, we can insert (30) into (36). Taking additionally (28) into account gives immediately $\dot{\eta} \|\mathbf{t}_T\| \ge 0$ for the mechanical dissipation of energy into the thermal field. Insertion of the remaining parts of (36) into (42) shows after some straightforward calculations, that the second law of thermodynamics is valid for the thermomechanical contact problem at hand. The resulting weak forms and the constitutive contact laws are summarised in Table 1.

**3. Isogeometric discretisation**

To perform the spatial discretisation for the numerical solution of the variational formulation in (43) and (44), each do-
main is subdivided into a finite set of non-overlapping elements $e \in \mathbb{E}$

$$
\mathcal{B}_0 \approx \mathcal{B}_0^h = \bigcup_e \mathcal{B}_{0,e}, \quad \forall e \in \mathbb{E}. \tag{48}
$$

For the numerical solution of the coupled problem, a finite element framework is considered. Therefore, we introduce finite
dimensional approximations of the deformation and its variation, such that

$$
\boldsymbol{\varphi}^h = \sum_{A \in \omega} \mathbf{R}^A \boldsymbol{q}_A, \quad \delta \boldsymbol{\varphi}^h = \sum_{A \in \omega} \mathbf{R}^A \delta \boldsymbol{q}_A. \tag{49}
$$

Here, $A \in \omega = \{1, \dots, n_{\text{node}}\}$ and $\mathbf{R}^A(\mathbf{X}) : \mathcal{B}_0 \to \mathbb{R}$ are global shape functions. In particular, we make use of NURBS based shape functions

$$
\mathbf{R}^A = \mathbf{R}_{p,q,r}^{i,j,t}(\xi) = \frac{N_{i,p}(\xi) M_{i,q}(\eta) L_{t,r}(\zeta) w_{i,j,t}}{\sum_{i=1}^{n} \sum_{j=1}^{m} \sum_{l=1}^{l} N_{i,p}(\xi) M_{j,q}(\eta) L_{t,r}(\zeta) w_{i,j,l}}, \quad (50)
$$

where *n*, *m*, *l* denote the number of control points along each parametric direction. In addition, *p*, *q*, *r* denotes the order of the non-rational B-Splines *N*, *M* and *L*, recursively defined as follows

$$
N_{i,p} = \frac{\xi_i - \xi_i}{\xi_{i+p} - \xi_i} N_{i,p-1}(\xi) + \frac{\xi_{i+p+1} - \xi_i}{\xi_{i+p+1} - \xi_{i+1}} N_{i+1,p-1}(\xi), \quad (51)
$$
beginning with

$$N_{i,0}(\xi) = \begin{cases} 1 & \text{if } \xi_i \le \xi < \xi_{i+1} \\ 0 & \text{otherwise} \end{cases}. \qquad (52)$$

The definition for M and L follows analogously. Furthermore, $w_{i,j,t}$ are NURBS weights, for details see Piegl & Tiller [28]. The shape functions are associated with a net of control points $q_A$, such that a geometrical map $\tilde{\mathcal{G}}: \tilde{B}_0 \to B_0$ can be defined to link the parameter and the physical space

$$\boldsymbol{\varphi}^h := \tilde{\mathcal{G}}(\xi) = \boldsymbol{R}_{p,q,r}^{i,j,t}(\xi) \boldsymbol{q}_{i,j,t}, \qquad (53)$$

see da Veiga et al. [5]. The finite element mesh is now defined via knot vectors, which subdivide the parameter space into elements.

The discretisation of the thermal field is given as follows

$$\theta^h = \sum_{A \in \omega} R^A \Theta_A, \quad \delta\theta^h = \sum_{A \in \omega} R^A \delta\Theta_A, \qquad (54)$$

such that we obtain a conforming discretisation of both fields.

### 3.1. Semi-discrete formulation of the coupled problem

For the ease of exposition, we discretise first the weak form in (8) as well as in (9) and add the contact contributions later on. The semi-discrete balance of linear momentum reads

$$\delta \boldsymbol{q}_A \cdot \left[ \boldsymbol{M}^{AB} \ddot{\boldsymbol{q}}_B + \int_{B_0} \boldsymbol{S}^{AB} \mathrm{d}\boldsymbol{V} \boldsymbol{q}_B \right] = \delta \boldsymbol{q}_A \cdot \left[ \boldsymbol{F}^{A,\mathrm{ext}} \right], \qquad (55)$$

whereas the semi-discrete energy balance equation is given by

$$\delta \Theta_A \left[ \Gamma^{AB} \Theta_B - \int_{B_0} K^{AB} \mathrm{d}V \Theta_B \right] = \delta \Theta_A \left[ Q^{A,\mathrm{ext}} \right]. \qquad (56)$$

Here, the consistent mass matrix

$$\boldsymbol{M}^{AB} = \int_{B_0} \rho_0 R^A R^B \mathrm{d}V, \qquad (57)$$

as well as the internal and external contributions

$$\int_{B_0} S^{AB} \mathrm{d}V = \int_{B_0} \nabla_X (R^A) \cdot \boldsymbol{S}^h \nabla_X (R^B) \mathrm{d}V, \qquad (58)$$

$$\boldsymbol{F}^{A,\mathrm{ext}} = \int_{B_0} R^A \boldsymbol{B} \mathrm{d}V + \int_{\partial B_0^h} R^A \boldsymbol{T} \mathrm{d}A \qquad (59)$$

are defined in the usual way. Note that the semi-discrete second Piola-Kirchhoff stress tensor is defined by $\boldsymbol{S}^h = 2\partial\Psi(\boldsymbol{C}^h, \theta^h)/\partial\boldsymbol{C}^h$ where $\boldsymbol{C}^h = \boldsymbol{q}_A \cdot \boldsymbol{q}_B \nabla_X(R^A) \otimes \nabla_X(R^B)$. Analogously, we obtain for the thermal contributions

$$\Gamma^{AB} = \int_{B_0} \dot{\eta}^h R^A R^B \mathrm{d}V, \qquad (60)$$

where the discrete entropy is given by $\eta^h = -\partial\Psi(\boldsymbol{C}^h, \theta^h)/\partial\theta^h$. Furthermore, the second term in (56) reads

$$\int_{B_0} K^{AB} \mathrm{d}V = \int_{B_0} \nabla_X (R^A) \cdot \boldsymbol{K}(\boldsymbol{C}^h, \theta^h) \nabla_X (R^B) \mathrm{d}V \qquad (61)$$

and the corresponding right hand side is given by

$$Q^{A,\mathrm{ext}} = \int_{B_0} R^A \tilde{R} \mathrm{d}V + \int_{\partial B_0^h} R^A \tilde{Q} \mathrm{d}A. \qquad (62)$$

### 3.2. Mortar method

To construct a Mortar based approach for the thermomechanical contact interface, the variational statements in (21) and (22) have to be discretised as well. Therefore, the space of admissible test functions for the discrete Lagrange multiplier field is introduced as

$$\mathcal{M}^h = \left\{ \delta \boldsymbol{t}^{(1),h} \in \mathcal{L}^2 \left( \partial B_0^{(1),c} \cap \partial B_0^{(2),c} \right) \right\}. \qquad (63)$$
To deal with non-conforming discretisations and NURBS of arbitrary order, we evaluate the geometrical map $\tilde{\mathcal{F}}(\xi)$ at the boundary of the parameter space, corresponding to the contact boundary in physical space to obtain a set of nodes $\tilde{\omega}^{(1)} = [\tilde{\boldsymbol{q}}_1, \dots, \tilde{\boldsymbol{q}}_{n_{\text{surf}}}]$ on surface (1), where $n_{\text{surf}}$ corresponds to the number of physical nodes on the surface geometry of $\mathcal{B}^{(1)}$. As shown in Hesch & Betsch [14], a linear interpolation of the Lagrange multipliers and their variations

$$ \mathbf{t}^{(1),h} = \sum_{A \in \tilde{\omega}^{(1)}} N^A \lambda_A, \quad \delta \mathbf{t}^{(1),h} = \sum_{A \in \tilde{\omega}^{(1)}} N^A \delta \lambda_A, \qquad (64) $$

where $N^A : \partial \mathcal{B}_0^{(1),c} \to \mathbb{R}$ are $(n-1)$ dimensional linear Lagrangian shape functions associated with nodes $A \in \tilde{\omega}^{(1)}$ in physical space, is sufficient. Note that this last statement holds even if we evaluate the nodal contributions $\lambda_A$ via penalty methods. Again, decomposing the tractions $\lambda_A$ into normal $\lambda_{N,A}$ and tangential $\lambda_{T,A}$ components and subsequently inserting them in the contact contributions within (43) yields

$$ G_{\varphi}^{c,h} = \lambda_{N,A} \mathbf{n} \cdot [\mathbf{n}^{AB} \delta \boldsymbol{q}_{B}^{(1)} - \mathbf{n}^{AC} \delta \boldsymbol{q}_{C}^{(2)}] + \lambda_{T,A} \cdot (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) [\mathbf{n}^{AB} \delta \boldsymbol{q}_{B}^{(1)} - \mathbf{n}^{AC} \delta \boldsymbol{q}_{C}^{(2)}], \qquad (65) $$

whereas the discrete contact contributions in (44) reads

$$ G_{\varphi}^{c,h} = -\delta \Theta_A^{(1)} \{ \gamma_1 \lambda_{T,B} \cdot (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) [\mathbf{m}^{ABC} \dot{\boldsymbol{q}}_C^{(1),s} - \mathbf{m}^{ABD} \dot{\boldsymbol{q}}_D^{(2),s}] - k_\theta |\lambda_{N,B}| [\mathbf{m}^{ABC} \Theta_C^{(1)} - \mathbf{m}^{ABD} \Theta_D^{(2)}] \} \\ - \delta \Theta_A^{(2)} \{ \gamma_2 \lambda_{T,B} \cdot (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) [\mathbf{m}^{ABC} \dot{\boldsymbol{q}}_C^{(1),s} - \mathbf{m}^{ABD} \dot{\boldsymbol{q}}_D^{(2),s}] + k_\theta |\lambda_{N,B}| [\mathbf{m}^{ABC} \Theta_C^{(1)} - \mathbf{m}^{ABD} \Theta_D^{(2)}] \}. \qquad (66) $$

Here, the Mortar integrals are defined as follows

$$ \begin{aligned} n^{AB} &= \int_{\partial \mathcal{B}_0^{(1),c}} N^A(\xi^{(1)}) R^B(\xi^{(1)}) dA, \\ n^{AC} &= \int_{\partial \mathcal{B}_0^{(1),c}} N^A(\xi^{(1)}) R^C(\xi^{(2)}) dA \end{aligned} \qquad (67) $$

and the triple Mortar integrals are given by

$$ \begin{aligned} m^{ABC} &= \int_{\partial \mathcal{B}_0^{(1),c}} R^A(\xi^{(1)}) N^B(\xi^{(1)}) R^C(\xi^{(1)}) dA, \\ m^{ABD} &= \int_{\partial \mathcal{B}_0^{(1),c}} R^A(\xi^{(1)}) N^B(\xi^{(1)}) R^D(\xi^{(2)}) dA. \end{aligned} \qquad (68) $$

The evaluation of the triple Mortar integrals $\bar{m}$ follows analogously using $R^A(\xi^{(2)})$ as first shape function in (68). The contact conditions required for the evaluation of the corresponding contact forces and the heat transfer of the semidiscrete system are summarised in Table 2.

For the numerical evaluation of the mortar integrals, a suitable segmentation process is necessary, subdividing the surface of the parameter space on both sides in triangles. We focus on the evaluation of the standard Mortar integrals, since the triple Mortar integrals follows analogously, extended by additional shape functions $R^A$. The main goal of this construction is to provide a common parametrisation of both surfaces to apply a Gauss quadrature. Therefore, a linear transformation based on bilinear, triangular shape functions $M^K$ is introduced

$$ \xi^{(i),h}(\boldsymbol{\eta}) = \sum_{K=1}^{3} M^K(\boldsymbol{\eta}) \xi_K^{(i)}, \qquad (69) $$

such that we obtain the segment contributions of the Mortar integrals

$$ \begin{aligned} n^{K\beta} &= \int_{\Delta} N^K(\xi^{(1),h}(\boldsymbol{\eta})) R^\beta(\xi^{(1),h}(\boldsymbol{\eta})) J_{\text{seg}} d\eta, \\ n^{K\zeta} &= \int_{\Delta} N^K(\xi^{(1),h}(\boldsymbol{\eta})) R^\zeta(\xi^{(2),h}(\boldsymbol{\eta})) J_{\text{seg}} d\eta, \end{aligned} \qquad (70) $$

which have to be assembled into the global system, see Hesch & Betsch [11,13] and the references therein. Note that the triangle symbols in Eq. (70) represent the triangle reference element. Finally, the Jacobian of each segment is required

$$ J_{\text{seg}} = \| \mathbf{A}_1(\xi^{(1),h}(\boldsymbol{\eta})) \times \mathbf{A}_2(\xi^{(1),h}(\boldsymbol{\eta})) \| \det(D\xi(\boldsymbol{\eta})), \qquad (71) $$

where $\mathbf{A}_\alpha(\xi) = R_\alpha^A(\xi) \boldsymbol{q}_A$ denote the tangential vectors in the reference configuration.

To construct the segments by calculating the nodes $\xi_K^{(j)}$ within the parameter space, a stable, but tedious method has been proposed in Hesch & Betsch [14]. Here, we propose a major simplification of the segmentation algorithm for isogeometric discretised bodies, which allows the usage of standard segmentation libraries. Therefore, we reuse the already introduced set of nodes $\tilde{\omega}^{(1)}$ on the physical contact boundary of $\mathcal{B}^{(1)}$ and apply a linear Lagrangian discretisation in analogy to the linear interpolation of the Lagrange multiplier field in (64). Based on this artificially discretised, virtual surface and a corresponding
<caption>**Table 2**
Semi-discrete form of the coupled thermomechanical contact problem.</caption>

$$
(1) \text{ Mechanical field} \qquad \sum_{i=1}^{2} \left( \delta \boldsymbol{q}_{A}^{(i)} \cdot \left[ M^{AB} \dot{\boldsymbol{q}}_{B}^{(i)} + \int_{\boldsymbol{g}_{0}} S^{AB} \mathrm{d}V \boldsymbol{q}_{B}^{(i)} - \boldsymbol{F}^{(i)A} \mathrm{ext} \right] \right) + \lambda_{N,A} \boldsymbol{n} \cdot \left[ n^{AB} \delta \boldsymbol{q}_{B}^{(1)} - n^{AC} \delta \boldsymbol{q}_{C}^{(2)} \right] + \lambda_{T,A} \cdot (\boldsymbol{I} - \boldsymbol{n} \otimes \boldsymbol{n}) \left[ n^{AB} \delta \boldsymbol{q}_{B}^{(1)} - n^{AC} \delta \boldsymbol{q}_{C}^{(2)} \right] = 0 \qquad (72)
$$

$$
(2) \text{ Thermal field} \qquad \sum_{i=1}^{2} \left( \delta\Theta_A^{(i)} \left[ \Gamma^{AB}\Theta_B^{(i)} - \int_{\Theta_B} K^{AB}\mathrm{d}V\Theta_B^{(i)} - Q^{(i)A,\mathrm{ext}} \right] - \delta\Theta_A^{(1)} \left\{ \gamma_1 \lambda_{T,B} \cdot (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) [m^{ABC} \dot{\mathbf{q}}_C^{(1),S} - m^{ABD} \dot{\mathbf{q}}_D^{(2),S}] - k_\theta |\lambda_{N,B}| [m^{ABC} \Theta_C^{(1)} - m^{ABD} \Theta_D^{(2)}] \right\} \right. \\
\qquad \left. - \delta\Theta_A^{(2)} \left\{ \gamma_2 \lambda_{T,B} \cdot (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) [\bar{m}^{ABC} \dot{\mathbf{q}}_C^{(1),S} - \bar{m}^{ABD} \dot{\mathbf{q}}_D^{(2),S}] + k_\theta |\lambda_{N,B}| [\bar{m}^{ABC} \Theta_C^{(1)} - \bar{m}^{ABD} \Theta_D^{(2)}] \right\} = 0 \qquad (73)
$$

$$
(3) \text{ Semi-discrete contact conditions}
$$
$$
\begin{align}
& \bullet \text{ Normal contact} \tag{74} \\
& \phantom{\bullet \text{ Normal contact}} g_N^A \le 0, \quad \lambda_{N,A} \ge 0, \quad \lambda_{N,A} g_N^A = 0, \quad g_N^A = \mathbf{n} \cdot [\mathbf{n}^{AB} \mathbf{q}_B^{(1)} - \mathbf{n}^{AC} \mathbf{q}_C^{(2)}] \nonumber
\end{align}
$$

$$
\bullet \text{ Tangential contact} \qquad (75) \\
\bar{\phi}_{c,A} := \|\lambda_{T,A}\| - \mu|\lambda_{N,A}| \le 0, \quad \dot{\zeta}^A \ge 0, \quad \bar{\phi}_{c,A} \dot{\zeta}^A = 0, \quad \dot{\mathbf{g}}_T^S A = \dot{\zeta}^A \frac{\lambda_{T,A}}{\|\lambda_{T,A}\|}
$$

$$
\text{with} \tag{76}
$$
$$
\dot{\mathbf{g}}_T^S A = (\mathbf{I} - \mathbf{n} \otimes \mathbf{n}) [\eta^{AB} \dot{\mathbf{q}}_B^{(1),S} - \eta^{AC} \dot{\mathbf{q}}_C^{(2),S}]
$$

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Virtual segmentation surface (yellow) on an arbitrary curved NURBS geometry (blue). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)</figcaption></figure>
virtual discretisation of contact boundary of $B^{(2)}$, we execute the standard segmentation procedure, see Fig. 1 for a virtual segmentation surface on an arbitrary curved NURBS surface.

**Remark 2.** We use this simplification only for the segmentation process, i.e. for the calculation of the local coordinates $\xi_K^{(i)}$ and not for the evaluation of the Mortar integrals. The speed up is tremendous, furthermore, the definition of the segment corners is not as strict as the definition of the closest point projection used within the KTS method, see Temizer [38]. More importantly, the proposed method yields a set of valid, non-overlapping segments, we only shift small area contributions from one segment to the next.

**Remark 3.** Due to the presence of potentially higher order NURBS, the chosen Gauss quadrature has to be adjusted carefully. A standard 4 point integration is not sufficient for e.g. quadratic NURBS, especially for the evaluation of the triple Mortar integrals.

**Remark 4.** The variation of the normal gap $\delta g_N$ for an initially curved surface may include additionally the variation of the normal vector, which is dropped for convenience in Eq. (65) and the following. This is sometimes done in computational contact mechanics based on the mortar method (see e.g. Laursen et al. [19]). In contrast to NTS and KTS methods, where the variation of the normal vector drops out (see e.g. Wriggers [40]), the normal components for Mortar methods are assembled from all involved segments, such that the corresponding variation does not vanish in all situations. Leaving the variation out of the balance equation may have considerable impact to the angular momentum at the contact interface. This has already been observed in Puso and Laursen [29, footnote p. 606] and Puso and Laursen [30, footNote 1 and 2 p. 4896].

## 4. Temporal discretisation

In a final step, the semi-discrete equations of motion have to be discretised in time in order to obtain a set of non-linear algebraic equations to be solved via a Newton method. Consider a sequence of times $t_0, \dots, t_n, t_{n+1}, \dots, T$, where $(\bullet)_n$ and $(\bullet)_{n+1}$ denote the value of a given physical quantity at time $t_n$ and $t_{n+1}$, respectively. Assume the state variables at $t_n$, given by $(\mathbf{q}_{A,n}, \Theta_{A,n})$ are known and the time step size $\Delta t = t_{n+1} - t_n$ is given. Then, the algorithmic approximation to the weak form (55) is defined by

$$ \delta \mathbf{q}_A \cdot \left[ \frac{2}{\Delta t^2} M^{AB} (\mathbf{q}_{B,n+1} - \mathbf{q}_{B,n} - \Delta t \mathbf{v}_{B,n}) + \int_{B_0} \nabla_x (R^A) \cdot \mathbf{S}_{n,n+1}^h \nabla_x (R^B) \, dV \mathbf{q}_B \right] = \delta \mathbf{q}_A \cdot \left[ \mathbf{F}_{n+1/2}^{A, \text{ext}} \right], \quad (77) $$

whereas the approximation of (56) is given by

$$ \delta \Theta_A \left[ \Gamma_{n,n+1}^{AB} \Theta_{B,n+1/2} - \int_{B_0} K_{n,n+1}^{AB} dV \Theta_{B,n+1/2} \right] = \delta \Theta_A \left[ Q_{n+1/2}^{A, \text{ext}} \right]. \quad (78) $$

Here, $\Gamma_{n,n+1}^{AB}$ and $K_{n,n+1}^{AB}$ are standard mid-point approximations of (60) and (61), cf. Hesch & Betsch [12] for details. Moreover, the algorithmic version of the second Piola–Kirchhoff stress tensor reads

$$ \mathbf{S}_{n,n+1}^{h} = 2 \frac{\partial \Psi(\mathbf{C}_{n+1/2}^{h}, \Theta_{n+1/2}^{h})}{\partial \mathbf{C}^{h}}, \quad \mathbf{C}_{n+1/2}^{h} = \frac{1}{2} (\mathbf{C}_{n+1}^{h} + \mathbf{C}_{n}^{h}). \quad (79) $$

The used evaluation of the right Cauchy–Green tensor $\mathbf{C}_{n+1/2}^{h}$ is equivalent to the energy–momentum method proposed by Simo and Tarnow [34], if applied to a St.Venant–Kirchhoff material. This approach is easy to implement and remarkably stable for multi-field problems.

### Normal contact conditions

Following the approach outlined in Hesch & Betsch [13], we evaluate the mortar integrals at time $t_n$, whereas the La-grange multipliers $\lambda_{N,A,n+1}$ remain constant throughout the time step. The Karush–Kuhn–Tucker conditions (74) are implemented via an active set strategy and the weak form reads

$$ G_{c,N}^{\phi,h} = \lambda_{N,A,n+1} \mathbf{n}_{n+1/2} \cdot \left[ \mathbf{n}_{n}^{AB} \delta \mathbf{q}_{B}^{(1)} - \mathbf{n}_{n}^{AC} \delta \mathbf{q}_{C}^{(2)} \right]. \quad (80) $$

Note that we evaluate the corresponding constraints at time $t_{n+1}$.

### Tangential contact conditions

The tangential contributions to the weak form follow analogously

$$ G_{c,T}^{\phi,h} = \lambda_{T,A,n+1/2} \cdot (\mathbf{I} - \mathbf{n}_{n+1/2} \otimes \mathbf{n}_{n+1/2}) \left[ \mathbf{n}_{n}^{AB} \delta \mathbf{q}_{B}^{(1)} - \mathbf{n}_{n}^{AC} \delta \mathbf{q}_{C}^{(2)} \right]. \quad (81) $$
In particular, we employ a trial state-return map strategy to determine the Coulomb frictional traction. The trial state is given by

$$ \lambda_{T,A,n+1}^{\text{trial}} = \lambda_{T,A,n} + \epsilon_T (\mathbf{I} - \mathbf{n}_{n+1/2} \otimes \mathbf{n}_{n+1/2}) [\mathbf{n}_n^{AB} \Delta \mathbf{q}_B^{(1)} - \mathbf{n}_n^{AC} \Delta \mathbf{q}_C^{(2)}]. \quad (82) $$

Here, $\Delta \mathbf{q}_B^{(1)} = \mathbf{q}_{B,n+1}^{(1)} - \mathbf{q}_{B,n}^{(1)}$ and $\Delta \mathbf{q}_C^{(2)} = \mathbf{q}_{C,n+1}^{(2)} - \mathbf{q}_{C,n}^{(2)}$. To prevent large errors in angular momentum, a modification of the form

$$ \lambda_{T,A,n+1}^{\text{trial}} = \lambda_{T,A,n} + \epsilon_T (\mathbf{I} - \mathbf{n}_{n+1/2} \otimes \mathbf{n}_{n+1/2}) \left[ (\mathbf{n}_{n+1}^{AB} - \mathbf{n}_n^{AB}) \mathbf{q}_{B,n+1/2}^{(1)} - (\mathbf{n}_{n+1}^{AC} - \mathbf{n}_n^{AC}) \mathbf{q}_{C,n+1/2}^{(2)} \right] \quad (83) $$

is often used. For both approaches, the slip function can be obtained for each node A separately

$$ \hat{\phi}_{c,A,n+1} = \| \lambda_{T,A,n+1}^{\text{trial}} \| - \mu |\lambda_{N,A,n+1}| \quad (84) $$

and we obtain the frictional tractions

$$ \lambda_{T,A,n+1} = \begin{cases} \lambda_{T,A,n+1}^{\text{trial}}, & \text{if } \hat{\phi}_{c,A,n+1} \le 0, \\ \mu |\lambda_{N,A,n+1}| \frac{\lambda_{T,A,n+1}^{\text{trial}}}{\|\lambda_{T,A,n+1}^{\text{trial}}\|}, & \text{elseif } \hat{\phi}_{c,A,n+1} > 0. \end{cases} \quad (85) $$

Once the frictional tractions at time $t_{n+1}$ are calculated, the corresponding midpoint approximation reads

$$ \lambda_{T,A,n+1/2} = \frac{1}{2} (\lambda_{T,A,n+1} + \lambda_{T,A,n}), \quad (86) $$

where we follow the arguments outlined in Armero & Petöcz [1].

### Thermal contact conditions

Eventually, the thermal contributions to the weak form are given by

$$ G_c^{\theta,h} = -\delta\Theta_A^{(1)} \left\{ \begin{aligned} & \gamma_1 \lambda_{T,B,n+1} \cdot (\mathbf{I} - \mathbf{n}_{n+1/2} \otimes \mathbf{n}_{n+1/2}) \left[ m_n^{ABC} \frac{\Delta\mathbf{q}_C^{(1),s}}{\Delta t} - \bar{m}_n^{ABD} \frac{\Delta\mathbf{q}_D^{(2),s}}{\Delta t} \right] - k_0 |\lambda_{N,B,n+1}| \left[ m_n^{ABC} \Theta_{C,n+1/2}^{(1)} - \bar{m}_n^{ABD} \Theta_{D,n+1/2}^{(2)} \right] \\ & - \delta\Theta_A^{(2)} \left\{ \gamma_2 \lambda_{T,B,n+1} \cdot (\mathbf{I} - \mathbf{n}_{n+1/2} \otimes \mathbf{n}_{n+1/2}) \left[ \bar{m}_n^{ABC} \frac{\Delta\mathbf{q}_C^{(1),s}}{\Delta t} - \bar{m}_n^{ABD} \frac{\Delta\mathbf{q}_D^{(2),s}}{\Delta t} \right] + k_0 |\lambda_{N,B,n+1}| \left[ \bar{m}_n^{ABC} \Theta_{C,n+1/2}^{(1)} - \bar{m}_n^{ABD} \Theta_{D,n+1/2}^{(2)} \right] \end{aligned} \right\}. \quad (87) $$

assuming that (82) is used. The modified form in (83) can be applied using $\frac{\Delta m^{ABC}}{\Delta t}$ and $\frac{\Delta m^{ABD}}{\Delta t}$ instead of $\frac{\Delta \mathbf{q}_C^{(1),s}}{\Delta t}$ and $\frac{\Delta \mathbf{q}_D^{(2),s}}{\Delta t}$ within the last expression.

**Remark.** The proposed second order accurate time integration scheme with stabilisation via a modified evaluation of the second Piola-Kirchhoff stress tensor can be replaced by a more simple implicit Euler approach, to which we refer to as endpoint rule. Therefore, we evaluate all corresponding terms at $t_{n+1}$, which will be used in the subsequent section for comparison.

## 5. Numerical examples

To demonstrate the accuracy and performance of the proposed methodology we investigate several quasi-static and transient numerical examples. In particular, we investigate first a modified patch test to outline the capabilities of the thermomechanical Mortar method for NURBS. In the second example, we focus our investigation on the frictional response, whereas our primary interest in the third example relies on the thermal response, although all example utilise the fully coupled Mortar method. Finally, a transient impact simulation demonstrates the accuracy of the frictional Mortar contact algorithms in terms of the first and second law of thermodynamics.

A non-linear and fully coupled Ogden material model with the associated strain energy density function

$$ \Psi(\lambda_1, \lambda_2, \lambda_3, \theta) = \sum_{A=1}^{3} \sum_{p=1}^{3} \frac{\mu_p(\theta_0) \frac{\theta}{\theta_0} (\bar{\lambda}_A^{\alpha_p} - 1) + \kappa(\theta_0) \frac{\theta}{\theta_0} \beta^{-2} (\beta \ln(J) + J^{-\beta} - 1) - 3\alpha_0 \kappa(\theta_0) \gamma^{-1} (J^{-\gamma} - 1) (\theta - \theta_0)}{c_0(\theta - \theta_0 - \theta \ln(\theta/\theta_0))}, \quad (88) $$

(see Table 3 for the material data and Holzapfel and Simo [15] for details on the function), as well as a thermomechanical Neo-Hookean material model, defined as

$$ \Psi(\mathbf{C}, \theta) = \frac{\bar{\mu}}{2}(\text{tr}(\mathbf{C}) - 3 - 2 \ln(J)) + \frac{\lambda}{2}[\ln^2(J) + (J-1)^2] - 3\alpha_0\lambda(\theta - \theta_0)(J^{-1}\ln(J) + J^{-1}) + c_0(\theta - \theta_0 - \theta \ln(\theta/\theta_0)), \quad (89) $$
<table><caption>**Table 3**
Material properties.</caption><tr><td>Ogden model</td><td>μ<sub>1</sub> = 6.30 · 10<sup>5</sup><br/>μ<sub>2</sub> = 0.012 · 10<sup>5</sup><br/>μ<sub>3</sub> = -0.10 · 10<sup>5</sup></td><td>α<sub>1</sub> = 1.3<br/>α<sub>2</sub> = 5.0<br/>α<sub>3</sub> = -2.0</td></tr><tr><td>Heat capacity</td><td>c<sub>0</sub> = 1830</td><td></td></tr><tr><td>Density</td><td>ρ<sub>0</sub> = 950</td><td></td></tr><tr><td>Linear expansion coefficient</td><td>α<sub>0</sub> = 22.333 · 10<sup>-5</sup></td><td></td></tr><tr><td>Bulk modulus</td><td>κ(θ<sub>0</sub>) = 2.0 · 10<sup>8</sup></td><td></td></tr><tr><td>Empirical coefficients</td><td>β = 9.0, γ = 2.50</td><td></td></tr><tr><td>Thermal conductivity</td><td>K<sub>0</sub>(θ<sub>0</sub>) = 0.15</td><td></td></tr><tr><td>Softening parameter</td><td>ω<sub>K</sub> = 0.004</td><td></td></tr></table>

are used in the examples. The strain energy density function in (88) is written in terms of principal stretches $\lambda_A$, where $\tilde{\lambda}_A = J^{-1/3}\lambda_A$ with $J = \lambda_1\lambda_2\tilde{\lambda}_3$. Furthermore, the thermal conductivity tensor within the Duhamel's law (5) is set to

$$ \dot{\mathbf{K}}(\mathbf{C}, \theta) = K_0(\theta)\mathbf{C}^{-1}, \qquad (90) $$

with

$$ K_0(\theta) = K_0(\theta_0)[1 - \omega_k(\theta - \theta_0)]. \qquad (91) $$

## 5.1. Patch test

First, we investigate a modified patch test with flat as well as curved interfaces using the Ogden material model. As shown in Fig. 2, two independently meshed blocks are tied together via the proposed Mortar method, assuming that both bodies remain in contact without tangential sliding and without thermal resistance at the interface. A Neumann boundary is applied to the top surface, whereas the bottom surface is clamped, such that the body can expand in tangential direction. The upper block consists of 3 × 3 × 3 quadratic NURBS elements while the lower block comprises of 8 elements, using quadratic NURBS as well. Thus, altogether 189 nodes with 756 unknowns controls the discrete geometry as well as the thermal fields To obtain a quasi-static solution, the density is set to zero, i.e. only $\Gamma^{AB}$ within the thermal field depends on the differentiation with respect to time. Fig. 3 (left) shows the von Mises stress distribution of the solution for a uniform pressure field of $\sigma = 10^6$, applied to the upper block. The initial temperature is set to 293.15. Due to the applied linear expansion coefficient, the body heats up, as shown in Fig. 3 (right). The solution shows a nearly perfect and uniform temperature distribution for the proposed linear approximation of the Lagrange multiplier field.

For further investigations, a highly curved interface is applied. Note that quadratic NURBS are able to represent the geometry exactly, independent of the chosen number of elements, such that both interfaces are geometrically conform and errors we obtain are discretisation errors of the Mortar interface. Again, a uniform stress field with $\sigma = 2.5 \times 10^5$ has been applied to the top surface. The error in the stress field decreases and converges to the correct, uniform solution with higher element numbers at the interface, see Fig. 4. As already shown in Hesch and Betsch [11], the results of the Mortar method for the given curved interface are in general superior to results of collocation type methods.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Reference configuration with flat (left) and curved (right) interface.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Stress distribution (left) and temperature distribution (right).</figcaption></figure>

Since we assume that the error in the stress field is of geometrical nature, we investigate the scalar valued thermal field again. Therefore, Neumann boundary condition as well as the thermal resistance within the interface are set to zero. A uniform, stressfree configuration with initial temperature of $θ_0 = 323.15$ is predefined for the upper block, whereas an initial temperature of $θ_0 = 293.15$ is applied to the lower block. The temperature fields at the top and bottom surfaces are fixed. Fig. 5 (left) shows the solution after a single time step with $Δt = 0.05$. Due to the given heat capacity and the thermal conductivity, heat flows across the interface, until the solution converges to the static solution, as shown in Fig. 5 (right) after a total simulation time of $t = 10^4$.

## 5.2. Rotating contact of two elastic bodies

The next example consists of an intender and a block. For both bodies the material behavior is assumed to be governed by the afore introduced Neo-Hookean model with Lamé parameters corresponding to a Young's Moduls of $E = 10$ and a Poisson's ratio of $v = 0.3$. Furthermore, the heat capacity is $c_0 = 300$ and the expansion coefficient $α_0 = 0.05$ and the reference temperatures of both bodies are set to 293.15. For the contact interface, the friction coefficient is set to $μ = 0.5$, the efficiencies to $γ_1 = γ_2 = 0.5$, the tangential penalty parameter $ε_T = 2 × 10^3$ and the constant heat exchange parameter to $k_θ = 0.2$.

The system itself remains quasistatic, i.e. the density of the bodies in contact is again set to zero. The initial configuration consists of $8 × 8 × 8$ elements for the upper and $9 × 9 × 9$ elements for the lower block, see Fig. 6. The bottom surface of the lower block is fixed in space. The top surface of the intender (upper block) moves downwards and starts afterwards with a

<figure><img src="image_5.png" /><figcaption>Fig. 4. Stress distribution for different discretisations.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Thermal distribution after the first step (left) and converged quasi-static solution (right).</figcaption></figure>

rotational movement. The problem is similar to the one shown in Temizer et al. [39,37]. In Fig. 7 the deformed configuration after the first phase of the movement is shown along with the segmentation of both surfaces.

Overall 4 quasi-static time steps of size 0.05 for the first phase and 40 time steps of size 0.02 have been used for the second phase. In Fig. 8 the side and top view after the rotational movement with a frictional coefficient of $\mu = 0.5$ is shown.

Next, different friction coefficients with $\mu \in [0.25, 0.5, 1]$ are applied and the resulting twisting angular momentum acting on the upper, twisted body is plotted over time, see Fig. 9.

At last in Fig. 10 a cut through both blocks after the rotational movement with $\mu = 0.5$ is shown. Note that the temperature distribution depends in this particular example mainly on the linear expansion coefficient. The different contributions will be shown in detail in the subsequent example.

### 5.3. Ironing

In this third example, we focus on the thermal contributions at the interface, using a setup similar to the ironing problem described in Puso and Laursen [30]. The geometry in its reference configuration is as shown in Fig. 11, where the upper block of size $1 \times 1 \times 1$ consists of 27 quadratic NURBS elements and the lower block of size $9 \times 4 \times 3$ consists of 1680 quadratic

<figure><img src="image_7.png" /><figcaption>Fig. 6. Reference configuration of the rotating problem.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Deformed configuration and segmentation.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Final deformed configuration, side and top view.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Twisting momentum plotted over time for different friction coefficients.</figcaption></figure>

NURBS elements. The material behaviour is again assumed to be governed by a Neo-Hookean model, using $\bar{\mu} = 5000/13$, $\lambda = 7500/13$ as Lamé parameters for the upper block and $\bar{\mu} = 50/13$, $\lambda = 75/13$ for the lower block. The values correspond to a Young's modulus of $E = 10$ and $E = 1000$, respectively and to a Poisson's ratio of $\nu = 0.3$ for both blocks.

The setup of the thermal field is given as follows: The initial temperature of the upper block is set to $\theta_0 = 298.15$ and the heat capacity to $c_0 = 1.83 \cdot 10^4$ in order to avoid a fast cool-down, while the corresponding settings for the lower block are given by $\theta_0 = 293.15$ and $c_0 = 1.83 \cdot 10^3$. The thermal conductivity $K_0 = 0.55$ and the expansion coefficient $\alpha_0 = 22.333 \cdot 10^{-5}$ are set for both blocks.

For the contact interface, the effusivities are $\gamma_1 = \gamma_2 = 0.5$, the heat exchange coefficient is $k_\theta = 3000$ and the friction coefficient is $\mu = 0.2$. Accordingly, frictional dissipation is taken into account.

The lower surface of the lower block is fixed, whereas the movement of the upper surface of the upper block is predefined. In particular, the block is first moved downwards and pressed into the lower block and then moved sidewards, such that it slides over the upper surface of the lower block. In Fig. 12, the temperature distribution at $t = 12$ is shown. On the right, a cut at the symmetry plane of both blocks is shown along with the marking of a specific node. The temperature at this marked position on the mesh is plotted over time in Fig. 13.
<figure><img src="image_11.png" /><figcaption>Fig. 10. Temperature distribution after the rotational movement.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Reference configuration of the ironing example.</figcaption></figure>

Due to the frictional contact forces, a small bow wave is running at the front of the upper block with a compression and an expansion part. In connection with the expansion coefficient, this leads to the small temperature wave between $t = 8$ and $t = 10$. Afterwards, the marked node on the mesh gets in contact with the upper block, which leads to a drastic temperature increase. Note that this example is quite difficult, since we obtain high stresses at the corners of the block and, furthermore, high curvatures for the given coarse mesh.

In the final setup for this problem, we focus on the frictional contributions. Therefore, the thermal expansion coefficient is set to zero and the initial temperature at 293.15 is assumed to be equal for both bodies, i.e. the only heat source is the frictional dissipation, see Fig. 14 for the resulting temperature distribution.

## 5.4. Impact problem

This last example deals with a fully coupled, transient thermoelastic impact problem. The reference configuration is shown in Fig. 15. Quadratic NURBS are used for the discretisation of both bodies in contact. The lower, larger body of size $2 \times 2 \times 1$ consists of 125 elements, whereas the upper, smaller body of size $1 \times 1 \times 1$ consists of 27 elements. The previously introduced Ogden model has been modified as follows: The density of the smaller body is set to 100 such that the body re-bound from the lower body whereas the initial velocity is set to $[3, 3, -2]$. Furthermore, the linear expansion coefficient has been set to zero for both bodies, such that we can investigate effects of the thermomechanical contact interface in detail.

The friction coefficient is set to $\mu = 0.5$, the effusivities to $\gamma_1 = \gamma_2 = 0.5$, the tangential penalty parameter $\epsilon_T = 10^8$ and the constant heat exchange parameter to $k_0 = 0.2$. The initial temperature of the lower block is 293.15, whereas the
<figure><img src="image_13.png" /><figcaption>Fig. 12. Temperature distribution at $t = 12$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Temperature at the marked position over time.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Temperature distribution at time $t = 9$ for frictional sliding.</figcaption></figure>

temperature of the upper block is 343.15. In Fig. 16 the configurations at time $t = 0.05$ and $t = 0.2$ are depicted. Due to the frictional impact, the upper body starts to rotate.

The von Mises stresses and temperature distribution at the bottom surface of the upper block just after the initial impact at time $t = 0.02$ is given in Fig. 17. Note that the upper block cools down at its bottom surface due to the contact
<figure><img src="image_16.png" /><figcaption>Fig. 15. Reference configuration of the thermoelastic impact problem.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Configuration at time $t = 0.05$ (left) and at time $t = 0.2$ (right).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. Von Mises stresses (left) and temperature (right) distribution.</figcaption></figure>

with the lower body, which has a significantly lower initial temperature. This effect exceeds the amount of heat generated by the frictional dissipation. Furthermore, the stress field has its maximum values at the front of the contact surface, which is correctly represented in the temperature distribution, since the heat exchange across the interface
<figure><img src="image_19.png" /><figcaption>Fig. 18. Different energies for midpoint and endpoint rule plotted over time.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. Changes in total energy for the midpoint and endpoint rule plotted over time.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. Changes in the first component in total angular momentum plotted over time.</figcaption></figure>

depends directly on the contact forces, see (35). Within the initial impact phase after the first few time steps both surfaces are geometrically in a perfect flat contact situation, such that for a frictionless scheme the stress distribution would be uniform.
<figure><img src="image_22.png" /><figcaption>Fig. 21. Convergence plot, average of the absolute change in total energy plotted over the time step size.</figcaption></figure>

Next, the inner energy, the mechanical energy (i.e. the kinetic plus the Helmholtz energy) as well as the total energy (i.e. the kinetic plus the inner energy) are plotted over time for both, the proposed midpoint type scheme as well as an endpoint scheme, see Fig. 18. In contrast to the endpoint rule, the proposed midpoint type scheme conserves the total energy very well. For details on a logarithmic scale, see Fig. 19, where the black line represents the Newton tolerance. The difference between the mechanical and the total energy represents the dissipated energy, transferred to the thermal field. Furthermore, as requested by the second law of thermodynamics, dissipation is always non-negative.

The changes of the first component in total angular momentum are plotted over time in Fig. 20. As shown, the midpoint type approximation conserves angular momentum algorithmically, except for the impact phase. During impact, the changes in angular momentum are considerably smaller compared with the implicit Euler scheme. As before, the black line represents the Newton tolerance. Additionally, results for the proposed midpoint evaluation without the variation of the normal vector (see Remark 3 in Section 3) are plotted, labelled as Midpoint 2 in Fig. 20. As can be seen, the differences are marginal.

A constant time step size of $\Delta t = 0.0005$ has been applied for the shown diagrams. In Fig. 21 a convergence plot for different time step sizes is depicted. Accordingly, the average of the absolute change in total energy has been plotted over the time step size. As can be seen, the endpoint rule is first, the proposed midpoint approximation second order accurate.

## 6. Conclusions

The proposed linear approximation of the Lagrange multiplier field delivers a flexible framework for the construction of Mortar contact algorithms for Isogeometric Analysis (IGA) with NURBS. It was demonstrated that this flexibility is preserved even for complex frictional impact problems. The IGA framework additionally benefits from a simplified Mortar segmentation procedure, leading to significantly higher computational efficiency while maintaining satisfactory accuracy for contact problems.

The extension to thermomechanical contact problems follows in a straightforward manner through the introduction of triple Mortar integrals. This novel formulation allows for modeling the energy transfer between the contact surfaces and helps incorporate the effect of the frictional dissipation on the thermal field, leading to a fully coupled thermomechanically consistent frictional Mortar contact formulation for NURBS of arbitrary order. In this formulation, the IGA surface geometry is endowed with a contact discretisation through an interface formulation which constitutes a link to established numerical techniques for Mortar methods. The flexibility of this formulation renders it a convenient starting point for the construction of similar thermomechanical Mortar contact algorithms for advanced IGA discretisations based on hierarchical refinement [33] and T-splines [7].

## Acknowledgements

Support for this research was provided by the Deutsche Forschungsgemeinschaft (DFG) under Grant HE 5943/3-1. This support is gratefully acknowledged.

## References

[1] F. Armero, E. Petöcz, A new dissipative time-stepping algorithm for frictional contact problems: formulation and analysis, Comput. Methods Appl. Mech. Eng. 179 (1999) 159–178.

[2] M. Bahrami, J.R. Culham, M.M. Yovanovich, G.E. Schneider, Review of thermal joint resistance models for nonconforming rough surfaces, Appl. Mech. Rev. 59 (2006) 1–12.
[3] M. Bahrami, M.M. Yovanovich, E.E. Marotta, Thermal joint resistance of polymer-metal rough interfaces, J. Electron. Packaging 128 (2006) 23-29.

[4] J.A. Cottrell, T.J.R. Hughes, Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, Wiley, New York, 2009.

[5] L.B. da Veiga, D. Cho, L.F. Pavarino, S. Scacchi, Overlapping Schwarz methods for isogeometric analysis, SIAM J. Numer. Anal. 50 (2012) 1394-1416.

[6] C. Agelet de Saracibar, Numerical analysis of coupled thermomechanical frictional contact problems. Computational model and applications, Arch. Computat. Methods Eng. 5 (1999) 243-301.

[7] R. Dimitri, L.D. Lorenzis, M.A. Scott, P. Wriggers, R.L. Taylor, G. Zavarise, Isogeometric large deformation frictionless contact using T-splines, Comput. Methods Appl. Mech. Eng. 269 (2014) 394-414.

[8] M. Franke, C. Hesch, P. Betsch, An augmentation technique for large deformation frictional contact problems, Int. J. Numer. Methods Eng. 94 (2013) 513-534.

[9] J.J. Fuller, E.E. Marotta, Thermal contact conductance of metal/polymer joints: an analytical and experimental investigation, J. Thermophys. Heat Transf. 15 (2001) 228-238.

[10] M. Gross, P. Betsch, Energy-momentum consistent finite element discretization of dynamic finite viscoelasticity, Int. J. Numer. Methods Eng. 81 (2010) 1341-1386.

[11] C. Hesch, P. Betsch, Transient 3D domain decomposition problems: frame-indifferent mortar constraints and conserving integration, Int. J. Numer. Methods Eng. 82 (2010) 329-358.

[12] C. Hesch, P. Betsch, Energy-momentum consistent algorithms for dynamic thermomechanical problems – application to mortar domain decomposition problems, Int. J. Numer. Methods Eng. 86 (2011) 1277-1302.

[13] C. Hesch, P. Betsch, Transient three-dimensional contact problems: mortar method. Mixed methods and conserving integration, Comput. Mech. 48 (2011) 461-475.

[14] C. Hesch, P. Betsch, Isogeometric analysis and domain decomposition methods, Comput. Methods Appl. Mech. Eng. 213 (2012) 104-112.

[15] G.A. Holzapfel, J.C. Simo, Entropy elasticity of isotropic rubber-like solids at finite strains, Comput. Methods Appl. Mech. Eng. 132 (1996) 17-44.

[16] S. Hüeber, B.I. Wohlmuth, Thermo-mechanical contact problems on non-matching meshes, Comput. Methods Appl. Mech. Eng. 198 (2009) 1338-1350.

[17] J.-Y. Kim, S.-K. Youn, Isogeometric contact analysis using mortar method, Int. J. Numer. Methods Eng. 89 (2012) 1559-1581.

[18] T.A. Laursen, On the development of thermodynamically consistent algorithms for thermomechanical frictional contact, Comput. Methods Appl. Mech. Eng. 177 (1999) 273-287.

[19] T.A. Laursen, M.A. Puso, J.D. Sanders, Mortar contact formulations for deformable-deformable contact: past contributions and new extensions for enriched and embedded interface formulations, Comput. Methods Appl. Mech. Eng. 205-208 (2012) 3-15.

[20] L.D. Lorenzis, I. Temizer, P. Wriggers, G. Zavarise, A large deformation frictional contact formulation using NURBS-based isogeometric analysis, Int. J. Numer. Methods Eng. 87 (2011) 1278-1300.

[21] L.D. Lorenzis, P. Wriggers, G. Zavarise, A mortar formulation for 3D large deformation contact using NURBS-based isogeometric analysis and the augmented Lagrangian method, Comput. Mech. 49 (2012) 1-20.

[22] J. Lu, Isogeometric contact analysis: geometric basis and formulation of frictionless contact, Comput. Methods Appl. Mech. Eng. 200 (2011) 726-741.

[23] C.V. Madhusudana, Thermal Contact Conductance, Springer, New York, 1996.

[24] M.E. Matzen, T. Cichosz, M. Bischoff, A point to segment contact formulation for isogeometric, nurbs based finite elements, Comput. Methods Appl. Mech. Eng. 255 (2013) 27-39.

[25] C. Miehe, Entropic thermoelasticity at finite strains. Aspects of the formulation and numerical implementation, Comput. Methods Appl. Mech. Eng. 120 (1995) 243-269.

[26] V.G. Oancea, T.A. Laursen, A finite element formulation of thermomechanical rate-dependent frictional sliding, Int. J. Numer. Methods Eng. 40 (1997) 4275-4311.

[27] B.N.J. Persson, B. Lorenz, A.I. Volokitin, Heat transfer between elastic solids with randomly rough surfaces, Eur. Phys. J. E 31 (2010) 3-24.

[28] L. Piegl, W. Tiller, The NURBS Book, second ed., Springer, 2010.

[29] M.A. Puso, T.A. Laursen, A mortar segment-to-segment contact method for large deformation solid mechanics, Comput. Methods Appl. Mech. Eng. 193 (2004) 601-629.

[30] M.A. Puso, T.A. Laursen, A mortar segment-to-segment frictional contact method for large deformations, Comput. Methods Appl. Mech. Eng. 193 (2004) 4891-4913.

[31] M.A. Puso, T.A. Laursen, J.M. Solberg, A segment-to-segment mortar contact method for quadratic elements and large deformations, Comput. Methods Appl. Mech. Eng. 197 (2008) 555-566.

[32] S. Reese, S. Govindjee, Theoretical and numerical aspects in the thermo-viscoelastic material behaviour of rubber-like polymers, Mech. Time-Dependent Mater. 1 (1998) 357-396.

[33] D. Schillinger, L. Dedè, M.A. Scott, J.A. Evans, M.J. Borden, E. Rank, T.J.R. Hughes, An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces, Comput. Methods Appl. Mech. Eng. 249-252 (2012) 116-150. Higher Order Finite Element and Isogeometric Methods.

[34] J.C. Simo, N. Tarnow, The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics, Z. angew. Math. Phys. (ZAMP) 43 (1992) 757-792.

[35] N. Strömberg, L. Johansson, A. Klarbring, Derivation and analysis of a generalized standard model for contact, friction and wear, Int. J. Solids Struct. 33 (1996) 1817-1836.

[36] I. Temizer, A mixed formulation of mortar-based frictionless contact, Comput. Methods Appl. Mech. Eng. 223-224 (2012) 173-185.

[37] I. Temizer, A mixed formulation of mortar-based contact with friction, Comput. Methods Appl. Mech. Eng. 255 (2013) 183-195.

[38] I. Temizer, P. Wriggers, T.J.R. Hughes, Contact treatment in isogeometric analysis with NURBS, Comput. Methods Appl. Mech. Eng. 200 (2011) 1100-1112.

[39] I. Temizer, P. Wriggers, T.J.R. Hughes, Three-dimensional mortar-based frictional contact treatment in isogeometric analysis with NURBS, Comput. Methods Appl. Mech. Eng. 209-212 (2012) 115-128.

[40] P. Wriggers, Computational Contact Mechanics, second ed., Springer-Verlag, 2006.

[41] B. Yang, T.A. Laursen, A large deformation mortar formulation of self contact with finite sliding, Comput. Methods Appl. Mech. Eng. 197 (2008) 756-772.

[42] G. Zavarise, P. Wriggers, E. Stein, B.A. Schrefler, Real contact mechanisms and finite element formulation – a coupled thermomechanical approach, Int. J. Numer. Methods Eng. 35 (1992) 767-785.