Janine C Mergel
email: janine.mergel@rwth-aachen.de

Julien Scheibert

Roger A Sauer
email: sauer@aices.rwth-aachen.de

Contact with coupled adhesion and friction: Computational framework, applications, and new insights

Keywords: van der Waals interactions, computational contact mechanics, nonlinear finite element methods, peeling, elastomer contact

. We here develop a corresponding finite element framework that captures 3D finitestrain contact of two deformable bodies. This framework is suitable to investigate sliding friction even under tensile normal loads. First, we demonstrate the capabilities of our finite element model using both 2D and 3D test cases, which range from compliant tapes to structures with high stiffness, and include deformable-rigid and deformable-deformable contact. We then provide new results on the onset of sliding of smooth elastomer-glass interfaces, a setup that couples nonlinear material behavior, adhesion, and large frictional stresses. Our simulations not only agree well with both experimental and theoretical findings, they also provide new insights into the current debate on the shear-induced reduction of the contact area in elastomeric contact.

Introduction

Soft materials like compliant tapes, elastomers, and biological adhesive pads (appearing e.g. in insects and lizards) play a major role in a large variety of dry, solid contact. This kind of contact usually features large adhesive (or tensile) stresses, large frictional stresses, and large deformations, simultaneously. Present computational contact models do not appropriately capture all these three features at once. However, this would be desirable in order to reproduce and interpret a wide range of experimental observations for such systems, including gecko pads [START_REF] Autumn | Evidence for van der Waals adhesion in gecko setae[END_REF], tape peeling [START_REF] De Zotti | Bending to kinetic energy transfer in adhesive peel front microinstability[END_REF], or rubber friction [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF].

Dry friction is often described with the classical Amontons-Coulomb law of friction

F t = µ F n , F n > 0, ( 1
)
where µ is a coefficient of friction relating the sliding friction force F t to the normal load F n . However, in applications dominated by adhesion, this friction force is often found to be proportional to the real contact area, A real , i.e., the total area of small microasperities in actual contact [START_REF] Carpick | Scratching the surface: Fundamental investigations of tribology with atomic force microscopy[END_REF][START_REF] Degrandi-Contraires | Sliding friction at soft micropatterned elastomer interfaces[END_REF][START_REF] Yashima | Normal contact and friction of rubber with model randomly rough surfaces[END_REF][START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF]. In this case,

F t = τ 0 A real , (2)
where τ 0 denotes a frictional shear strength that depends on both materials of the interacting surfaces. It was thus natural to propose the following intermediate law [START_REF] Derjaguin | Molekulartheorie der äußeren Reibung[END_REF][START_REF] Bowden | Mechanism of metallic friction[END_REF]:

F t = µ F n + τ 0 A real , (3)
which interpolates between the two limit cases of Eqs. ( 1) and ( 2). In the following, Eq. ( 3) will be referred to as the "extended" Amontons' law of friction. According to [START_REF] Ruths | Surface forces and nanorheology of molecularly thin films[END_REF], it is suitable to describe the friction force between dry surfaces sliding over each other in the presence of adhesion. As shown experimentally by [START_REF] Homola | Fundamental experimental studies in tribology: The transition from "interfacial" friction of undamaged molecularly smooth surfaces to "normal" friction with wear[END_REF] and others, it depends on the specific application which of these terms has stronger influence. The transition between both terms is also discussed in [START_REF] Berman | Amontons' law at the molecular level[END_REF], [START_REF] Gao | Frictional forces and Amontons' law: From the molecular to the macroscopic scale[END_REF], and [START_REF] Jagota | Adhesion, friction, and compliance of bio-mimetic and bio-inspired structured interfaces[END_REF]. Sometimes using a slightly different notation, relation (3) was considered in the context of microtribology [START_REF] Briscoe | A study of the friction and adhesion of polyethylene-terephthalate monofilaments[END_REF], molecular dynamics [START_REF] Sivebaek | Frictional properties of confined polymers[END_REF][START_REF] Mo | Friction laws at the nanoscale[END_REF], or adhesion and friction of biologic and bio-mimetic systems [START_REF] Zeng | Frictional adhesion of patterned surfaces and implications for gecko and biomimetic systems[END_REF][START_REF] Hill | Orientation angle and the adhesion of single gecko setae[END_REF]. It was further used by [START_REF] Tabor | Friction -the present state of our understanding[END_REF] to state a pressure-dependent, effective friction coefficient, which was then incorporated into a computational model by [START_REF] Wriggers | Finite element formulation of large deformation impact-contact problems with friction[END_REF]. Moreover, Eq. ( 3) is also known as Mohr-Coulomb criterion in soil mechanics. For a comprehensive review of tribological models, in general, we refer to [START_REF] Vakis | Modeling and simulation in tribology across scales: An overview[END_REF]. Note that Eqs.

(1) to (3) refer to the total forces F n and F t applied at the interface. Locally, the normal contact stress may vary between tension and compression within the same macroscopic contact area (see e.g. measurements by [START_REF] Eason | Stress distribution and contact area measurements of a gecko toe using a high-resolution tactile sensor[END_REF] for gecko toes). It is thus desirable to formulate the two macroscopic friction laws (2) and (3) in terms of local contact tractions, resulting in two general continuum models for dry adhesion and friction. This is exactly what was recently proposed in [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF]:

1) Model DI: In this model, the local sliding resistance is constant everywhere along the contact interface, with a value equal to τ 0 . This yields Eq. ( 2) at the macroscale. The local sliding resistance itself is independent of the (generally varying) contact pressure, which in turn is expressed as function of the (small, but non-vanishing) local normal distance between the contacting surfaces. The model is thus called DI, for Distance-Independent. 2) Model EA: This corresponds to a local version of Eq. (3), and is denoted EA for Extended Amontons. Its local sliding resistance is the sum of a constant term (like in model DI) and a term that linearly depends on the (local) normal pressure between the surfaces.

As will be seen, both models are capable of capturing friction even for zero or tensile contact pressures. This capability is an important feature of the class of models investigated here, and is motivated by e.g. soft and compliant bio-adhesive pads, which are able to generate friction forces under tensile normal loads.

The first aim of the present paper is to formulate a computational framework for adhesive-friction models in general, and for the two above-mentioned models in particular. For the sake of self-consistency, we first briefly review models DI and EA as well as the underlying assumptions. We then derive the equations that are necessary to implement these models into a nonlinear 3D finite element (FE) formulation based on large-deformation continuum mechanics. Such formulations go back to [START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems[END_REF], and were extended subsequently to consider e.g. wear [START_REF] Strömberg | Derivation and analysis of a generalized standard model for contact, friction and wear[END_REF], irreversible adhesion [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF], reversible adhesion [START_REF] Sauer | A contact mechanics model for quasi-continua[END_REF], and multiscale contact [START_REF] Wriggers | Multi-scale approach for frictional contact of elastomers on rough rigid surfaces[END_REF]. Regarding non-adhesive 3D frictional contact, important computational advances were made in the context of surface smoothing [START_REF] Padmanabhan | A framework for development of surface smoothing procedures in large deformation frictional contact analysis[END_REF][START_REF] Krstulović-Opara | A C 1 -continuous formulation for 3D finite deformation frictional contact[END_REF], mortar methods [START_REF] Puso | A mortar segment-to-segment frictional contact method for large deformations[END_REF][START_REF] Gitterle | Finite deformation frictional mortar contact using a semi-smooth Newton method with consistent linearization[END_REF][START_REF] Dittmann | Isogeometric analysis and thermomechanical mortar contact problems[END_REF], moving cone formulations [START_REF] Wriggers | The moving friction cone approach for three-dimensional contact simulations[END_REF], isogeometric analysis [START_REF] De Lorenzis | A large deformation frictional contact formulation using NURBSbased isogeometric analysis[END_REF][START_REF] Temizer | Three-dimensional mortar-based frictional contact treatment in isogeometric analysis with NURBS[END_REF], and unbiased friction algorithms [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF]. The latter work is used as a basis for the present formulation. A literature survey on more broadly related models (including e.g. cohesive zone models) is given in [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF].

The second aim of this paper is to provide various examples of application for our framework. The first three especially serve to illustrate the wide range of systems that can be studied using this framework, including various shapes and dimensionalities. We then provide a detailed study of a fourth example, which is currently a matter of scientific debate in the literature: the onset of sliding of Hertzian elastomer-glass contact [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF]Menga et al., 2018[START_REF] Menga | Corrigendum to "Do uniform tangential interfacial stresses enhance adhesion?[END_REF][START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF][START_REF] Khajeh Salehani | Modeling adhesive contacts under mixed-mode loading[END_REF]Papangelo and Ciavarella, 2019;[START_REF] Mcmeeking | The interaction of frictional slip and adhesion for a stiff sphere on a compliant substrate[END_REF][START_REF] Wang | Sphere and cylinder contact mechanics during slip[END_REF][START_REF] Lengiewicz | Finite deformations govern the anisotropic shear-induced area reduction of soft elastic materials[END_REF]. We here present new results that agree well with recent experimental [START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] and theoretical results [START_REF] Chen | Non-slipping JKR model for transversely isotropic materials[END_REF]Papangelo et al., 2019). These results lead us to suggest that the shear-induced reduction of the contact area, discussed in those references, may exist even in the absence of adhesion as well as for compressible materials.

The remainder of this paper is structured as follows. In Section 2 we state our general computational framework, and outline models DI and EA. Section 3 contains the algorithmic treatment of adhesive friction as well as the resulting finite element formulation. In Section 4 we illustrate the validity and applicability of our computational models by discussing 2D and 3D applications ranging from soft and compliant tapes to rather stiff structures. A detailed study of the onset of sliding for smooth elastomer-glass contact follows in Section 5. Section 6 concludes this paper.

Continuum modeling of dry adhesion and friction

The models for adhesive friction discussed here are expressed as functions of the (generally varying) local normal distance between two surfaces. We thus first provide the fundamentals to describe this distance mathematically. For this purpose, we introduce a co-variant description for both surfaces, which will be described in the following. As commonly done in large-deformation continuum mechanics, we use uppercase letters for variables in the reference configuration of a body (denoted B 0 ), and small letters for variables in the current configuration (B).

Contact kinematics

For a certain material point on the contact surface of one body, x k ∈ ∂ c B k (k = 1, 2), we need to determine its closest projection point, x p ∈ ∂ c B ( = 2, 1), that minimizes the distance between x k and the surface ∂ c B of the neighboring body (Fig. 1). To this end, we assume that ∂ c B can be parametrized by two convective coordinates ξ = ξ 1 , ξ 2 such that x p = x (ξ p ). At this (still unknown) point x p , ∂ c B is characterized by its co-variant and contra-variant tangent vectors, a p α and a α p (α = 1, 2), and by its surface normal n p . These vectors are defined as

a p α := ∂x (ξ) ∂ξ α ξ = ξ p , (4)
a α p := a αβ p a p β , a αβ p = a p αβ -1 , a p αβ = a p α • a p β , (5)
n p := a p 1 × a p 2 a p 1 × a p 2 , ( 6
)
where × denotes the cross product, and the definition of a α p in Eq. ( 5) contains a summation over β = 1, 2. The coordinates ξ p of the projection point x p can then be determined by solving the (generally nonlinear) equations [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF]. Once this projection point is found, one can define a normal gap vector, g n , and a scalar normal gap, g n , as

(x p -x k ) • a p α = 0, α = 1, 2; (7) see e.g. Appendix B of
g n := x k -x p , g n := g n • n p . ( 8
)
Figure 1: Closest projection point, xp = x (ξ p ), of x k onto surface ∂cB , and illustration of the tangent vectors a p α (α = 1, 2), surface normal np, and normal gap vector g n .

In order to distinguish between tangential sticking and sliding, we (for now only conceptually) introduce a tangential gap vector, g t , which is decomposed into 1) a reversible ("elastic") part, ∆g e , associated with a non-vanishing, tangential stiffness of the interface during sticking, and 2) an irreversible ("inelastic") part, g s , due to local sliding,

g t = ∆g e + g s . (9)
The concept of an elastic, tangential gap is used for the algorithmic treatment of friction, and explained in detail in Section 3.1.

Modeling of normal tractions

To allow for both tensile (adhesive) and compressive (repulsive) normal tractions at the contact interface between B k and B , we use a model based on an integrated Lennard-Jones (LJ) potential [START_REF] Sauer | A contact mechanics model for quasi-continua[END_REF][START_REF] Sauer | Formulation and analysis of a three-dimensional finite element implementation for adhesive contact at the nanoscale[END_REF]. According to this model, the contact traction (force per area) at point

x k ∈ ∂ c B k , due to B , is given by T n,k = θ k J T n (g n ) n p , T n (g n ) = A H 2πr 3 0 1 45 r 0 g n 9 - 1 3 r 0 g n 3 . ( 10
)
Note that this contact law is formulated with respect to the undeformed reference configuration. The function T n is shown in Fig. 2, together with several characteristic parameters.

Characteristic parameters:

• Equilibrium distance: g eq = r 0 / 6 √ 15

• Max. adhesive traction:

T max = √ 5 A H 9π r 3 0 • Location of -T max : g max = r 0 / 6 √ 5
• Work of adhesion: [START_REF] Sauer | Formulation and analysis of a three-dimensional finite element implementation for adhesive contact at the nanoscale[END_REF] for frictionless adhesion and repulsion as well as characteristic parameters; T 0 = A H / 2πr 3 0 ; the blue dashed line shows a regularization of Tn based on linear extrapolation at gn = geq [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF].

W adh = 3 √ 15 A H 16π r 2 0 Figure 2: Normal contact traction in the model of
We here consider a surface-to-surface contact formulation in which we iterate over both contact surfaces (see also Section 3.2). The index k (or ) thus switches from 1 to 2 (or vice versa). According to Eq. ( 10), T n,k depends on the normal vector and normal gap from Eqs. ( 6) and ( 8) as well as on the following quantities: the Hamaker constant A H , the characteristic length r 0 in the Lennard-Jones potential, the current volume change J around x p ∈ ∂ c B , and a scalar θ k involving the inclination and deformation of the two bodies. In the following we use two assumptions discussed in full detail in [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF]: θ k ≈ 1 and J ≈ J c , where J c is the local stretch of the surface ∂ c B , given by

J c := j a J A , j a := det a p αβ , J A := det A p αβ . (11)
The entries of a p αβ and A p αβ (α, β = 1, 2) are determined from Eq. ( 5), inserting the co-variant tangent vectors from the current and the reference configurations, a p α and A p α = ∂X (ξ)/∂ξ α | ξ = ξ p , respectively. Fig. 2 further shows a regularization of T n based on linear extrapolation for small normal gaps, g n < g reg . This regularization is also shown in the following figures (Figs. 3 and4), and increases the robustness of the computational model by preventing infinitely large slope magnitudes. We apply it in all our examples of Sections 4 and 5, using g reg = g eq .

Due to the strong repulsion in T n (Eq. ( 10) and Fig. 2), the two contacting surfaces are mainly (without regularization: always) separated by a small, but positive, normal gap. This directly affects the definition of the real contact area. Although this area is not explicitly required for our computational framework, it is important for the interpretation of the results. We thus define a gap, g area , to determine whether or not two neighboring surface points belong to the contact area. It is reasonable to set g area somewhere between g eq (below which the normal traction is compressive) and g max (above which tensile tractions decay). We found the sensitivity of the contact area to be large for g area ≈ g eq , and minimal at g max . In this work, we thus consider the real area of contact as that part of the interface that is separated by less than g area = g max .

Modeling of tangential tractions

We now present a general framework for combined adhesion, repulsion, and sliding friction. As specific examples we outline the two models DI and EA mentioned in the introduction. The equations stated in the remaining section, as well as in the algorithm description (Section 3.1), all refer to friction between a surface point x k ∈ ∂ c B k (k = 1, 2) and the neighboring surface ∂ c B ( = 2, 1). For a shorter notation, we omit the indices k and in the following.

In our framework we assume that the vector of tangential tractions, due to sticking and sliding friction, depends on both the normal gap from Eq. ( 8) and the tangential slip vector introduced in Eq. ( 9), i.e. t t = t t (g n , g t ). After defining a non-negative function t slide (g n ) for the (magnitude of the) frictional resistance during sliding, we consider friction laws of the general form t t (g n , g t ) < t slide (g n ) during sticking, = t slide (g n ) during sliding.

In Eq. ( 12) we assume that, for equal g n , the tangential traction required to initiate sliding (after sticking) is equal to the traction in the final sliding state. The motivation, validity, and restrictions of this assumption are addressed in detail in [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF].

Model DI: Distance-independent sliding friction in the contact area

The first model, denoted DI (for Distance-Independent), is the local extension of Eq. ( 2). It includes a constant sliding threshold for normal gaps smaller than a certain cutoff distance g cut :

t slide (g n ) = τ DI , g n ≤ g cut , 0, g n > g cut , (13)
where τ DI > 0 is a constant interface parameter. Note that the adhesive-friction model DI is stated in the current configuration of the bodies, which means that the sliding threshold τ DI refers to the actual area of the contact interface. To overcome the discontinuity in Eq. ( 13) at g n = g cut , for our computational framework we use the following regularization,

t slide (g n ) = τ DI 1 + e kDI(gn-gcut) , (14)
where k DI > 0 is a sufficiently large regularization parameter for the transition between τ DI and zero. Its inverse, 1/k DI , can be regarded as a characteristic decay length. Both the original and the regularized versions of model DI are shown in Fig. 3(a), together with the adhesion model from Section 2.2. here, T 0 = A H / 2πr 3 0 ; the regularization of Tn is discussed in Section 2.2; for εt see Section 3.1; the asterisks mark the distances gmax (where |Tn| = Tmax, yellow), gcut (purple), and 1.5 gcut (light blue). Fig. 3(b) shows the resulting friction law, i.e., the traction-separation law in tangential direction for arbitrary, but fixed, normal distances. In this friction law, we regularize the transition between tangential sticking ( t t < t slide , see Eq. ( 12)) and sliding ( t t = t slide ) by introducing a penalty parameter ε t . This parameter represents the finite tangential stiffness of the interface during sticking, and allows for a small, but reversible, tangential motion (∆g e in Eq. ( 9)). Further details are provided in Section 3.1.

The sliding traction of model DI, τ DI , is independent from the normal gap, and thus from the normal traction. If the same value is chosen for both g cut and g area (see Section 2.2), then the sliding threshold will be equal to τ DI everywhere within the real contact area A real . Under such a condition, model DI yields Eq. (2) at the macroscopic scale, with τ DI = τ 0 . In this work, we always use model DI like this, with g cut = g area = g max . For future comparison with model EA (Section 4), we further introduce a friction parameter2 µ DI := τ DI / T max to relate the sliding traction to the maximum adhesive traction from Fig. 2.

A constant sliding resistance has also been considered in previous studies, e.g. in the context of sliding graphene sheets [START_REF] Deng | Adhesion-dependent negative friction coefficient on chemically modified graphite at the nanoscale[END_REF], adhesive (and small-deformation) contact between a rigid sphere and an elastic half-space (Menga et al., 2018[START_REF] Menga | Corrigendum to "Do uniform tangential interfacial stresses enhance adhesion?[END_REF][START_REF] Wang | Sphere and cylinder contact mechanics during slip[END_REF], or for non-adhesive friction of hyperelastic materials [START_REF] Lengiewicz | Finite deformations govern the anisotropic shear-induced area reduction of soft elastic materials[END_REF]. The strength of model DI is that it combines adhesion and friction with large deformations and arbitrary-shaped geometries. It includes the model of [START_REF] Deng | Adhesion-dependent negative friction coefficient on chemically modified graphite at the nanoscale[END_REF] as a special case with g cut = g eq , while providing an additional regularization to overcome the discontinuity at g cut . Note that (the unregularized version of) model DI is not to be confused with the classical Maugis-Dugdale model. That model includes a constant normal traction for pure debonding, while the blue curve in Fig. 3(a) concerns the tangential traction for sliding friction.

Model EA: Extended Amontons' law in local form

The second adhesive-friction law, model EA (for Extended Amontons) is intended as a local form of Eq. (3). Unlike model DI, this model explicitly depends on the (local) normal contact tractions of our adhesion model, and thus on T n (g n ) from Eq. ( 10). Let us first choose a cutoff distance somewhere between the equilibrium distance g eq of T n and the location g max of -T max (see Fig. 2):

g cut = s cut g max + (1 -s cut ) g eq , s cut ∈ [0, 1]. ( 15
)
We then consider a sliding threshold that is proportional to the shifted curve T n (g n ) + |T n (g cut )|, where the function value T n (g cut ) is either zero or negative (because g cut ≥ g eq ):

T slide (g n ) =    µ EA J c T n (g n ) -T n (g cut ) , g n < g cut , 0, g n ≥ g cut . (16)
The scalar µ EA > 0 denotes the friction parameter for model EA (see also footnote 2). The first term in Eq. ( 16) (top) expresses the dependence on the normal traction that is formulated, by choice, w.r.t. the reference configuration. If desired, both the normal and tangential tractions can be mapped to the current configuration using t n da = T n dA. The second term corresponds to an additional constant representing the effect of adhesion. Since g cut ≤ g max (Eq. ( 15)), the sliding resistance in Eq. ( 16) is always positive. Fig. 4 illustrates model EA for the parameter values s cut = 0, 0.5, and 1. In analogy to Fig. 3, the left-hand side of Fig. 4 depicts the relation between the normal gap, g n , and the sliding threshold, T slide , while the right-hand side shows the resulting friction laws. Note that the function T slide (g n ) is C 1 -continuous only if s cut = 1 (Fig. 4(e)). In any other case, a kink occurs at g n = g cut , which requires a special algorithmic treatment in the solution strategy (Section 3.3).

For g cut = g eq (i.e., s cut = 0, see Figs. 4(a) and 4(b)), model EA provides a frictional resistance (i.e., T slide > 0) only for positive, and thus compressive, normal tractions. This corresponds to the local equivalent of the classical Coulomb-Amontons law for non-adhesive contact (Eq. ( 1)). It also follows the approach of many cohesive zone models [START_REF] Chaboche | Numerical analysis of composite systems by using interphase/interface models[END_REF][START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF][START_REF] Del Piero | A unified model for adhesive interfaces with damage, viscosity, and friction[END_REF], which consider both local tension and compression, but include sliding friction only where the normal traction is compressive. Note that, in particular for adhesion-dominated setups the model behavior may become very sensitive in this case. This is demonstrated in [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF] for the setup investigated in Section 4.1. It also agrees with the finding stated in Section 2.2, namely that the real contact area is very sensitive to g area for g area ≈ g eq .

As soon as g cut > g eq , a frictional resistance can arise also in regions where the normal traction is slightly tensile (up to the value T n (g cut )). Note that, since in the present work we use g area = g max (Section 2.2), and since g cut ≤ g max , the tangential tractions are always zero outside the real contact area.

Computational framework

This section contains the algorithmic treatment of adhesive friction and its implementation into a nonlinear finite element code. We also address the validity and restrictions of our framework.

Algorithmic treatment of adhesive friction

In computational friction algorithms, sticking and sliding are often realized by means of a predictorcorrector approach similar to that used for elastoplasticity, see, e.g., [START_REF] Simo | Computational Inelasticity[END_REF]: One first predicts tangential sticking, checks whether the criterion for sliding is satisfied, and if so, applies a return map to determine the sliding distance.

The following algorithm is formulated in the current configuration, i.e., in terms of the variables t t and t slide . However, as shown in [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF], it is also possible to analogously state the algorithm with respect to T t and T slide , the variables with which model EA is formulated in Section 2.3.2. Let us first define the sliding criterion

f s (t t , t slide ) := t t -t slide , (17)
which must satisfy [START_REF] Wriggers | Computational Contact Mechanics[END_REF] f s < 0 for sticking, = 0 for sliding. Figure 4: Model EA: Extended Amontons' law in local form illustrated for µ EA = 0.2, J c ≡ 1, and three different values of scut; here, T 0 = A H / 2πr 3 0 ; the regularization of Tn is discussed in Section 2.2; for εt see Section 3.1; the colored asterisks mark the distances 0.5 r 0 (green), geq (where Tn = 0, orange), and gmax (yellow).

These two cases are illustrated in Fig. 5(a) for model DI. Here, the colored surface marks f s = 0, for which t t must be equal to the sliding threshold t slide from Eq. ( 14). As a consequence, this surface corresponds to nothing else but t slide revolved around the g n -axis, cf. During sliding we must enforce f s (t t , t slide ) = 0, and determine both the direction and the magnitude of the resulting tangential gap. Since this projection is carried out at a certain (pseudo-)time step, the normal gap from Eq. ( 8), g n , is considered to be arbitrary but fixed. As shown in [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF], both models EA and DI fulfill convexity (i) of the domain of feasible tractions t t satisfying f s (t t , t slide ) ≤ 0, and (ii) of the function f s within this domain. These properties ensure that the mapping performed as a corrective step is unique. Let us now define the non-negative dissipation during sliding,

D s (t t ; Lg s ) := t t • Lg s , (19)
which represents the energy loss per time and area. The term

Lg s := ξα s a s α , a s α := ∂x (ξ) ∂ξ α ξ = ξ s , (20)
denotes the relative slip velocity, i.e., the Lie derivative of the slip vector g s associated with sliding (Eq. ( 9), see also [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF]). Note that the dot denotes the material time derivative of ξ α s (keeping material point X p fixed). According to the principle of maximum dissipation [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF], for a given slip velocity Lg s , the actual (physically true) traction t t maximizes the dissipation among all feasible tractions t * t . By keeping Lg s arbitrary but fixed, and by introducing the Lagrange multiplier γ ≥ 0, we formulate the following constrained minimization problem

∂L t * t , γ; Lg s ∂t * t t * t = tt = 0 (21) with Lagrangian L t * t , γ; Lg s := -D s t * t ; Lg s + γ • f s t * t , t slide , (22)
and the Karush-Kuhn-Tucker (KKT) conditions for optimality,

f s ≤ 0, γ ≥ 0, f s • γ = 0. ( 23
)
From Eq. ( 21) we then obtain the evolution law for the slip velocity during sliding

Lg s = γ n t , n t := ∂f s ∂t t = t t t t . (24)
According to Eq. ( 24), the multiplier γ corresponds to the magnitude of the tangential slip velocity. Since during sliding f s = 0, the traction satisfies

t t = t slide n t , (25)
as required. The time derivatives ξα s in Lg s can finally be determined by inserting Eq. ( 24) into Eq. ( 20), and contracting the result with the tangent vectors a α s , giving

ξα s = γ n t • a α s , α = 1, 2. ( 26
)
The following procedure corresponds to the friction algorithm proposed by [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF], but extended to a general sliding threshold t slide . For the derivation of the resulting equations we refer to the above-mentioned paper. The key idea is to decompose the convective coordinates, ξ α p , of the projection point into components that are either related to irreversible sliding, ξ α s , or to reversible ("elastic") sticking, ∆ξ α e (cf. Eq. ( 9)). As described in detail in the following, for sliding the traction is determined such that it satisfies Eq. ( 18); see also Fig. 5.

Given all quantities at a pseudo-time step t n , one proceeds to the next step t n+1 as follows:

1) Assume tangential sticking and compute a corresponding trial traction (Fig. 5)

t n+1 trial = ε t x n+1 ξ n+1 p -x n+1 (ξ n s ) . (27)
The large, but finite, penalty parameter ε t is used for regularization, i.e., to allow for a small, reversible tangential gap, given by x n+1 ξ n+1 p -x n+1 (ξ n s ). ε t is visualized in Fig. 3(b) for model DI, and e.g. in Fig. 4(d) for model EA.

2) Insert t n+1 trial into the criterion

f s t n+1 trial , t n+1 slide = t n+1 trial -t n+1 slide ( 28
)
to check for sticking or sliding: (a) If f s < 0, the point x p is sticking tangentially; thus

t n+1 t = t n+1 trial , ξ α n+1 s = ξ α n s , α = 1, 2. ( 29
)
(b) If f s ≥ 0, the point x p is sliding, and an additional correction is required to determine the actual traction t n+1 t such that it satisfies f s = 0. Regarding Fig. 5, in this case the trial traction lies outside the domain enclosed by the colored surface, and must be mapped back onto this surface. First compute the incremental update of the Lagrange multiplier in Eq. ( 24),

∆γ n+1 = f s t n+1 trial , t n+1 slide ε t , (30)
and then determine the updated coordinates ξ α n+1 s (α = 1, 2) associated with sliding as well as the corresponding tangential traction,

ξ α n+1 s = ξ α n s + ∆γ n+1 n n+1 t • a α n+1 p , n n+1 t = t n+1 trial t n+1 trial , (31)
t n+1 t = t n+1 trial -ε t ∆γ n+1 n n+1 t . (32)
Eq. ( 31) is based on the approximation a α n+1 s ≈ a α n+1 p (see [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF])3 , and on the discretization of Eq. ( 26) by means of the implicit Euler method. The traction t n+1 t in Eq. (32) satisfies Eq. ( 25). In the case of frictionless sliding (i.e., t n+1 slide = 0), one simply sets t n+1 t = 0 and ξ α n+1 s = ξ α n+1 p .

Contact contributions to the finite element equations

In order to discretize the corresponding model equations in space, we use a standard Galerkin finite element (FE) method. For the sake of brevity we only discuss the contact contributions here. To this end, we introduce a global FE contact force vector f c (u), which depends on the global displacement vector u for all finite element nodes. Following the sign convention by [START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF], we now decompose f c (u) into the normal force f n due to adhesion and repulsion and the tangential force f t due to friction,

f c (u) := f n (u) -f t (u). ( 33
)
The vectors f n and f t are obtained by assembling the contributions from each finite element on the contact surfaces. For contact between two deformable bodies, we consider the two-half-pass approach by Sauer andDe Lorenzis (2013, 2015). We thus iterate over the elements of both contact surfaces (k = 1, 2), and compute the elemental contribution to the normal forces from

f e n,k = - Γ e c0k N T k T n,k dA k . (34)
Here, Γ e c0k ⊂ ∂ c B h 0k , T n,k is given by Eq. ( 10), and the array N k contains the nodal shape functions associated with the current element. The tangential forces are determined similarly,

f e t,k = - Γ e ck N T k t t,k da k = - Γ e c0k N T k T t,k dA k , (35)
depending on whether the model is defined in the current (model DI, Eq. ( 14)) or the reference configuration (model EA, Eq. ( 16)).

Since the vector f c is nonlinear with respect to u, we linearize it to solve the governing equations with Newton's method. To this end, we derive the tangent matrix associated with f c , denoted k c := ∂f c /∂u. Like the global contact forces f n and f t , k c is assembled from the contributions of each single surface element, see Appendix A. Finally, we evaluate the integrals in Eqs. ( 34) and ( 35) and in Appendix A by means of Gaussian quadrature.

Active set strategy

As seen in Fig. 4(a), 4(c), and 4(e), the curve T slide (g n ) of model EA is C 1 -continuous at the cutoff distance g cut only if g cut = g max (i.e., s cut = 1). Otherwise, it is C 0 -continuous at g cut . Moreover, we must generally distinguish between discrete states of contact like tangential sticking and sliding. Following [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF], we here consider an active set strategy. To this end, we slightly modify the two active sets that are often used for unilateral and frictional contact: (i) "(compressive) contact vs. no contact" and (ii) "sticking vs. sliding". Since the function T n (g n ) in Eq. ( 10) is continuous anyway (and thus does not require any active set at all), we simply replace active set (i) by "frictional vs. frictionless sliding". Note that it is also reasonable to define certain cutoff values for which both the normal and tangential contact stresses are considered to be negligible. This on one hand increases the efficiency, and on the other hand circumvents computational problems in the sliding algorithm or closest point projection (due to large separations of the bodies).

For changes in the (discrete) contact states, the FE force vector f c (u) actually becomes non-differentiable, cf. Appendix A.2.1 and Appendix A.2.2. If the applied load step is too large, or if large parts of the contact surface change their state at once, this may lead to a phenomenon called cycling, or also zig-zagging [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. That is, the solution algorithm alternates between two states in consecutive iterations. This can be generally overcome e.g. by means of semi-smooth Newton methods. Instead, we here pursue a strategy that is based on a temporary freezing of the active sets; see also [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. Accordingly, we keep the active sets fixed at the beginning of the iterations, until the residual falls below a certain tolerance. We then assume that the changes in the state are sufficiently small to prevent cycling. Further improvement could be made by adjusting the load step adaptively when the body starts or stops (full) sliding, in particular when it changes its direction of sliding.

A special (and careful) treatment may be required for model EA with s cut = 0, i.e., g cut = g eq (see also the end of Section 2.3.2). In this specific case, the switch from a non-zero to zero frictional resistance is directly located at the equilibrium distance, g eq , where the normal traction is zero (Fig. 4(a)). Thus, for setups in which large parts of the contact area are close to this distance (like those in Sections 4.1 and 4.3), the tangential traction can then alternate back and forth between the states of frictional and frictionless sliding. This may either require unnecessarily small load steps, or even lead to a failure of the freezing strategy. As demonstrated in [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF], however, this issue is successfully overcome by choosing a cutoff distance g cut that is slightly smaller than g eq . For such a case we recommend considering Eq. ( 15) with a parameter s cut that is negative but close to zero (e.g., s cut = -0.001). Thus, the contact pressure must first exceed a certain threshold to generate a frictional resistance.

Validity, restrictions, and possible extensions of our models

Let us now discuss the validity and restrictions of our models, in particular regarding computational aspects. For physical assumptions see [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF]. The general properties of the applied two-halfpass algorithm are discussed in Sauer andDe Lorenzis (2013, 2015). A possible extension to incorporate sticking friction larger than the sliding resistance t slide or T slide is briefly addressed in [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF].

Jump-off and jump-to contact

Strong adhesion combined with soft materials may cause sudden jump-off or jump-to contact, see e.g. [START_REF] Chaboche | Interface debonding models: a viscous regularization with a limited rate dependency[END_REF] and [START_REF] Sauer | An Atomic Interaction based Continuum Model for Computational Multiscale Contact Mechanics[END_REF]. In that case the (de-)bonding process inherently becomes physically unstable, and the assumption of quasi-static conditions (as considered in this work) is not valid anymore. According to [START_REF] Raous | Interface models coupling adhesion and friction[END_REF], this is caused by the non-convexity and the softening character (i.e., a decreasing slope with increasing distance) of traction-separation laws like that in Fig. 2. A similar effect can occur also for sliding friction [START_REF] Mróz | Contact friction models and stability problems[END_REF].

The issues mentioned above can be generally avoided by solving the model equations with an arc length or continuation method [START_REF] Sauer | An Atomic Interaction based Continuum Model for Computational Multiscale Contact Mechanics[END_REF]. Another possibility would be to account for a viscous (velocitydependent) regularization as done in several cohesive zone models involving friction [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF][START_REF] Chaboche | Interface debonding models: a viscous regularization with a limited rate dependency[END_REF][START_REF] Del Piero | A unified model for adhesive interfaces with damage, viscosity, and friction[END_REF]. This approach is reasonable because viscous effects often occur not only in the bulk material but also at the interface. Besides, it would also be possible to account for inertial effects and to consider the problem to be dynamic. Since the study of dynamic and viscous effects lies outside the scope of this paper, in the following application examples we focus on the attachment or detachment process itself, but not on sudden jump-to or jump-off contact, respectively.

Going from the nano-to the macroscale

As shown in Fig. 2, the normal contact traction T n from the adhesion model by [START_REF] Sauer | Formulation and analysis of a three-dimensional finite element implementation for adhesive contact at the nanoscale[END_REF] decays rapidly for increasing distances g n . Depending on the involved materials, T n typically acts only within the range of several (tens of) nanometers. In addition, if g n approaches zero, the traction increases to infinity while its slope approaches minus infinity. These issues can cause numerical problems on one hand, and require a fine spatial resolution (finite element size and applied displacements) on the other hand.

We thus recommend the following strategy to avoid such problems. First, a sufficiently accurate surface discretization is crucial to reduce artificial oscillations [START_REF] Sauer | Enriched contact finite elements for stable peeling computations[END_REF], which could cause convergence problems in the numerical solution procedure. In Sections 4 and 5 we thus use specially surface-enriched finite elements based either on Hermite [START_REF] Sauer | Enriched contact finite elements for stable peeling computations[END_REF] or Non-Uniform Rational B-Spline (NURBS) shape functions (Corbett andSauer, 2014, 2015). Second, we recommend regularizing the normal tractions as discussed in Section 2.2. Third, it is generally possible to calibrate the model parameters from Section 2.2 (A H , T max , and W adh ) such that they match experimental data. This leads to a larger length parameter r 0 , which automatically regularizes the curve T n in Eq. (10). For further comments and restrictions see [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF].

Application examples

With the three following examples we illustrate the properties and qualitative behavior of our two models, and demonstrate the generality of our method. We consider not only different dimensions (2D and 3D), but also different contact types (rigid/deformable and deformable/deformable). A fourth, and more detailed, example follows in Section 5.

2D peeling of a thick strip from a rigid surface

We first consider a massless, beam-like strip that initially adheres to a flat and rigid substrate (Fig. 6). Along this interface, the normal traction is initially zero, so that the gap between the strip and the substrate is equal to the equilibrium distance (g n = g eq ) from Fig. 2. The strip is then peeled off the substrate, moving the midpoint of its right boundary by the vertical displacement u. Such a setup is representative e.g. for tape peeling. To avoid shear locking due to bending, we discretize the bulk of the strip with 320 × 8 quadratic (Q2) finite elements. The contact surface is further enriched by quadratic NURBS shape functions (Q2N2.1 elements) [START_REF] Corbett | NURBS-enriched contact finite elements[END_REF]. In this example, as well as in the following subsection, we consider the Neo-Hookean material model stated in [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF]. For the geometry, material, and contact parameters see Tab. 1. All results, shown here and in the following, are normalized by Young's modulus E and the unit length L 0 .

E ν A H r 0 L 0 2 GPa 0.2 10 -19 J 0.4 nm 1 nm
Table 1: 2D strip peeling: Model parameters.

4.1.1. Contact tractions Fig. 6 illustrates the peeling of the strip for model EA, g cut = g max , and a small friction parameter (µ EA = 0.01). Figs. 7(a) and 7(b) show the normal and tangential contact tractions directly before and after the onset of peeling (u = 5 L 0 and 10 L 0 , respectively), as well as for two of the configurations shown in Fig. 6. As seen in Fig. 7(a), for small u (like 5 L 0 ), adhesion induces large tensile normal tractions at the right boundary of the adhesive zone. When these tractions reach T max (at u ≈ 7 L 0 ), first detachment occurs, and a peeling front nucleates. This front then propagates to the left, visible as a sharp peak in Fig. 7(a). Some distance ahead of the peeling front, the normal tractions are repulsive due to the finite bending stiffness of the strip. Together with the peeling front, also a sliding front nucleates and propagates to the left (Fig. 7(b)). This front marks the boundary between the region still sticking to the substrate and the region already sliding. Note that a sliding region ahead of the peeling front is indeed observed in experiments (see e.g. [START_REF] Newby | Effect of interfacial slippage on viscoelastic adhesion[END_REF]). When the sliding front has reached the left boundary of the strip, the (remaining part of the) contact interface is fully sliding. As expected, for model EA the tangential traction distribution looks similar to the normal traction distribution within the sliding region. When repeating the simulation from Fig. 7 with model DI, a friction parameter of µ DI = µ EA = 0.01, and the same cutoff distance (g cut = g max ), we find the following behavior (not shown): The normal tractions are almost identical to those from model EA, while the tangential tractions remain constant at τ DI = µ DI T max within the entire sliding region. In addition, the sliding front evolves very similarly in both models.

Peeling forces

Let us now examine the horizontal and vertical pull-off forces for different friction parameters, focusing on model EA (Figs. 8(a) and 8(b)).

With increasing friction parameter, large parts of the strip initially remain sticking to the substrate, leading to both large horizontal and vertical forces. As illustrated in Figs. 8(a) and 8(b), the force curves slightly change for a different value for g cut (sufficiently larger than g eq , see Section 2.3.2). In analogy to the contact tractions, using model DI instead (with µ DI = µ EA and g cut = g max ) does not cause any significant changes in the forces. The good agreement between models EA and DI is generally expected if (i) the surface stretches are small (such that J ck , J c ≈ 1, and thus t t,k ≈ T t,k ), (ii) the contacting interfaces are mainly separated by the equilibrium distance g eq , and (iii) if the model parameters are chosen as

g eq < g cut ≤ g max , µ DI := µ EA |T n (g cut )| T max . ( 36
)
Note that models EA and DI will behave differently from each other if major parts of the contact area are in tension or compression, or if the surface stretch has a strong influence. This is discussed in the following two examples. Further deviation may be caused by the regularization parameter k DI of Eq. ( 14). For model DI the regularized traction law in Eq. ( 14) is C 1 -continuous for arbitrary cutoff distances g cut . In contrast, Eq. ( 16) for model EA has a kink at g cut if g cut = g max (see e.g. Fig. 4(c)). Nevertheless, due to our active set and freezing strategies (Section 3.3), our algorithm works robustly for all investigated cases.

To understand both the shape of the force curves in Figs. 8(a) and 8(b) and the evolution of their maximum when varying the friction parameter µ EA , we next examine the evolution of the peeling and sliding fronts during peeling. Fig. 8(c) hence shows the nominal length L peel of the part that is peeled off the substrate already. It further illustrates the nominal sliding length, L slide , which we define as the initial contact length minus the nominal length of the currently sticking part. The sliding front starts propagating directly after imposing a small, vertical displacement. Larger friction parameters result in later propagation, because the local tangential tractions must be larger to overcome the sliding threshold, T slide (Eq. ( 16)). In contrast, the peeling front nucleates (i.e., L peel > 0) only after a finite vertical displacement, u ≈ 20 L 0 , where the precise value slightly varies with µ EA (Figs. 8(b) and 8(c)). At the nucleation of the peeling front, the vertical force decreases afterwards for small µ EA , but continues to increase (although with a smaller slope) for the two largest values of µ EA .

In the range of µ EA considered here, we can distinguish between two types of evolution for L peel : (i) For µ EA 0.02, the sliding front reaches the boundary of the strip (L slide = 150 L 0 ) before the peeling front nucleates (L peel > 0). In this case, L peel increases almost linearly, and the strip detaches slightly earlier for increasing µ EA . (ii) For µ EA 0.02, the sliding length has not reached the limit of 150 L 0 before the peeling front nucleates. In this case, L peel initially evolves independently from µ EA . Then, when the entire strip starts sliding, L peel suddenly propagates with a different slope, which abruptly causes the horizontal and vertical forces to decrease again (Figs. 8(a) and 8(b)).

In the field of adhesive peeling, one is often interested in the fracture energy that is required to separate a unit area from the interface [START_REF] Creton | Fracture and adhesion of soft materials: A review[END_REF]. To evaluate this energy, we first compute the difference between the external work done during peeling and the internal energy due to elastic deformation ( * per out-of-plane width): ∆Π * = Π * ext -Π * int . The fracture energy Γ * , shown in Fig. 8(d), then corresponds to the derivative of ∆Π * with respect to the peeling length, Γ * = ∂(∆Π * )/∂L peel . Again, the behavior is found to be different for friction parameters smaller and larger than about 0.02. Below this value, Γ * is essentially independent from the peeling length. It is minimal for the frictionless case, for which Γ * approaches W adh from Fig. 2. This is expected because there is no dissipation due to friction, so that the fracture energy corresponds to the work of adhesion per unit area. For increasing values of µ EA , Γ * then increases, because an additional energy is necessary to overcome frictional dissipation within the sliding region, as discussed in [START_REF] Lu | Effect of interfacial slippage in peel test: Theoretical model[END_REF]. For friction parameters beyond 0.02, the fracture energy strongly increases (by a factor of 2 to 5) at full sliding (L slide = 150 L 0 ). Afterwards, large sliding distances are associated with each increment of u, inducing a much larger frictional dissipation, and thus fracture energy.

Our results suggest that the evolution of the peeling forces and the fracture energy can significantly depend on both the friction parameter and the finite size of the strip. For the two highest friction parameters, the forces evolve very similarly at the beginning, because the sliding region is still small compared to the initial contact length. This behavior is more related to the behavior expected for peeling of adhesive tapes, for which the friction parameter is probably considerably larger. Interestingly, our results in Fig. 8(d) predict that for small peeling lengths (about 20 -60 L 0 ), the fracture energy reaches its maximum for a finite value of µ EA ≈ 0.02. Such a result is potentially of practical interest for the design of optimal adhesives. This could also be tested experimentally, although we anticipate that the viscoelastic nature of real adhesives, not accounted for here, may affect the fracture behavior significantly.

2D contact of two deformable cylinders

In this second example, we consider friction between two Hertzian solids. This setup is particularly important for the field of tribology, because such solids are frequently used to represent either micro-contact within rough surfaces or single contact between smooth solids. Although Hertzian geometries have been investigated in numerous studies (see [START_REF] Vakis | Modeling and simulation in tribology across scales: An overview[END_REF]), their combination with adhesion, friction, finite deformations, and hyperelasticity (as done here) is much rarer.

We consider two identical, deformable cylinders sliding along each other in plane strain (see Fig. 9(a) and Tab. 2). Their vertical positions are fixed such that their nominal overlap is 25 % of their radius, R. We then impose a horizontal displacement u from 0 to 2 R. For each cylinder we use approximately 1, 500 elements and Hermite enrichment [START_REF] Sauer | Enriched contact finite elements for stable peeling computations[END_REF] at the contact surfaces. As discussed in Section 2.2, the volume change at the surfaces is approximated by J ≈ J c in Eq. ( 10). Since in this example (and for our choice of parameters) J c has only little influence [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF], we proceed with J c ≈ 1. Figs. 9(b) to 9(d) show the two cylinders sliding across each other for three different adhesion/friction models. Although the (two-half-pass) friction algorithm does not balance the contact traction explicitly [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF], for each deformation state, the stress fields of both cylinders are pointsymmetric up to machine precision (not shown).

E ν A H r 0 R 1 GPa 0.3 2.54 • 10 -20 J 0.4 nm 40 nm
Fig. 9(b) corresponds to the frictionless, i.e. purely adhesive, case. The corresponding horizontal and vertical forces (per out-of-plane width), acting on the lower cylinder, are shown as black curves in Figs. 10(a) and 10(b). Even in the frictionless case, there exists a finite horizontal force between the cylinders, caused by the geometrical setup. This force is positive (negative) for u < R (u > R), where contact occurs on the trailing (leading) sides of the two asperities (see also Fig. 9(a)). Thus, this setup cannot be reduced to an shapes of the two corresponding curves are rather similar, suggesting that the evolution of the horizontal force mainly results from the sum of the tractions, rather than from their distribution, which differs in the two models. The main deviation can be seen just before detachment of the cylinders, where model DI systematically predicts larger horizontal forces as well as a more abrupt force drop. The results stated above should be put more into perspective. If the cylinders were driven by a horizontal load instead of a displacement, a snap-through instability would occur at the maximum of the F *

x curves, marked with asterisks in Figs. 10(a) to 10(c). Note that in real tribological tests, finite stiffness springs are used to drive the system (see also [START_REF] Papangelo | On the effect of the loading apparatus stiffness on the equilibrium and stability of soft adhesive contacts under shear loads[END_REF]), so that snap-through may occur even when the springs themselves are displacement-driven. Thus, in practice the duration of contact may be much smaller than predicted in Figs. 10(a) to 10(c). Finally note that increasing friction reduces both the possibility of negative forces and the tendency for a snap-through before the (midpoints of the) cylinders have moved across each other (u = R).

The interaction between two cylinders, as considered here, has recently gained renewed attention in the context of wear modeling. In particular, Molinari and co-workers performed atomistic simulations to investigate the conditions for formation of debris when two circular asperities collide tangentially. They found that, depending on the size of the formed junction, such asperities can either deform plastically, or create debris via crack propagation [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF]. In the latter case, a crack nucleates at the trailing corner of each asperity, where it connects to its underlying bulk (Fig. 9(a)). While the contact geometry is similar to ours, the constitutive laws used by [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF] are considerably different. They especially consider a perfectly adhering interface, at which particles coming into contact remain attached to each other as they would be in the bulk. We suggest that our friction models could be used to investigate the effect of a contact interface with a more realistic behavior, including finite adhesion and friction. Note that our simulation setup could be readily modified to include also an elastic bulk below the asperity, as done in [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF].

To illustrate this, Fig. 10(d) shows the evolution of the normal stress at the (left) trailing corner of the bottom cylinder, where the crack nucleation starts in [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF]. Here, the cases from Figs. 9(b) to 9(d) are considered. For high friction, as in case 9(c) with model EA, the stress is always tensile (positive), with a maximum reached when the two cylinders pass each other (u ≈ R). Assuming a critical stress at which de-bonding between the asperity and its rigid base initiates, our simulations may allow to investigate for which parameters and geometrical properties this de-bonding would occur. The outcome of such an investigation may not be trivial, because both the stress profile along the base (see the inset in Fig. 10(d)) and the stress evolution at the trailing corner can differ significantly for different models.

3D peeling of a tape-like membrane

We now study 3D peeling of a soft and compliant tape (Fig. 11(a) and Tab. 3), which is modeled as a membrane (Sauer et al., 2014). Like other pure membrane models, this formulation only captures in-plane stresses, but does not account for any bending moments.

E • T ν A H r 0 L 0 0.02 N/m 0.3 10 -20 J 0.4 nm 1 nm
Table 3: 3D peeling of a membrane: Model parameters; T denotes the thickness of the tape.

Initially, the entire tape is pre-stretched and attached to a rigid and motionless substrate. The (uniform) pre-stretch of the tape is chosen as λ = 1.001. This is sufficient to avoid any local in-plane compression during peeling (which cannot be resisted by a pure membrane). The frictional sticking to the substrate prevents the tape from directly shrinking back to its un-stretched configuration. Afterwards, we start peeling off one side of the tape by prescribing a displacement along the angle θ = 45 • . Note that this angle corresponds to the direction of motion, and not to the angle between the tape and the substrate at the peeling front.

Since (pure) membranes do not have any bending stiffness, the normal gap from Eq. ( 8) between the tape and the plane always satisfies g n ≥ g eq during peeling (in contrast to the behavior observed in Section 4.1). From Fig. 2 then follows that the normal contact traction is purely adhesive/tensile within the entire contact area. As a consequence, this particular setup cannot be investigated with friction models that yield a tangential resistance only under compression. Such models include both model EA with s cut = 0 (corresponding to the classical Amontons-Coulomb law for friction) and the model of [START_REF] Deng | Adhesion-dependent negative friction coefficient on chemically modified graphite at the nanoscale[END_REF]. Apart from that, even though it would be possible to model tangential sticking by means of a cohesive zone model, the membrane would slide without any frictional resistance after debonding.

The membrane is described by an incompressible Neo-Hookean material model (see Sauer et al. (2014)) with the parameters given in Tab. 3. To exploit the symmetry of the tape shown in Fig. 11(a), we discretize only one half of it by means of 540 quadratic NURBS (N2.1) elements (Sauer et al., 2014), and apply suitable boundary conditions along the center line. Fig. 11(b) shows the peeling process once for a large and once for a small friction parameter, leading either to strong sticking or immediate sliding. As can be seen in Figs. 11(c) and 11(d), in both cases the forces in the vertical direction are very close to each other for models EA and DI. The horizontal forces, however, are similar only for small friction. This is caused by the flexibility of the tape (see also Fig. 11(b) top), for which the surface stretch J ck is not negligible anymore. It also explains the differences between the behavior observed here and the results from Section 4.1.2. Note that the kink in Figs. 11(c) and 11(d) appears at the onset of full sliding, also as observed in Section 4.1.2.

As this example demonstrates, we are able to model adhesive friction even for very soft structures exhibiting a small or negligible bending stiffness. These include adhesive tapes as well as adhesive pads of insects. Even though the bending stiffness of such adhesive pads is sufficiently large at small length scales (in order to prevent self-sticking and resulting entanglement), at larger length scales, the bending stiffness becomes negligible, as it scales with thickness cubed. In summary, the proposed computational model is able to describe sliding friction during peeling also for a negligible bending stiffness, which would not be possible with conventional friction models or cohesive zone models.

Onset of frictional sliding for elastomer-like contact

Let us finally investigate adhesive friction between a smooth elastomer cap and a rigid substrate, as considered experimentally e.g. in [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF] and [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF], theoretically e.g. in Papangelo and Ciavarella (2019), and numerically considering linear elasticity [START_REF] Khajeh Salehani | Modeling adhesive contacts under mixed-mode loading[END_REF]. To this end, we bring a rigid plate into contact with a deformable, cylindrical cap by applying a constant normal force F n (Fig. 12(a)). Due to adhesion between the cap and the plate, this force can be compressive, zero, or even slightly tensile. Keeping F n constant, we then move the rigid plate horizontally, while keeping the lower boundary of the cap fixed. As demonstrated in a short preliminary study [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF], model DI is well-suited to investigate this specific setup. We here consider 2D plane-strain conditions, nonlinear and nearly incompressible Neo-Hookean material behavior [START_REF] Bonet | Nonlinear Continuum Mechanics for Finite Element Analysis[END_REF], and the model parameters from Tab. 4. The cap itself consists of 42,300 linear (Q1) elements with quadratic NURBS enrichment (Q1N2.1 elements) on the surface [START_REF] Corbett | NURBS-enriched contact finite elements[END_REF]. In order to prevent volumetric locking arising from the near-incompressibility, we consider reduced integration for the hydrostatic (also denoted dilatational) part of the stress tensor. See also [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF] for further details. To facilitate a comparison with 3D experimental and theoretical results, in the following we will denote the (unspecified) out-of-plane width of the cap by W .

E

Adhesive friction between cap and substrate

Fig. 12(b) illustrates the contact interface deformation for sliding at F n = 0, including the stress distribution in the bulk of the cap as well as the corresponding contact tractions. The negative peaks in the normal traction indicate adhesive stresses at the two contact edges. Fig. 13(a) shows the tangential force versus the prescribed displacement for different normal forces, while Fig. 13(b) depicts the corresponding contact length L c (see also Fig. 12(a)). The combination of both, i.e., the contact length in dependence of the tangential force, is shown in Fig. 13(c). We find that for all investigated normal loads (compressive, zero, and tensile), the contact length is a decreasing function of the tangential force. By nature of model DI, this decrease stops when all points within the contact have reached the frictional shear strength µ DI T max (see the squares). This explains why all curves in Fig. 13(c) end at the dotted line. For the two shortest initial contact lengths (black solid lines in Figs. 13(a) to 13(c)), Newton's method applied in the (quasi-static) computation stops converging when the contact length decreases down to a certain, common value. A likely explanation is that for these cases, the tensile normal force is sufficiently large so that the cap snaps from the substrate before full tangential sliding. This behavior can also be observed in experiments [START_REF] Waters | Mode-mixity-dependent adhesive contact of a sphere on a plane surface[END_REF][START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF]. To test our hypothesis, we consider the generalized plane-strain JKR model by [START_REF] Chen | Non-slipping JKR model for transversely isotropic materials[END_REF], and compute the critical length L detach for which the cap detaches from the substrate. Inserting the parameters from Tab. 4 results in L detach = 4.42 L 0 , which is very close to the smallest possible lengths observed in our simulations (see the dashed black lines and the black crosses in Figs. 13(b) and 13(c)). This good agreement suggests that the lower limit for the contact length is indeed related to a (physically unstable) contact separation.

Let us now compare our results with experiments on shearing of smooth elastomer spheres on glass plates. As explained in [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF], these comparisons are only qualitative due to length scale differences between computations (µm range, related to the adhesion model by [START_REF] Sauer | Formulation and analysis of a three-dimensional finite element implementation for adhesive contact at the nanoscale[END_REF]), and experiments (mm range). In such experiments, the initially circular contact area shrinks to an ellipse during the onset of sliding [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF][START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF]. [START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] demonstrated experimentally that (i) the width of the contact area perpendicular to the direction of sliding nearly remains constant; and (ii) the length L c parallel to the direction of sliding is well-captured by the quadratic fit

L c (F t ) = L c,0 -ξ F 2 t . (37)
In Eq. ( 37), L c,0 is the initial contact length (at F t = 0), and ξ is a parameter that is independent from F t , but dependent on the applied normal load, or equivalently the initial contact length.

As the dash-dotted lines in Fig. 13(c) indicate, Eq. ( 37) is also well-suited to fit the curves from our numerical results, using the least-squares method. Although the caps in the experiments by [START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] and those in our simulations have different dimensions (spherical in the experiments vs. cylindrical in the simulations), we observe consistent behavior. This motivates us investigating whether the dependence of the parameter ξ on the initial contact length/area is similar, too. [START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] observed for their experiments ξ(A c,0 ) ∼ A -2 c,0 , with A c,0 being their (circular) initial contact area. Regarding the initial contact length L c,0 , this translates with A c,0 ∼ L 2 c,0 into ξ(L c,0 ) ∼ L -4 c,0 . For comparison, Fig. 13(d) shows the parameter ξ fitted for our numerical results in dependence of L c,0 . As can be seen, these data are also linear with logarithmic scales. Our simulations yield ξ(L c,0 ) ∼ L β c,0 with β = -3.71, which is in rather good agreement with the -4 inferred from [START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF].

Note that the experimental results by [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF]) also match well with the fracture-mechanics models by Papangelo and Ciavarella (2019) and Papangelo et al. (2019). Those models describe the evolution of the contact area of adhesive sphere-plane contact under shear, assuming either a circular or an elliptic contact area, respectively. For an initially circular area and for small F t , Papangelo and Ciavarella (2019) predict the relation A c (F t ) = A c,0 -α F 2 t with α(A c,0 ) ∼ A -5/4 c,0 . This equation can be related to the parameter ξ from Eq. (37) as follows: Motivated by their experiments, [START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] assume that the contact area shrinks only along the shear direction and has an elliptic shape afterwards. With these two assumptions, the authors show that the exponents of α and ξ differ by the value -1/2, so that one expects ξ(A c,0 ) ∼ A -7/4 c,0 . This, in turn, translates into ξ(L c,0 ) ∼ L -7/2 c,0 , which gives an exponent close to the -3.71 fitted for our numerical results (Fig. 13(d)). In summary, our exponent β lies in between the two values observed experimentally and determined theoretically.

Motivated by this good agreement, we investigate another finding by Papangelo et al. (2019): The authors conclude that the contact area does not always decrease as F 2 t , but rather as F η t , where η is an exponent depending also on the normal load. η is found to be close to 2 only for sufficiently large normal loads (such as those applied in [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF]), while increasing with smaller loads. To test whether our simulations produce a similar behavior, in Fig. 14(a) we plot the evolution of (1 -L c / L c,0 ) as a function of F t . We here consider the adhesive cases from Fig. 13, except those for which the contact area vanishes before full sliding. For each normal force, the curves are linear in logarithmic scales, indicating that indeed (1 -L c / L c,0 ) ∼ F η t . The least-square fitting of η to the numerical results is shown as dashed lines in Fig. 14(a). Fig. 14(b) shows the evolution of the fitted exponent η as a function of the normal force F n . Although η monotonously decreases with increasing F n , it is close to 2 for all considered normal loads, explaining the proper match of the quadratic fits in Fig. 13(c). Overall, both plots in Fig. 14 are in good qualitative agreement with those of Fig. 10 in Papangelo et al. (2019).

To summarize, even though we have different contact dimensions here, our 2D simulations agree very well with both experiments [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] and theoretical models (Papangelo and Ciavarella, 2019; Papangelo et al., 2019) for adhesive friction of smooth elastomer spheres and glass plates.

Compressible vs. nearly incompressible material

Since most experiments on the shear-induced reduction of the contact area include elastomers [START_REF] Savkoor | The effect of tangential force on the contact of elastic solids in adhesion[END_REF][START_REF] Waters | Mode-mixity-dependent adhesive contact of a sphere on a plane surface[END_REF][START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF], all simulations from Figs. 13 and 14 are based on nearly incompressible material. To test whether our findings also extend to compressible materials, we repeat the simulations discussed in the previous section with ν = 0.4 (instead of 0.49). These new results are shown in the supplementary material as Fig. S1. Below, we just summarize both common points and differences with respect to the nearly incompressible case.

Overall, the qualitative behavior does not change, meaning that the main features of the onset of sliding of Hertzian contact is not specific to incompressible materials. Only small quantitative differences can be observed. In the compressible case: (i) both the initial contact length and the final tangential forces are slightly larger; (ii) the exponent β is -3.54 rather than -3.71; and (iii) exponent η approaches slightly smaller values, so that the quadratic fits of the L c /F t curves are slightly less accurate.

Adhesive vs. non-adhesive friction

In the linear elastic, fracture-based models of Papangelo and Ciavarella (2019) and Papangelo et al. (2019), the presence of adhesion is a necessary condition for a reduction of the contact area. This could suggest that the same was true also for our simulations in Section 5.1, which include finite deformations and adhesion/friction described explicitly rather than being lumped in a phenomenological mode-mixity function. With our numerical model, we are able to test this hypothesis by considering five additional cases without adhesion. To this end, we simply set the normal traction T n in Eq. ( 10) to zero for g n > g eq . The normal forces are chosen such that the initial contact lengths are close to five adhesive cases (while the latter require a considerably smaller normal load).

As the five gray curves in Fig. 13(c) show, the contact length decreases under tangential shear even in the absence of adhesion. This observation alone implies that, in contrast to the prevailing view in the literature, adhesion is not necessary to generate a shear-induced decrease of the contact length. The curves still seem to be well-fitted by the quadratic decay (37), but the decrease in L c is lower than for the adhesive cases. As a consequence, the reduction parameter ξ also becomes smaller for the adhesionless cases (Fig. 13(d)), indicating that adhesion enhances the amount of area reduction.

Assuming that ξ still follows a power law of the form ξ(L c,0 ) ∼ L β c,0 , its exponent is found to be β = -2.42 (instead of -3.71 with adhesion). This difference suggests that adhesion influences the area reduction especially for small initial contact lengths or, equivalently, for small normal loads. The enhancement due to adhesion, expressed by the factor (ξ adh -ξ non-adh )/ξ non-adh , is e.g. +163 % for L c,0 ≈ 8.5 L 0 , but only +24 % for L c,0 ≈ 15 L 0 (Fig. 13(d)). For this specific setup, we thus expect adhesion to play a negligible role for an initial contact wider than approximately 17 L 0 . Finding a general prediction for L c,0 beyond which adhesion is not relevant anymore is an interesting question left for future studies.

If we (in theory) considered contact between an incompressible, linearly elastic material and a rigid halfspace, we would expect the normal and tangential displacements in the bulk to remain uncoupled [START_REF] Johnson | Contact Mechanics[END_REF]. In that case, the tangential shear would not have an influence on the contact length. From this we can conclude that in the adhesionless case, the reduction of the contact length is expected to mainly result from the nonlinear deformation of the cap. This assumption has recently been confirmed in 3D simulations for adhesionless friction between a hyperelastic sphere and a rigid plane under high normal loads [START_REF] Lengiewicz | Finite deformations govern the anisotropic shear-induced area reduction of soft elastic materials[END_REF]. That paper identifies various elementary mechanisms that all contribute to changes in the contact area. Investigating those mechanisms with our model, and comparing their influence with the 3D case, is again a promising topic for future work, but beyond the scope of this paper.

Conclusion

This work provides a computational framework to incorporate theoretical models for local adhesion and friction into 3D, large-deformation simulations. These models include the two theoretical continuum models DI and EA, which were recently proposed by [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF]. Such models are applicable to natural and technical systems, in which friction is either considerably influenced or purely dominated by adhesive effects. Up to a certain distance of the adhering surfaces, both models DI and EA capture friction even for zero or tensile contact tractions. This distinguishes them from existing approaches, and is motivated by soft bio-adhesive pads that are able to generate friction even under tensile normal loads.

While [START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF] contains the motivation of the models DI and EA, e.g. in terms of new experimental results, in the first part of this article we derive the corresponding model equations necessary for a 3D, large-deformation finite element (FE) formulation. We then describe the algorithmic treatment of sticking and sliding friction, using the unbiased friction algorithm by [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF]. Finally, we state the resulting contact FE forces as well as their tangent matrices for two deformable solids in 3D.

In the second part, we investigate various application examples to illustrate the capabilities of our simulation framework, as well as the physical properties of models DI and EA. The considered systems include flexible tapes with a negligible bending stiffness, for which sliding friction cannot be modeled with existing approaches like cohesive zone or other friction models. Overall, although qualitative, our results provide interesting perspectives for future studies, e.g. in the field of wear, or for the design of adhesives. It would also be promising to apply our peeling strip model to horizontal peeling, as reported experimentally for elastomers e.g. in [START_REF] Ponce | Effect of friction on the peeling test at zero-degrees[END_REF]. Besides, our models are suitable to computationally investigate peeling of gecko spatulae [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF]Gouravaraju et al., 2020a,b), or "friction hairs" at the feet of insects [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF]. They can be further combined with a geometrically exact beam formulation [START_REF] Sauer | A geometrically exact finite beam element formulation for thin film adhesion and debonding[END_REF], as done in [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF].

In the third part of this article, we finally study in more detail an example that is inspired by the onset of sliding of Hertzian elastomer-glass contact. To this end, we investigate friction of a soft cap and a rigid plate to compare our results both with experimental findings [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF] and theoretical investigations [START_REF] Chen | Non-slipping JKR model for transversely isotropic materials[END_REF]Papangelo and Ciavarella, 2019;Papangelo et al., 2019). As these results demonstrate, our computational model is able to capture remarkably well the qualitative behavior of smooth elastomer-glass interfaces. Our findings suggest that the shear-induced area reduction of such interfaces is not specific to (nearly) incompressible materials. They also suggest that adhesion is not a necessary ingredient to model this phenomenon, but enhances an effect due to finite, non-linear material deformations. These results contribute to the current debate on the shear-induced area reduction in elastomer contact [START_REF] Sahli | Evolution of real contact area under shear and the value of static friction of soft materials[END_REF][START_REF] Sahli | Shear-induced anisotropy in rough elastomer contact[END_REF]Menga et al., 2018[START_REF] Menga | Corrigendum to "Do uniform tangential interfacial stresses enhance adhesion?[END_REF][START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF]Papangelo and Ciavarella, 2019;Papangelo et al., 2019;[START_REF] Khajeh Salehani | Modeling adhesive contacts under mixed-mode loading[END_REF][START_REF] Scheibert | Onset of sliding of elastomer multicontacts: Failure of a model of independent asperities to match experiments[END_REF][START_REF] Mcmeeking | The interaction of frictional slip and adhesion for a stiff sphere on a compliant substrate[END_REF][START_REF] Wang | Sphere and cylinder contact mechanics during slip[END_REF][START_REF] Lengiewicz | Finite deformations govern the anisotropic shear-induced area reduction of soft elastic materials[END_REF].

Appendix A. Tangent matrices for adhesive and frictional contact

This section summarizes the tangent matrices that are required for the linearization of the contact terms in the governing equations (Section 3.2). For details on the derivation of these matrices we refer to Sauer andDe Lorenzis (2013, 2015) and [START_REF] Mergel | Advanced Computational Models for the Analysis of Adhesive Friction[END_REF].

As becomes apparent in Fig. 1, the normal gap g n , and hence also the elemental contact force f e c,k = f e n,k -f e t,k , depend on the surfaces of both interacting bodies, k, = 1, 2, k = . For linearization we thus need the following two tangent matrices The scalar J e c is the surface stretch at projection point x p : J e c = a p 1 × a p 2 / A p 1 × A p 2 (see also Eq. ( 11)). The derivative of J e c with respect to x k can be computed from In Eq. (A.9), N ,β contains the partial derivatives of the nodal shape functions N with respect to coordinate ξ β .

  (a) Sliding threshold for varying normal gap gn.(b) Friction law in 2D for different but fixed gn.

Figure 3 :

 3 Figure3: Model DI: Constant sliding traction within the cutoff distance gcut for τ DI = 0.5 Tmax, gcut = 1.1 r 0 , and k DI = 80/r 0 ; here, T 0 = A H / 2πr 3 0 ; the regularization of Tn is discussed in Section 2.2; for εt see Section 3.1; the asterisks mark the distances gmax (where |Tn| = Tmax, yellow), gcut (purple), and 1.5 gcut (light blue).

  (a) Sliding threshold for varying gn; scut = 0. (b) 2D friction law for different but fixed gn; scut = 0. (c) Sliding threshold for varying gn; scut = 0.5. (d) 2D friction law for different but fixed gn; scut = 0.5. (e) Sliding threshold for varying gn; scut = 1. (f) 2D friction law for different but fixed gn; scut = 1.

  Figs. 3(a) and 5(a). In analogy, Fig.5(b) indicates sticking and sliding for model EA with g cut = g max . The remaining quantities in Fig.5will be explained in the following.

  Model EA with gcut = gmax (i.e., scut = 1).

Figure 5 :

 5 Figure 5: Illustration of the return map performed in the corrector step for tangential sliding.

Figure 6 :

 6 Figure 6: 2D strip peeling: Setup (left) and three peeling states for model EA with gcut = gmax and µ EA = 0.01 (right); the colors show the first invariant of the Cauchy stress, tr σ, divided by E.

  (a) Normal traction for model EA. (b) Tangential traction for model EA.

Figure 7 :

 7 Figure 7: 2D strip peeling: Traction at the contact area (in dependence of the horizontal position X normalized by L 0 ) for model EA with µ EA = 0.01 and gcut = gmax; the dash-dotted lines in Fig. (b) indicate sliding.

  (a) Horizontal force * vs. displacement. (b) Vertical force * vs. displacement. (c) Peeling and sliding lengths vs. displacement. (d) Fracture energy * vs. peeling length.

Figure 8 :

 8 Figure 8: 2D strip peeling: (a) & (b) Horizontal and vertical forces * ( * : per out-of-plane width); solid colored lines: model EA with gcut = gmax; dashed colored lines: model EA with gcut = (geq + gmax)/2; solid gray lines: model DI with µ DI = µ EA and gcut = gmax; (c) nominal peeling length L peel (solid) and nominal sliding length L slide (dash-dotted); (d) fracture energy * ; the white dots in Figs. (a) -(d) mark the configurations plotted in Fig. 6; the filled squares indicate the onset of full sliding rounded to next integer value of u.

  (a) Horizontal force * for model EA. (b) Horizontal force * for model DI up to µ DI = 0.9. (c) Comparison of models EA and DI. (d) Normal stress at trailing corner.

Figure 10 :

 10 Figure 10: 2D contact of two deformable cylinders: (a) -(c) Horizontal forces * ( * : per out-of-plane width) acting on each body for models EA and DI, both with gcut = gmax; asterisks indicate the horizontal displacement at which snap through would occur in a horizontally load-controlled setup; (d) normal stress σyy / E at the (left) trailing corner of the bottom cylinder for the cases from Fig. 9(b) -9(d); the inset in Fig. (d) shows σyy / E along the base of the lower cylinder for a displacement of u = 0.5 R; the white dots in Figs. (a) -(d) mark the configurations shown in Fig. 9(b) -9(d), respectively.

  (a) Problem setup and boundary conditions.

  Detachment at θ = 45 • dominated by either sticking (top) or sliding (bottom) for model EA and gcut = gmax; the colors show the logarithmic area stretch ln Jc; the figure shows both halves of the tape. (c) Sticking-dominated peeling: µ EA = µ DI = 0.05. (d) Sliding-dominated peeling: µ EA = µ DI = 0.01.

Figure 11 :

 11 Figure 11: 3D peeling of a membrane: (a) Setup; (b) detachment at θ = 45 • for strong and weak friction; (c) & (d) forces in horizontal (dashed lines) and vertical (solid lines) directions for models EA and DI, the dots mark the configurations of Fig. (b).

  (a) Problem setup. (b) Sliding under zero normal load.

Figure 12 :

 12 Figure 12: Adhesive friction of a 2D cap: (a) Setup height of 10 L 0 refers to the undeformed cap before contact); (b) contact tractions (top panel) as well as deformation and stress distribution (bottom panel) at the interface during full sliding for Fn = 0 with model DI, gcut = gmax, and µ DI = 1; the vertical axis of the deformation plot is stretched by a factor of 2; the colors show the first invariant of the Cauchy stress, tr σ / E; the dashed line indicates the contour of the cap before horizontal sliding.

  (a) Tangential force vs. horizontal displacement. (b) Contact length vs. horizontal displacement. (c) Contact length vs. tangential force. (d) Parameter ξ vs. initial contact length.

Figure 13 :

 13 Figure 13: Adhesive friction of a 2D cap: (a) & (b) Tangential force and contact length, obtained for model DI (gcut = gmax, µ DI = 1) and different normal forces; the two solid black lines indicate contact vanishing under sufficient shear; (c) contact length in dependence of the tangential force; the five gray lines correspond to non-adhesive friction (see Section 5.3); (d) parameter ξ from Eq. (37) fitted for all the cases from Fig. (c) with (non-vanishing) adhesive and non-adhesive contact.

  (a) Contact length vs. tangential force (logarithmic). (b) Parameter η from Fig. (a) vs. normal force.

Figure 14 :

 14 Figure 14: Adhesive friction of a 2D cap: (a) Logarithmic representation of the contact length in dependence of the tangential force as considered in Papangelo et al. (2019); the figure shows the adhesive cases from Fig. 13 with stable (non-vanishing) contact; (b) parameter η (fitted from (1 -Lc / L c,0 ) ∼ F η t , and illustrated in Fig. (a) as dotted lines) versus the normal force.

  u e k contains the nodal displacements for element Γ e ck , and u e is the displacement vector of all elements Γ e c that are affected by Γ e ck . In analogy to Eq. (33), the tangent matrices in Eq. (A.1) can be decomposed into the contributions k e c,kk = k e n,kk -k e t,kk and k e c,k = k e n,k -k e t,k .Appendix A.1. Tangent matrices for normal (adhesive and repulsive) contact For normal contact (Section 2.2) with θ k ≡ 1 and J ≈ J c , the first matrix is given by k e n,kk = -

  (A.5), 1 denotes the identity tensor, a p α and n p are given by Eqs. (4) and (6), and c αβ p . The second tangent matrix, containing the derivatives with respect to the neighboring nodes u e ,

  (a) Tangential force vs. horizontal displacement. (b) Contact length vs. horizontal displacement. (c) Contact length vs. tangential force. (d) Parameter ξ vs. initial contact length. (e) Contact length vs. tangential force (logarithmic). (f) Parameter η from Fig. (e) vs. normal force.

Figure S1 :

 S1 FigureS1: Supplementary material -adhesive friction of a 2D cap: Test cases of Sect. 5.1 (panels analogous to Figs.13 and 14), but for compressible material (ν = 0.4 instead of ν = 0.49).

Table 2 :

 2 2D contact of two deformable cylinders: Model parameters.

Table 4 :

 4 Adhesive friction of a 2D cap: Model parameters.

		ν	T max	W adh	L 0
	2 MPa	0.49	0.33 MPa	0.027 J/m 2	1 µm

Preprint submitted to the Journal of the Mechanics and Physics of Solids December 1, 2020

To avoid any confusion with the classical coefficient of friction (defined as the ratio between the frictional resistance and the pressure), for our two models we use the terminology "friction parameter" instead.

This approximation simplifies the linearization of the computational formulation considerably, while introducing a negligible error proportional to 1/εt.

Acknowledgements

This work was initially supported by the German Research Foundation (DFG) under grants SA1822/5-1 and GSC111. JS acknowledges support from the Institut Carnot Ingénierie @ Lyon and from LABEX MANUTECH-SISE (ANR-10-LABX-0075) of the Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). We further thank Dr. Thang X. Duong (RWTH Aachen University) and Antonio Papangelo, Ph.D. (Politecnico di Bari) for helpful comments.

equivalent cylinder-plane contact (as commonly done for two cylinders [START_REF] Johnson | Contact Mechanics[END_REF]), because this would not capture any horizontal forces. Figs. 9(c) and 9(d) show analog snapshots for models EA and DI with a large friction parameter (µ EA = µ DI = 0.9). In contrast to the previous example (where models EA and DI where almost undistinguishable), the deformations of the cylinders considerably differ for both models, see e.g. u = R. We thus investigate also the horizontal forces (Figs. 10(a) and 10(b)), while varying both µ EA and µ DI . For µ EA = µ DI , model EA systematically produces considerably larger forces than model DI. This is explained by the fact that the normal contact tractions are strongly compressible (unlike in the strip example), in particular when the two cylinders pass the apex of each other. Such large normal tractions automatically increase the sliding threshold T slide in model EA, according to Eq. ( 16). In contrast, model DI depends on the normal traction only implicitly, via a changing contact area. Compared to the horizontal force, the vertical force (not shown) changes less with increasing µ EA or µ DI , especially for µ DI ≤ 0.9.

In compression-dominated problems like this setup, the friction parameters µ EA and µ DI cannot be related to each other via Eq. ( 36) anymore. We thus qualitatively compare cases that provide similar maximum horizontal forces instead. One of these cases corresponds e.g. to µ EA = 0.5 and µ DI = 3 (Fig. 10(c)). The Appendix A.2. Tangent matrices for tangential (sticking and sliding) contact

For a friction model stated in the reference configuration (like model EA), the associated tangent matrices have a form that is very similar to Eqs. (A.2) and (A.8):

In contrast, for a model in the current configuration (like model DI), we must additionally linearize the surface stretch appearing in da k = J e ck dA k ; this leads to a second term in k e t,kk ,

∂u e da k , (A.12) see e.g. [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF]. The partial derivatives appearing in Eqs. (A.10) to (A.12) are specified in the following for both sticking and sliding friction. Unless stated otherwise, all quantities are evaluated at the current pseudo-time step, t n+1 ; see Section 3.1. Like in that section, we first discuss the case for which the traction is defined in the current configuration.

Appendix A.2.1. Sticking friction During sticking, the tangential traction corresponds to the trial value t trial ; see Eqs. ( 27) and ( 29). The derivatives of this trial value with respect to both x k and u e are given by

Regarding a model in the reference configuration, one obtains the same two expressions for ∂T trial /∂x k and ∂T trial /∂u e , but with a penalty parameter that refers to the reference area, dA k .

Appendix A.2.2. Sliding friction During sliding, the contact traction satisfies Eq. ( 25). After introducing The surface stretch J e c and its derivatives are given in Eqs. ( 11