Je Souhaite Remercier Tout D'abord Mon Directeur De Thèse

M Zhi-Qiang

Professeur Feng

Encadrant De Thèse

M Yu Cong

M Vladislav

A Yas- Trebov

J ' Ajoute

qui ont toujours su entretenir une

Chapter 1

General introduction 1.1 Background

Dynamic frictional contact involving recoverable interface adhesion represents a frequent phenomenon of contact. Widely observed in nature, with examples on both the macroscopic scale such as biological sticky pads of lizards and insects [4], and the microscopic scale such as cell to cell contact [5,[START_REF] Di Cio | Cell sensing of physical properties at the nanoscale: Mechanisms and control of cell adhesion and phenotype[END_REF], recoverable adhesive contact has attracted significant attention in research. In the adhesive contact, the adhesion effect is caused by intermolecular forces, such as Van Der Waals forces and hydrogen bond, see Figure 1.1. There are numerous examples of intermolecular forces in real life: adhesives contain polyvinyl alcohol and water, which fuse to form hydrogen bonds, resulting in strong adhesion [START_REF] Kinloch | Adhesion and adhesives: science and technology[END_REF]; Gecko's feet has nearly five hundred thousand keratinous hairs or setae, which allow them to crawl on walls due to van der Waals forces [START_REF] Autumn | Adhesive force of a single gecko foot-hair[END_REF], see Figure 1.2.

Then, interface adhesion anisotropy is another topical issue, which have been investigated both experimentally and theoretically [START_REF] Gao | Shape insensitive optimal adhesion of nanoscale fibrillar structures[END_REF][START_REF] Gao | Mechanics of hierarchical adhesion structures of geckos[END_REF][START_REF] Yao | Mechanics of robust and releasable adhesion in biology: Bottom-up designed hierarchical structures of gecko[END_REF] by researchers in the area of biomechanics, leading to numerous applications, such as bio-mimetic adhesive materials [START_REF] Gorb | Biomimetic mushroom-shaped fibrillar adhesive microstructure[END_REF][START_REF] Meng | Tree frog adhesion biomimetics: opportunities for the development of new, smart adhesives that adhere under wet conditions[END_REF], as Figure 1.3. Some of them incorporate anisotropic interface properties of ad- hesion [START_REF] Beisl | Confined evaporation-induced selfassembly of colloidal lignin particles for anisotropic adhesion[END_REF][START_REF] Tardy | Exploiting supramolecular interactions from polymeric colloids for strong anisotropic adhesion between solid surfaces[END_REF][START_REF] Jin | Biomimetic bidirectional switchable adhesive inspired by the gecko[END_REF][START_REF] Afferrante | The ultratough peeling of elastic tapes from viscoelastic substrates[END_REF][START_REF] Carbone | Origin of the superior adhesive performance of mushroom-shaped microstructured surfaces[END_REF][START_REF] Das | Stick-slip friction of gecko-mimetic flaps on smooth and rough surfaces[END_REF][START_REF] Xue | Stick-slip friction of pdms surfaces for bioinspired adhesives[END_REF]. In the area of anisotropic friction modelling, recent contributions have led to numerous orthotropic interface models. We cite in particular the development of orthotropic slip functions [START_REF] Mróz | An anisotropic friction and wear model[END_REF][START_REF] Zmitrowicz | A theoretical model of anisotropic dry friction[END_REF][START_REF] Zmitrowicz | Mathematical descriptions of anisotropic friction[END_REF], orthotropic dry interface model [START_REF] He | Anisotropic dry friction between two orthotropic surfaces undergoing large displacements[END_REF] and elasto-plastic interface model [START_REF] Buczkowski | Elasto-plastic interface model for 3d-frictional orthotropic contact problems[END_REF]. Konyukhov et al. proposed a series of contributions which implement anisotropic interface adhesion based on covariant description of the interface kinematics [START_REF] Konyukhov | Covariant description of contact interfaces considering anisotropy for adhesion and friction: Part 1. formulation and analysis of the computational model[END_REF][START_REF] Konyukhov | Covariant description of contact interfaces considering anisotropy for adhesion and friction.: Part 2. linearization, finite element implementation and numerical analysis of the model[END_REF][START_REF] Michaloudis | An interface finite element based on a frictional contact formulation with an associative plasticity model for the tangential interaction[END_REF]. Moreover, adhesive properties are often observed on the surface of soft materials [START_REF] Dispersyn | Assessment of hyperelastic material models for the application of adhesive point-fixings between glass and metal[END_REF][START_REF] Dias | Development of a constitutive hyperelastic material law for numerical simulations of adhesive steelglass connections using structural silicone[END_REF], therefore many efforts have been made to investigate the link between material properties and adhesion, such as rubber [START_REF] Roberts | Looking at rubber adhesion[END_REF][START_REF] Raphael | Rubber-rubber adhesion with connector molecules[END_REF][START_REF] Tiwari | The effect of surface roughness and viscoelasticity on rubber adhesion[END_REF], biological soft tissues [START_REF] Cozens | Comparative adhesion of chemically and physically crosslinked poly(acrylic acid)-based hydrogels to soft tissues[END_REF][START_REF] Ferraris | Surface structuring by electron beam for improved soft tissues adhesion and reduced bacterial contamination on Ti-grade 2[END_REF][START_REF] Tabatabaei | Osteo-mucosal engineered construct: In situ adhesion of hard-soft tissues[END_REF][START_REF] Taguchi | Bonding of soft tissues using a novel tissue adhesive consisting of a citric acid derivative and collagen[END_REF].

In the area of numerical modelling, despite the efforts exerted over the last decade to develop adhesive contact algorithms [START_REF] Bazrafshan | On the role of adhesion and roughness in stick-slip transition at the contact of two bodies: A numerical study[END_REF][START_REF] Liprandi | A theoretical-numerical model for the peeling of elastic membranes[END_REF][START_REF] Mergel | Continuum contact models for coupled adhesion and friction[END_REF], modelling realistic 3D problems of contact and friction with recoverable adhesion involving both bonding and de-bonding between soft matters, is still a challenging topic [START_REF] Cocou | A dynamic unilateral contact problem with adhesion and friction in viscoelasticity[END_REF][START_REF] Salehani | Modeling adhesive contacts under mixedmode loading[END_REF], as well as, literature on modelling schemes accounting for interface adhesion anisotropy is still in initial state [START_REF] Kato | A model of anisotropic adhesion for dynamic locomotion control[END_REF][START_REF] Kato | Anisotropic adhesion model for translational and rotational motion[END_REF]. Achieving such models requires first, an appropriate description of the contact laws in both normal and tangential directions with a proper account for reversibility of the 3D interface adhesion, and a robust and stable resolution algorithm, that can deal with the computational difficulties inherent to the problem non-linearities. Secondly, adhesion anisotropy should also be considered to meet requirements of real applications. It has been demonstrated that numerous factors during the bonding process may affect the final state of adhesion anisotropy [START_REF] Liu | Controlled adhesion anisotropy between two rectangular grooved surfaces[END_REF]. In this regard, both bonding and de-bonding processes including the interplay between adhesion forces and the state of damage of the adhesive bonds must be taken care of. Therefore, properly modelling the bonding and de-bonding processes becomes one of our main focuses. In the numerical modelling of surface adhesion associated with biological soft tissues, the schemes of soft tissue structure [START_REF] Fougeron | Combining freehand ultrasound-based indentation and inverse finite element modeling for the identification of hyperelastic material properties of thigh soft tissues[END_REF][START_REF] Hou | A new model of soft tissue with constraints for interactive surgical simulation[END_REF][START_REF] Marino | Stress and strain localization in stretched collagenous tissues via a multiscale modelling approach[END_REF] and interfacial adhesion [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF][START_REF] Maugis | Adhesion of spheres: The JKR-DMT transition using a Dugdale model[END_REF][START_REF] Cocou | A dynamic unilateral contact problem with adhesion and friction in viscoelasticity[END_REF] have been respectively significant advances, however, little effort has been exerted to propose a modelling scheme for biological soft tissues with surface adhesion [START_REF] Leng | Modeling of experimental atherosclerotic plaque delamination[END_REF].

Such modelling schemes, in addition to the previously described adhesive contact model, should also involve a precise construction of structural continuum constitutive models of soft tissue that incorporate information about the tissue morphology allowing for the investigation of the interrelation between structure and function in response to mechanical loading [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF].

100 μm 1.2 Literature review on related subjects

In order to construct a numerical model capable of dealing with dynamic three-dimensional contact problems with adhesion and friction, we adopt the bi-potential contact approach [START_REF] De Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF] and the Raous-Cangémi-Cocou (RCC) adhesive model [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF] respectively, which can thus be used to simulate adhesive contact problems with different hyperelastic materials.

Raous-Cangémi-Cocou adhesive model

Concerning the adhesive interface law, a number of models have been developed over the past decades. The most prominent ones include Johnson-Kendall-Roberts (JKR) [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF],

Maugis-Dugdale (MD) [START_REF] Maugis | Adhesion of spheres: The JKR-DMT transition using a Dugdale model[END_REF], and Greenwood and Johnson models [START_REF] Greenwood | The mechanics of adhesion of viscoelastic solids[END_REF]. These models, proposed as early as the 1970s and considered as reference in the area ever since, provided the theoretical basis for the contact and friction modelling of adhesive interfaces. Although these early models are limited to simple, normal-load scenarios, they inspired numerous subsequent researches that offer extension to mixed-load schemes involving normal and tangential loads. We also note a few recent achievements based on finite element continuum contact models that incorporate mixed mode constitutive interface laws [START_REF] Huang | A mechanical model for the adhesive contact with local sliding induced by a tangential force[END_REF][START_REF] Salehani | Modeling adhesive contacts under mixedmode loading[END_REF][START_REF] Mergel | Contact with coupled adhesion and friction: Computational framework, applications, and new insights[END_REF],

all of which can describe the reaction of adhesive interfaces under complex load involving tension and shear. In order to simulate complex interface behaviours with reversible adhesion, we adopt the Raous-Cangémi-Cocou (RCC) model [START_REF] Raous | Quasistatic Signorini problem with Coulomb friction and coupling to adhesion[END_REF][START_REF] Raous | Friction and adhesion[END_REF][START_REF] Raous | Interface models coupling adhesion and friction[END_REF][START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF], which over the years has confirmed its robustness in dealing with adhesive frictional contact. The RCC model incorporates a complete set of interface law involving friction and reversible adhesion. It describes the strength of interface adhesion by prescribing an intensity parameter β [START_REF] Fremond | Contact with adhesion[END_REF]. Varying between 0 and 1, as shown in Figure 1.4, β characterizes the damage level of the interface adhesive bonds (0 refers to the state of complete de-bonding, 1 refers to complete bonding), and subsequently describes the reversible de-bonding and bonding process as function of the geometrical configuration of the contact interface [START_REF] Cocou | A dynamic unilateral contact problem with adhesion and friction in viscoelasticity[END_REF][START_REF] Raous | Recoverable adhesion and friction[END_REF]. From thermodynamic point of view, β derives from a free surface energy and a surface dissipation pseudo-potential.

In this regard, the RCC interface model can be considered as a particular case of the unified adhesion interface model given in [START_REF] Del Piero | A unified model for adhesive interfaces with damage, viscosity, and friction[END_REF], which is similar to the Generalized Standard Material (GSM) [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF] for material modelling.

Bi-potential method

The second aspect that requires attention is the severe non-linearities inherent in contact dynamics. In addition, the non-smooth and multivalued nature of the adhesive interface law gives rise to further computational difficulty. It is therefore necessary to apply robust, and stable algorithms to ensure iteration convergence, solution accuracy with balanced efficiency. A large number of algorithms for the modelling of contact problems by the finite element method have been presented in the literature. See, for example, the monographs by Wriggers [START_REF] Wriggers | Computational Contact Mechanics[END_REF] and Yastrebov [START_REF] Yastrebov | Numerical Methods in Contact Mechanics[END_REF], and the references therein. General computational methods for numerical treatment of contact constraints include penalty method [START_REF] Luenberger | Penalty and barrier methods[END_REF], Lagrangian multiplier method [START_REF] Bertsekas | Constrained Optimization and Lagrange Multiplier Methods[END_REF] and augmented Lagrangian method [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF]. The classical penalty function method is a common algorithm for solving constrained optimization problems.

However, contact boundary conditions and friction laws represent significant numerical difficulty, then it is tricky for the user to choose appropriate penalty factors. The method may become unstable with numerical oscillations when the system approaches the state of contact. In contrary, the Augmented Lagrangian Method is a convenient variant that overcomes the aforementioned disadvantages of the penalty method. The Augmented Lagrangian Method was first introduced to deal with constrained minimization problems.

Since friction problems are not a minimization problem, the method has been extended by Alart and Curnier [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF], Simo and Laursen [START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF] to suit for problems of contact and friction.

We propose to use the bi-potential theory which was also developed based on the augmented Lagrangian method, and in the first place, to solve contact problems in the context of what is called implicit standard materials (ISM) [START_REF] De Saxcé | New inequality and functional for contact with friction: The implicit standard material approach[END_REF][START_REF] De Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF]. Compared to two previous approaches, the bi-potential method couples the two variational inequalities of the with one unique inequality. The frictional contact problem is treated in a reduced system by means of a predictor-corrector solution algorithm in this method, where the corrector can be analytically found with respect to the three possible contact statuses: τ ∈ K µ (sticking), τ ∈ K * µ (no contact) and τ ∈ 3 -(K µ ∪ K * µ ), where K * µ is the polar cone of K µ , see Figure 1.5. Introduced in the 1990s, the approach has been recently extended to problems involving hyperelatic or elastic-to-plastic contact [START_REF] Feng | An elasto-plastic contact model applied to nanoindentation[END_REF][START_REF] Peng | Automatic contact detection between rope fibers[END_REF][START_REF] Zhou | A computational strategy for the modeling of elasto-plastic materials under impact loadings[END_REF] with interface wear [START_REF] Ning | Uzawa algorithm to solve elastic and elastic-plastic fretting wear problems within the bipotential framework[END_REF][START_REF] Ning | A Newton-like algorithm to solve contact and wear problems with pressure-dependent friction coefficients[END_REF]. In the area of adhesive contact modelling, the bi-potential theory has been recently applied to solve 2D interface adhesion between elastic materials [START_REF] Terfaya | A bipotential method coupling contact, friction and adhesion[END_REF].

Hyperelastic materials

Hyperelastic material refers to a constitutive model of an ideal elastic material in which the stress-strain relationship is expressed as a strain energy density function. Linear elastic models do not properly explain the observed material behaviour for all deformable materials, such as rubber, which is non-linearly elastic, isotropic, incompressible. Hyperelasticity provides a means of modelling the stress-strain behaviour of such materials [START_REF] Muhr | Modeling the stress-strain behavior of rubber[END_REF]. The Saint Venant-Kirchhoff model is one of the simplest hyperelastic material model, which is an extension of the geometrically linear elastic material to the hyperelastic regime. In addition, there are numerous hyperelastic material models, such as Neo-Hookean model [START_REF] Ogden | Non-linear elastic deformations[END_REF],

Mooney-Rivlin model [START_REF] Mooney | A theory of large elastic deformation[END_REF], Ogden model [START_REF] Ogden | Large deformation isotropic elasticity-on the correlation of theory and experiment for incompressible rubberlike solids[END_REF], Gent model [START_REF] Gent | A new constitutive relation for rubber[END_REF] etc. As well as, some constitutive models to simulate the anisotropic hyperelasticity of biological soft tissues, such as Holzapfel-Gasser-Ogden (HGO) model [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF][START_REF] Eberlein | An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies[END_REF]. More recently, Cai et al. proposed new polyconvex constitutive models of soft tissues [START_REF] Cai | A simple polyconvex strain energy density with new invariants for modeling four-fiber family biomaterials[END_REF][START_REF] Cai | Integrity basis of polyconvex invariants for modeling hyperelastic orthotropic materials-application to the mechanical response of passive ventricular myocardium[END_REF].

In this thesis, we adopt the Blatz-Ko hyperelastic model, which is used to model compressible foam-type polyurethane rubbers [START_REF] Blatz | Application of finite elastic theory to the deformation of rubbery materials[END_REF], to investigate the adhesive contact problems between soft matters, this material model will always be present in the content of Chapters 2 and 3. In Chapter 4, we select HGO model to simulate the surface adhesion of collagenous biological soft tissues. The collagen fibers behaviour lead to the anisotropy [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF], which is loaded in tension and buckled under compression [START_REF] Holzapfel | Comparison of a multi-layer structural model for arterial walls with a fung-type model, and issues of material stability[END_REF]. Therefore, fibers arrangement has a significant effect on the mechanical behaviour of soft tissues. Meanwhile, the matrix of soft tissues behaves in an isotropic manner, hence the energy densities of collagenous soft tissues contain isotropic and anisotropic parts [START_REF] Weiss | Finite element implementation of 124 incompressible, transversely isotropic hyperelasticity[END_REF][START_REF] Balzani | A polyconvex framework for soft biological tissues. adjustment to experimental data[END_REF], and each collagen fiber family has an independent anisotropic energy density. Following the formulation of HGO model, it has been recently extended to problems involving modelling of atherosclerotic plaque delamination [START_REF] Leng | Modeling of experimental atherosclerotic plaque delamination[END_REF] and brain tissue [START_REF] Wu | Explicit modeling of white matter axonal fiber tracts in a finite element brain model[END_REF], fibers arrangement effect [START_REF] Peyraut | A closed form solution for the uniaxial tension test of biological soft tissues[END_REF], and different hyperelastic material comparison of matrix [START_REF] Chen | Anisotropic hyperelastic behavior of soft biological tissues[END_REF]. In our work, we extend the HGO model combining with Yeoh hyperelastic constitutive law [START_REF] Yeoh | Some Forms of the Strain Energy Function for Rubber[END_REF][START_REF] Renaud | The Yeoh model applied to the modeling of large deformation contact/impact problems[END_REF] to represent the anisotropic hyperelastic behaviour of soft tissues.

Dissertation structure

This thesis is organized in the following manner:

In the current Chapter, we first present the background of the thesis, explaining the current research results and the application directions of adhesive contact in the related field. We then present, in the field of numerical simulation, the current challenges in constructing the adhesive contact model capable of dealing with related quasi-industrial problems. Furthermore, to overcome this challenge, the contact method, the adhesive model and the material model chosen for this thesis are all illustrated.

In Chapter 2, an extended, ready-to-implement 3D model for quasi-industrial problems of contact with friction and recoverable interface adhesion between soft material is formulated using the Raous-Cangémi-Cocou (RCC) interface model and a bi-potential based resolution method. According to the RCC description, the recoverable adhesive interface behaviour derives from a free surface energy and a surface dissipation pseudo-potential.

The obtained interface law describes both the de-bonding process of adhesive links due to tangential and normal interface deformation, and reversely, the bonding process that takes place when two surfaces approach closely enough. We then propose an associated formulation coupling 3D extended interface law and Blatz-Ko hyperelastic material, that enables modelling large deformations of foam type soft matters under conditions of contact and friction with recoverable adhesion. In the end, the subsequent local contact nonlinear equations are solved using a Newton-like algorithm within the bi-potential framework.

Numerical examples are performed to demonstrate the capacity of the proposed approach.

Based on the formulation of Chapter 2, in order to deal with more complex contact interface situations, an orthotropic adhesion model is proposed in Chapter 3 to solve adhesive contact problems with orthotropic interface properties between hyperelastic bodies.

The model proposes a straightforward description of interface adhesion with orthotropic adhesion stiffness, whose components are conveniently expressed according to the local orthogonal Cartesian coordinate system. Based on this description, a set of extended unilateral and tangential contact laws has been formulated. Furthermore, we use an element-wise scalar parameter β to characterize the strength of interface adhesive bonds, and the effects of damage. Therefore, complete cycles of bonding and de-bonding of adhesive links with the account for orthotropic interface effects can be modelled. The proposed model has been tested on cases involving both tangential and unilateral contact kinematics. The test cases allowed emergence of orthotropic interface effects between elastomer bodies involving hyperelasticity. Meanwhile, the model can be implemented with minimum effort, and provides inspiration for the modelling of adhesive interface effects in areas of applications such as biomechanics.

In Chapter 4, a numerical formulation is proposed based on the Holzapfel-Gasser-Ogden In Chapter 6, we conclude the thesis and give perspectives on future research.

Chapter 2

A bi-potential contact formulation for recoverable adhesion between soft bodies based on the RCC interface model

Introduction

In this chapter, a 3D numerical formulation for contact problems with friction and recoverable interface adhesion between soft material is formulated using the Raous-Cangémi-Cocou (RCC) interface model and a bi-potential based resolution method. The recoverable adhesive interface law describes the bonding and debonding processes due to tangential and normal interfacial deformation. We then incorporate a three-dimensional extended interface law and Blatz-Ko hyperelastic materials for modelling large deformations of foam-like soft materials under contact conditions with friction and recoverable adhesion. Finally, the subsequent local contact nonlinear equations are solved using a Newton-like algorithm within the bi-potential framework.

In the following, in Section 2.2, after a brief description of the contact kinematics, the RCC model is described, next, the complete framework of the adhesive contact law is constructed, which includes extended formulations of Signorini contact law and Coulomb friction rules. Then we present its implementation within the bi-potential framework, and provide the formulation of the Blatz-Ko hyperelastic material. In Section 2.3, the complete finite element formulation of the problem, including the resolution algorithm, is provided.

To validate the framework, we present numerical examples in Section 2.4. In the end, a few concluding remarks are drawn in Section 2.5.

Problem setting

Contact kinematics

We describe in this section the geometric definitions and notations related to the contact kinematics. Let's consider two deformable bodies B 1 and B 2 coming into contact with N c contact points. Each body is discretized with finite elements with nodal positions represented by X 1 (for B 1 ) and X 2 (for B 2 ) defined in the global coordinate system.

Contact points belonging to B 1 are denoted by P α 1 (α = 1, 2, ...N c ), and accordingly P α 2 .

Positions of P α 1 and P α 2 can be written using an interpolation matrix B 1 (accordingly B 2 )

as:

X(P α 1 ) = B 1 X 1 , X(P α 2 ) = B 2 X 2 . (2.1)
We consider on each P α 1 a local orthogonal coordinate system, formed by T 1 , T 2 and N, representing respectively the tangential, and normal direction vectors defined with respect to the global coordinates. Therefore, P α 1 can be seen as the projection point of P α 2 on B 1 .

We can build the relative position between P α 1 and P α 2 by

X α = X(P α 2 ) -X(P α 1 ), (2.2)
with X(P α 1 ) and X(P α 2 ) the position vectors of P α 1 and P α 2 in the global coordinates. We can then introduce x α , the local relative position vector of the contact point α, by projecting X α in the system (T 1 , T 2 , N):

x α =            x α t 1 = T T 1 X α x α t 2 = T T 2 X α x α n = N T X α            . (2.3)
We can thus express the local position vector x α as function of the global vector X:

B 1 B 2 P 1 P 2 g N T 1 T 2 Figure 2.1: Contact kinematics x α = H α X α , (2.4)
where H α is the transition matrix obtained by combining Eqs.(2.1,2.2,2.3). Similar relations can be determined with respect to contact forces. The local gap vector between two contact points can be derived from the incremental form of Eq.(2.4):

x α i+1 = H α ∆X α i + g α , (2.5)
with g α = (0, 0, g α ) T , the initial gap vector.

Then, let's denote the local and global contact force vectors with respectively r α and R α . By writing the virtual work

(r α ) T δx α = (R α ) T δX α , (2.6)
we obtain the relation between contact force vectors expressed in local and global coordinate systems:

R α = H T α r α . (2.7)
Here, due to the presence of adhesion on the contact interface, contact reaction r α is composed of the cumulative effects due to both dry contact and the interface adhesion, hence

r α = rα + rα , (2.8)
in which we use r to denote contact reactions associated with unilateral contact and tangential friction, and r contact forces due to interface adhesion. Note that the above relation can be projected to the local coordinate system according to the normal and tangential directions:

     r α n = rα n + rα n r α t = rα t + rα t .
(2.9)

We now assemble all the N c contact points based on Eqs.(2.4,2.5,2.7), we obtain the following geometric and kinematic relations:

     x = H∆X + g R = H T r , (2.10)
with

H =       H 1 . . . H Nc       , x =       x 1 . . . x Nc       , r =       r1 + r1 . . . rNc + rNc       , g =       g 1 . . . g Nc       . (2.11)

RCC model for recoverable adhesion

We use in this work the RCC model to describe the effect of recoverable adhesion between contact surfaces introduced by Raous et al. [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF], this model accounts for unilateral contact, friction and adhesion, based on an energy description of the contact interface, involving a free surface energy Ψ written as:

Ψ (x t , x n , β) = C n 2 x 2 n β 2 + C t 2 x t 2 β 2 -wβ + + (x n ) + Q (β) , (2.12)
and a pseudo-potential of the surface dissipation Φ:

Φ( ẋt , x n , β) = µ|r n -C n x n β 2 | ẋt + b 2 | β| 2 . (2.13)
In these expressions, β is a scalar parameter measuring the intensity of adhesion [START_REF] Fremond | Contact with adhesion[END_REF],

with β ∈ [0, 1]. Specifically, β = 0 represents no adhesion, β = 1 indicates perfect adhesion.

Therefore, any β ∈ (0, 1) refers to partial adhesion between contact surfaces. Other parameters in Eqs.(2.12,2.13) include: C t and C n : parameters characterizing the initial adhesive stiffness when adhesion is complete, w: decohesion energy threshold, : indicator function that assures unilateral contact (x n 0), and meaningful values of the degree of adhesion.

The subscript

Q indicates Q = {η | 0 η 1}, µ: friction coefficient, b: surface viscosity.
Deriving the surface free energy Eq.(2.12), we obtain the expression of the normal force of adhesion:

r ad n = C n x n β 2 , (2.14)
and the tangential force of adhesion:

r ad t = C t x t β 2 .
(2.15)

Both adhesion forces are dependant on the degree of adhesion β. Then deriving energy functions Eq.(2.12) and Eq.(2.13) with respect to β and β yields the incremental expression of β which gives its evolution in time:

           b β ≥ 0 with β = 0 b β = w -(C n x 2 n + C t x t 2 )β with 0 < β < 1 b β ≤ w -(C n x 2 n + C t x t 2 ) with β = 1 .
(2.16)

In Eq.(2.16), we can see that two components may impact the variation of β: the decohesion energy w and the elastic energy of the interface. When interface elastic energy prevails, β becomes negative, leading to decreasing β. Otherwise, β is positive, then β increases. We can view this adhesive model as a special spring system whose elasticity incorporates damage and self-recoverable behaviours. In this regard, the value of β can be seen as the degree of damage of the spring, whose stiffness is adjustable based on β.

Therefore, the decrease of the degree of adhesion β corresponds to the process of spring damage and breaking. Inversely, it can be seen as a recovering process of the spring stiffness.

Adhesive contact law and friction rule Modified Signorini law with adhesion

We recall the unilateral contact law, also called Signorini law, which for classical dry contact is characterized by conditions of non-penetration and non adhesion. By using rα n to denote local normal contact force on the point α due to dry contact, and the contact distance x n , we have

     x α n = ∆x α n + g = 0, rα n 0 x α n = ∆x α n + g > 0, rα n = 0 ⇒ x α n rα n = 0 .
(2.17)

The first relation eliminates geometric penetration between contact surfaces. The second inequality indicates the absence of adhesion forces between dry contact surfaces once they are separated. For adhesive contact, since contact forces result from both the effects of dry contact and adhesion, the classical conditions of unilateral contact should be modified by considering Eq.(2.8), hence

     x α n = 0, r α n -rα n 0 x α n > 0, r α n = rα n ⇒ x α n (r α n -rα n ) = 0 . (2.18)
Here, normal adhesive forces rα n are zeros with surfaces in contact. They will appear when contact surfaces start to separate (the second relation), and rα n will tend to maintain the contact surfaces together. By considering Eq.(2.14), a modified Signorini condition with account for adhesion writes By assuming perfect adhesion (β = 1) at x n = 0, the state of interface adhesion that evolves with x n can be distinguished by three major phases:

     x α n = 0, r α n -C n x α n β 2 0 x α n > 0, r α n = C n x α n β 2 ⇒ x α n (r α n -C n x α n β 2 ) = 0 . ( 2
(i) Fully bonded adhesion: Adhesion bonds remain undamaged (β = 1). In this phase, elastic energy due to x n does not exceed the decohesion threshold w. Hence, linear relationship dominates the adhesion force vs. displacement curve (light green area in Figure 2.2).

(ii) Adhesion with damage: This phase is highlighted by the light cyan area in Figure 2.2. In this phase, β decreases as the decohesion energy w is overpassed. Damage starts to accumulate on adhesion bonds. Adhesion force rα n = C n x α n β 2 continues to increase briefly with x n , before it decreases under the effect of the decreasing quadratic term β 2 , that represents the effect of damage to the interface adhesion.

(iii) Separation: Contact surfaces are separated due to broken adhesion bonds. β sig-nificantly decreases during the process. According to Eq.(2.16), the decreasing β only tends towards zero without exactly reaching zero. It is therefore convenient to consider a limit of β, that we denote by ε, below which the adhesion bonds can be considered as completely broken. In practice, we adopt an empirical ε = 10 -4 which is associated with adhesion forces on the 10 -9 N magnitude according to our test in Section 2.4.

-4 -3.5 -3 -2.5 - With β on the order of 10 -4 , we observe negligible adhesion forces on the 10 -9 N order.

Modified Coulomb friction rule with adhesion

Classically, friction problems are studied using

r α t = -µ(r α n -rα n ) x α t x α t + rα t ∀ x α t = 0 (sliding) , (2.21)
in which rα t , the adhesive tangential force on contact point α can be calculated by considering Eq.(2.15):

rα t = -C t x α t β 2 , (2.22)
and in the normal direction, contact forces are

     r α n -rα n = rα n ∀ x α n = 0 (unseparated) r α n -rα n = 0 ∀ x α n > 0 (separated) .
( By assuming perfect adhesion (β = 1) at x t = 0, the state of interface adhesion that evolves with x t can be distinguished, similar to the normal scenario described in the previous section, by three major phases: (i) fully bonded adhesion, (ii) adhesion with damage, and (iii) separation. Here, since both the slip vector x α t and the tangential adhesion force

x t x lim x ti C t x ti β 2 -C t x ti β 2 r t1
Complete separation Denbonding process

β = 1 1 > β > ε ε > β

Complete contact law with adhesion

By combining the modified Signorini law and Coulomb rule, we obtain the complete contact law with the account for interface adhesion as follows:

Separation : x α n > 0, r α = rα Sticking :

x α n = 0 and x α t = 0, r α = rα Sliding :

x α n = 0 and x α t > 0, r α n = rα n r α t = -µr α n x α t x α t -C t x α t β 2 , (2.24)
in which rα n refers to the normal contact force on point α when surfaces are in contact. In the Sticking situation, since no relative motion occurs, adhesive forces are absent, contact force vector rα lies in the classical Coulomb cone K µ , defined by

K µ = {r α ∈ R 3 | r α n 0, r α t -µr α n 0}. (2.25)
However, with the appearance of relative motion, either following the normal direction (Separation case), or the tangential direction (Sliding case), the contact force vector r α exceeds the boundary of the classical Coulomb cone K µ due to the adhesive forces rα .

Contrary to the classical Coulomb model for dry friction, the resultant contact force r α

will not remain on the boundary of the Coulomb Cone since the relation between r α t and r α n is no longer linear, but subject to variations due to evolving β, x α t and x α n . We cannot conclude an explicit expression relating r to x. In the work of Terfaya et al.

[77], the adhesion is directly incorporated into the bipotential [START_REF] De Saxcé | The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms[END_REF]. We have adopted a different approach where the progression of adhesion is solved at the resolution level by the augmented Lagrangian method, which offers as accurate results [START_REF] Hu | A bi-potential contact formulation for recoverable adhesion between soft bodies based on the RCC interface model[END_REF].

Contact law within the bipotential method

Based on augmented Lagrangian method, the bi-potential method has been developed to deal with contact and friction problems using a reduced system and a predictor-corrector Uzawa algorithm. For unilateral frictional contact, compared to classical methods that requires resolution of two minimum problems or variational inequalities: the first for unilateral contact and the second for friction, the bi-potential resolution unifies unilateral contact and friction, thus requires one single, unique inequality. From the perspective of contact geometry relations, the bi-potential algorithm can be attributed to the category of "node-to-segment" (NTS) contact algorithms. Comparative algorithms include sequential multi-pass NTS approaches, and more recently, the improved virtual-slave-node-to-segment (VTS) approach [START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF], which guarantees accurate assessment of contact interface pressure requiring only a single-pass scheme. Comparison of the presented bi-potential method with other contact algorithms is provided in Appendix A.

The bipotential function and inequality of contact law is as follows:

b c (-x α , r α ) = - (-x α n ) + Ku (r α ) + µr α n || -x α t || (2.26) b c (-x α , r α ) -b c (-x α , r α ) ≥ -x α • (r α -r α ), ∀ r α ∈ K µ , (2.27)
where is indicator function.

-and K µ represent respectively the negative real numbers and Coulomb cone.

The indicator functions become null when the variables -x α and r α comply with the restraining conditions.

We multiply both sides of the inequality (2.27) a parameter ρ, which is used to ensure numerical convergence, and substitude (2.26) into (2.27):

ρµ(r α n -r α n )|| -x α t || + [r α -(r α -ρx α )] • (r α -r α ) ≥ 0 . (2.28)
Taking into account the decomposition x = x t + x n n, the following inequality has to be satisfied:

(r α -r * α ) • (r α -r α ) ≥ 0, ∀ r α ∈ K µ , (2.29)
where the modified augmented contact force r * α is defined by:

r * α = r α -ρ(x + µ|| -x α t ||n) , (2.30)
r α is the projection of r * α onto the closed convex Coulomb cone:

r α = P roj(r * α , K u ) . (2.31)
According to the three different contact states, the projection procedure becomes:

if µ||r * α t || < -r * α n then r α = 0 separating elseif ||r * α t || ≤ µr * α n then r α = r * α sticking else r α = r α * -( ||r α * t ||-µr α * n 1+µ 2
)(

r α * t ||r α * t || + µn) sliding .
(2.32)

Blatz-Ko hyperelastic model for soft materials

Blatz-Ko hyperelastic model [START_REF] Blatz | Application of finite elastic theory to the deformation of rubbery materials[END_REF] is widely used to describe behaviours of compressible foam type soft materials. In practical situations, such materials undergo large deformations. To deal with the geometrical transformation with large deformation, we use the deformation gradient tensor F for the soft bodies in contact:

F = I + ∇u, (2.33)
where I is the unity tensor and u the displacement vector. The right Cauchy-Green deformation tensor C is defined as C = F T F, and the Green-Lagrangian strain tensor

E = 1 2 (C -I).
In the case of hyperelastic law, there exists a strain energy density function W which is a scale function of one of the strain tensors, whose derivative with respect to a strain component determines the corresponding stress component. This can be expressed by

S = 2 ∂W ∂C , (2.34)
where S is the second Piola-Kirchhoff stress tensor. In the particular case of isotropic hyperelasticity [START_REF] Ciarlet | Unilateral problems in nonlinear, three-dimensional elasticity[END_REF], Eq.(2.34) can be written by

S = 2 I 3 ∂W ∂I 3 C -1 + ∂W ∂I 1 + I 1 ∂W ∂I 2 I - ∂W ∂I 2 C , (2.35)
where I i denotes the three invariants of the right Cauchy-Green deformation tensor C:

I 1 = C ii ; I 2 = (I 2 1 -C ij C ij )/2; I 3 = det(C).
(2.36)

The Blatz-Ko strain energy density function is given as follows:

W = G 2
I 2 I 3 + 2 I 3 -5 , (2.37)
where G is the shear modulus. By deriving the energy density (2.37) with respect to the three invariants, we obtain

∂W ∂I 1 = 0; ∂W ∂I 2 = G 2I 3 ; ∂W ∂I 3 = G 2 - I 2 I 2 3 + 1 √ I 3 . (2.38)
Reporting the result in the second Piola-Kirchhoff stress tensor (2.35) gives

S = G JC -1 -C -2 , (2.39)
where J = det(F), the Cauchy stress tensor σ is calculated from the second Piola-Kirchhoff stress tensor as follows:

σ = 1 J FSF T . (2.40)
Eq.2.41 can also be written as:

S(E) = G J(2E + I) -1 -(2E + I) -2 , (2.41)
in order to construct the tangential stiffness matrix K through the finite element nonlinear structural analysis, we need to determine the stress-strain tangent tensor D:

D ijkl = ∂S ∂E =G -2J(2E + I) -1 ik (2E + I) -1 lj + J(2E + I) -1 lk (2E + I) -1 ij + 2 (2E + I) -1 ik (2E + I) -2 lj + (2E + I) -2 ik (2E + I) -1 lj .
(2.42)

Numerical implementation 2.3.1 Finite element formulation of the nonlinear problem

Since contact between soft bodies involves treatment of nonlinear kinematic relations and hyperelastic constitutive models (Section 2.2.5), we formulate the nonlinear finite element problem within the framework of large deformations. In this chapter, we use Green-Lagrangian strain tensor E which comprises both linear and nonlinear terms, as function of nodal displacements u:

E = B L + 1 2 B N L (u) u, (2.43)
where B L is the matrix relating the linear strain term to nodal displacements, and

B N L (u),
relates the nonlinear strain term to nodal displacements. From Eq.(2.43), the incremental form of the strain-displacement relationship can be written as:

δE = B L + B N L (u) δu. (2.44)
Using the principle of virtual displacement, we can write the virtual work δU of the problem as:

δU = δu T Mü + δu T A u + V 0 δE T S dV -δu T F ext -δu T R = 0, ( 2
δU = δu T Mü + δu T A u + δu T V 0 B L + B N L (u) T S dV -δu T F ext -δu T R = 0. (2.47)
We can identify in Eq.(2.47) the vector of internal force:

F int = V 0 B L + B N L (u) T SdV. (2.48)
Since δu is arbitrary, a set of nonlinear equations can be obtained as

Mü + A u + F int -F ext -R = 0. (2.49)
It is noted that the stiffness effect is taken into account by the internal force vector F int .

Eq.(2.49) can be transformed into

M ü = F + R, where F = F ext -F int -A u, (2.50)
with the initial conditions at t = 0 u = u0 and u = u 0 .

(2.51)

Taking the derivative of F int with respect to the nodal displacements u gives the tangent stiffness matrix as

K = ∂F int ∂u = V 0 B L + B N L (u) T ∂S ∂u + ∂B T N L (u) ∂u S dV. (2.52)
In addition, by considering Eqs.(2.44, 2.41), the tangent stiffness matrix can be written as the sum of the elastic stiffness matrix K e , the geometric stiffness (or initial stress stiffness) matrix K σ and the initial displacement stiffness matrix K u :

K = K e + K σ + K u , (2.53)
with

K e = V 0 B T L DB L dV K σ = V 0 ∂B T N L ∂u S dV K u = V 0 B T L DB N L + B T N L DB L + B T N L DB N L dV.
(2.54)

Numerical integration algorithm

Now we need to integrate Eq.(2.50) between consecutive time configuration t and t + ∆t. The Newmark method is the most common method which is based on a second order algorithm. However, higher order approximation does not necessarily mean better accuracy and may even be redundant in impact problems. In cases presenting sudden change of contact conditions (impact, release of contact), we observe discontinuous velocity and acceleration, which lead to excessive regularity constraints that may cause serious errors. For this reason, we use the method of Non-Smooth Contact Dynamics (NSCD) [START_REF] Jean | The non-smooth contact dynamics method[END_REF] involving a first order time stepping algorithm. Implementation of this algorithm for adhesion problems has been investigated in [START_REF] Raous | Unilateral contact, friction and adhesion: 3D cracks in composite materials[END_REF]. Based on NSCD, Eq.(2.50) can be transformed into:

M d u = F dt + R dt . (2.55)
This algorithm is based on the following approximations:

t+∆t t M d u = M ut+∆t -ut (2.56) t+∆t t F dt = ∆t (1 -ξ) F t + ξ F t+∆t (2.57) t+∆t t R dt = ∆t R t+∆t (2.58) u t+∆t -u t = ∆t (1 -θ) ut + θ ut+∆t , (2.59)
where 0 ≤ ξ ≤ 1; 0 ≤ θ ≤ 1. In the iterative solution procedure, all the values at time t + ∆t are replaced by the values of the current iteration i + 1; for example,

F t+∆t = F i+1 .
A standard approximation of F i+1 gives

F i+1 = F i int + ∂F ∂u (u i+1 -u i ) + ∂F ∂ u ( ui+1 -ui ) = F i int -K i ∆u -A i ∆ u . (2.60)
Finally, we obtain the recursive form of (2.55) in terms of displacements:

Ki ∆u = Fi + Fi acc + R i+1 u i+1 = u i + ∆u , (2.61)
where the so-called effective terms are given by

Ki = ξ K i + ξ θ ∆t A i + 1 θ ∆t 2 M i (2.62) Fi acc = - 1 θ∆t 2 M i (u i -u t -∆t ut ) (2.63) Fi = (1 -ξ) F t int + F t ext + ξ F i int + F t+∆t ext . (2.64)
At the end of each time step, the velocity is updated by

ut+∆t = 1 - 1 θ ut + 1 θ ∆t (u t+∆t -u t ) . (2.65)
By setting θ = 1 2 , this scheme is then called the implicit trapezoidal rule and it is equivalent to the Tamma -Namburu method in which the acceleration need not be computed [START_REF] Tamma | A robust self-starting explicit computational methodology for structural dynamic applications: Architecture and representations[END_REF].

It is noted that Eq.(2.61) is strongly non-linear, because of large rotations and large displacements of solid, for instance in multibody contact/impact problems. Besides, as mentioned above, the constitutive law of contact with friction is usually represented by inequalities and the contact potential is even non differentiable. Instead of solving this equation in consideration of all nonlinearities at the same time, Feng [START_REF] Feng | 2D or 3D frictional contact algorithms and applications in a large deformation context[END_REF] has proposed a solution strategy which consists in separating the nonlinearities in order to overcome the complexity of calculation and to improve the numerical stability. As ∆u and R are both unknown, Eq.(2.61) cannot be directly solved. First, the vector R is determined by the bi-potential method and the adhesive model in a reduced system, which only concerns contact nodes. Then, the vector ∆u can be computed in the whole structure, using adhesive contact reactions as external loading.

The iterative solution procedure involving contact modeling is written as

Numerical results

The algorithm presented above has been implemented within the in-house finite element code FER/Contact. In this section, four numerical examples based on contact simulations are presented to show normal and tangential behaviours of the adhesive contact interface under unidirectional and mixed loading conditions.

Indentation on adhesive hyperelastic material

The adhesion effect is usually most significant in the normal direction. In order to clearly show the evolution of β during the complete process of bonding and de-bonding, the first example simulates the normal adhesive contact between an elastic semi-sphere and a hyperelastic block, shown in Figure 2.6(a). The density of two bodies is: ρ = 2500 kg/m 3 (sphere indentor); ρ = 1000 kg/m 3 (block). Blatz-Ko hyperelastic material model is used and the shear modulus G for the indentor and the hyperelastic block are respectively 2.6 × 10 10 Pa and 16 × 10 6 Pa. In this case, the sphere indentor behaves as a rigid body compared to the block.

A time dependent displacement is prescribed on the upper surface of the semi-sphere, so that a complete cycle of indentation is performed in 4 seconds. The sequence of β evolution is thus distinctive on each point.

Rolling adhesion of a hyperelastic wheel

This example investigates the rolling adhesion of a hyperelastic wheel confined between 2 rigid plates. As shown in Figure 2.8(a), the upper and lower plates exert compression on the hyperelastic wheel, and slide simultaneously in opposite directions, driving the wheel in rotation under both the effects of interface adhesion and friction. As a result of the interface adhesion, the rotating wheel presents inclined, asymmetrical geometry during rotation. Since the rotating motion is cyclic, any point belonging to the wheel's adhesive surface will cyclically go through bonding and de-bonding process. The recoverability of interface adhesion is thus accounted for. The present case follows the next loading sequence: the upper plate first descends vertically for 5 × 10 -4 m at the velocity of 0.1 m/s, exerting slight compression on the wheel. Then, still on the upper plate, we prescribe a sliding motion at the velocity of 1 m/s so as to drive the compressed wheel in rotation.

We investigate the effect of material properties on the adhesion by testing shows the state of deformation and Von Mises stress distribution of the wheel during its rotation. Due to the combined effects of the interface friction, which exerts pure tangential force on the wheel, and the interface adhesion, which results in both normal and tangential forces on the wheel surface, the rotating wheel deforms into inclined, asymmetrical geometry. This morphology is the result of the adhesion force (attraction) that appears at the separation (de-bonding) between the plate and the wheel. In case adhesion is absent and under the exclusive effect of friction, the wheel will not present inclined shape during rotation, but remain in the configuration of

Adhesive friction between a hyperelastic plate and a deformable semi-cylinder

In this example, we investigate the adhesive friction of a hyperelastic plate that slides on top of a deformable semi-cylinder, as shown in Figure 2.10(a). Both the plate and the semi-cylinder are modelled by Blatz-Ko hyperelastic material, based on the same material property with shear modulus G = 10 MPa. Radius of the cylinder is 5 mm, and the plate thickness H = 2 mm. The plate is sufficiently long so as to ensure contact between the plate and the cylinder during the simulation. While the plate is allowed to slide horizontally, the bottom surface of the semi-cylinder is fixed. The simulated scenario involves 2 stages.

On the first stage, the upper plate descends for 1 mm to exert a slight compression on the cylinder. Then on the second stage, a lateral displacement is prescribed on the plate at a constant velocity. As a result of the combined effect of friction and adhesion, the cylinder is dragged to deform, and we investigate the interfacial behaviour during the process. In particular, by varying the descent velocity of the first stage, we modulate the total time of compression before sliding, during which bonding process takes place. This will have impact on the final adhesion level (characterized by β) before de-bonding starts at the onset of the sliding stage. In order to explore the influence of the adhesion level β on the subsequent adhesive friction behaviour, we set up 5 groups of cases with for each group a different descent velocity (summarized by whereas the most rapid descent (blue curve) did not allow enough time for the formation of complete bonding. In this case, de-bonding was already initiated after β reached 0.2.

We then investigate the influence of friction coefficient µ on the de-bonding behaviour, by prescribing varying friction coefficients µ while considering the same plate descent velocity.

We report in

3D frictional adhesive twisting

We investigate the evolution of interface behaviours of a 3D twist tribosystem (Figure 2.12) under the combined effect of adhesion and friction. The system is composed of an elastomer block that slides on a rigid surface under twisting load. The elastomer block is as can be seen in the frame t = 0.009 s, both cases present consistent configuration, with synchronized twists and similar distribution of normal forces.

We then explore the effect of friction coefficients on the combined adhesive-frictional interface behaviour, which includes interface forces and adhesion strength. Using 3 sets of friction coefficients µ = 0.2, 0.4 and 0.6, we carry out simulations based on the previous adhesive tribosystem, and investigate, for each friction coefficient, the evolution of normal and tangential reactions, and the intensity of adhesion (represented by β). are observed on cases with higher friction coefficients. The effect of higher friction coefficients on the de-bonding process can be further confirmed by investigating the evolution of β, as shown in Figure 2.15, in which light yellow colour indicates the state of complete bonding of interface adhesives. We observe directly that higher friction coefficients significantly delay the onset of be-bonding process.

We also investigate the evolution of tangential forces on the 3 sets of simulations. In 0.4 / t = 0.007 s, and µ = 0.6 / t = 0.008 s, then, on frames µ = 0.2 / t = 0.007 s, µ = 0.4 / t = 0.008 s, and µ = 0.6 / t = 0.009 s, and so on. This observation also results from the effect of higher friction coefficients on delaying the onset of be-bonding process, and subsequently the appearance of every distribution pattern of tangential forces.

Chronologically, at the beginning of loading, tangential forces are most significant on the outskirts of the contact area since linear velocity is higher. However, for the same reason, this is also where the onset of de-bonding initiates and propagates towards the centre area.

Consequently, the peak of tangential forces is observed as an evolving circular band, whose

Concluding remarks

In this chapter, an extended 3D formulation for quasi-industrial problems of adhesive contact with recoverable interface between soft materials under large deformation is im- Future extensions of this chapter include, for example, the account for anisotropic interface behaviours and the effect of interface fatigue which is a common phenomenon in adhesive applications involving cyclic loads.

plemented
Chapter 3

A bi-potential contact formulation of orthotropic adhesion between soft bodies

Introduction

In this chapter, we propose an orthotropic adhesion modelling between soft materials based on the adhesive contact model of Chapter 2. We extend, in the RCC model, derivatives of the free surface energy that yields a straightforward description of the interface adhesion orthotropy, whose stiffness components are conveniently expressed according to the local orthogonal Cartesian coordinate system. A set of extended unilateral and tangential contact rules incorporating the interface adhesion orthotropy is then formulated. We still combine the orthotropic interface law with Blatz-Ko hyperelastic materials for modelling large deformation contact problems with friction and orthotropic adhesion [START_REF] Hu | A bi-potential contact formulation of orthotropic adhesion between soft bodies[END_REF].

The remainder of the Chapter is organized in the following manner: in Section 3.2, we firstly redescribe the contact kinematics with a covariant description, and present the orthotropic extension of RCC adhesive model of contact. Then, the complete contact law is characterised, which includes an extended formulation of the unilateral and tangential rules of contact involving interface adhesion orthotropy. To validate the framework, we present numerical examples in Section 3.3. In the end, a few concluding remarks are drawn in Section 3.4.

Problem setting 3.2.1 Contact kinematics

We describe in this section the geometric definitions and notations related to the contact kinematics. Let's consider two deformable bodies B α , α = 1, 2 coming into contact. Deformation of the two bodies is represented by ϕ α , as shown in Figure 3.1, which maps the initial configuration to positions of the current configuration. We assume that contact occurs at the boundaries ϕ(Γ α c ) in the current configuration where Γ α c ⊂ ∂B α are possible contact surfaces of bodies B α . By assuming that the contact boundary describes, at least locally, a convex region, we can relate to P 2 and P 1 via the minimum distance problem [START_REF] Wriggers | Contact kinematics[END_REF]:

B 2 B 1 φ 2 φ 1 X 2 X 1 φ(B 2 ) φ(B 1 ) φ(X 2 ) = φ(X 1 ) Γ c 2 Γ c 1 φ( ) Γ c 1 φ( ) Γ c 2
d(ξ 1 , ξ 2 ) = x 2 -x 1 (ξ) , (3.1)
where x 2 and x 1 are the position vectors of two points in the global Cartesian coordinate system xyz, ξ = (ξ 1 , ξ 2 ) denotes the parametrization of the boundary ϕ(Γ 1 c ) via convective coordinates [START_REF] Wriggers | Contact constraints within coupled thermomechanical analysis-a finite element model[END_REF][START_REF] Laursen | A continuum-based finite element formulation for the implicit solution of multibody, large deformation-frictional contact problems[END_REF][START_REF] Schweizerhof | Covariant description for frictional contact problems[END_REF][START_REF] Konyukhov | On the solvability of closest point projection procedures in contact analysis: Analysis and solution strategy for surfaces of arbitrary geometry[END_REF]. d(ξ 1 , ξ 2 ) can be used to define the gap between two bodies.

In order to make Eq.(3.1) valid, x 1 needs to satisfy the following condition:

∂ ∂ξ α d(ξ 1 , ξ 2 ) = x 2 -x 1 (ξ 1 , ξ 2 ) x 2 -x 1 (ξ 1 , ξ 2 ) • x 1 ,α (ξ 1 , ξ 2 ) = 0, with α = 1, 2, (3.2)
where x 1 ,α (ξ 1 , ξ 2 ) is the tangent vector a α . With the cross product of tangent vectors, the normal vector n can be defined:

n = a 1 × a 2 ||a 1 × a 2 || , (3.3)
once the first term of Eq.(3.2) has the same direction as n, which proves that the current position of P 1 is the projection of P 2 on B 1 . Therefore, the normal relative displacement or gap g n is:

g n = (x 2 -x 1 ) • n. (3.4)
In the tangential sliding, the path of P 2 on the contact surface of B 1 is unknown, we only know the relative velocity vector. Therefore, the path of P 2 needs to be obtained by integrating over its relative velocity. The increment of tangential relative displacement, as shown in Figure 3.3, is:

dg t = a α dξ α , (3.5)
with dξ α = ξα , the tangential relative displacement can be computed as:

g t = t t 0 ξα dt a α , (3.6)
where t 0 is the initial time and t represent the current time. From Eq.(3.6), in order to obtain g t , we need to first calculate ξα by the following relation:

∂ ∂t [x 2 -x 1 (ξ 1 , ξ 2 )]•a α = [v 2 -v 1 -a β ξβ ]•a α +[x 2 -x 1 (ξ 1 , ξ 2 )]• ȧα = 0, with α, β = 1, 2, (3.7)
where v α = ẋα . We have ȧα = v α ,α + x α ,αβ ξβ , Eq.(3.7) can be developed as an expression containing ξβ :

(a αβ -g n b αβ ) ξβ = [v 2 -v 1 ] • a α + g n n • v α ,α , (3.8) n ξ 1 ξ 2 a 1 a 2 dg T φ( ) Γ c 1 Figure 3.3: Increment of tangential path. with            g n n = x 2 -x 1 a αβ = a α • a β b αβ = x α ,αβ • n , (3.9)
where a αβ and b αβ represent respectively the metric tensor and curvature tensor. Substituting ξβ from Eq.(3.8) into Eq.(3.6), we can solve the tangential slip g t .

Then the contact force vector r is defined as a covariant vector, which is expressed via the contravariant basis surface vectors a α and n:

r = r t + r n = r α t a α + r n n, α = 1, 2, (3.10)
where r t and r n are respectively tangential and normal component of contact force vector.

Let's denote the local and global contact force vectors with respectively r and R. The relation between contact force vectors expressed in local and global coordinate systems writes:

R = H T r, (3.11)
where H is the transition matrix. Here, due to the presence of adhesion on the contact interface, contact reaction r is composed of the cumulative effects due to both dry contact and the interface adhesion, hence r = r + r, (3.12) in which we use r to denote dry contact reactions, and r contact forces due to interface adhesion. Note that the above relation can be projected to the local coordinate system according to the normal and tangential directions:

     r n = rn + rn r t = rt + rt . (3.13)
Generalization of the RCC model is described in Cartesian coordinates.

RCC contact model with adhesion orthotropy

We develop an improved RCC contact model to describe the effect of orthotropic adhesion between contact surfaces. Introduced by Raous et al. [START_REF] Raous | A consistent model coupling adhesion, friction, and unilateral contact[END_REF], the original RCC model accounts for unilateral contact, friction and adhesion, based on an energy description of the contact interface, involving a free surface energy Ψ and a pseudo-potential of the surface dissipation Φ. Here, energy expressions Ψ and Φ are formulated based on displacements that we project to the local system (a 1 , a 2 , n), leading to tangential and normal components g t1 , g t2 and g n :

Ψ (g t1 , g t2 , g n , β) = C n 2 g 2 n β 2 + C t1 2 g 2 t1 β 2 + C t2 2 g 2 t2 β 2 -wβ + + (g n ) + Q (β) , (3.14)
Φ( ġt , g n , β) = µ|r n -C n g n β 2 ||| ġt || + b 2 | β| 2 . (3.15)
In the above expressions, β is a scalar parameter that measures the intensity of adhesion [START_REF] Fremond | Contact with adhesion[END_REF], with β ∈ [0, 1]. Specifically, β = 0 indicates the absence of adhesion, β = 1 refers to perfect adhesion. Hence, any β ∈ (0, 1) implies partial adhesion between contact surfaces.

Other parameters in Eqs. Deriving the surface free energy Eq.(3.14), we obtain the expression of the normal force of adhesion:

r ad n = C n g n β 2 , (3.16)
and the tangential forces of adhesion:

     r ad t1 = C t1 g t1 β 2
r ad t2 = C t2 g t2 β 2 .

(3.17)

Both adhesion forces are dependent on the degree of adhesion β. Then deriving energy functions Eq.(3.14) and Eq.(3.15) with respect to β and β yields the incremental expression of β which gives its evolution in time:

           b β ≥ 0 with β = 0 b β = w -(C n g 2 n + C t1 g 2 t1 + C t2 g 2 t2 )β with 0 < β < 1 b β ≤ w -(C n g 2 n + C t1 g 2 t1 + C t2 g 2 t2 ) with β = 1 . (3.18)
In Eq.(3.18), we can see that two components may influence β: the decohesion energy w, and the elastic energy of the interface i=n,t1,t2 C i g 2 i . When interface elastic energy prevails, β becomes negative, which leads to decreasing β. Otherwise, β is positive, then β increases. We can view Eq.(3.14) as a modified penalty function method where both are a spring model with zero rest length, except that the two springs are stretched in opposite directions.

Modified Coulomb slip rule with orthotropic adhesion

Classically, tangential problems are studied using the Coulomb model which is characterized by a set of rate-independent slip rules. The original Coulomb model describes tangential force that evolves with normal forces:

     rα t µr α n ∀ g α t = 0 (sticking) rα t = -µr α n g α t g α t ∀ g α t = 0 (sliding) . (3.19)
Here, with the consideration of adhesion, both tangential and normal forces are supplemented by contributions due to interface adhesion as shown in Eq.(3.12), the above rules become

     r α t µr α n ∀ g α t = 0 (sticking) r α t = -µ(r α n -rα n ) g α t g α t + rα t ∀ g α t = 0 (sliding) , (3.20)
in which rα t , the adhesive tangential force on the contact point α can be calculated by considering Eq.(3.17), and the orthotropic adhesive stiffness parameters C t1 and C t2 defined in Eq.(3.14)

rα t = -C t1 g α t1 β 2 -C t2 g α t2 β 2 =      -C t1 g α t1 β 2 -C t2 g α t2 β 2      . ( 3

Complete separation

Denbonding process

β = 1 1 > β > ε ε > β Fully bonded adhesion r t2 r t1 C t1 g ti -C t1 g ti -C t2 g ti C t2 g ti r t1 r t2 r t g t g lim g ti Figure 3
.4: Modified Coulomb rule with adhesion: evolution of tangential adhesive forces and the level of damage vs. slip state of interface adhesion that evolves with g t can be distinguished, similar to the normal scenario described in the previous section, by three major phases: (i) fully bonded adhesion, (ii) adhesion with damage, and (iii) separation. Here, since both the slip vector g α t and the tangential adhesion force vector rα t lie in the local plane (a 1 , a 2 ), their projection in the local system gives rise to expressions of tangential displacement and forces according to axis a 1 and a 2 . Furthermore, in orthotropic adhesion, distinct adhesion stiffness parameters C t1 and C t2 can be defined in the two principal axes. Hence, the critical tangential forces are:

rcrit t1 = -C t1 g t1 β 2 and rcrit t2 = -C t2 g t2 β 2 . (3.22)
The two critical forces are at the extreme points of the tangential forces ellipse, given by the equation:

(r α t1 ) 2 (C t1 β 2 g t ) 2 + (r α t2 ) 2 (C t2 β 2 g t ) 2 = 1 . (3.23)
The ellipse intersects the x-axis at C t1 β 2 g t and -C t1 β 2 g t . It intersects the y-axis at C t2 β 2 g t and -C t2 β 2 g t . To represent the adhesion orthotropy, any vector of adhesion force can be indicated on the ellipse, pointing from its centre to one point on the periphery. Then for 1 0 m m 6 mm

x y z Vertical compression On this stage, tangential effects have been sufficiently accumulated, leading to initiation of damages to the adhesive bonds. As a result, β significantly decreases, especially on remote nodes with respect to the centre, on which β falls back to 0, indicating rupture of the adhesive bonds. We also find contours of β plotted on the contact surface in Figure 3.6, where the effect of adhesion orthotropy can be distinguished. Since the adhesion stiffness in x axis C tx is 10 times stronger than in y axis, significant resistance to interface sliding can be expected in x axis. Therefore, rupture of the adhesive bonds first appears on the upper and lower peripheries of the contact interface, and gradually propagates towards the centre area. Meanwhile, peripheral areas near the left and right edges remain adhered due to stronger adhesion stiffness C tx in x axis.

Similar effects of adhesion orthotropy can be observed in Figure 3.7, which shows the distribution of the Euclidean norm of tangential adhesive forces on the contact surface and its evolution with time. We note that within areas where de-bonding is initiated, particularly near the upper and lower edges, the adhesion forces decrease quickly to zero.

On the contrary, we observe important adhesion forces in areas near the left and right edges since the adhesion orthotropy results in stronger resistance to sliding motions along the x axis. In conformity with the contours of β given in Figure 3.6, distribution of the adhesion forces in Figure 3.7 reflects identical effects of adhesion orthotropy, demonstrating better resistance to tangential interface effects in x axis compared to y axis.

Orthotropic adhesion in shear sliding

We investigate in this example behaviours of orthotropic adhesion in a test scenario involving shear sliding along varying orientations. Similar experimental setup which demonstrates microstructure based orthotropic adhesion has been explored in [START_REF] Khaled | Anisotropic dry adhesive via cap defects[END_REF]. Here, we model the interface adhesion orthotropy by considering distinctive tangential adhesive stiffnesses C tx and C ty , in x and y axis. The tested system is composed of an elastomer cylinder that slides on a rigid surface under tangential load, which is oriented along varying orientations on each test. As shown in Figure 3 complete load cycle involving bonding and de-bonding. In Figure 3.9(c), we note indistinguishable time history of β during the stage of adhesion bonding. However, initiation of de-bonding does not take place simultaneously for all the cases. It arises first in the case of sliding along x axis, in which direction the adhesion stiffness is the lowest. For the same reason, this scenario also exhibits the lowest adhesion force at the onset of de-bonding process (blue curve in Figure 3.9(d)). Comparatively, with the sliding direction approaching y axis, stronger adhesion stiffness is involved. We observe accordingly retarded initiation of de-bonding, accompanied by increased adhesion forces (red, yellow and purple curves in

Orthotropic adhesive twisting

In this example, we investigate the evolution of interface behaviours of a 3D twist tribosystem (Figure 3.10) by considering both isotropic and orthotropic adhesions. The system is composed of an elastomer block that slides on a rigid surface under twisting load. The elastomer block is 3 mm high, and has a 10 × 10 mm square section. For the isotropic case, the tangential adhesive stiffness C t = 1 × 10 10 N.m -3 , and for the orthotropic case We begin by investigating the effect of interface adhesion by comparing cases with and without the interface adhesion orthotropy. Figures 3.11-3.12 compare respectively the evolution of adhesion damage parameter β, and the tangential adhesion forces Rt , between the isotropic and orthotropic cases during the twisting process. For each group of comparison, 5 frames of result are extracted in chronological order to represent the evolving twist process. This allows us to highlight for each time instant, differences between the isotropic and orthotropic cases in terms of β and Rt distributions. In Figure 3.11, we use dark red colour to indicate complete bonding of the interface adhesives. As we apply twist kinematics to the elastomer body, tangential interface effects start to appear on the contact interface. They become first noticeable on the outskirts of the contact area where interface sliding is most significant. Damage to the adhesive bonds is thus initiated with decreasing β emerging at the corners of the contact interface, where also the first de-bonded area is observed. Then with the increasing load, de-bonding propagates from the outskirt area towards the centre, whereas the bonded region gradually shrinks until complete disappearance. During the process, the bonded region appears within a round area in the isotropic case. However, when adhesion orthotropy is involved, since stronger resistance to de-bonding is encountered in the x axis where tangential adhesive stiffness is more significant, delayed de-bonding is observed following the x axis, leading to an elliptical bonded region. We also investigate the evolution of tangential forces on the same setup. In Figure 3.12, Euclidean norms of tangential forces are depicted, allowing us to observe the evolving intensity of tangential forces on the contact interface. Chronologically, at the beginning of load, tangential forces are most significant on the outskirts of the contact area since linear velocity is higher. This is also where de-bonding is initiated and propagates towards the centre. Consequently, the peak of tangential forces appears in the form of an evolving circular band, whose radius decreases with the twist load, until gradually disappears in the centre of rotation, leading to complete de-bonding of the interface adhesives. In the case of orthotropic adhesion, the circular band appears in the form of an ellipse since stronger tangential adhesive stiffness is involved in x axis, following which de-bonding requires more efforts. This observation is in accordance with the evolution of β during the simulation.

||R t || (kN)

Figure 3.12: Comparison between isotropic and orthotropic adhesive twisting: Evolution of tangential adhesion forces Rt in two cases during the debonding process. In each square area, the colour progresses from dark red to blue, which represents the variation of the Euclidean norm of tangential adhesion from maximum to zero.

Concluding remarks

In this chapter, we proposed an orthotropic adhesion model to deal with problems of adhesive contact with orthotropic interface properties between hyperelastic bodies. This model has been implemented within the bi-potential method, based on a set of extended unilateral and tangential contact laws. The behaviour of orthotropic adhesion is described by adhesion stiffness, whose components can be expressed according to the local coordinate system. In this model, the strength of interface adhesive bonds and the effect of interfa-cial damage are characterized by a scalar parameter β, therefore an entire bonding and debonding process of the adhesive links with the account for orthotropic interface effects can be modelled. The proposed approach has been tested on cases involving both tangential and unilateral contact kinematics, which allowed emergence of orthotropic interface effects between soft bodies. Owing to the straightforward description of the contact rules, the presented approach can be easily implemented. Therefore, immediate implementation of this orthotropic adhesion model within a third-party software can be suggested for direct application on real problems.

HGO hyperelastic model for biological soft tissues

Holzapfel-Gasser-Ogden(HGO) model is widely used to describe anisotropic hyperelastic behaviours of biological soft tissues [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF]. It is usually assumed that anisotropy is commonly attributed to the behaviour of collagen fibers [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF], while the ground substance is considered to be isotropic, energy densities modeling transversely isotropic and orthotropic soft tissues are split in isotropic and anisotropic components [START_REF] Weiss | Finite element implementation of 124 incompressible, transversely isotropic hyperelasticity[END_REF][START_REF] Balzani | A polyconvex framework for soft biological tissues. adjustment to experimental data[END_REF].

W = W iso + n a=1 W a ani . (4.1)
Each anisotropic energy density W a ani refers to a preferred direction of the material. The number of fiber families n is generally set to 1 model tissues as ligaments or tendons while it is set to 2 to represent the behaviour of arterial walls. For example,to model the embedded collagen fibers of soft biological arterial tissues, HGO constitutive law [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF][START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF] superposes two transversely isotropic energies with two distinct preferred directions a 1 and a 2 corresponding to two fiber families:

a 1 =            cosθ sinθ 0            , a 2 =            cosθ -sinθ 0            . ( 4.2)
The phenomenological angle θ represents the angle between the collagen fibers and the circumferential direction for a strip extracted, for example, from the media of artery, as shown in Figure 4.1. To deal with the geometrical transformation with large deformation, a 1 a 2 θ θ

J = det(F) = 1. (4.3)
According to the Zhang-Rychlewski's theorem [START_REF] Zhang | Structural tensors for anisotropic solids[END_REF], the condition of material symmetry is satisfied if structural tensors are additionally included in the strain energy density representation. Transversely isotropic densities can then be expressed with the three invariants I 1 , I 2 , and I 3 of the right Cauchy-Green deformation tensor and two additional mixed invariants J 4 and J 5 [START_REF] Spencer | Isotropic polynomial invariants and tensor functions[END_REF][START_REF] Boehler | Introduction to the invariant formulation of anisotropic constitutive equations[END_REF][START_REF] Zheng | Tensors which characterize anisotropies[END_REF].

I 1 = tr(C), I 2 = tr(I 3 C -T ), I 3 = det(C), J 4 = tr(CM), J 5 = tr(C 2 M), (4.4)
where M is the so-called structural tensor representing the transverse-isotropy group and referring to a preferred direction a of the material:

M = a ⊗ a. (4.5)
It is noted that Eqs. (4.4,4.6) give

J 4 = tr(F T Fa ⊗ a) = Fa 2 , (4.6)
where the double brackets represent the usual Euclidean norm. The square root of J 4

represents thus the stretch in the fiber direction. It can also be interpreted as the radial coordinate of Fa in a cylindrical coordinate system where the polar angle γ represents the deformed angle between the collagen fibers and the circumferential direction ( In the case of hyperelastic law, there exists a strain energy density function W which is a scale function of one of the strain tensors, whose derivative with respect to a strain component determines the corresponding stress component. In the particular case of anisotropic

S =2 I 3 ∂W ∂I 3 C -1 + ∂W ∂I 1 + I 1 ∂W ∂I 2 I - ∂W ∂I 2 C + ∂W ∂J 1 4 M a 1 + ∂W ∂J 2 4 M a 2 + ∂W ∂J 1 5 (CM a 1 + M a 1 C) + ∂W ∂J 2 5 (CM a 2 + M a 2 C) , (4.7)
To uncouple the deviatoric part to the dilatational part of the response, the volume preserving part F = J -1/3 F of the deformation is introduced [START_REF] Weiss | Finite element implementation of 124 incompressible, transversely isotropic hyperelasticity[END_REF]. The modified invariants related to Ĉ = FT F = J -2/3 C are expressed from Eq.(4.4) by

Î1 = I 1 J -2/3 , Î2 = I 2 J -4/3 , Ĵa 4 = J a 4 J -2/3 , Ĵa 5 = J a 5 J -4/3 . (4.8)
The exponential type HGO density adopted in this work uses these modified invariants as follows:

W = Ŵ ( Î1 , Ĵa 4 ) + W H (J), (4.9)
Generally, soft biological tissues are assumed to be incompressible. Eq.(4.10), which was proposed by Horgan and Saccomandi [START_REF] Horgan | Finite thermoelasticity with limiting chain extensibility[END_REF], represents a penalty term added to the finite element model to account for the incompressible behavior of the material. Here, d is the material incompressibility parameter. The initial bulk modulus K is defined as

K = 2/d. W H (J) = 1 d J 2 -1 2 -lnJ , (4.10)
Ŵ ( Î1 , Ĵa 4 ) = W iso ( Î1 ) + 2 a=1 W ani ( Ĵa 4 ), (4.11)
W a ani =        k 1 2k 2 e k 2 ( Ĵa 4 -1) 2 -1 if Ĵa 4 ≥ 1 0 if Ĵa 4 < 1 . (4.12)
This energy density is case-sensitive with respect to Ĵa 4 because the case of Ĵa 4 < 1 represents the shortening of the fibers which is assumed to generate no stress. The proof of convexity of Eq.(4.12) with respect to F is given in [START_REF] Schröder | A variational approach for materially stable anisotropic hyperelasticity[END_REF][START_REF] Balzani | A polyconvex framework for soft biological tissues. adjustment to experimental data[END_REF]. The noncollagenous matrix of the media is modeled by the Yeoh model, which describes isotropic incompressible rubber-like materials [START_REF] Yeoh | Some Forms of the Strain Energy Function for Rubber[END_REF][START_REF] Renaud | The Yeoh model applied to the modeling of large deformation contact/impact problems[END_REF]. The energy density is given by:

W iso ( Î1 ) = 3 i=1 C i0 ( Î1 -3) i , (4.13)
where C i0 are material constants. In our particular case, Eq.(4.7) is reduced to

S = 2 I 3 ∂W ∂I 3 C -1 + ∂W ∂I 1 I + ∂W ∂J 1 4 M a 1 + ∂W ∂J 2 4 M a 2 . (4.14)
By deriving the energy density W with respect to the invariants (I 1 , I 3 and J a 4 ), we obtain:

∂W ∂I 1 = I -1/3 3 dW iso d Î1 , ∂W ∂I 3 = - 1 3I 3 Î1 dW iso d Î1 + 2 a=1 Ĵa 4 dW ani d Ĵa 4 + 1 2d 1 - 1 J 2 , ∂W ∂J a 4 = I -1/3 3 dW ani d Ĵa 4 , (4.15)
with

dW iso d Î1 = C 10 + 2C 20 ( Î1 -3) + 3C 30 ( Î1 -3) 2 , dW ani d Ĵa 4 = k 1 ( Ĵa 4 -1)e k 2 ( Ĵa 4 -1) 2 .
(4.16)

Numerical results

This section presents numerical examples to demonstrate the effect of anisotropic hyperelastic behaviour of soft tissues on surface adhesion. The algorithms presented above have been implemented into the in-house finite element code FER/Contact.

Stretch-release test

We investigate the adhesive anisotropic hyperelastic behaviour of soft tissues submitted to a stretch-release test between two fixed and rigid plates, where initially a tiny gap exists between the soft tissue and the rigid plate, as shown in Figure 4.3. A displacement along the positive x-axis is constantly prescribed on the right surface of the soft tissue and its left surface is fixed, causing the soft tissue to be stretched in the axial direction and expanded in the radial direction due to collagen fibers. As two fiber families are parallel to the x -z plane in our case, the extension of the soft tissue is along the y-axis, therefore upper and lower surfaces of the soft tissue will come into contact with two rigid plates during the stretching process. After stretching the soft tissue to ensure perfect adhesion (β = 1) between its contact surfaces and two rigid plates, we release the loading on the right surface to observe the adhesion behaviour during the soft tissue rebound. The test scenario allows observing consecutively two phenomena: first, the bonding process on the adhesive interface that takes place when contact is set up, then, initiation of the debonding process on the contact interface where normal separation and tangential sliding of the adhesive interface occur due to the soft tissue rebound. We investigate how the de-bonding area evolves with the soft tissue rebound, and how the evolution is affected by the fiber arrangement. Characteristics of the system are described in the following.

The soft tissue is 10 mm long with a square section of 2 × 2 mm In Figure 4.4, solid lines represent the displacement of the selected node along the y-axis and dashed lines are the displacement along the z-axis. We observe that radial extension occurs from 35 • to 45 • , indicating that their fibers have been loaded. Although the displacement U y at 50 • is negative, its U y and U z curves gradually separate, indicating that its fibers have also been loaded. However, U y and U z of 60 • have always overlapped from t = 0 s to t = 0.007 s, hence it is still in isotropic state, in other words, its fibers have not yet been loaded. According to the structural assumptions in [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF][START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: Continuum basis, computational aspects and applications[END_REF][START_REF] Zulliger | A strain energy function for arteries accounting for wall composition and structure[END_REF], the collagen is embedded as two families of fibers in the soft tissue matrix, which are symmetrically distributed with respect to the tensile direction. The embedded collagen fibers need to be rotated almost to the direction of loading before they can carry the load. This results in a significant radial extension and thus an increase in the soft tissue thickness. Additionally, due to the incompressibility constraint, the width of soft tissue reduces. As we observed in After release, we observe that the adhesive bonds break in the same direction as the rebound. In the initial phase of soft tissue rebound, from t = 0.007 s to t = 0.009 s, the propagation of adhesive fracture at 40 • is faster than 35 • , as the results of

Normal peeling of soft tissue

The second example is a normal peeling case, in which we investigate the effect of anisotropic hyperelastic behaviour on the soft tissue debonding process. The initial state of this example is a perfect adhesion (β = 1) between the soft tissue and a rigid substrate. A displacement loading along the positive y-axis is applied on the right side surface of the soft tissue and its left side surface is held in place, as shown in Figure 4.10. Here, two fiber families are equally parallel in the x -z plane, therefore, when the soft tissue is stretched axially it will expand in the y direction. According to the problem setting, the effect of variations in collagen fiber angle on debonding rates of soft tissue will be the focus of our investigation in this example. Characteristics of the system are described in the follow- For this purpose, we perform 4 sets of tensile tests on the current model. The U y and U z curves obtained are shown in Figure 4.12. We observe a similar phenomenon to that in Figure 4.4, except that the soft tissue expansion is less pronounced in this model.

According to Eq.(3.18), a greater variation of displacement during debonding will accelerate β reduction. During debonding, the radial expansion prevents the increasing of the normal displacement due to its perpendicularity to the contact surface, while the concavity leads to a continuous variation in the tangential displacement. Since the radial expansion is much smaller than the internal concavity in the current model, the U z variation is the key factor affecting the debonding rate for different fiber angles. This also confirms the distribution of U z displacements in Figure 4.11(a): the greatest U z variation at 35 • is associated with the largest debonding area.

Concluding remarks

In this chapter, we incorporated the HGO+Yeoh constitutive law and an adhesive contact model implemented by the bi-potential method and the RCC model. This approach has been tested through an adhesive contact of soft tissue during a stretch-rebound process and a normal peeling test. We investigated the effect of anisotropic hyperelastic behaviour on surface debonding by varying collagen fiber angles. The results clearly show the deformation differences in the soft tissue and on the contact surfaces. Owing to the straightforward descriptions of structural continuum constitutive models of soft tissues and the adhesive contact rules, the presented approach can be easily implemented into a finite element program. Therefore, this adhesive anisotropic hyperelastic model can be suggested to deal with the related issues in biological and medical application fields.

Unilateral cyclic loading

The adhesion effect is usually most significant in the normal direction. In order to clearly show the adhesion degradation over several loading cycles, the first example simulates the normal adhesive contact between a hyperelastic block and a rigid plate, shown in process, and the first degradation of β also occurs in this process. From the graph, we observe that the last bond breakage level is approximately 0.9, therefore during the first bonding process, the maximum value of β in the partially recoverable case is approximately 0.45 with the set degradation factor λ = 0.6, and β max in perfectly recoverable case rises to 1 as we expected. In the second bonding process, the bond break of the partially recoverable case precedes the perfect case due to its smaller β, which produces less adhesive force. For the second debonding, we also observe that β max in partially recoverable case continues to degrade from the previous cycle, whereas β max in the perfect case remains the same. distribution of the wheel remains consistent at all times, however, the deformation in the partial case diminishes with each debonding due to the adhesive degradation. Finally, it approximates to the morphology at the end of vertical compression, indicating that the adhesion effect is very small at this moment. In the second bonding, β max of the partial case can only reach 0.6 due to the degradation factor λ = 0.4, which also causes it to precede the perfect case in the second debonding. Due to the previous accumulation, the bonding-debonding process for the subsequent partially recoverable case is overall advanced.

Concluding remarks

In this chapter, a numerical formulation for contact problems with friction and partially recoverable adhesion between soft materials under large deformation is implemented by using the RCC interface model and the bi-potential resolution method. The RCC model proposes a straightforward description of the interface adhesion based on a local scalar parameter β. Both normal and tangential effects are taken into account by the adhesive interface model, involving both the process of bonding and de-bonding of the interface links. For achieving the adhesive degradation, a relationship between the maximum adhesive intensity and the number of debonding is also described. We have combined the RCC adhesive model with Blatz-Ko hyperelasticity to account for frictional contact of foam type soft material structures with partial recoverable interface under conditions of large deformation. The proposed approach has been tested on cases involving both unilateral and mixed-directional contact kinematics in perfect recoverability case and partial recoverability case. This model can be an effective tool for dealing with more complex adhesive contact problems. complex contact interfaces, we have extended the adhesive contact model of Chapter 2 with adhesive orthotropy. In Chapter 3, we proposed an orthotropic adhesion model to deal with problems of adhesive contact between hyperelastic bodies. This model has been implemented within the same contact laws of Chapter 2. The behaviour of orthotropic adhesion is described by adhesion stiffness, whose components can be expressed according to the local coordinate system. In this model, the entire bonding and debonding process of the adhesive links with the account for orthotropic interface effects is modelled. The proposed approach has been tested on cases involving both tangential and unilateral contact kinematics, which allowed emergence of orthotropic interface effects between soft bodies.

In Chapter 4, we proposed an adhesive anisotropic hyperelastic formulation to modeling soft tissues with surface adhesion. Due to the fibre-reinforced structure of collagenous soft tissues, the HGO model was used to simulate the deformation anisotropy arising from the arrangement of collagen fibres in the soft tissue, and we chose the Yeoh material model to represent the isotropic matrix in the soft tissues. The adhesive contact model, which is the same contact constitutive law as Chapter 2, implemented by bi-potential method and RCC model, based on a set of extended unilateral and tangential contact laws. The proposed approach has been tested through an adhesive contact of soft tissue during a stretch-rebound process and normal peeling. Due to the radial extension of soft tissue caused by collagen fiber arrangement, the effect of material anisotropy on the adhesive behaviour of contact interface has been investigated. [START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF]. Here, "NTS" refers to "node-to-segment" contact; "AR" to the technique of area regularization; "ME" to moment equilibrium; "AL" to augmented Lagrangian and "VTS" to the " Virtual-slave-node-To-Segment" approach.

As has been extensively studied by Zavarise et al. in [START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF] and recalled in Fig. (A.2), classical NTS contact algorithms, especially those using one-pass approaches introduce significant errors to contact stresses evaluation on non-conforming meshes. To obtain acceptable behaviours using classical NTS description, it is necessary to implement twopass sequential schemes in conjunction with Lagrangian multiplier method, or, develop improved one-pass schemes, for example the VTS ("virtual-node-to-segment") method.

VTS method extends the classical NTS approach by considering additional virtual slave nodes on the slave surface, leading to augmented slave segments.

In . "NTS" refers to "node-to-segment" contact; "AR" to the technique of area regularization; "ME" to moment equilibrium; "AL" to augmented Lagrangian and "VTS" to the " Virtual-slave-node-To-Segment" approach. The comparison highlights our result ("NTS-AL") among existing established methods, based on results reported in [START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF]. Title: Modelling of adhesive frictional contact problems for soft matters Keywords: Dynamic contact; Recoverable adhesion; RCC model; Bi-potential method; Hyperelastic materials Abstract: Dynamic frictional contact and recoverable interfacial adhesion involving soft materials represent a frequent contact phenomenon. In the numerical modelling field, constructing of a model capable of addressing contact with friction and recoverable adhesion has always been a challenging topic. In this thesis, an extended, readyto-implement 3D model for quasi-industrial problems of contact with friction and recoverable interface adhesion between soft material is formulated using the Raous-Cangémi-Cocou (RCC) interface model and a bi-potential based resolution method. The RCC model proposes a straightforward description of the interface adhesion based on a local scalar parameter, both normal and tangential effects are taken into account by the adhesive interface model, involving both the process of bonding and de-bonding of the interface links. This adhesive model has been implemented within the bi-potential method, based on a set of extended unilateral and tangential contact laws. We combine the 3D extended adhesive interface model with different hyperelastic models to investigate large deformation contact problems under various adhesive interface conditions. Such as, Blatz-Ko material model for large deformation contact problems under isotropic and orthotropic adhesive interface conditions; The mechanical behaviour of biological soft tissues with surface adhesion is investigated by using the Holzapfel-Gasser-Ogden (HGO)+Yeoh anisotropic hyperelasticity model. To illustrate the capability of the implemented model, we set up various test cases in each chapter to explore adhesive contact in normal, tangential and mixed directional scenarios for different material models and interface conditions, which brings us closer to quasi-industrial modelling situations.

Figure 1 . 1 :

 11 Figure 1.1: (a) Schematic of hydrogen bonds between polyvinyl alcohol (PVA)[1]; (b) Schematic of Van Der Waals force

Figure 1 . 2 :

 12 Figure 1.2: Structural hierarchy of the gecko adhesive system[2]

Figure 1 . 3 :

 13 Figure1.3: Bio-mimetic mushroom-shaped fibrillar adhesive microstructure[3]

Figure 1 . 4 :

 14 Figure 1.4: Graphic representation of adhesive intensity β

Figure 1 . 5 :

 15 Figure 1.5: Coulomb's cone K µ and its polar cone K * µ

(

  HGO) model incorporating interfacial adhesion to investigate the effect of anisotropic hyperelastic behaviours of soft tissues on surface adhesion. The HGO anisotropic hyperelastic constitutive law is widely used to model collagen fiber reinforced biological soft tissues, its anisotropy arises from collagen fiber behaviour, and the matrix of soft tissues is isotropic in nature. The energy densities of collagenous soft tissues contain both isotropic and anisotropic components, each collagen fiber family has an independent anisotropic energy density. In this work, the Yeoh hyperelastic constitutive law is adopted to model noncollagenous matrix of soft tissues. We use the same adhesive contact constitutive law constructed in Chapter 2, and then introduce the HGO anisotropic hyperelastic model into the contact model. Numerical examples are performed to demonstrate the effect of material anisotropy on surface adhesion.In reality, the adhesive effect of adhesive tapes decreases with frequent contact. However, the contact model with recoverable adhesion and friction mentioned in Chapter 2 describes a perfectly reversible adhesion, which means that the adhesion intensity β can always achieve 1 with a sufficient contact time. In response to more complex adhesive contact problems, Raous et al. propose a relationship between the β maximum value and the bond breakage level in adhesive bonds, which regulates the amount of adhesive degradation through the product of a scalar parameter and the bond breakage level per debonding process. In Chapter 5, we incorporate this relationship into our adhesive contact law for modelling adhesive degradation under cyclic loading. Numerical examples are performed to demonstrate the effect of adhesive degradation by comparing perfectly recoverable cases with partially recoverable cases.

Figure 2 . 2 :

 22 Figure2.2: Modified Signorini law with adhesion: graphic representation of normal adhesion forces and the level of damage evolving with contact distance. An empirical limit of β, denoted by ε, is adopted. With β < ε (ε = 10 -4 ), adhesion bonds are considered broken.

9 log 10 Figure 2 . 3 :

 91023 Figure 2.3: Evolution of adhesion forces as function of β on logarithmic scale: with β on the order of 10 -4 , adhesion forces become negligible on the 10 -9 N order.

Figure 2 .

 2 Figure 2.3 shows the evolution of adhesion forces as function of β on logarithmic scale. Starting from Point O, adhesion force first increases with the gap between contact surfaces, then due to the surface debonding, very quickly drops to insignificant levels.

Figure 2 . 5 :-

 25 Determine mass matrix M and damping matrix C -Read the data -Time iteration -Determine F ext -Detect contact conditions in local frame -Newton iteration -Compute K and F int -Compute K and F -Solve KΔu= F -Compute r, r and β by bipotential and adhesive model in local frame, R=H T (r+r) -Solve KΔu= F+R -Actualize u=u+Δu -Check convergence criteria if not convergence -Compute velocity -Gather element nodal displacement -Compute stress and strains

Figure 2 . 5 :

 25 Figure 2.5: The iterative solution procedure

Figure 2 .

 2 6(b) is the load curve showing the displacement of the upper surface of the sphere.

Figure 2 . 6 :

 26 Figure 2.6: Indentation on a hyperelastic material with adhesive surface: (a) Problem set; (b) Loaded displacement on the upper surface of the sphere.

Figure 2 .

 2 Figure 2.7(e) lists three different sets of adhesive parameters used in the test, whose results are reported in Figure 2.7(c). Cocou et al. investigated similar scenarios and obtained concordant results [41].Figure 2.7(a) shows the evolution of β on 7 contact nodes in

Figure 2 .

 2 7(a) shows the evolution of β on 7 contact nodes in Case 1. On any contact point, its horizontal distance from the center point O determines the time duration of the contact process on this point, involving bonding and de-bonding.

Figure 2 .

 2 Figure 2.7(b) shows the normal adhesion force of the 7 contact nodes in Case 1. Similarly, the distance from the center point O determines the sequence of separation, which however does not influence the adhesion force at the moment of separation.Figure 2.7(c)

Figure 2 .

 2 7(c) shows evolution of β on the contact point A under 3 groups of different adhesive parameters.

Figure 2 .Figure 2 . 7 :

 227 Figure 2.7: Indentation on a hyperelastic material with adhesive surface: (a) Evolution of β calculated on 7 contact nodes of the hyperelastic block. At t = 0 s, the only contact point between the indentor and the block is point O, see Figure 2.6(a), and d represents the horizontal distance between each node and the center point O; (b) Evolution of the adhesion force Rn calculated on 7 contact nodes of the hyperelastic block; (c) Evolution of β with 3 different sets of adhesive parameters (w and C n ) on node A (d = 3.4 mm); (d) Evolution of the adhesion force Rn based on 3 different sets of adhesive parameters, calculated on node A; (e)Table of the tested adhesive parameter sets.

  3 different shear modulus G = 5 × 10 6 Pa, 1 × 10 7 Pa, and 1.5 × 10 7 Pa for the hyperelastic wheel. Concerning the interface properties, the following parameters are used: friction coefficient µ = 0.4. Note that setting non-zero friction here is important to drive the wheel to rotate. The wheel rotates consequently under the combined effects of interface friction and adhesion. Parameters for the adhesive are : w = 20 J.m -2 , C n = C t = 2 × 10 9 N.m -3 , and b = 0.1 N.s.m -1 .

Figure 2 .

 2 Figure 2.8(c) shows the morphology and Von Mises stress distribution of the hyperelastic wheel just following application of the compression by the upper plate.Figure 2.8(d)

  Figure 2.8(d)

  Figure 2.8(c).

Figure 2 .

 2 Figure 2.9(a) shows the evolution of the adhesion parameter β, and the adhesion forces calculated on 3 selected nodes as indicated in Figure 2.8(b) with shear modulus G = 1 ×

Figure 2 . 8 :

 28 Figure 2.8: Rolling adhesion and friction of a confined hyperelastic wheel: (a) Problem setup; (b) The initial state of hyperelastic wheel; (c) Shape and Von Mises stress distribution of the confined wheel just before rotation; (d) Deformed shape and Von Mises stress distribution of the hyperelastic wheel during its rotation.

Figure 2 .Figure 2 . 9 :

 229 Figure 2.9: Rolling adhesion and friction of a confined hyperelastic wheel: (a) Evolution of β calculated on 3 selected nodes on the wheel surface. Locations of the 3 nodes are indicated in Figure 2.8(b). In this case, shear modulus G = 1 × 10 7 Pa; (b) Evolution of β calculated on the first node with 3 sets of shear modulus (G = 5 × 10 6 Pa, 1 × 10 7 Pa, and 1.5 × 10 7 Pa); (c) Traction-separation curves of the first selected node with 3 sets of shear modulus; (d) Evolution of the normal adhesion force Rn with time on the first selected node under 3 sets of shear modulus G.

  Figure 2.11(a)). Then for each group, we test on 5 different friction coefficients µ, so as to investigate the combined effect of friction and adhesion on the tangential interface behaviour. The adhesive parameters used in the simulations are: w = 20 J.m -2 , C n = 2 × 10 9 N.m -3 and b = 0.1 N.s.m -1 .

:Figure 2 . 10 :

 210 Figure 2.10: Adhesive friction between a hyperelastic plate and a deformable semi-cylinder: (a) Problem setup; (b) Distribution of Von Mises stress of the substrate and the cylinder at the end of push down; (c) Distribution of Von Mises stress at the end of calculation.

Figure 2 .

 2 Figure 2.10(b) and (c) present the Von Mises stress distributions of the sliding system, respectively at the onset of sliding, and during the sliding process. We post-process the

Figure 2 .

 2 11(b) 5 simulations based on 5 values of µ ranging from 0 to 0.8. All the 5 cases consider the same plate descent velocity of 1000 mm/s (Case 1 of Figure 2.11(a)). Results in Figure 2.11(b) indicate the formation of stronger bond (higher β) on rougher surfaces (greater µ). This can be interpreted by the fact that a rougher surface (higher µ) delays the onset of sliding motion, according to the Coulomb friction model, which results in longer time for better bonding of adhesive links. Therefore, we observe a concordant trend on the curves reflecting tangential adhesion forces. With greater friction coefficient (Figure 2.11(d)), the onset of de-bonding is delayed, creating increased level of adhesion force. Then, lower descent velocity on the first stage (Figure 2.11(e)) also creates the effect of delaying the onset of de-bonding, permitting better bonding and more significant adhesion forces.

Figure 2 . 11 :

 211 Figure 2.11: Adhesive friction between a hyperelastic plate and a deformable semi-cylinder: (a) 5 cases with different time and velocity of push down. This kind of setting is to ensure that the substrate has the same displacement during the whole process; (b) Evolution of β of blue node (see Figure 2.10(a)) with different friction coefficients µ in Case 1; (c) β evolution of blue node with different descent velocities under µ = 0; (d) Evolution of the tangential adhesion force Rt of blue node with different friction coefficients µ in Case 1; (e) Evolution of the tangential adhesion force Rt of blue node with 5 cases under µ = 0.

3Figure 2 . 12 :

 212 Figure 2.12: 3D adhesive frictional twisting: Problem setup and loading sequence (Step 1, compression and adhesion process; Step 2, twisting and de-bonding process)

Figure 2 .

 2 Figure 2.13 compares the normal reaction forces of the two cases during the twist process.9 frames of results are extracted in chronological order to represent the evolving twist process. On each frame, we compare distributions of the normal force calculated on the contact interface between dry friction and adhesive friction. Since the combined motion

Figure 2 . 14 presentsFigure 2 . 13 :

 214213 Figure 2.13: 3D adhesive frictional twisting: Evolution of the interface normal forces Rn distribution during the twisting process. Comparison between dry friction (left) and adhesive friction (right) on each frame of time. For both cases, friction coefficient µ = 0.4 is used.

Figure 2 . 16 ,Figure 2 . 14 :

 216214 Figure 2.16, Euclidean norm of tangential forces are presented, allowing us to observe the evolving intensity of tangential forces on the contact interface. In this figure, similar distribution of tangential forces can be observed on frames µ = 0.2 / t = 0.006 s, µ =

Figure 2 . 15 :Figure 2 . 16 :

 215216 Figure 2.15: 3D adhesive frictional twisting: Evolution of the adhesion intensity during the twisting load, cacluated using 3 sets of friction coefficients µ = 0.2, 0.4 and 0.6.

  by using the RCC interface model and the bi-potential resolution method. The RCC model proposes a straightforward description of the interface adhesion based on a local scalar parameter, and enables coupling the effect of adhesion, friction and unilateral contact within a unified framework. Both normal and tangential effects are taken into account by the adhesive interface model, involving both the process of bonding and debonding of the interface links. We have combined the 3D extended RCC adhesive interface model with 3D Blatz-Ko hyperelasticity to account for frictional contact of foam type soft material structures with recoverable interface under conditions of large deformation. To illustrate the ability of the implemented model to deal with real problems, we have tested various 3D test cases involving normal, tangential, and mixed-directional scenarios of adhesive contact with/without adhesion recoverability, which is very close to quasi-industrial modelling situations.

Figure 3 . 1 : 2 Figure 3 . 2 :

 31232 Figure 3.1: Finite deformation contact

  (3.14,3.15) include: C t1 , C t2 and C n : parameters characterizing the initial adhesive stiffness when adhesion is complete, w: decohesion energy threshold, : indicator function that assures unilateral contact (g n 0), and meaningful values of the degree of adhesion. The subscript Q indicates Q = {η | 0 η 1}, µ: friction coefficient, b: surface viscosity.

. 21 )

 21 With the consideration of interface adhesion, tangential forces are contributed by two mechanisms. The first mechanism is comparable to static friction by the classical Coulomb model. It vanishes once slip occurs. The second, arising from the effects of interface adhesion and defined by Eq.(3.21), gives rise to adhesive tangential force rα t which emerges with surface slip. The obtained rule of tangential contact with interface adhesion (Eqs.(3.20,3.21)) can be graphically interpreted by Figure 3.4. By assuming perfect adhesion (β = 1) at g t = 0, the

Figure 3 . 5 :Figure 3 . 6 )

 3536 Figure 3.5: Orthotropic adhesion of a soft body under compression on a rigid surface

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Orthotropic adhesion under compression: (a) Evolution of β calculated on 6 nodes on the contact interface, along the diagonal from the centre to the periphery; (b) Evolution of β on the contact interface and variation in the shape of the contact surface in debonding process. In each square area, the colour progresses from dark red to blue, which represents the damage of the adhesive strength β from perfect adhesion (β = 1) to complete separation (β = 0).

  .8, the elastomer cylinder is 2 mm high, and has a radius of 5 mm. The elastomer is modelled by Blatz-Ko material with shear modulus G = 2.1 × 10 5 MPa. The adhesive interface parameters are: w = 100 J.m -2 , C tx = 5 × 10 9 N.m -3 , C ty = 1 × 10 10 N.m -3 and b = 0.1 N.s.m -1 . The simulation scenario involves 2 stages. On the first stage, we prescribe a slight compression on the elastomer by descending its upper surface by 0.1 mm after contact. The compression activates the bonding process which leads to complete bonding on the adhesive interface. On the second stage, a lateral motion at the velocity of 0.1 m/s is applied on the cylinder's upper surface. Under the tangential effect on the contact interface, de-bonding is initiated and progresses until the rupture of adhesive bonds, which allows the cylinder to slide on the support surface. A group of 10 tests have been performed. On each test, we align the lateral motion to a new direction whose angle with respect to x axis, θ, increases from 0 • to 90 • by increments of 10 • .

Figure 3 .Figure 3 . 8 :

 338 Figure 3.9 presents the evolution of adhesion parameters calculated on the centre node that belongs to the contact surface of the elastomer cylinder, for the 10 calculations performed with θ ranging from 0 • to 90 • . Positions of the centre node at the moment of adhesion rupture are reported in Figure 3.9(a). Blue circles represent results based on orthotropic adhesion properties with C tx = 0.5C ty . Red circles are obtained considering the assumption of isotropic adhesion. For the isotropic cases, all the red circles are arranged

Figure 3 Figure 3 . 9 :

 339 Figure 3.9(d)).

C tx = 5 ×Figure 3 . 10 :

 5310 Figure 3.10: Comparison between isotropic and orthotropic adhesive twisting: Problem setup and loading sequence (Step 1, compression and adhesion process; Step 2, twisting and de-bonding process).

tFigure 3 . 11 :

 311 Figure 3.11: Comparison between isotropic and orthotropic adhesive twisting: Evolution of the adhesion intensity β in isotropic case and orthotropic case during the debonding process and their shape variation of the contact surface. In each square area, the colour progresses from dark red to blue, which represents the damage of the adhesive strength β from perfect adhesion (β = 1) to complete separation (β = 0).

  t= 0.006 s t= 0.007 s t= 0.008 s t= 0.009 s t= 0.01 s anisotropy (C tx = 5C ty ) isotropy (C tx = C ty )

  This chapter describe a numerical formulation based on the Holzapfel-Gasser-Ogden (HGO) model incorporating interfacial adhesion to investigate the effect of anisotropic hyperelastic behaviours of soft tissues on surface adhesion. The HGO anisotropic hyperelastic constitutive law is widely used to model collagen fiber reinforced biological soft tissues, its anisotropy arises from collagen fiber behaviour, and the matrix of soft tissues is isotropic in nature. The energy densities of collagenous soft tissues contain both isotropic and anisotropic components, each collagen fiber family has an independent anisotropic energy density. In this work, the Yeoh hyperelastic constitutive law is adopted to model noncollagenous matrix of soft tissues. We use the same adhesive contact constitutive law constructed in Chapter 2, and then introduce the HGO anisotropic hyperelastic model into the contact model.In the follows: in Section 4.2, we provide the formulation of the anisotropic hyperelastic constitutive law incorporating Yeoh hyperelastic model. Then the complete finite element formulation, including the resolution algorithm, is provided. To validate the framework, we present numerical examples in Section 4.3. In the end, a few concluding remarks are drawn in Section 4.4.

Figure 4 . 1 :

 41 Figure 4.1: Angle θ between the collagen fibers of the artery wall and the circumferential direction

Figure 4 1 Figure 4 . 2 :

 4142 Figure 4.2: Cylindrical coordinate system

2 ,

 2 and it has an initial gap of 0.06 mm from the rigid plate. The material parameters used for Yeoh part of the HGO model correspond to a skin model[117]: C 10 = 26912.5 kPa, C 20 = 37606.5 kPa, C 30 = 41596.3 kPa, k 1 = 996.6 kPa and k 2 = 524.6. Adhesive interface parameters are: w = 100 J.m -2 , C n = C t = 1×10 10 N.m -3 and b = 0.1 N.s.m -1 . The soft tissue is stretched by 0.7 mm in 0.007 s before releasing load.Before exploring the adhesive anisotropic hyperelastic behaviour, we investigate the effect of fiber arrangements on radial extension intensity by a conventional tensile test of

Figure 4 . 3 :

 43 Figure 4.3: Problem setup: in this case, collagen fibers parallel to the x -z plane; The left surface of soft tissue is fixed and we apply a load along the positive direction of x-axis on the right surface, releasing after the soft tissue has been stretched to a certain length.

Figure 4 . 4 :

 44 Figure 4.4: Tensile test of soft tissue (no upper and lower rigid plates): lateral displacements U y (solid lines) and U z (dashed lines) of a selected node (blue point) with 5 different fiber angles (35 • , 40 • , 45 • , 50 • , 60 • ), which is on the edge of the upper surface as shown in figure.

Figure 4 . 4 ,

 44 Figure 4.4, smaller fiber angle θ means closer to the tensile direction, therefore the radial extension of 35• occurs fastest and has the largest displacement U y at t = 0.007 s, the extension intensity decreases with increasing fiber angle θ, then, in the z-axis direction, the degree of concavity decreases with increasing θ as well, due to the incompressibility of tissues.

Figure 4 .

 4 Figure 4.5(a) shows U y evolution of 3 sets of selected fiber angles in the x -y plane view. During the loading phase (from t = 0 s to t = 0.007 s), we observe that the soft tissue with 35 • fiber angle expands faster than 40 • , while the soft tissue with 60 • fiber angle is depressed inwards in the y-axis direction, which is consistent with the results in Figure4.4.At the end of the stretching (t = 0.007 s), majority of the contact surfaces of the 35 • and 40 • cases touch the rigid plate. After releasing the load (from t = 0.007 s to t = 0.015 s), we observe that the rebound rate in the y-axis direction is inversely proportional to the radial extension rate. The reason is that, during stretching process, the march of soft tissue in extension direction is blocked by rigid plates, hence the work done for extension is converted into elastic potential energy of matrix, which is accumulated in the contact area of the soft tissue. As shown in Figure4.4, the radial extension of 35 • case is greater than 40 • , therefore the soft tissue with 35 • fiber angle needs to release more potential energy during rebound, which causes its U y to decrease slowly. At t = 0.015 s, we observe that the U y of 35 • and 40 • cases is still visible, while the soft tissue with 60 • has almost completely rebounded in the y-axis direction, which is the result of adhesion hindering the interface separation. Due to combination of the adhesive effect and the difference of rebound rate, the radial extension area of 35 • case is larger than 40 • at t = 0.015 s. In order to visualize the evolution of displacements in the x, y and z directions throughout the whole load-release process, we choose the same observation node of Figure4.4 and plot its displacement curves at three selected fiber angles, as shown in Figure4.5(b). Comparing U x curves, 35 • and 40 • cases are perfectly coincided during loading process, however, U x of 60 • is slightly higher than other two angles due to the lack of interference from contact interface. In rebound phase, we can observe that the disparity in U x between contacted cases and uncontacted case is immediately apparent, this phenomenon is not clearly observed in Figure4.5(a). The displacement evolution in y-axis or z-axis directions is expected, which shows the same tendency as in Figure4.4 in loading process and the identical results as observed in Figure4.5(a) in rebound stage.

Figure 4 .

 4 Figure 4.6(a) shows the evolution of adhesion intensity β on upper contact interface with 2 different sets of fiber angles (35 • , 40 • ) in rebound process. Both contact areas of two fiber angles are perfectly adhered with the rigid plates prior to release loading.

Figure 4 .

 4 5(a), the U y of 40 • decay more rapidly than 35 • due to a larger displacement variation, which accelerates the adhesive damage of soft tissue in 40 • case. From t = 0.011 s, the adhesive damage rate of 35 • gradually catches up, and its damage propagation line becomes sharper compared to 40 • . At t = 0.015 s, the adhesive damaged area of 35 • has overtaken 40 • . This phenomenon seems to be contrary to U y variation of these two fiber angles, the reason is that, in addition to the displacement variation in y-axis, which affects normal separation, the tangential sliding in z direction of the interface also breaks the adhesive bond. During soft tissue rebound, the U z gradually increases from the edge to the center, hence the β at the edge of the contact surface is the first to decrease, which create the curvature of damage propagation lines. Due to incompressibility, the larger displacement extension in the y-axis leads to a deeper concavity in the z-axis, therefore the U z variation of 35 • is larger than 40 • , as we can observed in the curves of Figure4.5(b). Under the dual actions of U y and U z , the soft tissue with 35 • fiber angle has a larger adhesive fractured area, and its damage propagation line is sharper due to a greater U z . Figure 4.6(b) shows the β evolution for the same observation node in Figure 4.4 at fiber angles of 35 • and 40 • . At t = 0.0068s, β increases to 1, indicating the achievement of complete bonding (Zone 1 in

Figure 4 .

 4 Figure 4.6(b)) of the adhesive interface. Comparing the two curves, since the soft tissue with 35 • fiber angle is the first to come into contact with the rigid plate, it grows faster than 40 • to β = 1 in Zone 1. In Zone 2, we release the loading, the soft tissue starts to rebound. As the observation point is in the middle of the edge, it takes some time for the damage propagation line to advance to that position, therefore the β of observation node is still 1. As the rebound of soft tissue proceeds, the fracture of adhesive bond advances

Figure 4 .

 4 Figure 4.7 shows the evolution of adhesive force in x-axis (R adh-x ) on upper contact interface with 2 different sets of fiber angles (35• , 40 • ) in rebound process. We note that R adh-x slides from right to left, which is consistent with the rebound direction in x-axis.And the adhesive forces shrink along the edge towards the center. The reason is the same as the curvature variation of the damage propagation line, β decreases from the outside to inside due to the release direction of U z during rebound, causing the R adh-x to shrink inwards. In addition, we observe that R adh-x of 40 • advances faster until t = 0.012s. Due to a faster reduction of β at 40 • in the early stages of rebound, as observed in Figure4.6, its adhesive forces in x-axis move quickly. Then, R adh-x of 35 • catches up as the effect of U z in 35 • case gradually becomes apparent.

Figure 4 .

 4 Figure 4.8 shows the distribution of adhesive forces in y-axis (R adh-y ) on upper contact interface with 2 different sets of fiber angles (35 • , 40 • ) in rebound process. We observe that R adh-y is barely visible on the contact surface of 35 • case until t = 0.011s, as the U y variation of 35 • is very subtle in the early stages of rebound, as shown in Figure 4.5(a).After t = 0.011s, as the upper surface begins to drop, R adh-y of 35 • gradually become visible in the area where β is still strong. The soft tissue with 40 • fiber angle is much weaker in radial extension, which allows its adhesive forces in the y-axis to play a role throughout the rebound process. The advance directions of R adh-y at both fiber angles are also consistent with the damage propagation line.

Figure 4 .

 4 Figure 4.9 shows the distribution of adhesive forces in z-axis (R adh-z ) on upper contact interface with 2 different sets of fiber angles (35 • , 40 • ) in rebound process. The adhesive forces at both angles in the z-axis direction are evident from the beginning of the rebound.

Figure 4 . 5 :Figure 4 . 6 : 2 R

 45462 Figure 4.5: (a) U y displacement evolution of 3 different sets of fiber angles (35 • , 40 • , 60 • ) in the x -y plane view; (b) The whole displacements U x , U y and U z of the same contact node of Figure 4.4 with 3 selected fiber angles.

Figure 4 . 7 :

 47 Figure 4.7: Evolution of adhesive force R adh-x on upper contact interface with 2 different sets of fiber angles (35 • , 40 • ) in rebound process

Figure 4 . 8 :

 48 Figure 4.8: Evolution of adhesive force R adh-y on upper contact interface with 2 different sets of fiber angles (35 • , 40 • ) in rebound process

  ing. The soft tissue dimension is 10 × 3 × 0.5 mm. The HGO+Yeoh material parameters are: C 10 = 26912.5 kPa, C 20 = 37606.5 kPa, C 30 = 41596.3 kPa, k 1 = 996.6 kPa and k 2 = 524.6. Adhesive interface parameters are: w = 100 J.m -2 , C n = C t = 1 × 10 10 N.m -3 and b = 0.1 N.s.m -1 . The right surface of the soft tissue is elevated by 3 mm in 0.015 s.

Figure 4 . 10 :

 410 Figure 4.10: Problem setup: in this case, collagen fibers parallel to the x -z plane;The left surface of soft tissue is fixed and we apply a normal load on the right surface to peel off the soft tissue. In the initial state, the contact surface of the soft tissue adheres perfectly to the substrate (β = 1).

Figure 4 .

 4 Figure 4.11(a) shows the deformation shape and U z displacement distribution of the soft tissue at t = 0.015s for 4 different fiber angles of 35 • , 40 • , 45 • , 50 • . For a clearer comparison of the debonding between 4 cases, we overlapped the deformed shape contours of 4 sets of soft tissue, as shown in Figure 4.11(b) and Figure4.11(c) respectively. We observe that the debonding contact area decreases with increasing fiber angle. According to previous results, a smaller fiber angle allows the soft tissue to reveal the anisotropic behaviour more rapidly. The anisotropy results in soft tissue expansion along the y-axis and accelerated inward concavity in the z-axis. However, which variable directly accelerates the breakage of the interfacial adhesive bond in the current example?

Figure 4 .

 4 Figure 4.13(a) shows β evolution of the contact surface at four different sets of fiber angles. We observe the same phenomenon as in Figure 4.6, where the adhesive bonds break fastest at 35 • and the debonding contact area is the largest and the debonding speed decreases with increasing angle. The combined effect of the tangential and normal displacements makes the bond break faster at the edges of the contact surface than in the central region, which results in an arc-shaped adhesive fracture propagation line.Figure

Figure

  4.13(b) shows β evolution curves of the selected point at 4 different sets of fiber angles, where the observation point is located in the middle of the contact surface edge (the red node shown in the figure). The sequence of decreasing β in the graph also corroborates the tendency in Figure4.13(a).

Figure 4 .

 4 [START_REF] Beisl | Confined evaporation-induced selfassembly of colloidal lignin particles for anisotropic adhesion[END_REF] shows the distribution of adhesive forces along the y-axis R adh-y on the contact surface. We observe that the adhesive forces propagate from right to left and are mainly apparent around the adhesive fracture propagation line. The propulsion of adhesive force is also consistent with β evolution for different fiber angles, i.e. 35 • is the fastest and 50 • the slowest.

Figure 4 . 12 :

 412 Figure 4.12: Tensile test of soft tissue (no contact): lateral displacements U y (solid lines) and U z (dashed lines) of a selected node (blue point as shown in figure) with 4 different fiber angles (35 • , 40 • , 45 • , 50 • ).

Figure 4 . 13 : 93 tFigure 4 . 14 :

 41393414 Figure 4.13: (a) β evolution on the contact interface with 4 different sets of fiber angles; (b) The evolution curves of β calculated on the red selected node (as shown in figure) for 4 different sets of fiber angles

Figure 5 . 1 :

 51 Figure 5.1: Comparative schematic between complete recoverability and recoverability degradation in β evolution

Figure 5 .Figure 5 . 2 :

 552 Figure 5.2: Interfacial adhesion fatigue in unilateral contact of elastomer: (a) Problem set; (b) Loaded displacement on the upper surface of elastomer.

Figure 5 .

 5 Figure 5.3(a) shows β evolutions of observation point, as shown in figure, in two different recoverable adhesion scenarios, where the blue line represents a perfect recoverability and the red line represents a partially recoverability. Since the initial contact interface is at perfect adhesion (β = 1), two lines of β coincide completely during the first debonding

Figure 5 . 3 (Figure 5 . 3 :

 5353 Figure 5.3: Interfacial adhesion fatigue in unilateral contact of elastomer: (a) Evolution of β in perfectly recoverable case and in partially recoverable case (λ = 0.6); (b) Evolution of normal adhesive force Rn in perfectly recoverable case and in partially recoverable case (λ = 0.6).

Figure 5 . 5 :Figure 5 .

 555 Figure 5.5: Adhesive fatigue in rolling adhesion of a hyperelastic wheel: Deformed shape and Von Mises stress distribution of the hyperelastic wheel in perfectly recoverable case and partially recoverable case (λ = 0.4)

Figure 5 . 6 :

 56 Figure 5.6: Adhesive fatigue in rolling adhesion of a hyperelastic wheel: Evolution of β of yellow dot (see Figure 5.4) in perfectly recoverable case and partially recoverable case (λ = 0.4)

In Chapter 5 ,Figure A. 1 :

 51 Figure A.1: Magnified contact interface configuration with and without surface penetration: comparison of the present contact algorithm ("NTS-AL") to other algorithms based on results reported in[START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF]. Here, "NTS" refers to "node-to-segment" contact; "AR" to the technique of area regularization; "ME" to moment equilibrium; "AL" to augmented Lagrangian and "VTS" to the " Virtual-slave-node-To-Segment" approach.

  Figure A.2: Contact patch test: comparison of several contact algorithms regarding the interface normal stresses. "NTS" refers to "node-to-segment" contact; "AR" to the technique of area regularization; "ME" to moment equilibrium; "AL" to augmented Lagrangian and "VTS" to the " Virtual-slave-node-To-Segment" approach. The comparison highlights our result ("NTS-AL") among existing established methods, based on results reported in[START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF].

  contact area regularization ("AR"), and the improved VTS method proposed by the work of Zavarise et al. We observed satisfactory contact geometry in Fig.(A.1.f) and the same level of accuracy as VTS method in Fig.(A.2) which confirm the capacity of augmented Lagrangian methods in enforcing geometrical relations of contact surfaces and improving the computational accuracy. 113 Titre : Modélisation des problèmes de contact frictionnel adhésif pour les matières molles Mots clés : Contact dynamique; Adhérence récupérable; modèle RCC; Méthode bi-potentielle; Matériaux hyperélastiques Résumé : Le contact frictionnel dynamique et l'adhésion interfaciale récupérable impliquant des matériaux mous représentent un phénomène de contact fréquent. Dans le domaine de la modélisation numérique, la construction d'un modèle capable de traiter le contact avec frottement et l'adhésion récupérable a toujours été un sujet difficile. Dans cette thèse, un modèle 3D étendu, prêt à être implémenté pour les problèmes quasi-industriels de contact avec friction et d'adhésion récupérable entre matériaux mous est formulé en utilisant le modèle d'interface Raous-Cangémi-Cocou (RCC) et une méthode de résolution basée sur le bi-potentiel. Le modèle RCC propose une description simple de l'adhésion de l'interface basée sur un paramètre scalaire local. Les effets normaux et tangentiels sont pris en compte par le modèle d'interface adhésif, impliquant à la fois le processus de collage et de décollage des liens de l'interface. Ce modèle adhésif a été implémenté dans la méthode bi-potentielle, basée sur un ensemble de lois de contact unilatérales et tangentielles étendues. Nous combinons le modèle d'interface adhésive étendu 3D avec différents modèles hyperélastiques pour étudier les problèmes de contact à grande déformation dans diverses conditions d'interface adhésive. Par exemple, le modèle de matériau de Blatz-Ko pour les problèmes de contact à grande déformation dans des conditions d'interface adhésive isotrope et orthotrope ; le comportement mécanique des tissus mous biologiques avec adhésion de surface est étudié à l'aide du modèle d'hyperélasticité anisotrope Holzapfel-Gasser-Ogden (HGO)+Yeoh. Pour illustrer la capacité du modèle mis en oeuvre, nous avons mis en place divers cas d'essai dans chaque chapitre pour explorer le contact adhésif dans des scénarios de direction normale, tangentielle et mixte pour différents modèles de matériaux et conditions d'interface, ce qui nous rapproche de situations de modélisation quasi-industrielles.

  Maison du doctorat de l'Université Paris-Saclay 2 ème étage aile ouest, Ecole normale supérieure Paris-Saclay 4 avenue des Sciences, 91190 Gif sur Yvette, France

  lie in the local plane (T 1 , T 2 ), their projection in the local system gives rise to expressions of tangential displacement and forces according to axis T 1 and T 2 . In the case of isotropic tangential behaviour, the adhesion stiffness can be described by a unique parameter C t . Therefore, vectors of tangential forces lie on a circle of radius C t x t β 2 . For any given slip value x t , one distinct circle can be drawn, which graphically leads to a conic representation of the adhesion force by swiping x t from 0 to +∞ as shown in Figure2.4.

	r t2	r t
	Fully bonded adhesion
	r t1
	r t2
	Figure 2.4: Modified Coulomb rule with adhesion: evolution of tangential adhesive forces
	and the level of damage vs. slip
	vector rα

t

  .[START_REF] Liu | Controlled adhesion anisotropy between two rectangular grooved surfaces[END_REF] where the second Piola-Kirchhoff stress tensor S, in the case of Blatz-Ko material model is given in Section 2.2.5 by Eq.(2.41). The vector of contact reaction force R is expressed in the global coordinate system. It is obtained by considering Eqs.(2.7,2.8,2.10) and includes in particular contributions due to adhesion: r determined according to the contact and friction rules given in Section 2.2.3.Other notations in Eq.(2.45) include V 0 , volume of the initial configuration; F

	R = H T (r + r),	(2.46)
	with r and

ext , vector of external loads; M, mass matrix; A, damping matrix; u, vector of velocity, and ü, vector of acceleration. Substituting δE from Eq.(2.44) into Eq.(2.45) results in

  Table of the tested adhesive parameter sets.

Remerciements

List of Figures

any given slip value g t , one distinct ellipse can be drawn, which graphically gives a conic representation of the adhesion force by swiping g t from 0 to +∞ as shown in Figure 3.4.

Numerical results

The algorithm presented above has been implemented within the in-house finite element code FER/Contact. In this section, three numerical examples based on contact simulations are presented to show orthotropic behaviours of the adhesive contact interface.

Orthotropic adhesion under compression

In this first example, we investigate the orthotropic interface adhesion of a hyperelastic soft body submitted to compressive load against a rigid surface. As shown in Figure 3.5, a vertical displacement is constantly prescribed on the upper surface of the soft body, pressing it against a fixed, rigid plate. The test scenario allows observing consecutively two phenomena: first, the bonding process on the adhesive interface that takes place when contact is set up, then, initiation of the de-bonding process on the contact interface where sliding occurs due to compression induced section expansion of the soft body. We investigate how the de-bonding area evolves with the compressive load, and how the evolution is affected by the interface adhesion orthotropy. Characteristics of the system are described in the following. The soft body is 6 mm high with a square section of 10 × 10 mm. It is modelled by Blatz-Ko hyperelastic material with a shear modulus of G = 2.1 × 10 5 MPa. Adhesive interface parameters are: w = 100 J.m -2 , C tx = 1 × 10 11 N.m -3 , C ty = 1 × 10 10 N.m -3 and b = 0.1 N.s.m -1 . Therefore, interface adhesive behaviour is orthotropic, with adhesive stiffness along x direction significantly stronger than that along y direction. We suppose that the system does not involve initial adhesion on the interface (adhesion strength parameter β = 0 at time 0). Knowing that the displacement variation in the z-axis direction is from the boundary to the centre, therefore we observe that the R adh-z are visible in the upper and lower edge regions. As a resistance, the direction of the R adh-z is always opposite to the direction of movement in this area. Moreover, the action of R adh-z is tightened to the left following the propagation of adhesive bond break. Interfacial adhesion fatigue under cyclic loading

Introduction

In reality, the reversible adhesion varies depending on the contact surface, the contact environment and the contact frequency, e.g. the adhesion of a tape degrades with successive sticking and removing. The RCC model used in Chapter 2 presents a perfect reversible adhesion, which means that a perfect adhesion (β = 1) can be achieved with a sufficient contact time, even if the bonding-debonding process is repeated numerous times. For modelling the degradation of adhesive effect, Raous et al. [START_REF] Raous | Recoverable adhesion and friction[END_REF] propose a relationship between the maximum adhesive intensity and the number of debonding. In this chapter, we add this adhesive degradation into our adhesive contact model to investigate the contact problems with partial recoverable adhesion.

In the follows: in Section 5.2, we describe the relation of degradation and reconstruct the differential equation of β. Then, we present numerical examples in Section 5.3. In the end, a few concluding remarks are drawn in Section 5.4.

Recoverability degradation of RCC adhesive model

Based on the RCC adhesive model, Raous et al. [START_REF] Raous | Recoverable adhesion and friction[END_REF] propose an evolution equation to achieve partial recoverability by regulating the maximum value of β (β ∈ [0, β max ]):

where λ represents the degradation factor with λ ∈ [0, 1] and β 0 max ∈ [0, 1] is the initial maximum adhesive intensity. The function of [ ] -is used to retain the negative part in the square brackets, in Eq.5.1, the result in the square brackets is maintained when β is negative, conversely ( β ≥ 0), the result in the square brackets equals zero. Moreover, Eq.5.1 can be developed to the following form:

Therefore, β max is reduced only when adhesive bond breakage ( β < 0) occurs, and the degree of each degradation depends on the product of the bond breakage level (β t -β 0 )

and the degradation factor λ. The differential equation (Eq.2.16) of β can be written as:

(5.3)

Assuming that we set a 5 cycles of bonding-debonding process with β 0 max = 1 and λ = 0.25, which means that the maximum adhesive strength β max decreases by 25% after each cycle. The comparative schematic between complete recoverability and recoverability degradation is shown as Figure 5.1.

Rolling adhesion of a hyperelastic wheel

For this example we use the exact same modelling as the rolling example in Chapter 2

to compare the difference between the perfect recoverability and the partial recoverability. In regard to meeting the requirements of dealing with adhesive friction problems at perfectly recoverable cases with partially recoverable cases.

To the best knowledge of the author, adhesive contact is a highly interesting topic in the medical and biological fields, especially in numerical analysis, where efficient and practical numerical models are very limited, and the commonly used commercial software does not perfectly address all the needs of the researcher. As the adhesion effect occurs in numerous scenarios, a robust, efficient adhesive contact model therefore has a very promising application. In the following we give our perspectives on potential applications based on the proposed recoverable adhesive contact model:

• In the biological and bio-mimetic field, in addition to the biological soft tissues and bio-mimetic materials mentioned in this thesis, adhesive contact models have many other application scenarios, such as, the effect of carcinogenesis on the cells adsorption, the effect of surface adhesion on the use of hydrogels in biomedical area, and the adhesion effect of some new micro structures of bio-mimetic materials, etc.

• Automated tape placement of carbon fibre composites is an important process for large aerospace component manufacturing, which creates adhesion between contact interfaces by heating the resin of the composite strip, pressing and cooling. However, the proposed adhesive contact model in this thesis is insufficient for modelling this scenario. In order to achieve this aim, a multi-physics adhesive contact model, which incorporate temperature-dependent material properties and adhesive parameters, needs to be constructed.

109

Appendix A

Comparison of the bi-potential method with other contact algorithms

To solve the adhesive interface law between hyperelastic bodies, a contact algorithm based on bi-potential theory is used. This algorithm, according to its description of contact kinematics, can be attributed to the category of "node-to-segment" approaches and, with regard to the resolution technique that enforces the contact geometry, belongs to the class of augmented Lagrangian methods. Let us refer to the present contact algorithm with "NTS-AL" (meaning "node-to-segment" contact using augmented Lagrangian resolution), and compare it with other established contact algorithms using alternative schemes of contact kinematics and resolution. In this regard, we consider the widely adopted contact patch test introduced by Taylor and Papodopoulos [START_REF] Taylor | On a patch test for contact problems in two dimensions[END_REF] and compare our results with those reported in [START_REF] Zavarise | A modified node-to-segment algorithm passing the contact patch test[END_REF]. The contact patch test investigates the capacity of a contact algorithm to correctly evaluate the normal contact stresses on contact interface, regardless of its discretization.