Florian Meray

Thibaut Chaise

Anthony Gravouil

Pierre Depouhon

Bruno Descharrieres

Daniel Nélias
email: daniel.nelias@insa-lyon.fr

A novel SAM/X-FEM coupling approach for the simulation of 3D fatigue crack growth under rolling contact loading

Keywords: 3D Fatigue Crack Growth, Rolling Contact Fatigue, Numerical Modeling, eXtented-Finite Element Method, Semi-Analytical Method

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Nowadays, it is well established that Rolling Contact Fatigue (RCF) is the most common failure mode of mechanical components such as rolling element bearings or gears [START_REF] Sadeghi | A review of rolling contact fatigue[END_REF]. Indeed, the complex localized multi-axial and non-proportional stress field, generated by the cyclic contact loading, can locally affect the materials of the contacting bodies and subsequently initiate micro-cracks. As already stated in the literature [START_REF] Tallian | Failure atlas for Hertz contact machine elements[END_REF], two groups of RCF failure mechanisms can be distinguished according to the starting position of these micro-cracks, namely subsurface initiated contact fatigue and surface initiated contact fatigue. The former failure process involves the development of micro-cracks around sub-surface material defects including, among other, inclusions, heterogeneities, voids, etc. Owing to the constant improvement of the cleanliness of materials used in modern machinery, sub-surface damage is now mostly prevented. Consequently, the most observed RCF failure mode in engineering application originates from surface contact fatigue. As depicted by Nélias et al. [START_REF] Nélias | Role of inclusions, surface roughness and operating conditions on rolling contact fatigue[END_REF], this failure mode is charac-terized by the apparition of micro-cracks near surface geometrical imperfections like dents or scratches, which act as stress raisers. In bearing applications, the presence of solid particles transported through the contact zone by the lubricant is often designated as the root cause of these failures [START_REF] Nélias | Detrimental effects of debris dents on rolling contact fatigue[END_REF]. Near-surface plastic deformation, surface contamination, surface roughness, finishing marks, insufficient lubrication and working conditions have also been identified as potential contributors of the surface initiated fatigue.

Over the time, coalescence of the surface originated microcracks may produced a contact-scale crack, usually oriented in the opposite direction to the friction force and inclined at shallow angle from the surface (15-30 • ) [START_REF] Olver | The mechanism of rolling contact fatigue: An update[END_REF]. After a stable propagation stage in shear mode driven by the compressive contact loading [START_REF] Bold | Shear mode crack growth and rolling contact fatigue[END_REF], several crack growth scenarios are possible [START_REF] Dubourg | A predictive rolling contact fatigue crack growth model: Onset of branching, direction, and growth -Role of dry and lubricated conditions on crack patterns[END_REF] (see Fig. 1):

➀ The main crack can branch towards the surface causing the formation of a macro-crater related to the so-called spalling phenomenon, ➁ The main crack can promote a sub-surface propagation leading to the complete failure of the mechanical component.

The latter scenario is usually preceded by the initiation of secondary cracks that branch towards the surface. This promotes the release of particles and facilitates crack detection in industrial applications. It is important to emphasize that, besides stresses induced by the contact, structural stresses and residual stresses due to the hardening and fabrication processes can also have a significant influence on the crack growth, which demonstrates the complexity of the problem.

To enhance the reliability of rolling element bearings and gears, a great interest is being addressed by the aeronautic industry to tackle the complex crack behavior under RCF.

Although some experimental investigations have recently been conducted [START_REF] Le | Rolling contact fatigue crack propagation in nitrided alloyed steels[END_REF][START_REF] Rycerz | Propagation of surface initiated rolling contact fatigue cracks in bearing steel[END_REF], analytical and numerical methods based on the Linear Elastic Fracture Mechanics (LEFM) remain the most widely used tool in the literature for investigating the influence of the numerous parameters involved in the RCF crack growth problem. Keer and Bryant [START_REF] Keer | A pitting model for rolling contact fatigue[END_REF] were amongst the first to propose a 2D numerical model dedicated to the analysis of a surface breaking crack embedded in a half-space and loaded by a simple Hertzian contact. The distributed dislocation technique has been employed by these authors to compute the Stress Intensity Factors (SIFs) at the crack tip. Bower [START_REF] Bower | The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks[END_REF] has extended the model of Keer and Bryant to study the role of lubricant which enters into the crack during the passage of rolling elements. The body force method, described in details in Ref. [START_REF] Murakami | Analysis of surface crack propagation in lubricated rolling contact[END_REF], has been used by Kaneta and Murakami in several papers [START_REF] Kaneta | Mechanism of crack growth in lubricated rolling/sliding contact[END_REF][START_REF] Kaneta | Mechanism of surface crack growth in lubricated rolling/sliding spherical contact[END_REF][START_REF] Kaneta | Propagation of semi-elliptical surface cracks in lubricated rolling/sliding elliptical contacts[END_REF] to investigate the SIF variations for semi-elliptical and semi-circular planar cracks under three-dimensional Hertzian contact pressure ranging from spherical to line configuration. Integral formulations involving Green's functions [START_REF] Fletcher | A simple method of stress intensity factor calculation for inclined surface-breaking cracks with crack face friction under contact loading[END_REF][START_REF] Fletcher | Rapid method of stress intensity factor calculation for semi-elliptical surface breaking cracks under three-dimensional contact loading[END_REF] and Weight functions [START_REF] Beghini | An application of the weight function technique to inclined surface cracks under rolling contact fatigue, assessment and parametric analysis[END_REF] as well as the Boundary Element Method (BEM) [START_REF] Akama | Boundary element analysis of surface initiated rolling contact fatigue cracks in wheel/rail contact systems[END_REF] have also been employed to solve similar problems.

However, the models cited above are limited to idealized crack geometry -edge crack in 2D and semi-circular or semielliptical planar crack in 3D -embedded in an elastic half-space and loaded by a theoretical Hertzian contact pressure. Moreover, residual stresses are generally defined as a constant value that does not represent real in-depth variation. For this purpose, the Finite Element Method (FEM) has been widely used to add complexity in this kind of simulation, by considering, for example, the true geometry of the contacting bodies [START_REF] Bogdanski | Numerical modelling of a 3D rail RCF 'squat'-type crack under operating load[END_REF][START_REF] Bogdanski | Modelling the three-dimensional behaviour of shallow rolling contact fatigue cracks in rails[END_REF][START_REF] Deng | Analysis of surface crack growth under rolling contact fatigue in a linear contact[END_REF][START_REF] Masoudi Nejad | Numerical study on fatigue crack growth in railway wheels under the influence of residual stresses[END_REF], elasticplastic behavior [START_REF] Li | Analysis of fatigue crack propagation in rails under combined action of wheel load and residual stress[END_REF] or complex lubricant behaviors along the crack faces [START_REF] Bogdanski | Quasi-static and dynamic liquid solid interaction in 3D squat-type cracks[END_REF][START_REF] Dallago | The role of lubricating fluid pressurization and entrapment on the path of inclined edge cracks originated under rollingsliding contact fatigue: Numerical analyses vs. experimental evidences[END_REF].

Even though FEM models enhance the capabilities to predict realistic crack behavior, the modeling of 3D long crack growth has not been tackled and most of the studies are usually restricted to one loading cycle. Indeed, despite the recent numerical developments, FEM requires time-consuming and laborious operations of re-meshing and field projection when the crack growths, which limits its usefulness for this kind of analysis. In addition, both specific elements and suitable fine discretization are needed around the crack front to accurately compute the SIFs.

The use of the eXtended-Finite Element Method (X-FEM) is an attractive alternative to alleviate these issues. Initially introduced by the US Northwestern University in 1999 [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], X-FEM is now well recognized in fracture mechanics for the simulation of 3D crack propagation. By using the partition of unity method [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF], it consists in enriching the standard finite element basis with additional functions that capture the discontinuity and the singularity generated by the presence of the crack in the structure. Thanks to this enrichment technique, bulk mesh does not necessarily conform to the crack shape and consequently field projections and re-meshing operations are not needed during the crack propagation. Moreover, X-FEM achieves an accurate mechanical solution even with a relatively coarse mesh. It should be noted that the level-set technique is usually combined with X-FEM to facilitate the representation and the tracking of 3D non-planar crack into the structure [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets -Part I[END_REF][START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets -Part II[END_REF].

Over the past two decades, X-FEM has been widely developed by the scientific community to tackle a large range of engineering issues related to crack propagation. Among others, problems involving contact along the crack faces have been addressed. For instance, Dolbow et al. [START_REF] Dolbow | An extended finite element method for modeling crack growth with frictional contact[END_REF] and Géniaut et al. [START_REF] Geniaut | A stable 3D contact formulation using X-FEM[END_REF] have respectively proposed a 2D frictional contact formulation and a 3D unilateral contact formation for cracks submitted to a static compressive loading. Bi-dimensional quasi-static crack problems with frictional contact have also been treated by Ribeaucourt et al. [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LATIN method[END_REF] and Liu et al. [START_REF] Liu | A contact algorithm for frictional crack propagation with the extended finite element method[END_REF]. Pierres et al. have further developed a two-scale X-FEM strategy dedicated to 3D quasi-static crack problem with unilateral [START_REF] Pierres | A two-scale extended finite element method for modelling 3D crack growth with interfacial contact[END_REF] and frictional contact [START_REF] Pierres | 3D two scale X-FEM crack model with interfacial frictional contact: Application to fretting fatigue[END_REF][START_REF] Pierres | Experimental and numerical analysis of fretting crack formation based on 3D X-FEM frictional contact fatigue crack model[END_REF]. The formulation is based on a three-field weak formulation that allows the crack faces to be discretized independently from the structure, and thus precisely describes the complex sequences of closing, opening, sticking and sliding along the interface. This model has already shown its ability to investigate the crack behavior under fretting [START_REF] Baietto | Fretting fatigue crack growth simulation based on a combined experimental and XFEM strategy[END_REF] and rolling contact [START_REF] Mai | Numerical simulation of rolling contact fatigue crack growth in rails with the rail bending and the frictional contact[END_REF] loadings. Fretting issues have also been addressed by Giner et al. [START_REF] Sabsabi | Experimental fatigue testing of a fretting complete contact and numerical life correlation using X-FEM[END_REF][START_REF] Giner | Fretting fatigue life prediction using the extended finite element method[END_REF] with an implementation of X-FEM in the Abaqus software [START_REF] Giner | An Abaqus implementation of the extended finite element method[END_REF]. The latter considers the contact between the crack faces [START_REF] Giner | Crack face contact in X-FEM using a segment-to-segment approach[END_REF] by using a mortar approach.

Although X-FEM is well-suited for 3D long crack propagation, dealing with tribological loading remains a challenging task. Indeed, such problems require managing different scales (structure scale, contact scale and local interfacial scale) inducing complex and prohibitive meshes. In that respect, multi-scale strategies and coupling algorithms within the X-FEM framework have been developed to reduce computational efforts. For example, Guidault et al. [START_REF] Guidault | A multiscale extended finite element method for crack propagation[END_REF] have proposed a global-local approach using the domain decomposition concept, whereas Rannou [START_REF] Rannou | A local multigrid X-FEM strategy for 3-D crack propagation[END_REF] and Passieux [START_REF] Passieux | Local/global non-intrusive crack propagation simulation using a multigrid X-FEM solver[END_REF] have used the multi-grid technique. More recently, Gibert [START_REF] Gibert | A 3D automatic mesh refinement X-FEM approach for fatigue crack propagation[END_REF] has combined X-FEM with an automatic Adaptive Mesh Refinement (AMR) procedure to address multi-scale issues. A non-intrusive multi-scale approach that couples a 3D Generalized Finite Element Method (GFEM) and Abaqus has been presented by Li et al. [START_REF] Li | Non-intrusive coupling of a 3-D Generalized Finite Element Method and Abaqus for the multiscale analysis of localized defects and structural features[END_REF]. They have used the non-intrusive global-local concept introduced by Allix and coworkers [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF]. However, most of these studies do not consider the non-linear frictional contact arising along the crack faces due to numerical difficulties.

To overcome the current limitations, a novel approach is proposed in this work to efficiently simulate 3D RCF long crack growth. The strategy is based on the coupling of the quasi-static X-FEM formulation proposed by Pierres for frictional contacting crack [START_REF] Pierres | 3D two scale X-FEM crack model with interfacial frictional contact: Application to fretting fatigue[END_REF][START_REF] Pierres | Experimental and numerical analysis of fretting crack formation based on 3D X-FEM frictional contact fatigue crack model[END_REF] with the Semi-Analytical Method (SAM) [START_REF] Nogi | Influence of a Hard Surface Layer on the Limit of Elastic ContactPart I: Analysis Using a Real Surface Model[END_REF][START_REF] Polonsky | A Fast and Accurate Method for Numerical Analysis of Elastic Layered Contacts[END_REF]. SAM is a numerical technique that allows to perform robust and fast modeling of three-dimensional contact problems between two semi-infinite bodies. Derived from boundary integral formulations, it consists in decomposing a complex contact problem into a sum of simple problems for which an analytical closed-form solution is known. A numerical superposition of each elementary solution is then performed in a contact solver algorithm to obtain the solution of the mechanical problem. Owing to the use of FFT acceleration techniques, computational costs are drastically reduced by several orders of magnitude compared to FEM. SAM is also more precise for a given mesh size and the pre-processing stage is straightforward. Through numerous works, SAM has been greatly developed to address a wide range of issues including the resolution of normal and tangential contact between bodies having non-linear behaviors such as elastic-plastic behavior [START_REF] Jacq | Development of a three-dimensional semi-analytical elastic-plastic contact code[END_REF][START_REF] Wang | Numerical Simulation for Three Dimensional Elastic-Plastic Contact with Hardening Behavior[END_REF], thermo-elastic behavior [START_REF] Liu | A three-dimensional thermomechanical model of contact between non-conforming rough surfaces[END_REF], thermo-elastic-plastic behavior [START_REF] Boucly | Contact analyses for bodies with frictional heating and plastic behavior[END_REF], heterogeneous elastic behavior [START_REF] Zhou | A fast method for solving three-dimensional arbitrarily shaped inclusions in a half space[END_REF][START_REF] Zhou | Multiple 3D inhomogeneous inclusions in a half space under contact loading[END_REF][START_REF] Koumi | Contact analysis in the presence of an ellipsoidal inhomogeneity within a half space[END_REF], heterogeneous elastic-plastic behavior [START_REF] Amuzuga | Fully coupled resolution of heterogeneous elasticplastic contact problem[END_REF], heterogeneous viscoelastic behavior [START_REF] Koumi | Rolling contact of a rigid sphere/sliding of a spherical indenter upon a viscoelastic half-space containing an ellipsoidal inhomogeneity[END_REF], elasto-plastic layered behavior [START_REF] Chen | Modeling elasto-plastic indentation on layered materials using the equivalent inclusion method[END_REF], elastic-damageable behavior [START_REF] Beyer | A damage model for fretting contact between a sphere and a half space using semi-analytical method[END_REF] and more recently viscoelastic layered behavior [START_REF] Wallace | Three-dimensional rolling/sliding contact on a viscoelastic layered half-space[END_REF].

The proposed SAM/X-FEM coupling derives from a globallocal strategy [START_REF] Sun | A global-local finite element method suitable for parallel computations[END_REF][START_REF] Ransom | Global/local stress analysis of composite panels[END_REF] similar to the so-called submodeling technique. It allows to benefit from the advantages of SAM and X-FEM in a single model without any alteration of their existing implementation. The method consists firstly in performing a global SAM analysis to solve the contact problem between two semi-infinite uncracked bodies that may contain residual stresses. Then, a data transfer procedure is employed to determine boundary conditions for the subsequent local analysis by extracting and mapping the results from the global SAM analysis on the boundary of the local model. It permits the use of a fully independent discretization scheme for the global and local models while avoiding the introduction of numerical errors. Fi-nally, the local analysis of the 3D crack problem is performed only in the region of interest with the help of the aforementioned X-FEM formulation which accounts for the frictional contact between the crack faces. As a consequence, computational efforts are just focused on the local X-FEM sub-domain and thus CPU times are drastically reduced compared to the existing model. The implementation of the non-intrusive SAM/X-FEM method within a crack propagation algorithm provides a robust and accurate tool, particularly useful for engineering and industrial applications.

This paper is organized as follows. Firstly, the developed SAM/X-FEM coupling is introduced in Section 2. Specifically, the data transfer procedure between the global model and the local model is detailed. Then, the crack growth algorithm dedicated to the simulation of 3D RCF long crack propagation is exposed in Section 3. In Section 4 and Section 5, results concerning the validation of the SAM/X-FEM strategy are shown. Finally, Section 6 highlights the great potential of these developments on a practical example.

The SAM/X-FEM method

In this section, the novel global-local SAM/X-FEM coupling strategy dedicated to the resolution of 3D crack problem under contact loading is introduced. A particular attention is paid to the procedure for the transfer of the mechanical quantities between the models.

Coupling procedure

First of all, the reference problem of 3D crack under contact loading is recalled. Consider an elastic contact problem between two bodies Ω 1 and Ω 2 . The half-space assumption is made here, which means that the contact area must be small compared to the characteristic dimensions of the bodies and to their radii of curvature. The materials under consideration are homogeneous and isotropic. A linear elastic behavior is assumed as well as the small strain and displacement assumptions. A crack, denoted by Γ c , is located within the body Ω 1 . For practical purposes, a restriction of Ω 1 , denoted by Ω L , is introduced. It defines a local area surrounding the crack, called zone of interest. Γ represents the surface boundary of Ω L . n Γ is the outward unit normal to Γ. Due to the compressive multiaxial stresses induced by the contact, complex states of opening and contact may arise together along the crack faces. This reference problem is illustrated in Fig. 2.

As previously mentioned, the resolution of this problem is performed by a global-local approach [START_REF] Sun | A global-local finite element method suitable for parallel computations[END_REF][START_REF] Ransom | Global/local stress analysis of composite panels[END_REF] that can be seen as a submodeling technique. This multi-model method allows to efficiently examine the behavior of large structures that require accurate stress analysis in a critical localized area which exhibits a local feature such as geometric detail (crack, hole, etc.) or localized complex non-linear phenomenon (plasticity, etc.). First, it consists in carrying out a coarse time-saving analysis disregarding the local feature to access the global response of the entire structure. Then, the results of the global model are used to conduct a more refined independent analysis only in

Ω 1 Ω 2 Γ c Γ Ω L n Γ Fig. 2
Reference problem of 3D crack under contact loading.

the critical sub-domain. It gives accurate informations related to the presence of the local detail while saving computational time and resources compared to a fully refined analysis over the whole structure.

Based on this concept, a global-local SAM/X-FEM coupling procedure is proposed to solve the problem under consideration. In details, the developed solving scheme is decomposed into three main steps (see Fig. 3): ➀ Global Analysis. The resolution of the 3D contact problem between the semi-infinite bodies Ω 1 and Ω 2 is performed by means of the in-house SAM solver without considering the crack (see Appendix A). The solution of the global analysis is denoted by:

s G = u G , σ G , p, q (1)
where u G and σ G correspond to the displacement and stress fields in Ω 1 and Ω 2 , p and q are the pressure and shear distributions acting on the contact surface. Note that an initial residual stress field can be taken into account by means of Eq. (A.12).

➁ Data transfer. After defining the position of the boundaries of Ω L , the solution of the global analysis s G is used to define boundary conditions for the subsequent local analysis. More details about this step are exposed in Section 2.2.

➂ Local Analysis. The boundary value problem of crack is solved by using the implementation in the Cast3M software [START_REF]Cast3M, Finite element software developed by the french atomic energy center (cea)[END_REF] of the quasi-static X-FEM method accounting for contact and friction along the crack faces (see Appendix B). This analysis is performed only in the restriction Ω L of Ω 1 where the crack is located. The solution of the local analysis is denoted by:

s L = u L , σ L , w, t (2)
where u L and σ L correspond to the displacement and stress fields in Ω L , w and t are the interface displacement and traction fields on the crack faces Γ c .

Hence, the solution provided by the SAM/X-FEM procedure can be summarized as the substitution of the solution of the global analysis by the solution of the local analysis in the zone of interest Ω L :

s =        s G in (Ω 1 ∪ Ω 2 ) \Ω L s L in Ω L (3)
This coupling approach is non-intrusive because the SAM and X-FEM codes do not require any modification. Thus, the in-house SAM solver can be easily coupled to other industrial computational softwares by means of this procedure.

Moreover, it should be emphasized that the proposed coupling relies on a unidirectional data transfer from the global SAM model to the local X-FEM model since there is no iterative procedure between the 2 models. Therefore, the influence of the local X-FEM model on the global SAM model is neglected here which means that the presence of crack does not disturb the contact pressure and shear distribution. This assumption is discussed by Giner et al. [START_REF] Giner | Numerical modelling of crackcontact interaction in 2D incomplete fretting contacts using X-FEM[END_REF] for crack propagation under fretting contact conditions.

As the crack is disregarded in the global analysis, the oneway SAM/X-FEM coupling procedure implies that the mechanical field disturbance caused by the crack must vanish on the coupling interface Γ. Consequently, the size of the local model should be carefully chosen to ensure that the location of Γ does not affect the solution of the local analysis. For cracking issues, the local sub-domain must enclose at least the region affected by the presence of the crack, which is linked to the crack geometry and loading conditions. Due to the difficulties in quantifying this area in such a 3D complex problem, the user has to investigate the sensitivity of the local solution with respect to the location of the coupling interface Γ in order to extract the appropriate minimum size of the sub-model. In practice, the size of the local model in each direction (see Fig. 4) is defined with the ratios α i given by:

α 1 = a a + d 1 ( 4
)
α 2 = b b + d 2 (5)
α i represents the proportion of cracked surface within the local model. Therefore, the smaller α i is, the larger the size of the local model is. In the first instance, for a crack size approaching the contact dimensions, it is recommended to take α i equal to 0.5. Note that these ratios can be increased when the crack becomes significantly larger than the contact area. A sensitivity analysis about the position of the coupling interface Γ is carried out in this work (see Section 6).

Besides the position of the coupling interface, the efficiency of the proposed numerical scheme also depends on the accuracy of the boundary conditions prescribed on the local model. This point is discussed hereafter.

Ω 1 s G = (u G ,σ G ,p,q) 1 Global Contact Problem 3 Local Crack Problem 2 Data Transfer Ω 2 s G Γ Γ Γ s L = (u L ,σ L ,w,t) Ω L Γ c

Global-local data transfer

The data transfer step is a crucial issue of the proposed SAM/X-FEM coupling. As previously explained, this step provides boundary conditions for the local X-FEM analysis using the results of the global SAM analysis. The mechanical quantities usually transferred across the coupling interface Γ may be the displacement field and/or the stress field. However, as already stated, applying displacements as boundary conditions can induce inaccuracies in the solution when the stiffness of the sub-domain changes between the global and local analyses. For this reason, only the stresses are transferred here because crack causes a stiffness change in the local model. More precisely, the equilibrium condition must be satisfied along the coupling interface Γ:

σ L • n Γ + σ G • (-n Γ ) = 0 ⇔ t L + t G = 0 on Γ (6)
In a discrete point of view, the enforcement of this condition is straightforward if the global model and the local model share the same stress calculation points along the interface Γ. However, in the SAM/X-FEM framework, such a situation cannot be fulfilled since the methods are founded on completely different formalisms. Moreover, in order to accurately capture the local behavior in the zone of interest, the discretization of the local model needs to be freely refined, regardless of the discretization of the global model.

Therefore, to deal with non-conforming and independent discretization scheme for the global and local models, a robust interpolation procedure (see Fig. 5) is developed to access t L (defined on the local model) from t G (defined on the global model). In what follows, E Γ are the 2D elements of the surface mesh M Γ of the local model boundary Γ. Moreover, P Γ i defines the i-th integration point of E Γ and x Γ i is its location. 1. First, a virtual 3D hexahedral P1 finite element is built using the 8 SAM calculation points surrounding P Γ i . 2. By using the values of global SAM stresses at the nodes of the virtual element, a linear interpolation is performed to evaluate the stress tensor at P Γ i :

σ L (x Γ i ) = 8 j=1 N h, j (x Γ i )σ G j ( 7
)
where N h, j are the 3D linear interpolation functions of a standard hexahedral P1 finite element. 3. Then, the surface traction vector t L is evaluated at P Γ i as follows :

t L (x Γ i ) = σ L (x Γ i ) • n Γ (x Γ i ) (8)
As it permits the use of completely different discretization for the two models, this interpolation procedure makes the preprocessing easier and avoids complex meshing operations at the coupling interface.

To enforce the equilibrium (6) on Γ within a X-FEM framework, the traction vector t L must be converted into nodal forces. So the elementary nodal force vectors are computed by integrating t L over each surface element E Γ :

F e = E Γ N T e t L dS (9) Global SAM Discretization Local X-FEM Discretization SAM stress calculation point Virtual hexahedral P1 finite element σ 1 G σ 2 G σ 3 G σ 4 G σ 5 G σ 6 G σ 7 G σ 8 G Fig. 5
Schematic overview of the stress interpolation process between the global SAM model and the local X-FEM model.

where N e is the matrix of 2D linear shape functions of the element E Γ . This integration is numerically realized with a Gaussian quadrature rule. As a last step, the contributions of each elementary nodal force vector are summed to access the nodal force vector F that satisfies (6) over the entire surface boundary Γ of the local X-FEM model:

F = F e (10)
It may be added that, in a classical global-local approach, a particular attention is paid to the convergence of the local problem solution in relation to the quality of the imposed boundary conditions. Indeed, a too coarse discretization of the global model can induce inaccuracies that affect the local model through the boundary conditions [START_REF] Cormier | Aggressive submodelling of stress concentrations[END_REF]. However, in the SAM/X-FEM context, the results of the global SAM model are derived from known analytical solutions. As a consequence, the response of the local model is not noticeably sensitive to the discretization of the global model.

Nodal force vector balancing

In the proposed approach, only the force boundary condition F is used as input for the local analysis. It is therefore a pure force boundary value problem. The local model can be seen as a floating sub-structure as there are no Dirichlet boundary conditions to constrain its movements. From a mathematical point of view, the resolution of such a problem is possible only if the prescribed nodal force vector is orthogonal to the kernel of the stiffness matrix K of the local model:

F ∈ Im(K) = Ker(K) ⊥ (11)
This condition can also be expressed as follows:

R T F = 0 ( 12
)
where R = Ker(K) is the null space of K and represents the rigid-body modes of the local model. As indicated by Eq. ( 12), solving a pure force boundary value problem requires that the projection of the nodal force vector F onto the rigid body mode space R T is zero. Physically, it means that the applied nodal forces must ensure the self-equilibrium of the floating local model and so prevent rigid body displacements. In other words, the sum of forces and moments applied over Γ must be zero.

However, the interpolation and integration steps depicted in Section 2.2 can induce a residual unbalanced of the nodal force vector, which raises difficulties to reach convergence of the local X-FEM analysis. To tackle this difficulty, a projector denoted P is introduced. The multiplication of this projector P by the nodal force vector F yields a self-balanced nodal force vector F p satisfying the orthogonality condition [START_REF] Murakami | Analysis of surface crack propagation in lubricated rolling contact[END_REF]:

F p = PF =⇒ R T F p = 0 (13)
with:

P = I -R(R T R) -1 R T (14)
This projection eliminates from the initial nodal forces F the components causing the residual unbalance, namely the part of F that is not orthogonal to the rigid body mode space. It is worth to point out that similar projection steps are generally implemented within dual Domain Decomposition solvers such as FETI algorithm [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Hirschler | A dual domain decomposition algorithm for the analysis of non-conforming isogeometric Kirchhoff-Love shells[END_REF], which also handles floating substructures. F p is then used as input for the local X-FEM analysis.

3D crack propagation under RCF

This section deals with the computation of the SIFs as well as the implementation of the SAM/X-FEM coupling within an automatic algorithm for RCF crack growth simulation.

Stress intensity factor calculation

After solving the mechanical problem with the SAM/X-FEM approach, the SIFs along the crack front are extracted separately by means of the interaction integral concept. Firstly introduced in Refs. [START_REF] Stern | A contour integral computation of mixed mode stress intensity factors[END_REF][START_REF] Yau | A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity[END_REF], this method is an extension of the so-called J-integral. It relies on a superimposition of two distinct equilibrium states of an elastic cracked body: the actual state related to the solution of the local X-FEM analysis and the auxiliary state associated to the solution of an arbitrary crack problem. The interaction integral is generally expressed in a domain form [START_REF] Nakamura | Antisymmetrical 3-D stress field near the crack front of a thin elastic plate[END_REF][START_REF] Gosz | Domain integral formulation for stress intensity factor computation along curved three-dimensional interface cracks[END_REF][START_REF] Gosz | An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions[END_REF], which is more appropriate for 3D finite element computation. In the present work, the domain interaction integral is evaluated at each discrete point P j of the crack front as follows [START_REF] Pierres | Experimental and numerical analysis of fretting crack formation based on 3D X-FEM frictional contact fatigue crack model[END_REF]:

I (L,aux) = D σ L k j u aux k,i + σ aux k j u L k,i -σ L kl ε aux kl δ i j q i, j dV - Γ + c ∪Γ - c t L k u aux k,i q i dS (15)
Note that the second integral over Γ + c ∪ Γ - c in Eq. ( 15) accounts for the frictional contact arising along the crack faces [START_REF] Dolbow | An extended finite element method for modeling crack growth with frictional contact[END_REF][START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LATIN method[END_REF]. Here, the integration domain D is supported by the X-FEM mesh. More precisely, it is built by adding a user-defined number of element layers N D around the element containing the point P j [START_REF] Prabel | Propagation de fissures tridimensionnelles dans des matériaux inélastiques avec xfem dans cast3m[END_REF] (see Fig. 6). The vector field q represents the virtual crack extension and is tangent to the crack plane and normal to the crack front. The magnitude of q is defined by a sufficiently smooth scalar function, denoted by α(x), which is equal to 1 at P j and vanishes on the boundaries of D. The implementation of the domain interaction integral within the X-FEM framework is described in Ref. [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets -Part I[END_REF].

By choosing the auxiliary fields wisely, the SIFs in the actual state can be expressed as a function of the domain interaction integral. For instance, by introducing the pure mode I Westergaards asymptotic fields into Eq. ( 15), the mode I SIF K I at P j can be evaluated with the following equation:

K j I = E 2(1 -ν 2 )
I (L,mode I) C α(s)ds [START_REF] Fletcher | A simple method of stress intensity factor calculation for inclined surface-breaking cracks with crack face friction under contact loading[END_REF] where C corresponds to the crack front segment included in D. The plane strain assumption is assumed here. In a similar way, K II and K III at P j can be calculated as:

K j II = E 2(1 -ν 2 ) I (L,mode II) C α(s)ds (17) K j III = E 2(1 + ν) I (L,mode III) C α(s)ds (18)

Crack growth algorithm

To automatically handle the numerical propagation of RCF crack, the SAM/X-FEM coupling is integrated into an algorithm discussed here (see Fig. 7).

First of all, the problem is initialized by setting the contact and loading conditions, the materials, the residual stresses, the initial crack geometry, the number of element layers N D for the calculation of the SIFs, and also α i for the definition of the adequate size of the local model. Furthermore, as the position of the global model boundaries does not affect the results (see Appendix A), the size of this model is set to fit closely with the local model dimensions and the loading range.

Generally, for RCF issues, crack propagation is studied incrementally at the cycle scale because this is a relatively slow temporal phenomenon. In that respect, the current loading cycle N i is discretized into several time steps describing the passage of the rolling elements over the crack (see Fig. 8). The resolution of the mechanical problem is performed for each loading position by means of the quasi-static SAM/X-FEM coupling approach (see Section 2). From the results of this computation, the domain interaction integral is employed to evaluate the SIFs (K j I ,K j II ,K j III ) at each point P j of the crack front (see Section 3.1).

At the end of the simulated cycle, the SIF history is used together with criteria to predict the crack growth direction θ j along the crack front. More precisely, the Maximum Shear Stress Range (MSSR) criterion is adopted here. It assumes that crack propagates in shear mode along the direction which maximizes the range of the equivalent shear SIF:

max θ ∆K τ (θ) = max θ max t K τ (θ, t) -min t K τ (θ, t) (19)
with:

K τ (θ, t) = 1 2 cos θ 2 [K I (t) sin(θ) + K II (t)(3 cos(θ) -1)] (20)
Alternative criteria can be found in the literature, such as the Houlier and Pineau's criteria [START_REF] Dubourg | A predictive rolling contact fatigue crack growth model: Onset of branching, direction, and growth -Role of dry and lubricated conditions on crack patterns[END_REF][START_REF] Hourlier | Propagation of fatigue cracks under polymodal loading[END_REF] for instance. However, the analysis of these criteria is out of the scope of this paper. It is important to emphasize that the contribution of mode III is

N i = N i + ΔN i s G = (u G ,σ G ,p,q)
Resolution of the Local Crack Problem (X-FEM)

s L = (u L ,σ L ,w,t)
Computation of the LEFM quantities

K I , K II , K III j j j
Determination of the growth rate and direction θ j , (da / dN) j

Geometric Update

End of Simulation

End of Loading Cycle Yes Fig. 7 Algorithm for the modeling of 3D crack propagation under RCF.

Cycle scale

Loading cycle i

Loading cycle i+1

N i+1 = N i + ΔN i ΔN i N i Fig. 8
Temporal discretization for fatigue crack growth simulation.

not included in most of these criteria. Indeed, the impact of this mode on the crack growth direction is poorly studied in the literature, especially for RCF issues.

The crack growth rate (da/dN) j is also estimated along the crack front with the help of a Paris-type law:

da dN j = C ∆K j eq m ( 21
)
where C and m are material constants and ∆K eq is a combination of the SIF ranges based on the energy release rate:

∆K eq = ∆K 2 I + ∆K 2 II + 1 1 -ν ∆K 2 III ( 22
)
Due to the lack of data for multi-axial and non-proportional loading, the Paris law coefficients C and m determined for mode I are used. This is questionable because mode I kinetics and mixed mode kinetics (I+II+III in sense of G) are probably different.

By assuming that the crack growth rate is constant over a block of cycles ∆N i , the crack front can be updated by fixing the maximum crack advance ∆a max at the node where the crack growth rate is maximum [START_REF] Shi | Abaqus implementation of extended finite element method using a level set representation for three-dimensional fatigue crack growth and life predictions[END_REF][START_REF] Ren | Three dimensional crack propagation through mesh-based explicit representation for arbitrarily shaped cracks using the extended finite element method[END_REF]. The corresponding cycle jump ∆N i is estimated from the Paris-type law:

∆N i = ∆a max C(∆K max eq ) m (23)
The crack advances ∆a j at the other nodes are determined by keeping constant the calculated cycle jump ∆N i along the crack front:

∆a j =        ∆K j eq ∆K max eq        m ∆a max ( 24
)
Note that the user-defined parameter ∆a max must be carefully chosen in order to accurately simulate propagation while limiting computational time and power (see Section 6). Finally, as explained above, the crack geometry and the level-sets are updated by adding triangular elements between the former and the updated crack front (see Fig. 9). The size of the local model is also updated to ensure that the boundaries have no influence on the solution of the local analysis when the crack propagates, and then the global model dimensions are adapted to fit with the updated local model. This procedure is then repeated for the next loading cycle at N i + ∆N i (see Fig. 8) with the updated meshes and so on. The algorithm is stopped when the crack propagation is no longer influenced by the contact loading.

Validation

In this section, the validation of the proposed SAM/X-FEM coupling strategy is reported through comparison with results extracted from the literature. The variations of the SIFs published by Kaneta et al. in Ref. [START_REF] Kaneta | Mechanism of surface crack growth in lubricated rolling/sliding spherical contact[END_REF] are taken as references.

Crack Front at N i + ΔN i Crack Front at N i θ j e 1 e 2 e 3
Δa j P j Fig. 9 Update of the crack geometry by adding triangular elements between the front at N i and the front at N i + ∆N i

Description of the numerical model

The configuration studied here consists in a semi-circular crack located at the surface of a semi-infinite body. This crack is loaded by a moving compressive loading arising from the passage of a rolling ball on the half-space (see Fig. 10). The material properties are chosen as : Young's Modulus E = 210 GPa and Poisson's ratio ν = 0.3. In what follows, p 0 and c are the maximum contact pressure and the radius of the contact area respectively. The friction coefficient between the contacting bodies is denoted µ contact . As full sliding conditions are assumed in this work, the shear distribution q is simply given by q = µ contact • p. Friction occurs when a power is transmitted through the contact interface. The shear distribution acts to decelerate the faster surface (driver surface) and accelerate the slower surface (follower surface). When the crack is located on the driver surface, the direction of the shear distribution is opposite to the rolling direction. Thus, the direction of the shear distribution and the direction of the contact motion are the same. As the contact moves from the negative x (x < 0) to positive x (x > 0), the shear distribution is therefore positive, i.e. µ contact > 0. Conversely, when the crack is located on the follower surface, the direction of the shear distribution and the rolling direction are the same. Thus, the direction of the shear distribution is opposite to the direction of the contact motion and µ contact < 0.

a corresponds to the crack radius and β represents its inclination with respect to the half-space surface. In this section, the latter is equal to 45 • . γ indicates the position along the semi-circular crack front (0 • ≤ γ ≤ 180 • ). As a consequence, γ = 90 • coincides with the deepest point of the crack front. The parameter x c is used to define the distance between the crack mouth and the contact center.

Five different cases are investigated here. They are summarized in Tab. 1.

For all cases, the local model is a cuboid domain whose dimensions are defined using α 1 = α 2 = 0.25. The dimensions are chosen to be relatively large here to avoid inaccuracy related to the position of the coupling interface. The discretization of the local model is realized by using eight-nodes X-FEM elements with 64 integration points. In the more refined zone where the crack is located, the characteristic element size is equal to ∆x L = ∆y L = ∆z L = a/10, which is fine enough to accurately capture the solution of the local problem. This leads to a X-FEM mesh of 59489 nodes. The triangular elements of the crack mesh have the same size as the bulk discretization, namely ∆x crack = a/10 (210 nodes). N D = 2 is considered for the calculation of the SIFs.

Ranging from x c /c = -3 to x c /c = +3, the rolling contact cycle is decomposed into 61 time steps. Therefore, the imposed displacement increment along the x-direction at each time step is a/10. Furthermore, the global SAM model is discretized in 139 × 85 × 40 computation points spaced by ∆x G = ∆y G = ∆z G = c/10, which is sufficient to cover the local model and the loading range (see Fig. 11).

It is important to emphasize that the contact and friction between the crack faces are not considered here in order to comply with the results provided by Kaneta et al. for these cases. Note that the reference SIFs have been obtained with the body force method. Details about the method can be found in the paper of Murakami et al. [START_REF] Murakami | Analysis of surface crack propagation in lubricated rolling contact[END_REF].

In the following section, the dimensionless SIFs are used for the comparison. They are given as a function of the nominal SIFs:

K i = F i p 0 √ πa (25)
with i = [I, II].

Results and Discussion

Fig. 12(a) and (b) exhibit the variations of F I and F II at γ = 90 • during the passage of the contact over the crack for cases #1 and #2 (crack radius a = c). At first sight, correct agreement can be observed between the models in terms of F I and F II variations, regardless of the contact friction coefficient. In order to quantify more precisely the differences between the two models, the relative errors based on the dimensionless SIF range are computed as follows:

ε i = ∆F i -∆F re f i ∆F re f i ( 26
)
with i = [I, II].
As can be seen in Tab. 2, the errors for cases #1 and #2 are small ranging from -5.79% to 4.21%. Note that the differences are slightly more important for F I than for F II . The discrepancies can be attributed to the formalism of the body force method used by Kaneta that is extremely different from that of the domain integral method employed in the present work (see Section 3). Despite this, the comparison is well acceptable, especially for this kind of simulation.

It is also interesting to point out that negative F I are observed in Fig. 12(a) because the contact and friction along the crack faces are not considered in this study.

(a) (b)

Fig. 12 Variation of the dimensionless SIFs at γ = 90 • during the loading cycle for cases #1 and #2 (crack radius a = c). The dashed curves refer to the results extracted from the Kaneta's reference work [START_REF] Kaneta | Mechanism of surface crack growth in lubricated rolling/sliding spherical contact[END_REF] while the solid curves represent the results obtained with the SAM/X-FEM strategy. The black and blue curves correspond to µ contact = -0.1 (case #1) and µ contact = 0.1 (case #2) respectively. (a) mode I dimensionless SIF and (b) mode II dimensionless SIF.

As above, Fig. 13(a) and (b) show the variations of F I and F II at γ = 90 • during the passage of the contact over the crack for cases #3, #4 and #5 (crack radius a = c/2). Similar evolutions of F I and F II during the passage of the rolling contact are observed between the reference and the proposed model for all the tested friction coefficients. The relative errors for these three last cases range from -7.32% to 1.63% (see Tab. 2), which demonstrates again the excellent agreement between the two models.

All these observations allow to state that the SAM/X-FEM coupling is successfully validated. Furthermore, the present SAM/X-FEM calculations take only about 32 minutes on a laptop with the following characteristics: 8 × Intel CoreTM i7-Tab. 2 Differences in terms of dimensionless SIF range between the reference [START_REF] Kaneta | Mechanism of surface crack growth in lubricated rolling/sliding spherical contact[END_REF] and the SAM/X-FEM model. 7820HQ CPU @ 2.90GHz -31Go RAM. Although the resolution of the contact between the crack faces is not performed, the numerical cost remains relatively low, which proves the great potential of the proposed strategy to simulate in a fast and robust way the problem of crack under moving contact.

Case

Comparison with FEM

To demonstrate the efficiency of these developments, the SAM/X-FEM coupling is now compared with the standard finite element method (FEM). Contrary to the previous section, the interfacial contact between the crack faces is considered here.

Case under consideration

As in Section 4, a rolling contact between a ball and a half space is considered here (see Fig. 10). The ball, which has a radius r = 10 mm, is submitted to a load W = 123 N. Under these loading conditions, the maximum pressure is p 0 = 1468.5 MPa and the contact radius c is equal to 0.2 mm. The two contacting bodies are made of a 16NCD13 bearing steel whose Young's Modulus and Poisson's ratio are E = 210 GPa and ν = 0.3 respectively.

The initial surface breaking crack has a semi-circular shape of radius a = 0.2 mm. Its inclination with respect to the contact surface β is chosen to be equal to 15 • in order to comply with classical experimental observations. γ designates the position along the initial semi-circular crack front (0

• ≤ γ ≤ 180 • ).
The distance between the crack mouth and the contact center is denoted by the parameter x c (see Fig. 10).

A rigid-body displacement along the x direction is applied to the ball to represent the rolling cycle. 61 displacement increments of a/10 in the range x c /c = [-3; 3] are used to describe this moving loading. Frictional effects related to this tangential displacement are not regarded here, thus µ contact = 0.

It is important to emphasize that, in contrast to Section 4, the interfacial contact between the crack faces is taken into account. It allows to faithfully represent the complex states arising along the crack interfaces due to the passage of the compressive contact loading. The interfacial friction coefficient µ crack is equal to zero for this comparison (µ crack = 0).

SAM/X-FEM model

The SAM/X-FEM model used for this comparison is now described. As suggested in Section 2, the local model dimensions are defined using α i = 0.5. To evaluate the relevance of this choice, a convergence study is performed at the deepest point of the crack front (γ = 90 • ). More precisely, the evolution of the relative error based on ∆K II is plotted against α i for several µ crack (see Fig. 14). This error is given by:

ε II = ∆K II -∆K re f II ∆K re f II ( 27
)
where ∆K re f II is the reference value obtained with the parameter α i = 0.25 applied in the validation part.

As expected, it can be observed that the decrease of α i , i.e. the enlargement of the local model, leads to the convergence of ∆K II towards the reference value. Furthermore, for a given α i , the relative error ε II decreases when µ crack increases. Indeed, the augmentation of µ crack significantly reduces the sliding between the crack faces and the stress disturbance generated by the crack within the local model. Here, α i = 0.5 allows to benefit from a good compromise between accuracy and computation efforts with less than 5% of error on the ∆K II value for µ crack = 0. Note that similar results can be observed along the crack front for mode I and III SIFs.

The local model is discretized through the use of eight-node hexahedral X-FEM elements with 64 integration points. To ensure the convergence of the local analysis, the prescribed size of the X-FEM elements is in the same proportions as the one employed in Section 4, namely ∆x L = ∆y L = ∆z L = a/10 = 0.02 mm. This results in a X-FEM mesh of 16233 nodes. The crack interface is discretized with flat triangular elements of size ∆x crack = 0.02 mm (210 nodes). The computation of the SIFs is performed with N D = 2.

Finally, the global SAM model is decomposed into 103 × 43 × 15 computation points spaced by ∆x G = ∆y G = ∆z G = c/10 = 0.02 mm (see Fig. 15).

Local X-FEM Model

Global SAM Model Fig. 15 SAM/X-FEM model used for the comparison with FEM.

FEM model

The finite element model used as reference for this comparison is presented. As shown in Fig. 16, only the body con-taining the crack is modelled. Indeed, contrary to the SAM/X-FEM analysis, the contact problem between the two bodies is not solved in this FE calculation. Thus, a pressure field is applied on the free surface of the body to represent the contact. More precisely, a pressure distribution which represents the sphere/plane contact is prescribed as stress boundary conditions:

p(x, y) = p 0 1 - (x 2 + y 2 ) c 2 (28)
The rolling cycle is represented by a displacement of this pressure distribution along the x direction. In terms of displacement boundary conditions, the bottom face of the body is fully constrained.

Crack

Fig. 16 FEM model used for the comparison with the SAM/X-FEM method.

As the body is considered as a half-space, the edges of the finite element domain must be placed far enough away from the contact zone to avoid the introduction of errors in the solution. In practice, the dimensions of the domain should be at least 20 times larger than the characteristic dimension of the contact zone. Thus, L x = L y = L z are equal to 20 × 2c.

Then, remeshing techniques implemented in Z-cracks [START_REF]Z-cracks, simulation of 3d crack propagation[END_REF][START_REF] Chiaruttini | An adaptive algorithm for cohesive zone model and arbitrary crack propagation[END_REF][START_REF] Chiaruttini | Advanced remeshing techniques for complex 3d crack propagation[END_REF] are used to explicitly introduce the crack geometry into the mesh of the body. A structured mesh around the crack front is constructed to avoid oscillations in the spatial evolution of the SIFs. The size of the elements in this area is equal to 0.002 mm (10 times smaller than the size of the X-FEM elements in the SAM/X-FEM model) in order to capture the crack tip singularity with standard finite elements. To minimize the number of dof, the mesh of the structure is coarser away from the zone of interest containing the crack. This FE mesh is composed of 108250 nodes. The SIFs are extracted using a domain energy integral technique similar to the one employed in this work. Finally, the mechanical problem is solved by means of the SAMCEF industrial finite element solver.

Results and Discussion

In order to compare the two models, the variations of SIFs in relation with the position of the contact center are plotted at γ = 90 • (see Fig. 17). The results obtained with the SAM/X-FEM model are similar to those obtained with the FEM model. For exemple, the maximum value of ∆K II is 5.671 MPa.m 0.5 for the FEM model and 5.722 MPa.m 0.5 for the SAM/X-FEM model, i.e. a difference of 0.90%. The minimal value of ∆K II is 2.574 MPa.m 0.5 for the FEM model and 2.558 MPa.m 0.5 for the SAM/X-FEM model, i.e. a difference of 0.62%. Moreover, the K III do not exhibit oscillations and are equal to zero because γ = 90 • is the symmetry point of the case under consideration.

Then, the comparison is extended to the entire crack front. As the crack is subjected to the three cracking modes simultaneously, the combination of the SIF ranges based on the energy release rate ∆K eq [START_REF] Deng | Analysis of surface crack growth under rolling contact fatigue in a linear contact[END_REF] is computed for each point of the crack front. As shown Fig. 18, the ∆K eq obtained with the two models are comparable, although notable differences are observed at the two vertices of the crack front.

Tab. 3 gives the relative errors in terms of ∆K eq between the FEM model and the SAM/X-FEM model at γ = 45 • , 90 • and 90 • . These errors range from 0.14% at γ = 90 • to 5.03% at γ = 135 • .

Tab. 3

Differences in terms of ∆K eq between the FEM model and the SAM/X-FEM model.

ε γ=45 • (%) ε γ=90 • (%) ε γ=135 • (%)
4.90 0.14 5.03

Thus, these results demonstrate the validity of the developments when the contact between the crack faces is considered. In addition, the advantages in terms of calcultation time can be highlighted (see Tab. 4). It is important to recall that the contact between the sphere and the half-space is not resolved within the FEM analysis, whereas it is tackled with the SAM/X-FEM approach. Despite this, the use of the SAM/X-FEM coupling leads to a reduction of 86% in computation time. The calculations have been carried out with a laptop having the same characteristics as the one used in Section 4.

(b) (c) (a)

Fig. 17 Variation of ∆K eq along the crack front. The solid curve refers to the results extracted from the SAM/X-FEM model while the dashed curve represents the results obtained with the FEM model.

Tab. 4

Comparison of the calculation times of the SAM/X-FEM model and the FEM model.

Model

Calculation Time SAM/X-FEM 2170 s (∼ 36 min) FEM 15120 s (∼ 4 h 12 min)

A practical example

A numerical example demonstrating the efficiency of the development to simulate long crack growth under rolling contact condition is presented in this section. After describing the SAM/X-FEM model, results concerning the first loading cycle are exposed as well as results about the propagation of the crack over several cycles. The influence of residual stresses on crack behavior is investigated.

Description of the numerical model

For this practical example, the case studied is strictly identical to the one investigated in Section 5, i.e. a semi-circular crack inclined at β = 15 • under a rolling contact between a ball and a half space. The parameters of this case are detailed in Section 5.1. The only difference is that the friction between the crack faces is now taken into account. Indeed, a moderate interfacial friction coefficient crack = 0.1 assumed here to simulate the oil lubrication of the crack faces.

As the case under is same, the parameters of the SAM/X-FEM model are identical to those given in Section 5.2.

Concerning propagation aspects, the crack growth direction and rate are predicted with both the MSSR criterion [START_REF] Akama | Boundary element analysis of surface initiated rolling contact fatigue cracks in wheel/rail contact systems[END_REF] and the Paris-type law [START_REF] Bogdanski | Modelling the three-dimensional behaviour of shallow rolling contact fatigue cracks in rails[END_REF]. Due to the lack of data for multi-axial and non-proportional loading, the Paris law coefficients of the 16NCD13 steel established by Oni et al. [START_REF] Oni | A comparison of near-threshold fatigue crack propagation in two high strength steels[END_REF] with mode I fatigue tests are used: C = 1 • 10 -11 and m = 2.97.

In this analysis, a non-constant maximum crack advance is adopted. More precisely, it increases linearly with the maximum crack front depth: ∆a max is equal to ∆x L = 0.02 mm at the beginning of the simulation, and takes the value of 3 • ∆x L = 0.06 mm when the crack front reaches a sufficient depth corresponding to three times the depth of the maximum contact shear stress 3•z * (with z * ≈ 0.49c for ball/plane contact). As the stress gradient generated by the contact becomes small at this depth, this choice achieves an acceptable compromise between accuracy and computational efforts.

As described in Section 3, the local and global models dimensions are automatically updated to comply with the crack extension.

Finally, two cases are investigated: #1 corresponds to the case without residual stress, whereas the case #2 corresponds to the case with residual compressive stresses. For the latter case, only the components σ i xx and σ i yy are non-zero values. It is assumed that they decrease linearly with the depth and vanish at z = 1 mm (see Fig. 19).

Fig. 19

In-depth profile of the residual stresses σ i xx and σ i yy .

Analysis of the first loading cycle

Fig. 20 depicts the SIF evolutions along the crack front during the first loading cycle for case #1.

First, by looking at Fig. 20(a), it can be seen that K I does not exhibit negative value since the contact between the crack faces is considered here. Furthermore, due to high compressive stresses arising from the contact loading, the value of K I remains relatively low and does not exceed K max I = 0.70 MPa.m 0.5 . When the contact loading approaches the crack (x/c = [-3; -1]), a stable phase of slight increase is observed. It is followed by a sudden variation of K I in the interval x/c = [-1; -0.2]. As already stated by Ribeaucourt in Ref. [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LATIN method[END_REF], this phenomenon is attributed to the spread of a "bubble" zone towards the crack front, which enhances the opening of the crack. This partial opening zone appears at x/c = -1 when the leading edge of the contact reaches the crack mouth. For clarity, this particular opening mechanism is schematically represented in Fig. 21. To support this illustration, the crack opening displacement [w N ] is plotted in Fig. 22 for different positions of the contact center. Afterward, when the contact center passes over the crack (x/c >= 0), the bubble zone disappears and the crack remains closed (K I = 0 MPa.m 0.5 ) on a wide portion of the crack front. However, an increase of K I that initiates at x/c = 0.7 is captured around the two vertices of the crack. It progressively decreases and spreads to the bottom of the crack when the contact moves away. This change is induced by a traction zone produced by the axial contact stress σ xx which takes positive values around the contact area (see Fig. 23). Concerning mode II (see Fig. 20(b)), it can be noticed that K II is almost zero when the contact approaches the crack mouth (x/c = [-3; -1]). Then, K II exhibits a reversal behavior for x/c >= -1 with a sign change at x/c = 1.8. So the direction of the relative sliding displacement along the crack faces is reversed when the contact area moves over the crack. This behavior is caused by a combination of the sub-surface shear stresses generated by the contact loading, even though the orthogonal shear stress σ xz is the major contributor. It is worth noting that σ xz follows a reverse trend similar to the K II variation (see Fig. 23). The extreme values K min II = -2.28 MPa.m 0.5 and K max II = 3.85 MPa.m 0.5 are obtained at γ = 90 • for x/c = 0 and x/c = 1.8 respectively.

Contrary to K I and K II which are symmetrical around the crack front center (γ = 90 • ), mode III SIF gets an antisymmetric behavior (see Fig. 20(c)). In the same manner as mode II, Comparing the three modes in terms of magnitude, it can be inferred that mode II dominates the other modes. Moreover, mode I is clearly lower than mode II and mode III. These observations support the choice here to assume a shear crack growth mechanism for RCF crack.

The analyzed results indicate that the crack experiences a highly complex multi-axial and non-proportional behavior, typically encountered during a rolling contact cycle.

To highlight the influence of residual stresses, the variations of SIFs in relation with the position of the contact center are plotted at γ = 45 • for case #1 and case #2 (see Fig. 24). As shown in Fig. 24(a), the introduction of compressive stresses promotes crack closure and consequently reduces the K I values. Furthermore, although the trends remain the same, the residual state causes a shift in the mean value of the K II and K III variations (see Fig. 24(b) and(c)). In other words, it induces an initial relative displacement between the crack faces that does not have a noticeable impact on the subsequent behavior of the crack.

Fig. 25 illustrates the variation of the SIF ranges along the crack front for cases #1 and #2. As previously stated, mode II prevails over the other modes. The maximum value of ∆K II is 6.13 MPa.m 0.5 for case #1 and 6.32 MPa.m 0.5 for case #2. These values, as well as the minimum of ∆K III , are γ = 90 • . The maximum of ∆K III , that is located at opposition positions γ = 50; 130 • , is equal to 1.91 MPa.m 0.5 for case #1 and 1.95 MPa.m 0.5 for case #2. The consistency and the quality of the results exposed above demonstrate the efficiency of the SAM/X-FEM coupling approach to tackle such complex problems. Indeed, the developed method provides smooth informations along the crack front that depicts accurately the particular behavior of crack subjected to rolling contact conditions. Moreover, the numerical stability of this tool shows the robustness of the development. Besides this, the benefits in terms of CPU times can also be highlighted. The calculation of the first loading cycle of case #1 requires 34 minutes with a laptop having the same characteristics as the one used in Section 4. This is quite acceptable for such simulations which include the resolution of frictional contact between the contacting bodies and the crack faces. Note that the introduction of residual stresses has a minor impact on calculation efforts (38 minutes for case #2).

Analysis of the crack propagation over several cycles

Results about crack growth are given in this section. Note that simulations have been carried out until the maximum depth of crack reaches 3 • z * . This corresponds to 26 cycle jumps. Fig. 26 shows the progression of the cracks during the propagation. As expected, the cracks grow in the most stressed zone, mainly along the travel direction (x-direction). Furthermore, it can be observed that the MSSR criterion predicts a nearly coplanar propagation for case #1 (see Fig. 26 In order to compare the cases studied, the simulated final crack shapes are superimposed in Fig. 27. It can be seen that these shapes are almost the same. This is not surprising given that the MSSR criterion uses only the SIF ranges to determine the crack direction.

As a guide, a -N curves are plotted in Fig. 28. It is found that the cracks speed gradually decreases as they extend in the depth. At the end of the simulation, the maximum number of cycles is 6.15 • 10 5 for case #1 and 5.42 • 10 5 for case #2. Therefore, the presence of residual stresses accelerates the propagation. This is due to the discrepancies in terms of ∆K II (see Fig. 25) that become more pronounced when the crack propagates.

Note that the SIF range threshold value ∆K th has not been taken into account in this analysis. The consideration of this parameter could arrest the crack and reduce the crack front extension in the zones where the SIFs are small.

Case #1 and case #2 require approximately 4 days by using

x y z Fig. 27 Comparison of the simulated final crack shapes. The black and blue shapes correspond to case #1 (without residual stress) and case #2 (with residual compressive stresses) respectively.

Fig. 28

Crack length a evaluated in the xz plane (symmetrical plane of the crack) against the number of cycles N. The black and blue curves correspond to case #1 (without residual stress) and case #2 (with residual compressive stresses) respectively. resources of a more powerful computing center. Despite this, the numerical costs remain relatively low for such simulations. Furthermore, no intervention of the user has been performed during the propagation, which underlines the robustness of the crack growth algorithm.

Conclusion

The present work proposes a novel numerical method to simulate in a fast, robust and accurate way the complex behavior of 3D non-planar crack under rolling contact conditions. To this end, a non-intrusive SAM/X-FEM coupling based on a global-local strategy has been developed together with a pow-erful data transfer procedure that allows the use of independent discretization for the global and local models. Besides significantly shortening computational time and resources, the developments allow to reduce the pre-processing efforts, which is well suitable in an industrial context. It offers a serious alternative that overcomes the current numerical limitations of other methods.

A validation based on a comparison with results from the literature and from a FEM model has been successfully performed. Then, the SAM/X-FEM approach has been applied to a practical example that focuses on the propagation of a semicircular crack under a spherical contact. It demonstrates the great potential of this strategy to handle in a fully automatic manner the 3D RCF crack growth. The influence of residual stresses on crack propagation has also been highlighted.

Concerning the potential improvements of the SAM/X-FEM method, the next step consists in extending the proposed one way coupling to an iterative one in order to consider the influence of the local detail at the global scale. In a future work, the proposed model will be used within an extensive parametric study that aims to understand the behavior of cracks under various rolling contact conditions. where µ contact is the friction coefficient between the contacting bodies. Note that the in-house SAM code is also able to deal with partial sliding conditions [START_REF] Gallego | Modeling of Fretting Wear Under Gross Slip and Partial Slip Conditions[END_REF][START_REF] Gallego | A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II and III[END_REF].

After solving the contact problem, the stresses in bulk can be calculated by discretizing the volume into a 3D regular grid of N x × N y × N z points spaced by ∆x along the x-direction, ∆y along the y-direction and ∆z along the z-direction. Following the same approach as above, the elastic stresses σ e are obtained by superimposing the contribution of each surface patch of pressure:

σ e I J (n) = N m=1 C p I J (m, n) + µ contact C q I J (m, n) p(m) (A.11)
where the coefficients of influence C p I J and C q I J correspond to the elastic stresses at the point n due to a rectangular patch of uniform pressure and shear at the point m. The expressions of C p I J and C q I J are given for example in Ref. [START_REF] Liu | Studying contact stress fields caused by surface tractions with a discrete convolution and fast fourier transform algorithm[END_REF]. As these analytical solutions implicitly consider the geometry of the halfspaces, the calculation zone may be restricted only to the region of interest where the stresses are needed. This contributes to a significant reduction in computational time.

Moreover, an initial residual stress field, denoted by σ i , can be considered by adding it to the stresses σ e generated by the elastic contact (A.11). Thus the total stress tensor σ is expressed as follows:

σ I J (n) = σ e I J (n) + σ i I J (n) (A.12)
Appendix B. X-FEM for 3D non-planar frictional contacting crack

In this appendix, the 3D quasi-static X-FEM model proposed by Pierres [START_REF] Pierres | 3D two scale X-FEM crack model with interfacial frictional contact: Application to fretting fatigue[END_REF][START_REF] Pierres | Experimental and numerical analysis of fretting crack formation based on 3D X-FEM frictional contact fatigue crack model[END_REF] for frictional contacting crack is recalled. Only the basics of the method are presented here. Partial details about its implementation in the industrial finite element software Cast3M [START_REF]Cast3M, Finite element software developed by the french atomic energy center (cea)[END_REF] are also provided. The main specificity of this mixed formulation is the weak coupling between the variables of the two problems by means of the Lagrange multiplier field λ. It allows the use of independent discretization scheme for the structure and the interface. Thus, the crack mesh can be freely adjusted to accurately capture the complex local contact non-linearities on the interface. Moreover, the subproblems have their own constitutive law. At the structure scale, the variables (u, σ) follow the stressstrain law for homogeneous isotropic material having a linear elastic behavior: It is worthwhile to note that ∆[w T ] corresponds to the increment of sliding displacement between two time steps and µ crack is the interfacial friction coefficient.

div σ + f = 0 in Ω (B.1) σ • n = F on ∂ σ Ω (B.2) u = ū on ∂ u Ω (B.
σ = K : ε(u) in Ω (B.

Appendix B.2. X-FEM discretization and numerical resolution

In order to take into account the crack influence at the structure scale, the displacement u are discretized by the means of X-FEM [START_REF] Sukumar | Extended finite element method in computational fracture mechanics: a retrospective examination[END_REF]. The method uses the partition of unity concept [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] to add enrichment functions into the classical FEM approximation. The generalized Heaviside function H and the singular F j i are introduced to capture respectively the discontinuity generated along the crack and the asymptotic behavior close to the crack front. The X-FEM approximation of the displacement u at a given time t ∈ [0; T ] is: where u i , a u i and b j i define the standard, discontinuous and singular degrees of freedom, N i correspond to the 3D linear finite element shape functions, N nodes are the nodes of the mesh defining the whole structure Ω, N crack are the nodes of the elements completely cut by the crack and N f ront are the nodes of the elements where the crack front is located (see Fig. B.32). Thanks to this enrichment technique, the mesh of the structure can be defined independently of the crack shape. Therefore, remeshing and field interpolation are not needed during the crack propagation.

u(x,
An explicit/implicit approach [START_REF] Prabel | Propagation de fissures tridimensionnelles dans des matériaux inélastiques avec xfem dans cast3m[END_REF][START_REF] Fries | Crack propagation with the extended finite element method and a hybrid explicit-implicit crack description: XFEM with a hybrid explicit-implicit crack description[END_REF] is employed to describe the geometry of 3D non-planar crack. First, the crack surface is discretized with flat triangular elements separately from the structure mesh (see Fig. B.33). This explicit description facilitates the update of the crack during propagation. Indeed, the crack extension is achieved by simply adding new triangular elements along the crack front. Then, a set of two level-sets [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets -Part I[END_REF][START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets -Part II[END_REF] is computed analytically with the help of the crack mesh. This implicit description is useful to define the X-FEM enrichment, the polar coordinates (r, θ) and the local axes close to the crack front. The consideration of the frictional contact between the crack faces requires the discretization of the displacement field w, the traction field t and the Lagrange multiplier field λ. The nodes of the triangular crack mesh are used as geometric support. In the same spirit as X-FEM, a discontinuous enrichment is used to distinguish the interface quantities on Γ + c and Γ - c : After introducing the discretized unknowns (B.15-B.18) into the three-field weak formulation (B.6), the quasi-static problem of 3D non-planar frictional crack is solved by using a stabilized non-linear LArge Time INcrement (LATIN) solver [START_REF] Gravouil | Stabilized global-local X-FEM for 3D non-planar frictional crack using relevant meshes[END_REF]. This iterative scheme prevents numerical oscillations classically observed on the dual variables in mixed formulations. An optimization of the solver parameters is proposed by Trollé et al. [START_REF] Trollé | Optimization of a stabilized X-FEM formulation for frictional cracks[END_REF] to ensure convergence quality and numerical efficiency.

w(x,

Movement 1 2 Fig. 1

 21 Fig. 1 Schematic representation of the main RCF crack growth scenarios.

Fig. 3

 3 Fig. 3 Global-local strategy for the resolution of the 3D crack problem under contact loading.

Fig. 4

 4 Fig. 4 Parametrization of the local model dimensions.

  Fig. 6 Construction with N D = 1 of the integration domain D for the computation of the SIFs at P j .

  Fig. 10 Semi-circular crack under a rolling spherical contact. (a) 3D view and (b) 2D view in the xz plane.

  Fig. 11 SAM/X-FEM model used for the validation (a = c).

Fig. 13

 13 Fig. 13Variation of the dimensionless SIFs at γ = 90 • during the loading cycle for cases #3, #4 and #5 (crack radius a = c/2). The dashed curves refer to the results extracted from Kaneta's reference work[START_REF] Kaneta | Mechanism of surface crack growth in lubricated rolling/sliding spherical contact[END_REF] while the solid curves represent the results obtained with the SAM/X-FEM strategy. The black, blue and red curves correspond to µ contact = -0.3 (case #3), µ contact = 0 (case #4) and µ contact = 0.3 (case #5) respectively. (a) mode I dimensionless SIF and (b) mode II dimensionless SIF.

Fig. 14

 14 Fig. 14 Evolution of the relative error ε II against the dimensioning parameter α i for different µ crack .

  Fig. 17Variation of the SIFs at γ90 • during the loading cycle. The black and blue curves correspond to the SAM/X-FEM model and the FEM model respectively. (a) mode I SIF, (b) mode II SIF and (c) mode III SIF.

Fig. 18

 18 Fig. 18Variation of ∆K eq along the crack front. The solid curve refers to the results extracted from the SAM/X-FEM model while the dashed curve represents the results obtained with the FEM model.

Fig. 20

 20 Fig. 20 Variation of the SIFs along the crack front during the loading cycle for case#1. (a) mode I SIF, (b) mode II SIF and (c) mode III SIF.

Fig. 21

 21 Fig. 21Schematic representation of the formation and the spread of the "bubble" opening zone during the rolling cycle.

22

  Crack opening displacement [w N ] for different positions of center.

Fig. 23

 23 Fig. 23 Evolution of normalized stresses generated by the contact along the xdirection. The solid curve refers to σ xx at y = z = 0 while the dashed curve corresponds to σ xz at y = 0, z = z * .

Fig. 24

 24 Fig. 24 Variation of the SIFs at γ = 45 • during the loading cycle. The black and blue curves correspond to case #1 (without residual stress) and case #2 (with residual compressive stresses) respectively. (a) mode I SIF, (b) mode II SIF and (c) mode III SIF.

Fig. 25

 25 Fig. 25 Variation of the SIF ranges along the crack front. The solid curves refer to case #1 (without residual stress) while the dashed curves represent case #2 (with residual compressive stresses). The black, blue and red curves correspond to ∆K I , ∆K II and ∆K III respectively.

  Fig. 26 3D view of the crack progression. (a) case #1 (without residual stress) and (b) case #2 (with residual compressive stresses).

Appendix B. 1 .

 1 Two-scale strategy and three-field weak formulation for 3D quasi-static frictional crack problem Let's consider a linear elastic body Ω containing a crack Γ c . The crack body is submitted to surface forces F on ∂ σ Ω and displacements ū on ∂ u Ω. Moreover, volume forces f are applied in Ω. n is the outward unit normal to Ω while n c is the outward unit normal to Γ c . Contact and friction can occur along the crack faces Γ + c and Γ - c (Γ c = Γ + c ∪ Γ - c ). The cracked body problem (Ω ∪ Γ c ) is divided into two subproblems: a problem related to the structure scale Ω and a problem related to the crack scale Γ c (see Fig. B.31). Each problem has its own primal and dual variables: the displacement field u and the stress tensor σ for the problem of structure, and the interface displacement field w and the interface traction field t for the problem of crack. At a given time t ∈ [0; T ], the solution (u, σ) of the problem of structure must satisfy the static equilibrium (B.1), and the Dirichlet and Neumann boundary conditions (B.2,B.3):

3 )

 3 Eqs. (B.4,B.5) ensure the coupling between the subproblems. These equations enforce the continuity conditions on Γ c of the primal (u, w) and dual (σ, t) variables:u = w on Γ c (B.4) σ • n c = t on Γ c (B.5)Eqs. (B.1-B.5) consist in the strong formulation of the fracture problem with contact and friction along the crack faces. An equivalent three-field (u, w, λ) weak formulation can be obtained by applying the principle of the virtual works:0 = -Ω σ : ε(u * )dΩ + Ω f • u * dΩ + Ω F • u * dΩ + Γ c λ • u * dS + Γ c (tλ) • w * + Γ c (uw) • λ * dS ∀u * ∈ U 0 , ∀w * ∈ W, ∀λ * ∈ Λ, ∀t ∈ [0; T ] (B.6)with:u ∈ U, U = {u + regularity/u = ū on ∂ u Ω} (B.7) u * ∈ U 0 , U 0 = {u * + regularity/u = 0 on ∂ u Ω} (B.8) w ∈ W, w * ∈ W, W = {w + regularity} (B.9) λ ∈ Λ, λ * ∈ Λ, Λ = {λ + regularity} (B.10)

  [START_REF] Bower | The influence of crack face friction and trapped fluid on surface initiated rolling contact fatigue cracks[END_REF] with K the Hooke tensor.At the crack scale, the variables (w, t) obey to the contact unilateral law and the Coulomb's friction law:Open Zone: [w N ] > 0 → t + = t -= 0 Contact Zone: [w N ] = 0 → t + = -t - Stick Zone: t T < µ crack • t N → ∆[w T ] = 0 Slip Zone: t T = µ crack • t N → ∃γ > 0/∆[w T ] = -γt + T (B.12)where the interface quantities are defined in the local basis (n c , t c ) attached to the crack:w = w N • n c + w T • t c and t = t N • n c + t T • t c (B.13)

Fig. B. 31

 31 Fig. B.31 Strategy for the cracked body problem.

  Fig. B.32 X-FEM nodal enrichments for the bulk and the crack interface.

Fig. B. 33

 33 Fig. B.33 Explicit description of the 3D non-planar crack geometry.

  t) ⋍

		N i (x)u i (t) +	H(x)N i (x)a u i (t)
	i∈N nodes		i∈N crack
		4		(B.15)
	+	N i (x)	F j b j i (t)
	i∈N f ront	j=1

  t) ⋍

	X-FEM Mesh	Triangular Crack Mesh
		Ψ m (x)w m (t) + H(x)Ψ m (x)a w m (t) (B.16)
		m∈N inter f ace

t(x, t) ⋍ m∈N inter f ace

Ψ ′ m (x)t m (t) + H(x)Ψ ′ m (x)a t m (t) (B.17) λ(x, t) ⋍ m∈N inter f ace Ψ ′ m (x)λ m (t) + H(x)Ψ ′ m (x)a λ m (t) (B.18)

Funding

The authors would like to thanks Airbus Helicopters and ANRT for funding this research (CIFRE grant # 2018/1631).

Appendix A. Semi-Analytical Method

The basics of SAM are recalled in this appendix. Since linear elasticity is assumed in the present work, the in-house SAM solver discussed here is only addressed to the resolution of 3D elastic contact problems.

Appendix A.1. Contact problem formulation

Let's consider two non-conforming contacting bodies denoted by Ω 1 and Ω 2 (see Fig. A.29). Small strain and small displacement assumptions are assumed. The materials are isotropic and homogeneous with linear elastic behavior in bulk. The size of the contact area A c is relatively small with regard to the characteristic dimensions of the bodies. Thus, Ω 1 and Ω 2 can be regarded as half-spaces. Generally, the formulation of the normal contact consists in a system of equations and inequalities which describe the physics of the problem at the interface. This formulation is recalled bellow:

• The load balance. The applied load W is equal to the integration of the contact pressure p(x) over the contact area A c :

• The surface separation. The distance between the contacting bodies h(x) is given by the summation of the initial separation h i (x), the rigid body displacement δ z and the normal elastic surface displacement of both bodies

• The contact condition. The distance h(x) is always positive due to the non-interpenetration of contacting bodies:

Parameters of the elastic contact problem.

For elastic half-spaces, the displacement u z can be expressed as a function of the contact pressure distribution p using the Boussinesq solution [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF]:

The elementary solutions G p z , usually called Green's functions or influence coefficients, correspond to the displacement of the contact surface at the point x due to an unit load at the point ξ. It should be noted that Eq. (A.6) satisfies the equilibrium in bulk, the linear elastic behavior and the half-space boundary conditions.

Appendix A.2. Numerical resolution

To numerically solve the contact problem, the potential contact surface is discretized into a 2D regular grid of N = N x × N y points which are spaced by ∆x along the x-direction and ∆y along the y-direction. Each point n is associated to a centered elementary rectangular zone of size dS = ∆x × ∆y where the pressure is assumed to be constant. In this manner, the pressure distribution p(x) is discretized into several patches of uniform pressure p(n) (see Fig. A.30).

By using the linear elasticity theory, the effect of each patch of pressure can be superimposed in order to determine the elastic response of the contact interface. In that respect, Eq. (A.6) can be rewritten in the following discretized form:

p(n)

The influence coefficients K p z represent the effect of a rectangular patch of uniform pressure at the point m on the normal elastic displacement at the point n. Love [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF] gave closed-form solution of K p z by integrating the Boussinesq concentrated load solution over a rectangular area. Detailed expressions can be found in the well-known book of Johnson [START_REF] Johnson | Contact mechanics[END_REF].

Then, Eqs. (A.1-A.5) are discretized and the contact problem is expressed by the following system:

where N c and N e define respectively the points inside and outside the contact area A c (N = N c + N e ).

The system of equations and inequalities (A.8) is solved simultaneously by the means of an iterative solver based on the Conjugate Gradient Method (CGM) [START_REF] Polonsky | A numerical method for solving rough contact problems based on the multi-level multi-summation and conjugate gradient techniques[END_REF]. This numerical resolution is generally coupled with Fast Fourier Transform (FFT) techniques [START_REF] Ju | Spectral analysis of two-dimensional contact problems[END_REF] to transform the time-consuming convolution product (A.7) into a simple term by term product in the frequency domain: ûz = K p z p (A.9)

In practice, the Discrete Convolution and Fast Fourier Transform (DC-FFT) technique, as proposed by Liu et al. [START_REF] Liu | A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses[END_REF], is used to significantly reduce the computational time while keeping a good accuracy of the solution.

Concerning the tangential contact problem, full sliding contact conditions are assumed in this work. As a consequence, the shear distribution q is simply calculated as follows: .10)