A chemo-mechano-thermodynamical contact theory for adhesion, friction and (de)bonding reactions

This work presents a self-contained continuum formulation for coupled chemical, mechanical and thermal contact interactions. The formulation is very general and hence admits arbitrary geometry, deformation and material behavior. All model equations are derived rigorously from the balance laws of mass, momentum, energy and entropy in the framework of irreversible thermodynamics, thus exposing all the coupling present in the field equations and constitutive relations. In the process, the conjugated kinematic and kinetic variables for mechanical, thermal and chemical contact are identified, and the analogies between mechanical, thermal and chemical contact are highlighted. Particular focus is placed on the thermodynamics of chemical bonding distinguishing between exothermic and endothermic contact reactions. Distinction is also made between long-range, non-touching surface interactions and short-range, touching contact. For all constitutive relations examples are proposed and discussed comprehensively with particular focus on their coupling. Finally, three analytical test cases are presented that illustrate the thermo-chemo-mechanical contact coupling and are useful for verifying computational models. While the main novelty is the extension of existing contact formulations to chemical contact, the presented formulation also sheds new light on thermo-mechanical contact, since it is consistently derived from basic principles using only few assumptions.


1
identity tensor in R 3 a α := ∂x/∂ξ α , α ∈ {1, 2}; covariant tangent vectors at surface point x; unit: [1] or [m] a α contravariant tangent vectors at x; α ∈ {1, 2}; unit: [1] or [m −1 ] depending on ξ α a p α and a α p ; tangent vectors at x p ∈ S 1 b k prescribed, mass-specific body force at x k ∈ B k ; unit: [N/kg] B k set of points contained in body k ∂B k set of points on the surface of body k ∂ q B k ⊂ ∂B k ; boundary of body k where heat fluxq k is prescribed ∂ t B k ⊂ ∂B k ; boundary of body k where tractiont k is prescribed D k rate of deformation tensor at x k ∈ B k ; unit: [1/s] D i k inelastic part of D k da k area element on the current surface S k ; unit: [m 2 ] dv k volume element in the current body B k ; unit: [m 3 ] E k Green-Lagrange strain tensor at x k ∈ B k ; unit-free E e k elastic part of E k E i k inelastic part of E k E total energy of the two-body system; unit: [J] e c bond-specific contact enthalpy on surface S c (per bonding site); unit: [J] φ non-dimensional bonding state at x c ∈ S c ψ k mass-specific Helmholtz free bulk energy at x k ∈ B k ; unit [J/kg] Ψ k = ρ 0k ψ k ; Helmholtz free bulk energy density per undeformed volume of B k ; unit: [J/m 3 ] ψ 12 Helmholtz free interaction energy between x 1 ∈ S 1 and x 2 ∈ S 2 ; unit: [J] ψ c bond-specific Helmholtz free contact energy on surface S c (per bonding site); unit: [J] Ψ c = N c ψ c ; Helmholtz free contact energy density per undeformed area; unit: [J/m 2 ] F k deformation gradient at x k ∈ B k ; unit-free g := x 2 − x 1 ; gap vector between surface points x 1 ∈ S 1 and x 2 ∈ S 2 ; unit: [m] g n := g · n; normal gap g t := g − g n n = (a α ⊗ a α )g; tangential gap vector g t =ξ α p a p α ; Lie derivative of g t at x 1 = x p ; unit: [m/s] g e reversible (elastic) part of g associated with sticking contact g i irreversible (inelastic) part of g associated with sliding contact η e k mass-specific external bulk entropy production rate at x k ∈ B k ; unit: [J/(kg K s)] η i k mass-specific internal bulk entropy production rate at x k ∈ B k ; unit: [J/(kg K s)] η i c area-specific internal contact entropy production rate at x c ∈ S c ; unit: [J/(K m 2 s)] h heat transfer coefficient between body B 1 and body B 2 ; unit: [J/(K m 2 s)] h k heat transfer coefficient between body B k and an interfacial medium; unit: [J/(K m 2 s)] J k change of volume at x k ∈ B k ; unit-free J sk change of area at x k ∈ S k ; unit-free J c reference value for the area change; either J c = J s1 or J c = J s2 k ∈ {1, 2}; body index K kinetic energy of the two-body system; unit: [J] µ c chemical contact potential per bonding site at x c ∈ S c ; unit: [J] M c = n c µ c ; chemical contact potential per current area at x c ∈ S c ; unit: [J/m 2 ] M c = N c µ c ; chemical contact potential per reference area at x c ∈ S c ; unit: [J/m 2 ] n k current density of bonding sites at x k ∈ S k and time t > 0; unit: [1/m 2 ] N k = J sk n k ; initial density of bonding sites at x k ∈ S k ; unit: [1/m 2 ] n c reference value for the current bonding site density; either n c = n 1 or n c = n 2 N c = J c n c ; reference value for the initial bonding site density; either N c = N 1 or N c = N 2 n b k current density of bonded bonding sites at x k ∈ S k and time t > 0; unit: [1/m 2 ] n ub k current density of unbonded bonding sites at x k ∈ S k and time t > 0; unit: [1/m 2 ] n b := n b 1 in case n b 2 = n b 1 n k outward unit normal vector at x k ∈ S k n p := n 1 (ξ α p , t); outward unit normal vector at x p ∈ S 1 P k subset of B k or S k p c contact pressure, i.e. normal part of t c ; unit: [N/m 2 ] q k heat influx per current area; unit: [J/(m 2 s)] q k heat flux vector per current area; unit: [J/(m 2 s)] q k prescribed heat influx per current area at x k ∈ S k ; unit: [J/(m 2 s)] q c k contact heat influx per current area at x k ∈ S c ; unit: [J/(m 2 s)] q c m := (q c 1 + q c 2 )/2; mean contact heat influx into B 1 and B 2 q c t := (q c 1 − q c 2 )/2; transfer heat flux from body B 2 to body B 1 q k entropy influx per current area; unit: J/(K m 2 s)] q k entropy flux vector per current area; unit: [J/(K m 2 s)] q k prescribed entropy influx per current area at x k ∈ S k ; unit: [J/(K m 2 s)] q c k contact entropy influx per current area at x c ∈ S c ; unit: [J/(K m 2 s)] ρ k current mass density at x k ∈ B k and time t > 0; unit: [kg/m 3 ] ρ 0k = J k ρ k ; initial mass density at x k ∈ B k ; unit: [kg/m 3 ] r k prescribed, mass-specific heat source at x k ∈ B k ; unit: [J/(kg s)] R c bonding reaction rate per current area; unit: [1/(m 2 s)] σ k Cauchy stress tensor at x k ∈ B k ; unit: [N/m 2 ] S k second Piola-Kirchhoff stress tensor at x k ∈ S k ; unit: [N/m 2 ] s k mass-specific bulk entropy at x k ∈ B k ; unit: [J/(kg K)] if J c = J s2 t t tangential contact traction, i.e. tangential part of t c T k temperature at x k ∈ B k ; unit: [K] T c temperature of the interfacial contact medium at x ∈ S c ; unit: [K] [[T ]] := T 2 − T 1 ; temperature jump (or temperature gap) across the contact interface u k mass-specific internal bulk energy at x k ∈ B k ; unit: [J/kg] U k = ρ 0k u k ; internal bulk energy density per undeformed volume of B k ; unit: [J/m 3 ] u 12 internal interaction energy between x 1 ∈ S 1 and x 2 ∈ S 2 ; unit: [J]   x c current position of a material point at time t > 0 on contact surface S c ; unit: [m] x p = x 1 (ξ α p , t) ∈ S 1 ; closest projection point of the projection x 2 ∈ S 2 → S 1

Introduction
Many applications in science and technology involve coupled interactions at interfaces. An example is the mechanically sensitive chemical bonding appearing in adhesive joining, implant osseointegration and cell adhesion. Another example is the thermal heating arising in frictional contact, adhesive joining and electrical contacts. Further examples are the temperaturedependent mechanical contact conditions in melt-based production technologies, such as welding, soldering, casting and additive manufacturing. A general understanding and description of these examples requires a general theory that couples chemical, mechanical, thermal and electrical contact. Such a theory is developed here for the first three fields in the general framework of nonlinear continuum mechanics and irreversible thermodynamics. Two cases are considered in the theory: touching contact and non-touching interactions. While the former is dominating at large length scales, the latter provides a link to atomistic contact models.
There is a large literature body on coupled contact models. General approaches deal, however, with two-field and not with three-field contact coupling, such as is considered here. General thermo-mechanical contact models have appeared in the early nineties, starting with the work of Zavarise et al. (1992) that combined a frictionless contact model with heat conduction. Following that, Johansson and Klarbring (1993) presented the full coupling for linear thermo-elasticity. This was extended by Wriggers and Miehe (1994) to large deformations using an operator split technique for the coupling -a staggering scheme based on successively solving mechanical and thermal subproblems. Oancea and Laursen (1997) extended this to a monolithical coupling formulation and proposed general constitutive models for thermo-mechanical contact. Following these initial works, many thermo-mechanical contact studies have appeared, studying the extension to wear (Strömberg et al., 1996;Stupkiewicz and Mróz, 1999;Molinari et al., 2001), rough-surface contact (Willner, 1999) multiscale contact (Temizer and Wriggers, 2010;Temizer, 2014Temizer, , 2016 and adhesion (Dittmann et al., 2019). Beyond that, there have also been many advancements in the computational description of thermo-mechancial contact (Agelet de Saracibar, 1998;Laursen, 1999;Strömberg, 1999;Pantuso et al., 2000;Adam and Ponthot, 2002;Xing and Makinouchi, 2002;Bergman and Oldenburg, 2004;Rieger and Wriggers, 2004;Hüeber and Wohlmuth, 2009;Dittmann et al., 2014;Khoei et al., 2018;Seitz et al., 2019). General chemo-mechanical contact models go back even further -to the work of Derjaguin et al. (1975) that describes molecular, e.g. van der Waals, adhesion between an elastic sphere and a half-space. Argento et al. (1997) then extended this approach to a general surface formulation for interacting continua. The work was then generalized to a nonlinear continuum mechanical contact formulation by Sauer (2006) and Sauer and Li (2007a,b). Subsequently, the framework was applied to the study of cell adhesion (Zeng and Li, 2011), generalized to various surface interaction models (Sauer and De Lorenzis, 2013), and combined with sliding friction models (Mergel et al., 2019). In a strict sense, the chemo-mechanical contact models mentioned so far are not coupled models. Rather, the chemical surface interaction is described by distancedependent potentials. Hence, the resulting problem is a single-field problem that only depends on deformation. A second field variable representing the chemical contact state is not used. This is different to the debonding model of Frémond (1988). There a state variable is introduced in order to describe irreversible damage during debonding. It is essentially a phenomenological debonding model, where the bond degrades over time following a first order ordinary differential equation (ODE). Raous et al. (1999) extended the model to sliding contact and implemented it within a finite element formulation. Its finite element implementation is also discussed in Wriggers (2006) in the framework of large deformations. Subsequently the model has been extended to thermal effects (Bonetti et al., 2009), applied to multiscale contact (Wriggers and Reinelt, 2009) and generalized to various constitutive models (Del Piero and Raous, 2010), among others. Even though the Frémond model is a coupled two-field model, it only describes debonding and not bonding. It therefore does not provide a general link to chemical contact reactions. The adhesion and debonding models mentioned above are similar to cohesive zone models (CZMs). Those propose phenomenological traction-separation laws for debonding, that are often derived from a potential, like the seminal model of Xu and Needleman (1993), thus ensuring thermodynamic consistency. CZMs have been extended to thermo-mechanical debonding through the works of Hattiangadi and Siegmund (2004); Willam et al. (2004); Fagerström and Larsson (2008);Özdemir et al. (2010); Fleischhauer et al. (2013) and Esmaeili et al. (2016), among others. Recently, CZMs have also been coupled to a hydrogen diffusion model in order to study fatigue (Busto et al., 2017). CZMs usually have a damage/degradation part that sometimes follows from an evolution law, e.g. see Willam et al. (2004). This makes them very similar to the Frémond (1988) model. Like the Frémond model, CZMs have not yet been combined with chemical contact reactions, and so this aspect is still absent in general contact models. Adhesion models based on chemical bonding and debonding reactions have been developed by Bell and coworkers in the late seventies and early eighties in the context of cell adhesion (Bell, 1978;Bell et al., 1984). Many subsequent works have appeared based on these models, for example to study substrate adhesion (Hammer and Lauffenburger, 1987), strip peeling (Dembo et al., 1988), nanoparticle endocytosis (Decuzzi and Ferrari, 2007), sliding contact (Deshpande et al., 2008), cell nanoindentation (Zhang and Zhang, 2008), cell migration (Sarvestani and Jabbari, 2009), cell spreading (Sun et al., 2009), focal adhesion dynamics (Olberding et al., 2010) and substrate compliance (Huang et al., 2011), among others. Similar chemical bonding models have also been used to describe sticking and sliding friction, see Srinivasan and Walcott (2009). The Bell model is an ODE for the chemical reaction. Assuming separation of chemical and mechanical time scales, this can be simplified into an algebraic equation that describes chemical equilibrium (Evans, 1985). Even though Bell-like models have been combined with contact-induced deformations in several of the works mentioned above, the contact formulations that have been considered in those works are not general continuum mechanical contact formulations as will be considered here. Another related topic is the field of tribochemistry that is concerned with the growth of socalled tribofilms during sliding contact. Chemical evolution models are used to describe the tribofilms in the framework of elementary contact models, e.g. see Andersson et al. (2012) and Ghanbarzadeh et al. (2016). There is also a review article on how mechanical stresses can affect chemical reactions at the molecular scale (Kochhar et al., 2015). But none of these works use general continuum mechanical contact formulations. Chemical contact reactions can be described by a state φ(x, t) that follows from an evolution law of the kindφ = f (φ). Mathematically, they are thus similar to the description of contact ageing (Dieterich, 1978;Rice and Ruina, 1983;Ruina, 1983), contact wear (Strömberg et al., 1996;Strömberg, 1999;Stupkiewicz and Mróz, 1999) and contact debonding (Frémond, 1988;Raous et al., 1999). The thermodynamics is, however, very different.
None of the present chemo-mechanical models is a general contact model accounting for the general contact kinematics, balance laws and constitutive relations. This motivates the development of such a formulation here. To the best of our knowledge, it is the first general thermo-chemo-mechanical contact model that accounts for large-deformation contact and sliding, chemical bonding and debonding reactions, thermal contact and the full coupling of these fields. It is derived consistently from the general contact kinematics, balance laws and thermodynamics, introducing constitutive examples only at the end. Its generality serves as a basis for computational formulations and later extensions such as membrane contact or thermo-electrochemo-mechanical contact.
The novelties of the proposed formulation can be summarized as follows. The formulation • derives a self-contained, fully coupled chemical, mechanical and thermal contact model, • accounts for general non-linear deformations and material behavior, • consistently captures all the coupling present in the interfacial balance laws, • highlights the similarities between chemical, mechanical and thermal contact • obtains the general material-independent constitutive contact relations, • provides several interface material models derived from the 2. law of thermodynamics, • is illustrated by three elementary contact solutions.
Chemical reactions here are restricted to reactions across the contact interface. In this, two assumptions are made: (1) that the reaction rate is much higher than any sliding rate, and (2) that there is no tangential diffusion of bonds on the contact surface. Reactions inside the bodies and on their free surfaces, that can also occur without contact, are not considered in this study.
The remainder of this paper is organized as follows: Sec. 2 presents the generalized continuum kinematics characterizing the mechanical, chemical and thermal contact behavior. The kinematics are required in order to formulate the general balance laws that govern the coupled contact system in Sec. 3. Based on these laws the general constitutive relations are derived in Sec. 4. Sec. 5 then gives several coupled constitutive examples satisfying these relations. In order to illustrate these, Sec. 6 provides three analytical test cases for coupled contact. The paper concludes with Sec. 7.

Continuum contact kinematics
This section introduces the kinematic variables characterizing the mechanical, thermal and chemical interaction between two deforming bodies B k , k = 1, 2, following the established developments of continuum mechanics (Holzapfel, 2000) and contact mechanics (Laursen, 2002;Wriggers, 2006). The novelty here consists in their contrasting juxtaposition.
The primary field variables are the current mass density ρ k = ρ k (X k , t), current surface bonding site density n k = n k (X k , t), current position x k = x k (X k , t), current velocityẋ k =: v k = v k (X k , t) and current temperature T k = T k (X k , t) that are all functions of space and time.
Here, the spatial dependency is expressed through the initial position X k = x k t = 0 , and the dot denotes the material time derivativė The two bodies can come into contact and interact mechanically, thermally and chemically on their common contact surface S c := ∂B 1 ∩ ∂B 2 , see Fig. 1. In case of long-range interaction, like van der Waals adhesion, they may also interact when not in direct contact. In this case, interaction is considered to take place on the (non-touching) surfaces regions S k ⊂ ∂B k . 2 The deformation within each body is characterized by the deformation gradient a.
b. c. from which various strain measures can be derived, such as the Green-Lagrange strain tensor The deformation generally contains elastic and inelastic contributions that lead to the additive strain decomposition 3 An example for an inelastic strain is thermal expansion. From F k also follows the quantity which governs the local volume change dv k = J k dV k (6) at x k ∈ B k between undeformed and deformed configuration. Similarly, the local surface area change da k = J sk dA k at x k ∈ S k is governed by the quantity where det s (...) is the surface determinant on S k , e.g. see Sauer et al. (2019). Time dependent deformation is characterized by the symmetric velocity gradient also known as the rate of deformation tensor. The velocity also gives rise to the identitieṡ where div s (...) is the surface divergence on S k . Mechanical interaction between the two bodies is characterized by the gap vector that is defined for all pairs x 1 ∈ S 1 and x 2 ∈ S 2 . Two cases have to be distinguished: i. Both x 1 and x 2 are material points -a setting that is suitable for describing long-range interactions; or more generally: ii. One of the points is not necessarily a material point -a setting that arises when using a closest point projection suitable for short-range interactions. In the latter case, one surface, say S 1 , is designated the so-called master surface, while the other, S 2 , is designated the so-called slave surface (Hallquist et al., 1985). The master surface is then used to define normal and tangential contact contributions. Therefore the master point x 1 in (12) is determined as the closest point projection of slave point x 2 onto master surface S 1 . Parameterizing S 1 by the curvilinear surface coordinates ξ α , α ∈ {1, 2}, the closest projection point is given by where ξ α p is the value of ξ α that solves the minimum distance problem Here a α := ∂x 1 /∂ξ α (α = 1, 2) are tangent vectors of S 1 , which at x p are denoted a p α , i.e. a p α := a α (ξ α p , t). Likewise n p = n 1 (ξ α p , t) is the surface normal at x p . The closest point projection is illustrated in Fig. 2. The gap can be decomposed into its normal and tangential parts where the latter is zero at x p . Hence, for contact g = g n n p and g n = 0, so thaṫ g =ġ n n p .
Based on (12) and (13), the material time derivative of g can also be written aṡ where, in agreement with definition (1),ẋ 1 = ∂x 1 /∂t| ξ α p = fixed andξ α p = ∂ξ α p /∂t| X 2 = fixed . The latter partg (with summation implied over α and often denoted L v g t ) is equal to the Lie derivative of the tangential gap vector (Wriggers, 2006). Likewise (16) is equal to the Lie derivative of the normal gap vector.g t is zero when ξ α p is fixed (s.t.ξ α p = 0), i.e. when x 1 is a material point as in case i above. In case ii,g t is equal to the relative tangential velocity between the two surfaces. Thus, ξ α p = 0 also denotes the case of tangential sticking, whileξ α p = 0 characterizes tangential sliding. The sliding motion is irreversible and accumulates over time. Combining (16) and (17), one can identify the velocity jump and see that it is the Lie derivative of the total gap vector in case ii. Even during sticking, it can be advantageous (see remark 2.2) to allow for (small) motion that is reversible upon unloading. This leads to the tangential velocity decompositioṅ whereξ α e captures the reversible (elastic) motion during sticking, whileξ α i captures the irreversible (inelastic) motion during sliding. Plugging (20) into (19) yields whereg e :=ġ n n p +ξ α e a p α characterizes normal contact and tangential sticking, whileg i :=ξ α i a p α characterizes tangential sliding. Decomposition (21) is analogous to decomposition (4), making the contact kinematics analogous to the kinematics of the bodies. It contains Eq. (19) for the caseξ α e = 0. Expressions forġ n andξ α p can be found e.g. in Wriggers (2006). Thermal contact is characterized by the temperature jump between surface points x 1 and x 2 , also denoted as the thermal gap (Agelet de Saracibar, 1998), and the contact temperature T c that is associated with an interfacial medium, e.g. a lubricant or wear particles, see Fig. 3. It is assumed that this medium is very thin, such that T c is constant through the thickness of the medium and can only vary along the surface. In order to characterize chemical contact, the bonding state variable φ = φ(x, t) is introduced on the contact surface, as discussed in the following section. It varies between φ = 0 for no bonding and φ = 1 for full bonding, and can be associated with a chemical gap (e.g. defined as 1 − φ). Tab. 1 summarizes the kinematic contact variables and their corresponding kinetic counterparts that will be discussed in later sections.
contact heat influx q c thermal (medium) contact temperature T c contact entropy s c Table 1: Energy-conjugated contact pairs. Remark 2.1: Apart from the tangential sticking constraintξ α p = 0, there is also the normal contact constraint g n ≥ 0. Normal and tangential contact are usually treated separately in most friction algorithms, but there are also unified approaches directly based on the gap vector g, e.g. see Wriggers and Haraldsson (2003) and Duong and Sauer (2019). Such an approach is also taken in the following sections.
Remark 2.2: There are reasons for relaxing the contact constraints g n ≥ 0 andξ α p = 0. One is to use a penalty regularization, which is typically simpler to implement than the exact constraint enforcement. Another is to use the elastic gap g e to capture the elastic deformation of (microscale) surface asperities during contact. In both cases g e is nonzero (but typically small) during contact.
Remark 2.3: The designation into master and slave surfaces introduces a bias in the contact formulation that can affect the accuracy and robustness of computational methods. The bias can be removed if alternating master/slave designations are used, as is done for 3D friction in Jones and Papadopoulos (2006) and Sauer and De Lorenzis (2015).
Remark 2.4: As seen above, the Lie derivativesg andg e only contain the relative changes of the normal and tangential gap. They do not contain the basis changesṅ p andȧ p α . This makes the Lie derivative objective and a suitable quantity for the constitutive modeling (Laursen, 2002), see Sec. 4.1.

Balance laws
This section derives the chemical, mechanical and thermal balance laws for a generally coupled two-body system. The resulting equations turn out to be in the same form as the known relations for chemical reactions (Prigogine, 1961) and thermo-mechanical contact (Laursen, 2002). The novelty here lies in establishing their coupling and highlighting the similarities between the chemical, mechanical and thermal contact equations. The derivation is based on the following three mathematical ingredients. The first is Reynold's transport theorem for volume integrals, It follows from substituting (6) on the left, using the product rule, and then applying (10). Applied to surfaces, Reynold's transport theorem simply adapts to It follows from substituting (7) on the left, using the product rule, and then applying (11). The second ingredient is the divergence theorem, where n k is the outward normal vector of boundary ∂B k . The third ingredient is the localization theorem, that can be equally applied to surface integrals.
Some of the field quantities appearing in the following balance laws can be expressed per mass, bonding site, current volume, current area, reference volume or reference area. The notation distinguishing these options is summarized in Tab Not all combinations are required here: the reaction rate, chemical potential, contact tractions and contact heat flux are only defined on the contact surface, while the body force and heat source are only defined in the bulk.

Conservation of mass and bonding sites
Assuming no mass sources, the mass balance of each body is given by the statement Applying (23) and (26), this leads to the local balance laẇ Due to (10), this ODE is solved by ρ k = ρ 0k /J k , where ρ 0k := ρ k t=0 is the initial mass density.
In order to model chemical bonding, the interacting surfaces are considered to have a certain bonding site density n k (the number of bonding sites per current area) composed of bonded sites and unbonded sites, i.e.
The number of bonding sites is considered to be conserved, i.e.
due to (24) and (26). Due to (11), this ODE is solved by n k = N k /J sk , where N k := n k t=0 is the initial bonding site density. Two cases will be considered in the following, see  In the first case we assume n ub k = 0 such that n k = n b k . Then no further equation is needed for n b k . An example are van-der-Waals interactions described by the Lennard-Jones potential discussed in Sec. 5.3. In the second case, discussed in the remainder of this section, further equations are needed for n b k . Since, in this case, the total number of bonded sites is the same on both surfaces, we must have If the two surfaces are very close, S 1 ≈ S 2 can be replaced by a common reference surface S c . Eq. (32) then implies that n b 1 = n b 2 =: n b due to localization (since (32) is still true for any subregion P c ⊂ S c ). This case is considered in the following. The evolution of n b is described by the chemical reaction where AB denotes the bond, while A and B are its components on the two surfaces prior to bonding, see Fig. 4b. The bonding reaction is characterized by the reaction rate composed of the bonding reaction rate − → R c (forward reaction) and the debonding reaction rate ← − R c (backward reaction). An example for these is discussed in Sec. 5.7. Given R c , the balance law for the bonded sites (considering no surface diffusion) then follows as (Sahu et al., 2017) where v c := v 1 = v 2 is the common velocity required to enable chemical bonding. For the unbonded sites on the two surfaces, we have the two balance laws where still v 1 = v 2 (as long as n b > 0). Summing (36) and (38) leads to (31), and so one equation (for each k) is redundant. For convenience, the non dimensional phase field is introduced, where n c is the reference value for the bonding site density, either picked as n c = n 1 , in case surface 1 is taken as reference, or n c = n 2 , in case surface 2 is taken as reference. 5 Then (36) can be rewritten into using Eq. (31). Eq. (40) is the evolution law for the chemical contact state.
Remark 3.1: As we are focusing on solids, the derivation of (40) assumes no surface mobility (or diffusion) of the bonds. In the context of fluidic membranes, the surface mobility of bonds is usually accounted for, e.g. see Brochard-Wyart and de Gennes (2002) and Freund and Lin (2004).

Remark 3.2:
In the derivation leading up to Eq. (31), v 1 = v 2 has been used. But it suffices to assume that the time scale of chemical reactions is much smaller than the timescale of sliding, such that the two reacting surfaces can be assumed stationary w.r.t. one another.

Momentum balance
The linear momentum balance for the entire two-body system is given by whereb k andt k are prescribed body forces and surface tractions. The latter are prescribed on the Neumann boundary ∂ t B k ⊂ ∂B k that is disjoint from contact surface S k . The bonding events, described by Eq. (40), are not considered to affect the momentum of the system. If the two bodies are cut apart, the additional contact interaction traction t c k needs to be taken into account on the interface (see Fig. 1). The individual momentum balance for every part P k ⊂ B k (k = 1, 2) then reads where t k is the traction on surface ∂P k that contains the cases Applying (23), (25) and (26) leads to the corresponding local form where σ k is the Cauchy stress at x k that is defined by the formula Using (42), Eq. (41) implies This simply states that the net contact forces are in equilibrium. It is the corresponding statement to Eq. (32) for mechanical contact. If the two surfaces touch, i.e. S 1 ≈ S 2 =: S c and da 1 ≈ da 2 , (46) implies that Hence there is no jump in the contact traction. Such a jump only arises if an interface stress (e.g. surface tension) is considered.
The angular momentum balances for the individual bodies and the entire two-body system can be written down analogously to (42) and (41), respectively. As long as no distributed body and surface moments are considered, the angular momentum balance of each body has the well-known consequence σ T k = σ k (Chadwick, 1976), while global angular momentum balance implies analogously to (46). This statement is relevant for separated surfaces (S 1 = S 2 ), but in the case of touching surfaces (S 1 ≈ S 2 ), it leads to the already known traction equivalence (47) at common contact points x 1 = x 2 .
Remark 3.3: In (47) the master surface S 1 is taken as reference surface for defining the reference traction t c , which is commonly done in mechanical contact formulations (Wriggers, 2006). In case of chemical and thermal contact, however, the two possible choices n c = n 1 or n c = n 2 for the reference bonding site density introduced in (39) are maintained in this treatment.

Energy balance
The energy balance for the entire two-body system is where the total energy in the system, is given by the kinetic energy and the internal energy which accounts for the individual energies u k = u k (x k ) in B k and the contact energy U c . U c can describe long range surface interactions, as discussed in Remark 3.6, or it can correspond to the energy of a third medium residing in the contact interface, e.g. a thin film of lubricants or wear particles. In the latter case, assuming a sufficiently thin film, U c can be expressed as the surface integral where u c is the contact energy per bonding site on contact surface S c . 6 Further,r k andq k in (49) denote external heat sources in B k and external heat influxes on ∂ q B k ⊂ ∂B k , respectively. The Neumann boundary ∂ q B k is considered disjoint from contact surface S k . An external heat source on the interface is not required here. As will be seen later, the present setup already accounts for the heat from interfacial friction and interfacial reactions.
If the two bodies are cut apart, the additional contact heat influx q c k needs to be taken into account on the interface (see Fig. 1), leading to the individual energy balance for each body where e k := u k + v k · v k /2 and t k satisfies (43). Further, q k is the heat influx on surface ∂P k that contains the cases q k =q k on ∂ q B k , Introducing the heat flux vector q k at x k ∈ B k that is defined by Stokes formula and using (45), the divergence theorem (25) can be applied to obtain and where D k is the symmetric velocity gradient from (9). Using (23), (26) and (44), the local form of (54) thus becomes With this and Eq. (44), the combined balance statement (49) implieṡ which is the corresponding statement to Eqs. (32) and (46) for thermal contact. Herė follows from (53) for bonding site conservation (30). If the two surfaces are very close or even touch, i.e. S 1 ≈ S 2 =: S c and da 1 ≈ da 2 (i.e. for very thin interfacial media), this implies that (since S c can be considered an arbitrary subregion of the contact surface) due to (47). Eq. (62) states that the energy rateu c and the mechanical contact power t c ·(v 1 −v 2 ) cause the heat influxes q c k on S k . Eq. (62) is equivalent to Eq. (6.63) of Laursen (2002). 7 Remark 3.4: Using (21), the contact power becomes showing that it contains elastic and inelastic contributions associated with sticking and sliding, respectively. The first one is associated with energy stored in the contact interface (defined by (82), later). This energy vanishes for an exact enforcement of sticking (g e = 0). The second contribution causes dissipation, which is also zero in the case of sticking (sinceg i = 0) and in case of frictionless sliding (since t c ⊥g i ). In principle, the second contribution could also contain stored energy. In the context of plasticity such a stored energy occurs for hardening (Rosakis et al., 2000). But this is not considered here.
Remark 3.5: Note, that contrary to the contact tractions governed by (47), the contact heat flux according to (62) can have a jump across the interface. This has also been recently explored by Javili et al. (2014).
Remark 3.6: For long-range interactions between two surfaces, U c can be written as 8 U c := S1 S2 n 1 n 2 u 12 da 2 da 1 , where u 12 = u 12 (x 1 , x 2 ) is an interaction energy defined between points on the two surfaces. For bonding site conservation (30) this leads toU c = S1 S2 n 1 n 2u12 da 2 da 1 , in (60). Localization in the form of (62) is not possible for long-range interaction so that Eq. (60) then remains the only governing equation.

Entropy balance and the second law of thermodynamics
The entropy balance for the entire two-body system can be written as 6 where the total entropy in the system, 7 In Laursen (2002) the contact traction tc and heat fluxes q c k are defined with opposite sign, and (62) is written per unit reference area. Further, uc does not contain chemical energy in Laursen (2002). 8 e.g. by coarse-graining the molecular interactions across the two surfaces (Sauer and Li, 2007b) accounts for the individual entropies s k in B k and the contact entropy s c that is associated with an interfacial medium. Further, η e k and η i k are the external and internal entropy production rates in B k ,q k is the entropy influx on the heat flux boundary ∂ q B k , and η i c is the internal entropy production rate of the interface. An external entropy production rate is not needed for the interface, as long as no heat source is considered in the interface. Like u c in (53), s c and η i c are bond-specific. If the two bodies are cut apart, the additional contact entropy influxq c k needs to be taken into account on the interface, leading to the individual entropy balances (k = 1, 2) whereq k is the entropy influx on surface ∂P k that contains the cases Introducing the entropy fluxq k in body B k defined by theorems (23), (25) and (26) can be applied to (68) to give the local form Plugging this equation into the total entropy balance (66) implies This is the corresponding entropy statement to Eqs. (32), (46) and (60). If the two surfaces touch, i.e. S 1 ≈ S 2 =: S c , it implies that The second law of thermodynamics states that In the absence of chemical contact, (73) then becomes equivalent to Eq. (6.66) of Laursen (2002).

General constitution
This section derives the general constitutive equations for the two bodies and their contact interface -as they follow from the internal energy and the second law of thermodynamics. In general, the internal energy is a function of the mechanical, chemical and thermal state of the system that is characterized by the deformation, the bonding state and the entropy. Introducing the Helmholtz free energy, the thermal state can be characterized by the temperature instead of the entropy. The derivation uses the framework of general irreversible thermodynamics established in chemistry (Onsager, 1931a,b;Prigogine, 1961;de Groot and Mazur, 1984) and mechanics (Coleman and Noll, 1964) following the coupled thermo-mechanical treatment of Laursen (2002) and the coupled chemo-mechanical treatment of Sahu et al. (2017). The novelty here is the extension to coupled thermo-chemo-mechanical contact leading to the establishment of the interfacial thermo-chemo-mechanical energy balance in Eq. (108), which is then discussed in detail.

Thermodynamic potentials
For each B k , the Helmholtz free energy (per mass), is introduced as the chosen thermodynamic potential. Thus, Inserting (59), then leads to For the interface, the Helmholtz free energy (per bonding site) is introduced as where T c is the interface temperature introduced in Sec. 2 (see Fig. 1). Thus, Inserting (62), then leads to In this study, the Helmholtz free energy within B k is considered to take the functional form where E e k is the elastic Green-Lagrange strain tensor introduced by Eq. (4). For the interface, the Helmholtz free energy is considered to take the form where g e = g n n p +ξ α e a p α is the elastic part of the gap as introduced by Eq. (21). From Eqs. (81) and (82) In (84)g e =ġ n n p +ξ α e a p α appears instead ofġ e (cf. Laursen (2002), Eq. (6.72)), since ψ c should be objective and hence not depend on the surface basis, i.e. ∂ψ c /∂n p and ∂ψ c /∂a p α vanish. The following subsections proceed to derive the general constitutive equations for the two bodies and their interaction based on ψ k and ψ c .

Constitutive equations for the two individual bodies
Inserting (71) and (77) into (74.1), gives Since this is true for anyr k , q k , T k , v k andψ k , we can identify the relations for the external entropy source and the entropy flux. We are then left with the well know dissipation inequality which, in view of (4), (83) and is the second Piola-Kirchhoff stress, becomes Here Ψ k := J k ρ k ψ k = ρ k0 ψ k is the Helmholtz free energy per reference volume, and D i k is the inelastic part of D k that is related to E i k via (88). Since (89) is true for any E k and T k , we find the constitutive relations Inserting (83) and (90) back into (77) then leads to the entropy evolution equation Here, the term σ k : D i k + ρ krk can be understood as an effective heat source. The equations in (90) are the classical constitutive equations for thermo-mechanical bodies. They happen to be very similar to the contact equations obtained in the following section.
Remark 4.1: The above derivation considers elastic and inelastic strain rates following from (4), which are analogous counterparts to the elastic and inelastic velocity jump, see (21), required for a general contact description. On the other hand, a decomposition of stress σ k is not considered, which corresponds to a Maxwell-like model for the elastic and inelastic behavior. For other constitutive models, that are not considered here, elastic and inelastic stress contributions can appear and need to be separated.
Remark 4.2: Apart from constitutive equations (90) and the corresponding field equations (44) and (91), one also needs to determine the strain decomposition of (4). Hence another equation is needed, e.g. an evolution law for the inelastic strain E i k . A simple example is thermal expansion, where E i k follows from temperature T k . Likewise, an evolution law for the inelastic gap g i will be needed, e.g. see Wriggers (2006).

Constitutive equations for the contact interface
Next we examine the case that there is a common contact interface S 1 ≈ S 2 = S c . Inserting (73) and (80) into (74.2) and using (86.2), we find Further inserting (21) and (84), and introducing the nominal contact traction t c := J c t c , we find where Ψ c := J c n c ψ c = N c ψ c and where denotes the chemical potential associated with the interface reactions. Here, J c is the reference value for the area stretch, chosen analogously to n c in (39). Since inequality (93) is true for any T k ,g e ,g i and φ, we find the constitutive relations for the contact tractions, for the interfacial entropy, for the reaction rate R c (from (40)), and for the contact heat fluxes. Multiplying by T 1 , T 2 and T c (that are all positive), the last statement can be rewritten into Since this has to be satisfied for any q c k , setting either q c 2 = 0 or q c 1 = 0 yields the two separate conditions (k = 1, 2) for q c k . They are, for example, satisfied for the simple and well-known linear heat transfer law where the constant h k ≥ 0 is the heat transfer coefficient between body B k and the interfacial medium. Introducing the mean influx into B k , and the transfer flux from B 2 to B 1 (see Fig. 3), such that q c 1 = q c m + q c t and q c 2 = q c m − q c t , one can also rewrite inequality (99) into Since this also has to be true for q c m = 0, we find the further condition which is, for example, satisfied by where the constant h ≥ 0 is the heat transfer coefficient between bodies B 1 and B 2 . From (62) we can further find that the mean influx is given by which in view of (21), (79), (84), (95.1) and (96) becomes i.e. the mean heat influx is caused by mechanical dissipation (from friction), chemical dissipation (from reactions), and entropy changes at the interface. The three terms are composed of the conjugated pairs identified in table 1. While the first two terms are always positive (due to (95.2) and (97)) and thus lead to an influx of heat into the bodies, the third term can be positive or negative. It thus allows for a heat flux from the bodies into the interface, where the heat is stored as internal energy, see the example in Sec. 6.2.2. Defining the contact enthalpy per bonding site e c from n c e c := n c u c − t c · g , which is the logical extension of the classical enthalpy definition 9 , one finds based on the Lie derivative of (19). Thus the mean heat influx is proportional to the enthalpy change at constant nominal contact traction. Following the classical definition (Laidler, 1996), the contact reactions can thus be classified according to q c m    > 0 exothermic contact reaction, = 0 isothermic contact reaction, < 0 endothermic contact reaction, in case there is no heat coming from mechanical dissipation (t c ·g i = 0). Corresponding examples are given in Sec. 6.2. (108) is analogous to the bulk equation (91). While (91) describes the entropy evolution in each body, (108) describes the entropy evolution of the contact interface. In case of transfer law (101), this evolution law explicitly follows from (102) and (108) as

Remark 4.3: Interface equation
Like (91), it can be rewritten as an evolution law for the temperature using the entropy-temperature relationship, i.e. (96). The analogy between (91) and (108) is not complete, as (91) contains an explicit heat source, while (108) contains a chemical dissipation term. In principle, an explicit heat source can also be considered on the interface, while a chemical reaction can also be considered in the bulk. This would lead to a complete analogy between (108) and (91). In the absence of chemical reactions, (108) resorts to the thermo-mechanical energy balance found in older works, see Johansson and Klarbring (1993); Strömberg et al. (1996); Oancea and Laursen (1997).

Remark 4.4:
A special case is s c = 0 ∀ t (for example because there is no interfacial medium), see the example in Sec. 6.2.1. In this caseṡ c = 0, such that (108) becomes 2q c m = t c ·g i − n c µ cφ , which is not an evolution law for T c anymore, but just an expression for the energy influx into B 1 and B 2 .
Remark 4.5: According to Eq. (108), equal energy flows from the interface into bodies B 1 and B 2 . The transfer heat flux q c t then accounts for the possibility that the bodies heat up differently. Alternatively, as considered in the works of Agelet de Saracibar (1998), Pantuso et al. (2000) and Xing and Makinouchi (2002), factors can be proposed for the relative contributions going into the two bodies. Pantuso et al. (2000) also account for a loss of the interface heat to other bodies, like an interfacial gas.
Remark 4.6: The nominal traction t c := J c t c is not a physically attained traction. Only the traction t c = n c ∂ψ c /∂g e is. Since n c = n c (g e ), one can also write t c = ∂Ψ c /∂g e , for Ψ c := n c ψ c .
Remark 4.7: Multiplying (108) by the area change J c yields where S c := N c s c and M c := N c µ c are the contact entropy and chemical potential per reference area. Due to Ψ c := N c ψ c , (94) and (96), they follow directly from Remark 4.8: Alternatively to transfer laws for q c 1 and q c 2 , such as (101), one can also propose laws for q c t and q c m . Examples are given by (106) and for some h m ≥ 0, as they satisfy condition (104).

Remark 4.9:
In case q c m = 0 and relation (101) holds, T c explicitly follows from (102) as One can then find the well known relation h −1 = h −1 1 + h −1 2 from (103) and (106), which also has to be true for q c m = 0 in case all h k are constants. Remark 4.10: In case of perfect thermal contact T c = T 1 = T 2 , the fluxes q c 1 and q c 2 become the unknown Lagrange multipliers to the thermal contact constraints T 1 = T c and T 1 = T c . Alternatively, q c t and q c m can be used as the Lagrange multipliers to the constraints T 2 = T 1 and T c = (T 1 + T 2 )/2.

Constitutive equations for non-touching (non-dissipative) interactions
n case of non-touching (long-range) interactions, the two surfaces remain distinct, i.e. S 1 = S 2 . For simplification we consider that the surfaces have the uniform temperatures T 1 and T 2 , and that the space between between S 1 and S 2 has no mass, temperature or entropy, i.e. s c ≡ 0. Further, since the surfaces don't touch, no bonding can take place, i.e. φ ≡ 0. The Helmholtz free energy of the interface then is so thatu 12 =ψ 12 . The Helmholtz free energy is now considered to take the form such thatψ Inserting this into (60) yields Noting that v k and q c k are arbitrary and can be independently taken as zero, we find for the interaction tractions, and for the interaction heat fluxes. At the same time (72) and (74.2) yield If the temperature is constant, one can multiply this by T 1 T 2 and insert (123) to get which is the corresponding statement to (105) for long-range interactions. The traction laws in (122) are identical to those obtained from a variational principle (Sauer and De Lorenzis, 2013).

Constitutive examples for contact
This section lists examples for the preceding constitutive equations for contact and long-range interaction accounting for the coupling between mechanical, thermal and chemical fields. The examples are based on the material parameters defined in Tab. 3. While many of the examples are partially known, their thermo-chemo-mechanical coupling has not yet been explored in detail.

Contact potential
A simple example for the contact potential (per unit reference area) is the quadratic function where 10 where n = n1 is the surface normal of the master body symbol name unit E n , E t normal & tangential contact stiffness N/m 3 C c interfacial contact heat capacity at T 0 J/(K m 2 ) K c bond energy density J/m 2 T 0 reference temperature K t max n , t max t normal & tangential bond strength N/m 2 g 0 , ψ 0 reference distance & reference energy m, J µ 0 , µ static and dynamic friction coefficient 1 η dynamic interface viscosity N s/m 3 h heat transfer coefficient J/(K s m 2 ) C r reaction rate constant 1/(J s m 2 ) is the surface identity on S c . Here, E n and E t denote the normal and tangential contact stiffness (e.g. according to a penalty regularization), C c denotes the contact heat capacity (at T c = T 0 ) and K c denotes the bond energy. (126) corresponds to an extension of the models of Johansson and Klarbring (1993); Strömberg et al. (1996) and Oancea and Laursen (1997) to chemical bonding. The parameters E • , C c and K c are defined here per undeformed surface area. They can be constant or depend on the other contact state variables, i.e. E • = E • (T c , φ), C c = C c (g e , φ) and K c = K c (g e , T c ). If they are constant, (95.1) and (114) yield the contact traction, entropy and chemical potential (per reference area) If E • , C c and K c are not constant, further terms are generated from (95.1) and (114). An example is given in (162).
Remark 5.1: As noted in Remark 3.4, the first term in (126) is zero if the sticking constraint is enforced exactly. In that case, the elastic gap g e is zero, while E c approaches infinity. If there is no mass associated with the contact interface, its heat capacity C c , and hence the second term in (126), is also zero. On the other hand, non-zero g e can be used to capture the deformation of surface asperities during contact (see Remark 2.2), while non-zero C c can be used to capture the heat capacity of trapped wear debris (Johansson and Klarbring, 1993).
Remark 5.2: The last part in (126) corresponds to a classical surface energy. In the unbonded state, (φ = 0) the free surface energy is K c /2.
Remark 5.4: One can change the last term of (126) into such that the minimum energy state is at φ = φ eq . This case implies that chemical equilibrium is a balance of bonding and debonding reactions, like in the model of Bell (1978). An example for this is given in Sec. 5.7. Like K c , the parameter 0 ≤ φ eq ≤ 1 can be a function of g e and T c .
Remark 5.5: The contact traction t c = t c /J c can be decomposed into the contact pressure p c = −n · t c and tangential contact traction t t = i t c (such that t c = −p c n + t t ). From (127)-(129.1) and • := E • /J c thus follow p c = − n g n and t t = t g et , where g n and g et are the normal and tangential parts of g e following from (15).
Remark 5.6: The normal contact contribution in (126) and (129.1) is only active up to a debonding limit. Beyond that it becomes inactive, i.e. by setting E n = 0. This is discussed in the following section.

Adhesion / debonding limit
An adhesion or debonding limit implies a lower bound on the contact pressure p c , i.e.
where the tensile limit (or bond strength) t max n can depend on the contact state, i.e. t max According to this, the bond strength is increasing with φ and decreasing with T c . Once the limit is reached, sudden debonding occurs, resulting in φ = 0. 11 This makes debonding (mechanically) irreversible (unless t max n = 0). Fig. 5a shows a graphical representation of this.  Debonding can also be described in the context of cohesive zone models. An example is the exponential cohesive zone model where t max and g 0 are material parameters that can depend on T c and φ. Eq. (133) is an adaption of the model of Xu and Needleman (1993). It leads to the contact traction according to (95.1) and Remark 4.6. Eq. (134) replaces traction law (129.1) and its limit (131). It is illustrated in Fig. 5b. Traction law (134) is reversible unless t max is a function of φ that drops to zero beyond some g e . Model (133) is similar to adhesion models for non-touching contact discussed next.

Non-touching adhesion
An example for non-touching contact interactions according to Sec. 4.4 is the interaction potential ψ 12 (ḡ k ) = ψ 0 1 45 where ψ 0 and g 0 are constants and where either k = 1 and = 2 or k = 2 and = 1. From (121) and (122) then follows with − ∂ψ 12 ∂x k = ψ 0 g 0 10 45 These expressions are valid for general surface geometries and arbitrarily long-range interactions. For short-range interactions between locally flat surfaces, these expressions can be integrated analytically to give (Sauer and De Lorenzis, 2013) t c k (x k ) = 2π n 1 n 2 g n ψ 12 (g n ) n p , where g n is the normal gap between point x k and surface S and n p is the surface normal of S at x p .
Remark 5.7: The surface interaction potential (135) can be derived from the classical Lennard-Jones potential for volume interactions (Sauer and Li, 2007b;Sauer and Wriggers, 2009).
Remark 5.8: Example (135) is only valid for separated bodies (ḡ k > 0). Formulation (136) can however also be used for other potentials ψ 12 that admit penetrating bodies (with negative distances). Examples, such as a penalty-type contact formulation, are given in Sauer and De Lorenzis (2013). Computationally, (136) leads to the so-called two-half-pass contact algorithm, which is thus a thermodynamically consistent algorithm.
Remark 5.9: Eq. (136) also applies to the Coulomb potential for electrostatic interactions (Sauer and De Lorenzis, 2013). However, in order to account for the full electro-mechanical coupling, the present theory needs to be extended.
Remark 5.10: We note again that for non-touching contact, the tractions t c 1 and t c 2 generally only satisfy global contact equilibrium (46), but not local contact equilibrium (47).
Remark 5.11: Eq. (138) is a pure normal contact model that does not produce local tangential contact forces. Tangential contact forces only arise globally when t c k acts on rough surfaces.

Sticking limit
Similar to the debonding limit (131), a sticking limit implies a bound on the tangential traction where the limit value t max t is generally a function of contact pressure, temperature and bonding state, i.e. t max t = t max t (p c , T c , φ). It is reasonable to assume that it is monotonically increasing with p c and φ and decreasing with T c , e.g. where is a temperature-and bonding state-dependent coefficient of sticking friction based on the constants µ b 0 and µ ub 0 describing the limits for full bonding and full unbonding, respectively. Once the limit is reached, the bonds break (φ = 0) and tangential sliding occurs, as is discussed in the following section. Fig. 6 shows a graphical representation of this. Figure 6: Tangential contact behavior: Sticking, debonding and sliding according to (139)-(141) and (143). The latter two are irreversible.

Sliding models
The simplest viscous friction model satisfying (95.2) is the tangential traction model where η is the positive definite dynamic viscosity tensor and g n is the normal gap. In the isotropic case, η = η1. One of the simplest dry friction models satisfying (95.2) for p c ≥ 0 is the Amontons-Coulomb where µ is the positive definite coefficient tensor for sliding friction. In the isotropic case, µ = µ1. For adhesion, where p c ≥ −t max n , this extends to In the adhesion literature this extension is often attributed to Derjaguin (1934), whereas in soil mechanics it is usually referred to as Mohr-Coulomb's law. Note that µ t max n can be taken as a new constant. The application of (144) to coupled adhesion and friction in the context of nonlinear 3D elasticity has been recently considered by Mergel et al. (2019Mergel et al. ( , 2021. Remark 5.12: A transition model between dry (rate-independent) and viscous (rate-dependent) sliding friction is Stribeck's curve, e.g. see Gelinck and Schipper (2000).
Remark 5.13: The above sliding models can be temperature dependent. An example for a temperature dependent friction coefficient µ = µ(T c ) is given in Laursen (2002). Temperature dependent viscosity models are e.g. discussed in Roelands (1966) and Seeton (2006).
Remark 5.14: The present formulation assumes fixed bonding sites that break during sliding friction. Sliding friction is thus independent from φ. Bonding sites can, however, also be mobile and hence stay intact during sliding. Sliding friction may thus depend on φ. This requires a mobility model for the bonds, e.g. following the work of Freund and Lin (2004).
Remark 5.15: The friction coefficient can be considered dependent on the surface deformation, as has been done by Stupkiewicz (2001).
Remark 5.16: Friction coefficients that are dependent on the sliding velocity and (wear) state have been considered in the framework of rate and state friction models (Dieterich, 1978;Ruina, 1983). Those models are usually based on an additive decomposition of the shear traction t t instead of an additive slip decomposition as is used here.
Remark 5.17: Note that t t causes the mechanical dissipation D mech = t t ·g i that leads to heating of the contact bodies due to Eq. (108). The mechanical dissipation is D mech =g i · ηg i /g n for model (142) and D mech = p cgi · µg i / g i for model (143); see Sec. 6.1 for an example.

Heat transfer
The simplest heat transfer model satisfying condition (105) is the linear relationship already given in (106), where h > 0 is the heat transfer coefficient, also referred to as the thermal contact conductance. Eq. (106) is analogous to the mechanical flux model in (129.1). Similar to the mechanical case, h → ∞ implies T 2 − T 1 → 0. In general, h can be a function of the contact pressure, gap or bonding state, i.e. h = h(p c , g n , φ). Various models have been considered in the past. Those usually consider h to be additively split into a contribution coming from actual contact and contributions coming from radiation and convection across a small contact gap. During actual contact (g n = 0, p c > 0) a simple model is the power law dependency on the contact pressure, where h 0 , h c , H and q are positive constants (Wriggers and Miehe, 1994;Laschet et al., 2004). This model is a simplification of the more detailed model of Song and Yovanovich (1988) that accounts for microscopic contact roughness. Another model for h is the model by Mikić (1974).
A new model has also been proposed recently by Martins et al. (2016). The dependency of (145) on the bonding state can, for example, be taken as i.e. assuming that h 0 and h c increase monotonically with φ. Here h b 0 , h ub 0 , h b c and h ub c are model constants. Out of contact (g n > 0, p c = 0, φ = 0), the heat tranfer depends on the contact gap. Considering small g n , Laschet et al. (2004) propose an exponential decay of h with g n according to where and h gas = k gas (T c ) g n + g 0 (T c ) correspond to the heat transfer across the gap due to radiation and gas convection, respectively.
Here, C trans , c rad , ε 1 and ε 2 are constants, while k gas and g 0 depend on the contact temperature T c . We note that all the above models are consistent with the 2. law as long as h > 0.
Helmholtz free contact energy Ψ c = n c ψ c . 13 Under these conditions the temperature change within each body is governed by the energy balance that follows directly from integrating (59) over the reference volume, writing dV k = H k dA k on the left hand side and applying the divergence theorem on the right. Considering the bulk energy where c k is the heat capacity per unit mass, leads to due to (75) and (90). For quasi-static deformations then followsu k = c k T kṪk /T 0 . Inserting this into (154) then gives where C k := H k ρ k0 c k is the heat capacity per unit contact area satisfying C 1 = C 2 .

Sliding thermodynamics
The first test case illustrates thermo-mechanical coupling by calculating the temperature rise due to mechanical sliding. The test case is illustrated in Fig. 7a: Two blocks initially at mechanical rest and temperature T 0 are subjected to steady sliding with the relative sliding velocity magnitude v := g i . Following Strömberg et al. (1996) and Oancea and Laursen is used with constant and C c . Eq. (158) leads to the entropy given in (129.2). The mean influx, given by (108) and (113), then becomes 13 Jc = 1 leads to Uc = Uc, Sc = Sc, Ψc = Ψc, etc.
where D mech = ηv 2 /g n , according to model (142) and D mech = µpv according to model (143). During stationary sliding (v = const.) the mechanical deformation is time-independent. From (157) thus follows where τ is the time scale of the temperature rise. Integrating this from the initial condition T c (0) = T 0 leads to the temperature rise which is shown in Fig. 7b.

Bonding thermodynamics
The second test case illustrates thermo-chemical coupling by calculating the temperature change due to chemical bonding. The test case is illustrated in Fig. 8a: Two blocks initially unbonded and at temperature T 0 are bonding and changing temperature. Now, the contact energy is used, which, due to (114), leads to the chemical potential and entropy This in turn leads to the internal energy and the bonding ODEφ according to (40), (78) and (150). ODE (165) can only be solved with knowledge about T c = T c (t).

Exothermic bonding
Exothermic bonding occurs for K 2 = 0 (which implies U c = Ψ c ). Considering c r = const., the bonding ODE becomes independent of T c and can be integrated from the initial condition φ(0) = 0 to give φ(t) = 1 − e −t/τ , where τ = n 1 /(c r K 0 ) is the time scale of the exothermic bonding reaction. Solution (166) is shown in Fig. 8b. The mean influx, given in (108) and (113), then becomes It satisfies the exothermic condition q c m > 0 and is shown in Fig. 9a. From (157) thus follows which can be integrated from the initial condition T c (0) = T 0 to give the temperature rise T c (t) = T 0 1 + κ 1 − e −2t/τ , with κ : shown in Fig. 9b for various values of κ. b. Figure 9: Bonding thermodynamics: a. heat influx q c m (t) and b. temperature rise T c (t) during exothermic bonding
Integrating this from the initial condition T c (0) = T 0 gives the temperature drop where the parameter κ must be smaller than unity for the temperature to remain physical. According to (170), this now leads to q c m (t) = − K 2 2τ 1 − κ (1 − κ e −2t/τ ) 2 e −2 t/τ , and satisfies q c m < 0. Fig. 10 shows the energy outflux and temperature drop of the contact interface for the endothermic case. b. Figure 10: Bonding thermodynamics: a. heat influx q c m (t) and b. temperature drop T c (t) during endothermic bonding

Debonding thermodynamics
The last test case illustrates thermo-chemo-mechanical coupling by calculating the temperature change due to mechanical debonding. The test case is illustrated in Fig. 11a: Two blocks initially bonded and at temperature T 0 are pulled apart leading to debonding and rising temperature. Now the contact energy density is used. When the bond breaks, elastic strain energy is converted into surface energy and kinetic energy. The former corresponds to the stored bond energy (and is given by the second term in (174)). If viscosity is present, the latter then transforms into thermal energy. Waiting long enough, the kinetic energy K completely dissipates into heat. The temperature rise can then be calculated from the energy balance. Setting the energy (per contact area) before and after debonding equal, gives where U mech is the mechanical energy at debonding and C tot := C 1 + C 2 + C c is the total heat capacity of the system. Considering linear elasticity gives U mech = 1 2 (t max n ) 2 /k eff , where k eff = (H 1 /E 1 + H 2 /E 2 ) −1 is the effective stiffness of the two-body system based on height H k and Young's modulus E k . This then leads to the temperature change with the positive constants κ := K 0 /(T 0 C tot ) and t 2 0 := T 0 C tot k eff . It is shown in Fig. 11b for various values of κ and t max n /t 0 . Increasing the bond energy (represented by κ) leads to lower  final temperatures, while higher bond strengths t max n lead to higher final temperatures. As a consequence, the final temperature can either be higher or lower than the initial temperature. But note that t max n /t 0 > κ − 1 to ensure T c > 0. This implies that for large κ > 1 higher stresses are needed to break the bond.

Conclusion
This work has presented a unified continuum theory for coupled nonlinear thermo-chemomechanical contact as it follows from the fundamental balance laws and principles of irreversible thermodynamics. It highlights the analogies between the different physical field equations, and it discusses the coupling present in the balance laws and constitutive relations. Of particular importance is (108), the equation for the mean contact heat influx. It identifies how mechanical dissipation, chemical dissipation and interfacial entropy changes lead to interfacial heat generation. This is illustrated by analytically solved contact test cases for steady sliding, exothermic bonding, endothermic bonding and forced bond-breaking. The proposed theory applies to all contact problems characterized by single-field or coupled thermal, chemical and mechanical contact. There are several applications of particular interest to the present authors that are planned to be studied in future work. One is the pressuredependent curing thermodynamics of adhesives (Sain et al., 2018). A second is the study of the bonding thermodynamics of insect and lizard adhesion based on the viscolelastic multiscale adhesion model of Sauer (2010). A third is the local modeling and study of bond strength and failure of osseointegrated implants (Immel et al., 2020). There are also interesting applications that require an extension of the present theory. An example is electro-chemo-mechanical contact interactions in batteries. Therefore the extension to electrical contact is required, e.g. following the framework of Weißenfels and Wriggers (2010). Another example is contact and adhesion of droplets and cells. Therefore the extension to surface stresses, surface mobility of bonding sites, contact angles and appropriate bond reactions is needed, e.g. following the framework of Sauer (2016a) and Sahu et al. (2017).