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

HAL
open science

# Finite element treatment of two-dimensional thermoelastic wear problems

Niclas Strömberg

► To cite this version:

Niclas Strömberg. Finite element treatment of two-dimensional thermoelastic wear problems. Computer Methods in Applied Mechanics and Engineering, 1999, 177 (3-4), pp.441-455. 10.1016/S0045-7825(98)00392-2. hal-01574166

HAL Id: hal-01574166

https://hal.science/hal-01574166

Submitted on 17 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Finite element treatment of two-dimensional thermoelastic wear
problems

Niclas Strömberg

Division of Mechanics, Department of Mechanical Engineering, Linköping Institute of Technology, S-581 83 Linköping, Sweden

The present paper concerns the numerical treatment of thermoelastic wear problems. The governing equations of thermoelasticity coupled to Signorini contact, Coulomb's friction and Archard's wear are formulated as a system of discrete equations. This equation system is solved, using a Bouligand differentiable Newton method, for five problems of didactic nature.

1. Introduction

In the present paper, a Newton method for thermoelastic wear problems is suggested. The continuum thermodynamical wear model presented in [20] is considered for a two-dimensional thermoelastic body which is unilaterally constrained by a rigid flat tangential moving support. In addition, it is assumed that no heat can flow through the support. For this system, starting from the governing equations, a system of finite element equations, which can be associated with an augmented Lagrangian formulation (see e.g. [1,11]), is derived. This system of equations is numerically solved for five problems of didactic nature using a Bouligand differentiable Newton method given by Pang [16].

Augmented Lagrangian formulations of two- and three-dimensional elastic wear problems were solved in [21,22] using Pang's Newton method. Furthermore, an extensive comparison of Pang's algorithm with an interior point method and a method commercially available in ABAQUS was performed for two-dimensional friction problems in [4]; numerical tests showed that the Newton method is both effective and robust. In the present paper, thermal properties are added to the augmented Lagrangian formulations, discussed in the papers above, in order to achieve a robust and efficient algorithm for thermoelastic wear problems. The method achieved is a direct method without any splitting between mechanical and thermal subproblems (see e.g. [18]). Numerical tests have proven that the method works viably. However, a crucial detail for making the method effective is to eliminate the temperature by algebraical manipulations, such that the temperature is not explicitly involved in the Newton step, but rather appears as an inner variable in the algorithm. Otherwise, the temperature field will not be solved accurately, which in turn, of course, would affect the quality of the mechanical part of the solution.

During the past years, a certain degree of progress has been made concerning the finite element analysis of thermomechanical contact problems.¹ Zavarise et al. [29] formulated a finite element model for thermal contact conduction; frictional heat generation was briefly discussed using a frictional tangential stiffness term and

¹ Concerning the analysis of thermomechanical contact problems of simple geometries such as half-planes, see [2].
additional empirical coefficients. In [26], frictionless thermal contact conduction problems were considered using a penalty formulation, where the penalty terms depended on micro-phenomenological assumptions. Thermoelastic contact problems including both thermal contact conduction and frictional heat generation were formulated and solved in [9] by splitting the problem in one mechanical part and one thermal part; the mechanical problem was solved using a frictional contact algorithm suggested by Klarbring and Björkman [10]. Wriggers and Miehe [25] treated frictional contact problems including thermal contact conduction and frictional heat generation using a normal compliance law and a radial return approach for the friction. In [28], augmented Lagrangian techniques were used for thermomechanical frictional contact problems described with micro-scopical laws; unfortunately, as it seems (see also [27]), the theoretical contact formulation is for bilateral contact and not for the proper unilateral formulation; in addition, a simplification of Coulomb's law of friction is used and the frictional heat generation is excluded. Recently, finite element formulations of thermomechanical frictional sliding in large deformations were given by Laursen and Oancea [12,15].

In this paper, we focus on the finite element treatment of frictional heat generation. Therefore, thermal contact conduction is not considered; of course, it is straightforward to include heat transfer through the contact interface (see [20]). Actually, wear heat generation is also included in the thermomechanical model presented herein. However, for physically realistic wear coefficients, the wear heat generation can be neglected in comparison to the frictional heat generation. Still, the internal state variable measuring the wear gap has a major influence on the frictional heat generation and vice versa. This fact is illustrated by numerical examples. For instance, it is numerically shown that two frictional hot spots appearing in a frictional analysis are removed if wear also is included in the constitutive settings.

The contents of the present study is as follows: in Section 2, the governing equations for a unilaterally constrained thermoelastic body are given. In Section 3, the governing equations are put together to an initial boundary value problem, which in turn is discretized both in space and time in order to reach a system of finite element equations. The system of equations obtained is considered in Section 4, where the adoption of Pang's Newton algorithm also is discussed. In Section 5, five numerical problems are studied. The problems differ only in the constitutive settings. That is, five different combinations of elasticity, thermoelasticity, Signorini contact, Coulomb's friction and Archard's wear are investigated for a particular initial boundary value problem. Finally, in Section 6, concluding remarks are given.

## 2. Governing equations

Let us consider a two-dimensional body in a quasi-static small displacement setting, occupying a domain $\Omega \in \mathbb{R}^2$ with a piecewise smooth boundary $\partial\Omega$, see Fig. 1. The displacement field of the body is denoted by $u$,

<figure><img src="image_3.png" /><figcaption>Fig. 1. The two-dimensional body considered.</figcaption></figure>
the heat flux vector by $q$, the Cauchy stress by $\sigma$, the infinitesimal strain tensor by $\epsilon$ and the location of a point by $x$. Vectors and tensors are represented by components given in a Cartesian basis ($e_1, e_2$) which is defined in Fig. 1. Moreover, the summation convention is applied to repeated indices and Kronecker's delta $\delta_{ij}$ is used.

The boundary $\partial\Omega$ is divided in three disjoint parts $\Gamma_t$, $\Gamma_u$ and $\Gamma_c$. The body is subjected to prescribed tractions $t$ on $\Gamma_t$ and zero displacements on $\Gamma_u$. It is further assumed that the body is unilaterally constrained on $\Gamma_c$ by a rigid support in the direction $n_c$. No heat transfer can take place through the support. Furthermore, the support is flat ended with a normal $n_s = e_2$ and subjected to a tangential displacement $\xi = \xi e_1$, where $\xi = \xi(t)$ is a function of time $t$. An important consequence of the small displacement assumption is that $n_c \approx -n_s$. This fact is utilized in the work by setting $n_c = -e_2$.

## 2.1. Equilibrium equations and principles of thermodynamics

Neglecting inertial forces and body forces, the equilibrium equations for the system considered above read

$$ \frac{\partial \sigma_{ij}}{\partial x_j} = 0 \quad (i=1,2) \quad \text{in } \Omega, \qquad (1) $$

$$ \sigma_{ij} n_j = t_i \quad (i=1,2) \quad \text{on } \Gamma_t, \qquad (2) $$

$$ \sigma_{12} = p_t, \quad \sigma_{22} = -p_n \quad \text{on } \Gamma_c, \qquad (3) $$

where $n_i$ is the $i$-th component of the outward unit normal on $\Gamma_t$ and $p = p_t e_1 - p_n e_2$ represents the contact traction vector acting on the support. Two work conjugate measures belong to $p$. These are $u_n = -u_2$ and $w_t = u_1 - \xi$. That is, $u_n$ is the normal displacement and $w_t$ is the relative tangential displacement between the contact surface and the support.

The balance of energy of the body and the contact surface are given by

$$ \rho \dot{e} = \sigma_{ij} \dot{e}_{ij} - \frac{\partial q_i}{\partial x_i} \quad \text{in } \Omega, \qquad (4) $$

$$ \dot{E} = p_n \dot{u}_n + p_t \dot{w}_t - q_2 \quad \text{on } \Gamma_c, \qquad (5) $$

where $\rho$ is the mass density of the bulk material, $e$ is the specific energy of the body, $E$ is the density of the internal energy of the contact surface and a superimposed dot represents a time derivative.

Defining the absolute temperature $T$, the specific entropy of the body $s$ and the entropy density of the contact surface $S$, the second principle of thermodynamics is given by the following Clausius-Duhem inequalities:

$$ \rho \dot{\psi} \le \sigma_{ij} \dot{\epsilon}_{ij} - \rho s \dot{T} - \frac{q_i}{T} \frac{\partial T}{\partial x_i} \quad \text{in } \Omega, \qquad (6) $$

$$ \dot{\psi} \le p_n \dot{u}_n + p_t \dot{w}_t - s \dot{T} \quad \text{on } \Gamma_c, \qquad (7) $$

where $\psi = e - sT$ and $\dot{\psi} = E - sT$ are the Helmholtz free energies of the body and the contact interface, respectively.

## 2.2. Constitutive assumptions—body

Established constitutive assumptions of isotropic thermoelasticity are utilized for the body.<fnref n="2" /> That is, defining Lame's elasticity coefficients $\bar{\lambda}$ and $\bar{\mu}$, the specific heat capacity $c$, the thermal dilatation coefficient $\alpha$ and a reference temperature $T_0$, the following free energy is considered for the body:

$$ \psi = \frac{1}{2\rho} (\bar{\lambda}(\epsilon_{ii})^2 + 2\bar{\mu}\epsilon_{ij}\epsilon_{ij} - (6\bar{\lambda} + 4\bar{\mu})\alpha(T - T_0)\epsilon_{ii}) - \frac{c}{2T_0}(T - T_0)^2, \qquad (8) $$

from which the following state laws are derived per definition:

<footnote n="2">² The material given in this section is classical and can be found in textbooks (see e.g. [13,14]).</footnote>
$$
\sigma_{ij} = \rho \frac{\partial \psi}{\partial \epsilon_{ij}} = \tilde{\lambda} \epsilon_{kk} \delta_{ij} + 2 \tilde{\mu} \epsilon_{ij} - (3 \tilde{\lambda} + 2 \tilde{\mu}) \alpha (T - T_0) \delta_{ij} \quad (i, j = 1, 2), \qquad (9)
$$

$$
-s = \frac{\partial \psi}{\partial T} = -\frac{1}{\rho} (3\tilde{\lambda} + 2\tilde{\mu}) \alpha \epsilon_{ii} - \frac{c}{T_0} (T - T_0). \quad (10)
$$

By using $\dot{e} = \dot{\psi} + s\dot{T} + sT$, and (8)-(10) in (2.4) and (2.6), we obtain

$$
0 = \frac{\partial q_i}{\partial x_i} + (3\tilde{\lambda} + 2\tilde{\mu})\alpha T \dot{\epsilon}_{kk} + \frac{\rho c}{T_0} T \dot{T}, \quad (11)
$$

$$
- \frac{q_i}{T} \frac{\partial T}{\partial x_i} \ge 0. \tag{12}
$$

In agreement with this latter inequality, representing the second law of thermodynamics, Fourier's law of heat
diffusion is adopted. That is,

$$
q_i = -k_t \frac{\partial T}{\partial x_i} \quad (i = 1, 2), \tag{13}
$$

where $k_t$ is the thermal conductivity. Inserting (2.13) in the balance of energy represented by (11) yields

$$
k_t \frac{\partial}{\partial x_i} \left( \frac{\partial T}{\partial x_i} \right) = (3\tilde{\lambda} + 2\tilde{\mu}) \alpha T \dot{\epsilon}_{kk} + \frac{\rho c}{T_0} T \dot{T}. \quad (14)
$$

This equation is rewritten, by neglecting $\dot{\epsilon}_{kk}$ and putting $T = T_0$, to read

$$
k_t \frac{\partial}{\partial x_i} \left( \frac{\partial T}{\partial x_i} \right) = \rho c T, \quad (15)
$$

which is the standard heat conduction equation of uncoupled thermoelasticity.

2.3. Constitutive assumptions—contact surface

The constitutive assumptions of the contact surface are given by Signorini's contact conditions, Coulomb's law of friction and Archard's law of wear. These laws were derived from a generalized standard material, by defining a certain internal state variable $\omega$ measuring the wear gap, in [20]. Utilizing this internal state variable approach, the free energy is assumed to be given by

$$
\Psi = I_C(u_n, \omega), \quad C = \{(u_n, \omega) : u_n - \omega - g \le 0\}, \tag{16}
$$

where $I_C$ is the indicator function of the closed convex set $C$ (i.e. $I_C = 0$ if $(u_n, \omega) \in C$, $I_C = +\infty$ if $(u_n, \omega) \notin C$) and $g$ is the initial gap between the body and the support. In addition, the following state laws are defined and derived:

$$
(p_n, -W) \in \partial \Psi(u_n, \omega) = \begin{cases} N_C(u_n, \omega) & \text{if } (u_n, \omega) \in C \\ \emptyset & \text{otherwise,} \end{cases} \tag{17}
$$

$$
-S = \frac{\partial \Psi}{\partial T} = 0, \tag{18}
$$

where $\mathcal{W}$ is the force associated to $\omega$, called the wear driving force. $\partial\mathcal{W}$ stands for the subdifferential of $\mathcal{W}$ and

$$
N_C(u_n, \omega) = \{(p_n, -W') : 0 \ge p_n(u'_n - u_n) - W(\omega' - \omega) \quad \forall (u'_n, \omega') \in C\}
$$

defines the normal cone of C. The contact law given in (17) represents an extension of Signorini's contact
conditions that includes wear.

To evaluate the implication of the choice of free energy and state laws given in (16)-(18) in conjunction with
the balance of energy in (5) and the Clausius-Duhem inequality in (7), we consider evolutions such that the
left-hand time derivative and the right-hand time derivative of the state variables and the state functions exist
with possibly different values, meaning that the time derivative of these functions and variables may not
necessarily exist. That is, the state functions and the state variables are continuous in time, but not necessarily
smooth. Thus, a superimposed dot stands for either the left-hand time derivative or the right-hand time
derivative.

Then, the above assumption on the free energy implies that $\dot{\Psi}=0$ for $(u_n, \omega) \in C$; otherwise $\dot{\Psi}$ is not
defined. Besides, $(u_n, \omega) \notin C$ implies that $\dot{\Psi} = +\infty$ which indicates that such states are impossible. Moreover,
(18) implies that $\dot{S}=0$. Consequently, $\dot{E} = \dot{\Psi} + \dot{S}\dot{T} + \dot{S}\dot{T} = 0$ for thermodynamical admissible processes, i.e. for
$(u_n, \omega) \in C$. Furthermore, the state law given in (17) is equivalent to

$$
\mathcal{W} = p_n \ge 0, \quad u_n - \omega - g \le 0, \quad p_n (u_n - \omega - g) = 0. \tag{19}
$$

These conditions imply also that $p_n \dot{u}_n = p_n \dot{\omega}$. Putting $\dot{\Psi} = \dot{E} = S = 0$ and $p_n \dot{u}_n = p_n \dot{\omega}$ in (5) and (7) yield

$$
q_2 = p_n \dot{\omega} + p_t \dot{w}_t, \tag{20}
$$

$$
p_n \dot{\omega} + p_t \dot{w}_t \ge 0. \tag{21}
$$

These two relationships represent the balance of energy and the second principle of thermodynamics for the
contact interface. Furthermore, for later use, we remark that the contact conditions in (19) are optimality
conditions to the following variational inequality:

$$
p_n \ge 0 : \quad (u_n - \omega - g)(p'_n - p_n) \le 0 \quad \forall p'_n \ge 0. \tag{22}
$$

Next, Coulomb's law of friction and Archard's law of wear are introduced. Coulomb's law of friction is given
by the following principle of maximal dissipation:

$$
|p_t| \le \mu(p_n)_+ : \quad \dot{w}_t (p'_t - p_t) \le 0 \quad \forall |p'_t| \le \mu(p_n)_+, \tag{23}
$$

where $\mu$ is a constant friction coefficient, $|\cdot|$ stands for the absolute value and $(p_n)_+ = (p_n + |p_n|)/2$. Assuming that $k$ is a constant wear coefficient, Archard's law of wear is represented by

$$
\dot{\omega} = k(p_n)_+ |\dot{w}_t|. \tag{24}
$$

One might remark that the use of $(p_n)_+$ in (23) and (24), instead of just $p_n$, implies that these two
representations of Coulomb's and Archard's laws are defined also for $p_n < 0$ and not only for $p_n \ge 0$.

Let us evaluate the consequences of Coulomb's law and Archard's law with respect to the principles of thermodynamics given in (20) and (21). Taking $p'_t = 0$ in (23) implies that $p_t \dot{w}_t \ge 0$. Using this and (24), it is shown that the second law of thermodynamics represented by (21) is satisfied. Furthermore, (23) implies that $p_t = \mu(p_n)_+ \operatorname{sgn}(\dot{w}_t)$ when $\dot{w}_t \ne 0$. This together with (24) in (20) result in the following energy balance:

$$
q_2 = (k(p_n)_+^2 + \mu(p_n)_+) |\dot{w}_t|, \quad (25)
$$

which represents the friction and wear heat generation at the contact interface. Thus, all heat generated by
friction and wear flows into the body.

3. The system of discrete equations

Given proper initial conditions, a load history $t(t)$ and a displacement history $\xi(t)$ on a time interval $[0, \tau]$, the following initial boundary value problem is stated: Find $u: [0, \tau] \to \mathcal{V}$, $p_n: [0, \tau] \to K_n$, $p_t: [0, \tau] \to F(p_n)$ and $T: [0, \tau] \to \mathcal{T}$ such that for each $t \in [0, \tau]$:

$$
\int_{\Omega} E_{ijkl} \frac{\partial u_k}{\partial x_j} \frac{\partial v_i}{\partial x_j} dV - \int_{\Omega} G(T - T_0) \frac{\partial v_i}{\partial x_i} dV = \int_{\Gamma_t} t_i v_i dA - \int_{\Gamma_c} (p_t v_1 - p_n v_2) dA \quad \forall v \in \mathcal{V}, \quad (26)
$$

$$
\int_{\Omega} \rho c \ddot{T} \theta \, dV + \int_{\Omega} k_t \frac{\partial T}{\partial x_i} \frac{\partial \theta}{\partial x_i} \, dV = \int_{\Gamma_c} (k(p_n)_+^2 + \mu(p_n)_+) |\dot{w}_t| \theta \, dA \quad \forall \theta \in \mathcal{T}, \quad (27)
$$

where
$$E_{ijkl} = \tilde{\lambda}\delta_{ij}\delta_{kl} + 2\tilde{\mu}\delta_{ik}\delta_{jl},$$

$$G = (3\tilde{\lambda} + 2\tilde{\mu})\alpha,$$

$$\mathcal{V} = \{\mathbf{v} : \mathbf{v}(\mathbf{x}) = \mathbf{0}, \mathbf{x} \in \Gamma_u\},$$

$$K_n = \{p_n : p_n(\mathbf{x}) \ge 0, \mathbf{x} \in \Gamma_c\},$$

$$F(p_n) = \{p_t : |p_t(\mathbf{x})| \le \mu(p_n(\mathbf{x}))_+, \mathbf{x} \in \Gamma_c\}$$

and $\mathcal{V}$ is a function space of temperatures. In addition, weak forms of the tribological laws given in (22)–(24) are assumed to be valid on $\Gamma_c$ (see [21]). Eq. (26) is the principle of virtual work, equivalent to the equilibrium equations given in (1)–(3), with the thermoelastic properties in (9) inserted. Finally, Eq. (27) is the balances of energy for the body and the contact surface given in (15) and (25) put together using (13) and $\mathbf{q} \cdot \mathbf{n}_c = -q_2$.

### 3.1. Space discretization

The initial boundary value problem stated above is treated by introducing finite element approximations of $\mathcal{V}$ and $\mathcal{T}$ (see e.g. [7]), and evaluating the integrals over $\Gamma_c$ by an appropriate quadrature rule depending on the choice of finite elements. Let $f$ be a function over $\Gamma_c$, then the integrals are approximated by

$$\int_{\Gamma_c} f \, dA \approx \sum_{M \in \eta_c} I^M f(\mathbf{x}^M), \qquad (28)$$

where $I^M$ are the weighting factors and $\eta_c$ is the set of integration points $\mathbf{x}^M$ on $\Gamma_c$. When applying this formula, the contact boundary is divided into contact elements. If these elements coincide with the edges of the displacement elements, then it is possible to choose an integration rule such that the integration points will coincide with the nodal displacement points of the contact surface. For two-dimensional problems, examples of such rules are the trapezoidal quadrature rule for linear elements and Simpson's rule for quadratic elements. A study of different integration rules can be found in [8].

After performing the space discretization outlined above, one obtains the following discrete equilibrium equation from (26):

$$\bar{K}d - \hat{K}(T - T_0) = F - C_n^T P_n - C_t^T P_t, \qquad (29)$$

where $d = \{d_k^A\}$ represents the nodal displacements, $T = \{T^A\}$ is the nodal temperatures and $T_0 = \{T_0^A\}$ is the nodal reference temperatures at the nodal points $\mathbf{x}^A$. Furthermore, letting $N^A = N^A(\mathbf{x})$ represents the shape functions,

$$\bar{K} = [\bar{K}_{ik}^{BA}], \quad \bar{K}_{ik}^{BA} = \int_{\Omega} E_{ijkl} \frac{\partial N^A}{\partial x_l} \frac{\partial N^B}{\partial x_j} dV,$$

$$\hat{K} = [\hat{K}_i^{BA}], \quad \hat{K}_i^{BA} = \int_{\Omega} GN^A \frac{\partial N^B}{\partial x_i} dV,$$

$$F = \{F_t^B\}, \quad F_t^B = \int_{\Gamma_t} t_i N^B dA,$$

$$C_n = [C_{ni}^{MB}], \quad C_{ni}^{MB} = -N^B(\mathbf{x}^M)\delta_{2i},$$

$$C_t = [C_{ti}^{MB}], \quad C_{ti}^{MB} = N^B(\mathbf{x}^M)\delta_{1i},$$

$$P_n = \{P_n^M\}, \quad P_n^M = I^M p_n(\mathbf{x}^M),$$

$$P_t = \{P_t^M\}, \quad P_t^M = I^M p_t(\mathbf{x}^M).$$

The discretization of the balance of energy in (27) can be expressed by

$$M\dot{T} + OT = L(d, P_n), \qquad (30)$$

where
$$
M = [M^{BA}], \quad M^{BA} = \int_{\Omega} \rho c N^A N^B dV,
$$

$$
O = [O^{BA}], \quad O^{BA} = \int_{\Omega} k_{t} \frac{\partial N^{A}}{\partial x_{i}} \frac{\partial N^{B}}{\partial x_{i}} dV,
$$

$$
L(d, P_n) = \{L^B\}, \quad L^B = \sum_{M \in \eta_c} N^B(x^M) \left( \frac{k}{I^M} (P_n^M)_+^2 + \mu (P_n^M)_+ \right) |\dot{w}_t^M|
$$

and, in addition, $w_t^M = N^A(x^M)d_1^A - \xi$.

Furthermore, the discretization of the tribological laws can be obtained by assuming that (22)–(24) hold in every integration point $M \in \eta_c$, i.e.

$$
P_n^M \ge 0 : \quad (u_n^M - \omega^M - g^M)(P_n^{M'} - P_n^M) \le 0 \quad \forall P_n^{M'} \ge 0, \qquad (31)
$$

$$
|P_t^M| \le \mu(P_n^M)_+ : \quad \dot{w}_t^M (P_t^{M'} - P_t^M) \le 0 \quad \forall |P_t^{M'}| \le \mu(P_n^M)_+, \qquad (32)
$$

$$
\dot{\omega}^M = \frac{k}{I^M} (P_n^M)_+ |\dot{w}_t^M|. \qquad (33)
$$

Here, $u_n^M = -N^A(x^M)d_2^A$, $\omega^M = \omega(x^M)$ and $g^M = g(x^M)$. In [21], the above relationships were derived from weak forms of (22)–(24) by using the integration rule in (28). Finally, we rewrite (31) in the following manner for any $r > 0$ (see [11]): find $P_n^M \ge 0$ such that

$$
(P_n^M - (P_n^M + r(u_n^M - \omega^M - g^M))(P_n^{M'} - P_n^M) \ge 0 \quad \forall P_n^{M'} \ge 0, \qquad (34)
$$

implying that $P_n^M$ is the projection of $P_n^M + r(u_n^M - \omega^M - g^M)$ onto $K_n^h = \{P_n^M : P_n^M \ge 0\}$.

3.2. Time discretization

The time rates appearing in (30), (32) and (33) are approximated using a backward Euler discretization, i.e.

$$
\dot{T}(t_{k+1}) \approx \frac{T(t_{k+1}) - T(t_k)}{t_{k+1} - t_k}, \qquad (35)
$$

$$
\dot{w}_t^M (t_{k+1}) \approx \frac{w_t^M(t_{k+1}) - w_t^M(t_k)}{t_{k+1} - t_k}, \qquad (36)
$$

$$
\dot{\omega}^M(t_{k+1}) \approx \frac{\omega^M(t_{k+1}) - \omega^M(t_k)}{t_{k+1} - t_k}. \qquad (37)
$$

Defining $\bar{w}_t^M = w_t^M(t_{k+1}) - w_t^M(t_k)$, and inserting (35) and (36) in (30) yield

$$
(M + (t_{k+1} - t_k)O)T(t_{k+1}) = Q + \bar{L}(d(t_{k+1}), P_n(t_{k+1})), \qquad (38)
$$

where $Q = MT(t_k)$ and

$$
\bar{L}(d, P_n) = \left\{ \sum_{M \in \eta_c} N^B(x^M) \left( \frac{k}{I^M} (P_n^M)_+^2 + \mu (P_n^M)_+ \right) |\bar{w}_t^M| \right\}. \quad (39)
$$

By using (36) and the definition of $\bar{w}_t^M$, Coulomb's law in (32) can be rewritten, for any $r > 0$, as the following variational principle: find $|P_t^M| \le \mu(P_n^M)_+$ such that

$$
(P_t^M - (P_t^M + r\bar{w}_t^M))(P_t^{M'} - P_t^M) \ge 0, \quad \forall |P_t^{M'}| \le \mu(P_n^M)_+, \qquad (40)
$$

which shows that $P_t^M$ can be expressed as the projection of $P_t^M + r\bar{w}_t^M$ onto $F^h(P_n^M) = \{P_t^M : |P_t^M| \le \mu(P_n^M)_+\}$ (see [11]).

Finally, using (36) and (37), Archard's law in (33) is approximated by
$$
\omega^M(t_{k+1}) \approx \omega_0^M + \frac{k}{l^M} (P_n^M(t_{k+1}))_+ |\bar{w}_i^M|,
$$

(41)

where $\omega_o^M = \omega^M(t_k)$.

4. The algorithm

The initial boundary value problem defined in the beginning of Section 3 can now be treated by solving a
sequence of the discrete equations given in (29), (34), (38), (40) and (41) using Pang's algorithm presented
below. However, numerically tests have proven that the search direction involved in the algorithm can not be
solved accurately for this representation of the equations. Comparisons by analytical solutions given in [3] have
shown that the Newton step smoothens out the temperature gradients in a wrong way; in addition, the algorithm
converges slowly. This depends probably on the great differences in the numerical values of K, K̂ and
(M + (t_{k+1} - t_k)O), leading to a system of ill-condition equations. This problem is overcome by putting together
(29) and (38) to form

$$
\bar{K}d - \tilde{K}(Q + \bar{L}(d, P_n)) + \hat{K}T_0 = F - C_n^T P_n - C_t^T P_t,
$$

(42)

where $\tilde{K} = \hat{K}(M + (t_{k+1} - t_k)O)^{-1}$. In such manner, the temperature is eliminated from the Newton step and appears instead as an inner variable in the algorithm.

Thus, letting $y = (d, P_n, P_t)$, the following system of equations are solved for each time increment:

$$
H(y) = \begin{cases}
\bar{K}d - \tilde{K}(Q + \bar{L}(d, P_n)) + C_n^T P_n + C_t^T P_t + \hat{K}T_0 - F \\
- P_n + \Pi_n(d, P_n) \\
- P_t + \Pi_t(d, P_n, P_t)
\end{cases}
$$

(43)

where

$$
\Pi_n(d, P_n) = \left\{ \left( P_n^M + r \left( u_n^M - \omega_o^M - \frac{k}{l^M} (P_n^M)_+ |\bar{w}_i^M| - g^M \right) \right)_+ \right\}_+
$$

(44)

is the projection of $P_n^M + r(u_n^M - \omega^M - g^M)$ onto $K_n^h$ with the approximation of Archard's wear law given in (41) inserted and

$$
\Pi_t(d, P_n, P_t) = \begin{cases} P_t^M + r\bar{w}_t^M & \text{if } |P_t^M + r\bar{w}_t^M| \le \mu(P_n^M)_+ \\ \mu(P_n^M)_+ & \text{sgn}(P_t^M + r\bar{w}_t^M) \text{ otherwise} \end{cases}
$$

(45)

is the projection of $P_t^M + r\bar{w}_t^M$ onto $F^h(P_n^M)$. The equation system in (43) is an extension of the quasi-augmented Lagrangian formulations considered previously in [21,22]. These formulations are simply obtained by setting $\hat{K} = 0$ in the above formulation.

4.1. Pang's Newton algorithm

In [21], it was shown that, for $\hat{K} = 0$, the equation system in (43) is Bouligand differentiable, i.e. Lipschitz continuous and directionally differentiable. Obviously, this holds also for the case with $\hat{K} \neq 0$ as $\bar{L}(d, P_n)$ is Lipschitz continuous and directionally differentiable. For solving systems of Bouligand differentiable equations, Pang [16] suggested the following globally convergent algorithm (concerning convergence result of the algorithm, see also [4]):

ALGORITHM. Let $\beta, \gamma$ and $\epsilon$ be given scalars with $\beta \in (0, 1)$, $\gamma \in (0, 1/2)$ and $\epsilon$ small. Repeat the following steps for each time increment $k + 1$:

(0) Let $y^0$ be given from the previous time step $k$ and set $l = 0$.

(1) Find a search direction $z = (z_d, z_n, z_t)$ such that
$$ H(y') + H'(y'; z) = 0, \qquad (46) $$

where $H'(y'; z)$ is the directional derivative.

(2) Let $\alpha' = \beta^{m_l}$, where $m_l$ is the first non-negative integer $m$ for which the following decrease criterion holds:

$$ g(y' + \beta^m z) \le (1 - 2\gamma\beta^m)g(y'), \quad g(y) = \frac{1}{2} H^T(y)H(y). $$

(3) Set $y^{l+1} = y' + \alpha'z$.

(4) If $g(y^{l+1}) \le \epsilon$, then terminate with $y^{l+1}$ as an approximate zero of $H(y)$. Otherwise, replace $l$ by $l+1$ and return to step 1.

After convergence in each time step is achieved, the variables $T$, $w_t^M$ and $\omega^M$ are updated.

## 4.2. Approximation of directional derivatives

The time-consuming part of the algorithm is to find the search direction in (46). Generally, this system of equations is non-linear due to the non-linearity of $H'(y; z)$ in the second argument, which in turn depends on the non-differentiability of $H(y)$. Of course, $H'(y; z) = \nabla H(y)z$ at states where $H(y)$ is differentiable. The nonlinearity of (46) might limit the effectiveness of the algorithm. This drawback is overcome by introducing an appropriate approximation of $H'(y; z)$ such that (46) results in a system of linear equations (see [4,21,22]).

The non-differentiability of $H(y)$ concerns the functions $\tilde{L}(y) = \tilde{L}(d, P_n)$, $\Pi_n(y) = \Pi_n(d, P_n)$ and $\Pi_l(y) = \Pi_l(d, P_n, P_l)$; otherwise, $H'(y; z) = \nabla H(y)z$. In order to obtain a linear equation system in (46), the following approximations of $\tilde{L}'(y; z)$, $\Pi'_n(y; z)$ and $\Pi'_l(y; z)$ are introduced:

$$ \tilde{L}'(y; z) \approx \left\{ \sum_{M \in \eta_{c1}} N^B(x^M) \left( \left( \frac{2kP_n^M}{I^M} + \mu \right) |\bar{w}_t^M| z_n^M + \left( \frac{kP_n^M}{I^M} + \mu \right) \frac{P_n^M \bar{w}_t^M z_{dl}^M}{|\bar{w}_t^M|^2} \right) \right\}, $$

$$ \Pi'_n(y; z) \approx \begin{cases} 0, & M \in \eta_{c2} \\ z_n^M + r z_{dn}^M, & M \in \eta_{c3} \\ z_n^M + r \left( z_{dn}^M - \frac{k |\bar{w}_t^M|}{I^M} \left( z_n^M + \frac{P_n^M \bar{w}_t^M z_{dl}^M}{|\bar{w}_t^M|^2} \right) \right), & M \in \eta_{c4} \end{cases}, $$

$$ \Pi'_l(y; z) \approx \begin{cases} 0, & M \in \eta_{c5} \\ z_l^M + r z_{dl}^M, & M \in \eta_{c6} \\ \mu z_n^M \operatorname{sgn}(P_l^M + r \bar{w}_l^M), & M \in \eta_{c7} \end{cases}, $$

where

$$ \eta_{c1} = \{M : P_n^M |\bar{w}_l^M| > 0\}, $$

$$ \eta_{c2} = \{M : P_n^M + r(u_n^M - \omega_o^M - \frac{k}{I^M} (P_n^M)_+ |\bar{w}_l^M| - g^M) \le 0\}, $$

$$ \eta_{c3} = \{M : P_n^M |\bar{w}_l^M| \le 0, P_n^M + r(u_n^M - \omega_o^M - \frac{k}{I^M} (P_n^M)_+ |\bar{w}_l^M| - g^M) > 0\}, $$

$$ \eta_{c4} = \{M : P_n^M |\bar{w}_l^M| > 0, P_n^M + r(u_n^M - \omega_o^M - \frac{k}{I^M} P_n^M |\bar{w}_l^M| - g^M) > 0\}, $$

$$ \eta_{c5} = \{M : P_n^M \le 0\}, $$

$$ \eta_{c6} = \{M : P_n^M > 0, |P_l^M + r \bar{w}_l^M| \le \mu P_n^M\}, $$

$$ \eta_{c7} = \{M : P_n^M > 0, |P_l^M + r \bar{w}_l^M| > \mu P_n^M\} $$
and

$$z_{dn} = C_n z_d = \{z_{dn}^M\}, \quad z_{dt} = C_t z_d = \{z_{dt}^M\}, \quad z_n = \{z_n^M\} \text{ and } z_t = \{z_t^M\}.$$

## 5. Numerical examples

Five different problems are considered for a structure subjected to boundary conditions and a displacement history as given in Fig. 2. The problems differ only in the constitutive settings, i.e. the same boundary conditions and displacement history (except for the number of cycles) are used for all five problems. The dimensions of the structure are 3.0 × 1.0 cm². Furthermore, the structure is fixed at the top and unilaterally constrained at the bottom by a rigid support using a negative initial gap of $g = -0.1$ µm. In addition, the support is subjected to a cyclic tangential displacement history with an amplitude of 1.0 mm according to Fig. 2. In the context of wear, this type of boundary condition can be called a fretting condition. Although $\xi$ is small, this condition will result in global sliding fretting. It is a major principle difference between the contact stresses obtained for global sliding fretting compared to stick-slip fretting. This is discussed in more detail further on. The thermoelastic properties of the structure are approximated by using 42 × 28 bilinear quadrilateral elements.<fnref n="3" /> The potential contact surface at the bottom of the structure is approximated by 42 contact elements using the trapezoidal rule in such manner that the integration points coincide with the nodal displacement points of the contact surface. Concerning the time discretization, each cycle is treated using four time steps, i.e. $t_{k+1} = t_k + 0.01$ s.

For the system given above, the following problems, defined by five different constitutive settings, are considered: the first problem is pure elastic, meaning that no thermal properties are used, with $\mu = k = 0$, the second problem is pure elastic with $\mu = 0.3$ and $k = 0$, the third problem is pure elastic with $\mu = 0.3$ and $k = 1.0 \times 10^{-11}$ Pa⁻¹ (this is a typical wear coefficient for adhesive wear in steel assemblies), the fourth problem is thermoelastic with $\mu = 0.3$ and $k = 0$, and, finally, the fifth problem is thermoelastic with $\mu = 0.3$ and $k = 1.0 \times 10^{-11}$ Pa⁻¹.

The contact pressures of problem one ($\mu = k = 0$, no heat) and problem two ($\mu = 0.3$, $k = 0$, no heat) are compared in Fig. 3. The only difference between these two problems is the value of the coefficient of friction. Obviously, $p_n$ is constant in time for $\mu = 0$, but, for $\mu = 0.3$, $p_n$ depends on the sliding direction of the support. In Fig. 3, $p_n$ is plotted at states when $t = n \times 0.04$ s, where n is a positive integer. The solution will be the same at times when $\xi > 0$ and global sliding is developed. On the other hand, the reflected image of $p_n$ will appear when there is global sliding with $\xi < 0$. Furthermore, the dissipation is non-zero for $\mu = 0.3$ in contrast to the case for $\mu = 0$ when it is zero. Actually, the frictional dissipation is constant for each cycle. Thus, if one is interested in studying the dissipation for a large number of cycles such as the dissipation in frictional dampers, then the constitutive assumptions of problem two might not be sufficient due to the fact that real engineering materials eventually will be broken during a repeated frictional process and consequently the dissipation will eventually be equal to zero. This fact is included in the constitutive settings of problem three and five. For those

<figure><img src="image_4.png" /><figcaption>Fig. 2. The structure with the boundary conditions and the displacement history for the five problems considered.</figcaption></figure>

<footnote n="3">³ Young's modulus is 210 GPa, Poisson's ratio is 0.3, $\rho = 7800$ kg/m³, $\alpha = 12.0 \times 10^{-6}$ l/K, $c = 460$ J/kg·K and $k_t = 46$ W/m·K.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Comparison of the contact pressures for problem one ($\mu$=$k$=0, no heat) and problem two ($\mu$=0.3, $k$=0, no heat).</figcaption></figure>

problems the frictional dissipation turns toward zero after many cycles depending on the surface break down due to wear.

In Fig. 4, the contact pressure and the wear gap are plotted for problem three ($\mu$ = 0.3, $k$ = $1.0 \times 10^{-11}$ Pa$^{-1}$, no heat) after 30 cycles ($t$ = 1.2 s). At this instant, the contact pressure is approximately zero since the initial gap is almost worn away, i.e. $\omega \approx -g \approx 0.1$ µm over the contact surface. However, during the first cycle of this fretting process, the wear gap is approximately zero and the contact pressure is similar to the pressure obtained in problem two, see Fig. 3. Thus, although a physically realistic value on the wear coefficient is used, a large reduction of the contact pressure has been obtained within a few number of cycles. This depends not only on the global sliding conditions but also on the fact that the wear model used in this work does not include the transportation of wear debris in the contact interface but rather annihilates the wear debris. Of course, the rate of reduction in contact pressure can be decreased by diminishing the wear coefficient. Perhaps a better approach to obtain such a reduction would be to extend the formulation used in this work by a third body of wear debris, following the ideas of Godet [6]. Another approach in that direction is accomplished in problem five, where the rate of reduction in contact pressure is decreased by adding thermal properties to the constitutive settings of this problem.

One might notice that the result obtained above would be completely different if the fretting conditions would be chosen such that no global sliding would be developed. If the amplitude on $\xi$ would be taken sufficiently small such that no slip would be developed over the contact surface (global stick), then, of course, no wear would be initiated. However, in the intermediate with amplitudes on $\xi$ implying stick-slip, the contact pressure would turn toward zero in the slipping regions and the pressure would increase in the sticking regions (see [21,22]). In addition, large values and large gradients in the contact stresses would be developed in the regions between stick and slip. Concerning the field of fretting fatigue, such stress distributions might initiate cracks near the contact surface leading to fatigue failures of the structure. Similar observations have been found experimentally by using so-called fretting maps in e.g. [17,23].

<figure><img src="image_6.png" /><figcaption>Fig. 4. The contact pressure and the wear gap for problem three ($\mu$ = 0.3, $k$ = $1.0 \times 10^{-11}$ Pa$^{-1}$, no heat) plotted after 30 cycles.</figcaption></figure>
In Fig. 5, two frictional hot spots and two peaks in the contact pressure are plotted after 30 cycles for problem four ($\mu$ = 0.3, $k$ = 0, heat). After the first cycle the contact pressure is similar in shape to the pressure of problem two, see Fig. 3, but due to the thermal expansion induced by frictional heat generation the contact pressure changes rapidly from this shape during the evolution of cycles. The shape obtained in the temperature field and the contact pressure are very sensitive to the geometry and the boundary conditions. For instance, if the dimensions of the structure are changed to 1.0 x 1.0 cm², then only one frictional hot spot will appear. Furthermore, if the velocity of the support is decreased sufficiently, then no hot spots at all will appear. An obvious drawback of this constitutive model, concerning this type of evolution problems, is that no damages are induced due to the high temperatures and high contact stresses, neither in the bulk material nor in the contact surface. In real engineering materials such high temperatures and high stresses would result in damages of the material in some way, which in turn would reduce the amount of frictional heat generation and consequently result in lower temperatures and pressures than predicted here. This fact is included in the constitutive settings of problem five by adding Archard's law of wear to the constitutive laws used here. In that problem, the surface damages in form of wear imply that no frictional hot spots are developed.

In Fig. 6, the temperature field, the contact pressure and the wear gap are plotted after 210 cycles ($t$ = 8.4 s) for problem five ($\mu$ = 0.3, $k$ = 1.0 x 10⁻¹¹ Pa⁻¹, heat). At this instant, the solution has almost reached a steady state. The temperature field is approximately constant at 0.0987 K. Notice, in comparison to problem four, that no frictional hot spots are developed. Instead the increase in temperature is astonishingly low. This fact depends on the wear process, which decreases the frictional and wear heat generation. However, following the discussion on problem three, this reduction is probably too large due to the annihilation of wear debris. If the wear debris would be included in the model as a third body, then the frictional and wear heat generation would be higher, but not as high as in problem four. Although the temperature is very low, it has a major effect on the wear rate. After seven times more cycles compared to problem three, the contact pressure is still much higher here. This depends on the thermal expansion. Furthermore, the thermal expansion also increases the amount of wear; the wear gap is no longer constant and the maximum in wear gap obtained is more than 20% larger in comparison to problem three.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Two frictional hot spots and two peaks in the contact pressure are developed after 30 cycles for problem four ($\mu$ = 0.3, $k$ = 0, heat).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. The temperature field, the contact pressure and the wear gap plotted after 210 cycles for problem five ($\mu=0.3, k=1.0 \times 10^{-11} \text{ Pa}^{-1}$, heat). Notice that the scale on the temperature is given in milli-Kelvin.</figcaption></figure>

The temperature rise in fretting is a subject of considerable interest. Some researchers claim that the presence of $\alpha$-Fe$_2$O$_3$ in the wear debris must be taken as evidence that the interface temperature exceeds 300-600 K in fretting (see e.g. [24]). Obviously, such high temperatures are not predicted in problem five. Nevertheless, by neglecting the removal of material, as in problem four, even higher temperatures are predicted, see Fig. 5. However, measurements by other researchers indicate that the average temperature in contacting asperities subjected to fretting conditions is much lower; magnitudes of 20-30 K have been measured (see e.g. [5,19]). Such temperature increases are more comparable to the increases obtained by the wear model used in the present paper.

For fretting problems with stick-slip conditions, the frictional heat generation has a minor influence on the contact pressure and the wear gap using this wear model. Still, other initial and boundary conditions on the temperature and the heat flux vector can have a major effect on the contact pressure and the wear gap for stick-slip fretting problems. For instance, one might think of problems where external thermal boundary conditions induce stick-slip fretting. This type of problem can of course be treated with the algorithm suggested in this work.

Finally, some execution statistics for the five problems considered are provided. The number of Newton iterations per time increment (New./inc.) and the number of line-searches, i.e. Armijo iterations, per Newton iterations (Arm./New.) are given in Table 1. In order to get a fair comparison, the statistics are determined for

<table><caption>Table 1
Some execution statistics for the five problems considered</caption><thead><tr><td>Problem-constitutive setting</td><td>New./inc.</td><td>Arm./New.</td></tr></thead><tbody><tr><td>(1) μ = k = 0, no heat</td><td>1.0</td><td>1.0</td></tr><tr><td>(2) μ = 0.3, k = 0, no heat</td><td>1.0</td><td>1.0</td></tr><tr><td>(3) μ = 0.3, k = 1.0 × 10<sup>-11</sup> Pa<sup>-1</sup>, no heat</td><td>2.25</td><td>5.16</td></tr><tr><td>(4) μ = 0.3, k = 0, heat</td><td>2.7</td><td>3.96</td></tr><tr><td>(5) μ = 0.3, k = 1.0 × 10<sup>-11</sup> Pa<sup>-1</sup>, heat</td><td>1.725</td><td>3.07</td></tr></tbody></table>
the ten first cycles of each problem. Apparently, problem four is the most time-consuming to solve, and problem one and two are the fastest ones. The reason for problem two to have the same statistics as problem one is that there is only global sliding and no stick-slip developed in the solution of this problem as well as in the solutions of the other problems. If there would be stick-slip in the solutions, then the statistics would be higher for problem 2-5, at least the number of New./inc. would be higher.

## 6. Concluding remarks

In this work thermoelastic wear problems have been solved by a finite element method. The method is a direct method, meaning that no iterative strategy between a mechanical and a thermal part of the problem is used, but instead the finite element equations are solved using a Bouligand differentiable Newton method. The method obtained is effective and robust.

Numerically, it has been found that the internal state variable measuring the wear gap has a great influence on the temperature field generated by frictional heat. For a particular frictional heat generation problem, it is numerically shown that when this internal state variable is included in the analysis, then the increase in temperature is very low. On the contrary, if wear is excluded, then two frictional hot spots are developed for the same problem. However, although the temperature increase is very low for the wear problem, the frictional heat generation has a major effect on the wear rate.

The wear model utilized in this work does not include the transportation of wear debris, but instead the internal state variable measuring the wear gap rather annihilates the wear debris. From the numerical result discussed above, it is clear that the transportation of wear debris will affect the frictional heat generation which in turn will affect the wear rate. Therefore, for future plans, it would be of interest to extend this wear model by a third body of wear debris and to study its influence on the wear rate.

## Acknowledgments

This work was supported by the Swedish Research Council for Engineering Sciences (TFR). I wish to thank Anders Klarbring and Lars Johansson for fruitful discussions.

## References

[1] P. Alart and A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Methods Appl. Mech. Engrg. 92 (1991) 353-375.
[2] J.R. Barber and M. Cominou, Thermoelastic contact problems, in: R.B. Hetnarski, ed., Thermal Stresses III (North-Holland, Amsterdam, 1989).
[3] H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids (The Clarendon Press, Oxford, 1959).
[4] P.W. Christensen, A. Klarbring, J.S. Pang and N. Strömberg, Formulation and comparison of algorithms for frictional contact problems, Int. J. Numer. Methods Engrg., to appear.
[5] H.M. Ghasemi, M.J. Furey and C. Kadjas, Surface temperatures and fretting corrosion of steel under conditions of fretting contact, Wear 162-164 (1993) 357-369.
[6] M. Godet, The third body approach: a mechanical view of wear, Wear 100 (1984) 437-452.
[7] T.J.R. Hughes, The Finite Element Method, Linear Static and Dynamic Finite Element Analysis (Prentice-Hall, Inc., Englewood Cliffs, N.J., 1987).
[8] L. Johansson, A study on the discretization of a mixed formulation of the frictionless contact problem, M.Sc. Thesis, Linköping University, Sweden, 1986.
[9] L. Johansson and A. Klarbring, Thermoelastic frictional contact problems: modelling, finite element approximation and numerical realization, Comput. Methods Appl. Mech. Engrg. 105 (1993) 181-210.
[10] A. Klarbring and G. Björkman, A mathematical programming approach to contact problems with friction and varying contact surfaces, Comput. Struct. 30 (1988) 1185-1198.
[11] A. Klarbring, Mathematical programming and augmented Lagrangian methods for frictional contact problems, in: A. Curnier, ed., Proc. Contact Mechanics Int. Symp. (PPUR, 1992) 409-422.
[12] T.A. Laursen and V.G. Oancea, On the constitutive modeling and finite element computation of rate-dependent frictional sliding in large deformations, Comput. Methods Appl. Mech. Engrg. 143 (1997) 197-227.
[13] J. Lemaitre and J.-L. Chaboche, Mechanics of Solid Materials (Cambridge University Press, Cambridge, 1990).
[14] G.A. Maugin, The Thermomechanics of Plasticity and Fracture (Cambridge University Press, Cambridge, 1992).
[15] V.G. Oancea and T.A. Laursen, A finite element formulation of thermomechanical rate-dependent frictional sliding, Int. J. Numer. Methods Engrg. 40 (1997) 4275-4311.
[16] J.S. Pang, Newton's method for B-differentiable equations, Math. Operations Res. 15 (1990) 311-341.
[17] C. Petiot, L. Vincent, K. Dang Van, N. Maouche, J. Foulquier and B. Journet, An analysis of fretting-fatigue failure combined with numerical calculations to predict crack nucleation, Wear 181-183 (1995) 101-111.
[18] 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.
[19] E.S. Sproles Jr. and D.J. Duquette, Interface temperature measurements in the fretting of a medium carbon steel, Wear 47 (1978) 387-396.
[20] N. Strömberg, L. Johansson and A. Klarbring, Derivation and analysis of a generalized standard model for contact, friction and wear, Int. J. Solids Struct. 33 (1996) 1817-1836.
[21] N. Strömberg, An augmented Lagrangian method for fretting problems, Eur. J. Mech., A/Solids 16 (1997) 573-593.
[22] N. Strömberg, A Newton method for three-dimensional fretting problems, Int. J. Solids Struct. 36(14) (1999) 2075-2090.
[23] O. Vingsbo and S. Söderberg, On fretting maps, Wear 126 (1988) 131-147.
[24] R.B. Waterhouse, Fretting Corrosion (Pergamon Press, Oxford, 1972).
[25] P. Wriggers and C. Miehe, Contact constraints within coupled thermomechanical analysis—A finite element model, Comput. Methods Appl. Mech. Engrg. 113 (1994) 301-319.
[26] P. Wriggers and G. Zavarise, Thermomechanical contact—a rigorous but simple numerical approach, Comput. Struct. 46 (1993) 47-53.
[27] P. Wriggers and G. Zavarise, Application of augmented Lagrangian techniques for non-linear constitutive laws in contact interfaces, Comm. Numer. Methods Engrg. 9 (1993) 815-824.
[28] G. Zavarise, P. Wriggers and B.A. Schrefler, On augmented Lagrangian algorithms for thermomechanical contact problems with friction, Int. J. Numer. Methods Engrg. 38 (1995) 2929-2949.
[29] G. Zavarise, P. Wriggers, E. Stein and B.A. Schrefler, Real contact mechanisms and finite element formulation—a coupled thermomechanical approach, Int. J. Numer. Methods Engrg. 35 (1992) 767-785.