Valentine, Jérémy, Roozbeh, Tom, Tobias, Fatima Mauro Ramin

Jaehyun, Okan, Fabian Aurelia

Son, Thibault, Lucas, Manon Nicolas

Emma Birgitte

Étienne

Éric, Mathieu, Élias Quentin Yannick

Eugene Sarah

Amy Dimitri

Nicolas Filippo

Marc Gabriel

Estelle Nathalie Et Nelly

Alain Enfin

Soeur Cécile

Keywords: contact, plasticity, volume integral method, fast-Fourier transform, wear, up-scaling, rough surface, 1 Tzanakis et al., 2012. v contact, plasticité, méthode d'intégrale volumique, fast-Fourier transform, usure, transposition d'échelle, surface rugueuse, vi

I think the closing of the four-year long adventure that was my thesis can only be done by sharing the euphoria of success with the people who have made it possible. Foremost among them are my advisors Guillaume Anciaux and Jean-François Molinari. I could not have imagined a better environment for learning the craft of research than the one they have created.

With their encouragements and the challenges they put forward, I was able to better myself in every area my thesis touched on, from mathematics to mechanics, from programming to writing, from teaching to research. I absorbed so much from them I hope I was able to give back a little of what I learned along the way. I thank the jury members Marie Violay, Stéphanie Chaillat, Lars Pastewka and Brice Lecampion for giving me an opportunity to show that.

Perhaps the best quality of that research environment is the deep reservoir of collaborations that exists both in and out of the lab. I am grateful I was able to tap that great pool of people, who have been inspirational, helpful and friendly to me. I want to thank among them Marc Bonnet, whose tremendous help has allowed me to solidify the core of this thesis. I have enjoyed our discussions very much and hope to continue our collaboration on Fourier-based integral methods. I also want to acknowledge Jérémy Bleyer, whose expertise on optimization problems (among many other topics) amazes me, and who helped me understand the mathematical subtleties of contact, friction and plasticity.

The number of collaborations that I have had within the lab is too large to be recollected here. I have unforgettable relationships with all past and present members of the lab, both on a professional and personal level.

Abstract

All mechanical systems, naturally occurring or human-produced, are subjected to friction and wear at the interface of solid constituents. Large portions of energy dissipation and loss of material, in every-day life and industrial applications alike, are due to friction and wear. Mitigating their effects could save between 1% and 2% of the GDP of a developed country. 1 Some systems governed by friction and wear can have an even more important bearing on human lives, such as earthquakes nucleating from the sliding of tectonic faults. Despite the tremendous impact of tribological phenomena on society, their understanding has remained empirical, and to this day no predictive model has emerged. Interface processes such as friction and wear are difficult to investigate because of the large number of underlying physical phenomena (e.g. adhesion, fracture, etc.) and the difficulty of observing them at contact interfaces.

Although research endeavors into friction and wear have not produced predictive models, they have identified key components of tribological systems necessary to build such models. Central among them is the idea that solids may not be in contact across their apparent interface area, but instead a much smaller "true contact area." This true contact area is the result of the surfaces in contact being inevitably rough. In addition, contact pressures on roughness peaks are expected to cause plastic flow of material, drastically changing the properties of the contact interface, and the role it plays in tribological processes. Therefore, the aim of this PhD thesis is to develop tools for the modeling of elastic-plastic rough contact interfaces, and to study the applicability of knowledge of the contact state to the modeling of interface phenomena.

The first part of this objective is the development of a novel computational approach to volume integral methods, which are used to solve elastic-plastic rough surface contact. Volume integral methods have the advantage over the finite-element method in that they can represent exactly elastic constitutive behavior and semi-infinite bodies, which are commonly used in rough contact applications. This thesis develops a new fundamental solution used in a volume integral approach, which drastically improves computation times and required memory over previous approaches. Derived directly in the Fourier domain, this fundamental solution makes optimal use of the fast-Fourier transform while retaining the advantages of classical volume integral methods.

In the second part, this numerical approach is used to study the so called "Archard's wear coefficient", and to up-scale known micro-scale adhesive wear mechanisms to the macro-scale via rough contact simulations. These show that wear is an emergent process dependent on the interaction of micro-scale mechanisms: they demonstrate the role of plastic deformations in the crack nucleation process, and the necessity to look beyond the true contact area to Résumé Qu'ils soient naturels ou manufacturés, tous les systèmes mécaniques sont sujets au frottement et à l'usure aux interfaces entre leurs composants. Des pertes énergétiques et matérielles majeures, aussi bien au quotidien que dans les applications industrielles, sont dues au frottement et à l'usure. Les gains économiques que générerait un meilleur contrôle de ces effets sont estimés entre 1% et 2% du PIB d'un pays développé. 1 Plus important encore, certains systèmes soumis au frottement et à l'usure, comme les failles tectoniques, ont un impact sur des vies humaines. Malgré la considérable influence des phénomènes tribologiques sur nos sociétés, nous n'en avons qu'une compréhension empirique, sans capacité de prédiction. Les mécanismes physiques se déroulant aux interfaces, tels que le frottement et l'usure, sont complexes à analyser à cause des multiples processus dont ils sont issus (p. ex. l'adhésion, la rupture, etc.) et de la difficulté d'observation des interfaces de contact.

Malgré que la recherche n'ait pas produit de modèle prédictif, certains aspects clés des systèmes tribologiques ont néanmoins étés mis en évidence. Parmi eux, le fait que deux solides ne soient pas en contact sur toute leur interface, mais sur une aire réduite appelée « aire de contact, » joue un rôle central. Cette aire de contact résulte de l'inévitable rugosité des surfaces. De plus, les pressions de contact peuvent être suffisantes pour causer des déformations plastiques, changeant ainsi les propriétés de l'interface et son rôle des les processus tribologiques. L'objectif de cette thèse est donc de développer des outils pour la modélisation du contact entre solides élastoplastiques, et de transférer la connaissance de l'état de contact à des modèles de phénomènes d'interface.

La première partie de cet objectif comporte le développement d'une nouvelle approche numérique pour des méthodes d'intégrale volumique utilisées dans la résolution de problèmes de contact élastoplastique à surface rugueuse. Comparées aux éléments finis, les méthodes intégrales ont l'avantage d'exactement représenter le comportement élastique et les conditions aux limites infinies qui sont couramment utilisées pour des applications de contact de surfaces rugueuses. Cette thèse apporte une nouvelle solution fondamentale qui améliore considérablement le temps de calcul et le coût mémoire par rapport aux méthodes existantes. Cette solution fondamentale, dérivée dans le domaine de Fourier, fait un usage optimal de la fast-Fourier transform en conservant les avantages des méthodes intégrales volumiques classiques.

En seconde partie, cette approche est appliquée à l'étude de l'usure entre les échelles micro et macro à l'aide de simulations de contact élastoplastique. Ces dernières rendent explicite l'émergence de l'usure à partir de la rugosité des surfaces et des interactions entre des mécanismes à l'échelle micro. Le rôle prépondérant des déformations plastiques dans le processus de nucléation de fissure est démontré, ainsi que la nécessité d'aller au-delà de l'aire de contact Because the contact load must be carried over an area typically two orders of magnitude smaller than the apparent area of contact, the local pressures may cause plastic activity. The principal shear stress is maximum below the surface [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], so plastic deformations first develop below contacting points, then grow in volume and eventually reach the surface. It is common that disjoint surface contacts "share" a plastic region below the surface. These plastic strains modify significantly the stress state compared to an elastic material, which drastically changes the macroscopic picture of stress-driven phenomena such as wear. Note that this figure is the result of an actual elastic-plastic contact simulation. The original method used to obtain it is detailed in Chapter 2 of this manuscript. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.1. Generated self-affine rough surface. Color code indicates the asperity height. The "fractal" aspect of the roughness is represented by small asperities sitting on top of larger ones. . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.2. Histogram comparison of rough surface generation algorithms. Figure (a) shows the root-mean-square of heights, with the average value in dashed line and the colored zone containing ±σ around the average. The random phase algorithm by [START_REF] Wu | Simulation of Rough Surfaces with FFT[END_REF] gives consistent values around the analytical prediction, while the algorithm of [START_REF] Hu | Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis[END_REF] gives a much larger spread (approx. 10 times). Figure (b) shows the root-meansquare of slopes. Here the random phase algorithm does not reproduce the analytical value. This is due to the discrepancy between the continuous and discrete version of the PSD which is increased by the moment computation of eq. (1.47). This error decreases as the ratio L/λ l increases. . . . . . . . . 32 xi 1.3. Measures of the root-mean-square of slopes. Due to the finite difference approximation of the surface gradient, the geometric estimator is dependent on discretization and is therefore not a good measure of the root-mean-square of slopes. The spectral estimator computes the discrete version of eq. ( 1 (1950) solution and the proposed VIM.The displacements are shown along the x 1 = x 3 = 0 line. The agreement between the analytical and computed solutions is very good in the central part of the periodicity cell, while distortion induced by periodic conditions becomes apparent near its boundary.

2.3.

Normalized stresses due to a hydrostatic spherical inclusion in a halfspace. Comparison between the [START_REF] Mindlin | Thermoelastic Stress in the Semi-Infinite Solid[END_REF] solution and the proposed volume integral method. The stresses are shown along the x 3 axis. The red region is where the eigenstrain ε = αT I is imposed. One can observe a good agreement of the numerical approximation with the analytical solution outside the inclusion. Some Gibbs effect can be observed at the boundary of the inclusion, with spurious oscillations in the inclusion. This is due to the Fourier approximation of the discontinuous eigenstrain function.

2.4. Relative computation times for VIM and FEM. We compare the application of the operator N to an elastostatic FEM solve step (Cholesky factorization). The reference time T 0 is T 0 = 1.41 • 10 -3 s for the VIM and T 0 = 1.77 • 10 -2 s for the FEM. The scaling for large problem sizes agrees with the theoretical algorithmic complexities. For N = 2 21 the stiffness matrix factorization needed over 128 GB of memory, an amount two orders of magnitude larger than what is required for the VIM. . . . . . . . . . . . [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF]. Increase of the applied load beyond the initial yield shows that the pressure profile deviates from the elliptic Hertzian profile by flattening of the curve at the axis of symmetry, with a plateau whose extent increases with the load. The results of [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] however show oscillations of the pressure profile at high plastic loads which is not reproduced by our simulation. As there is, to our knowledge, no physical reason to these oscillations, they are likely due to the coarse discretization of the finite-element mesh they used in their study 2 .

Elastic perfectly-plastic Hertzian contact, comparison with

The simulation parameters are given in Appendix B.1.1 . . . . . . . . . . . 59 3.2. Elastoplastic Hertzian contact. Comparison between the elastic Hertz solution [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], a simulation using Algorithm 5 where the surface residual displacement is calculated by FEM, and a full VIM simulation.

Both simulations have identical surface discretization. As expected [START_REF] Johnson | An Experimental Determination of the Contact Stresses between Plastically Deformed Cylinders and Sphere[END_REF], the plastic pressure distribution deviates from the elliptical shape of the Hertzian distribution. The contact radius is larger in the plastic case.

The simulation parameters are given in Appendix B.1.2. . . . . . . . . . . . 60 3.3. Plastic strain norm |ε p | in elastoplastic Hertzian contact. Results obtained with Algorithm 5, with the plastic problem solved using our VIM approach (left) or a first order FEM (right). Note that in the FEM case, plastic deformations are piece-wise constant, whereas in the volume integral result they are interpolated between nodal values (cf. Section 2.3). Nonetheless, the solutions give similar plastic zone size and maximum plastic deformation norm. The simulation parameters are given in Appendix B.1.2. . . . . . . 60 4.1. Schematic representation of rough contact and corresponding wear mechanisms. A contact cluster forms a wear particle upon sliding if its area is larger than A * , otherwise the asperity in contact deforms upon sliding without breaking (hypothesis B.). This hypothesis brings asperity-level physics in the wear particle formation process through the critical cluster size A * . . . . 76 4.2. Evolution of contact area with increasing normal load. Simulation of forty realizations with λ l /λ s = 8 and fifty load steps in [0.001, 0.20] • W 0 . A * is taken as 2λ 2 s . (a) shows a rough surface sample. (b), (c) and (d) show the state of one realization at the loads indicated by dashed lines in (e). The video provided with the supplementary material shows the contact evolution for this realization. It can be seen from (b), (c) and (d) that the number and the size of clusters increase with W . (e) shows the combined effect of those two contact-area growth mechanisms on the total contact area A c . It also shows the size of the largest cluster, which increases dramatically when the growth of A c is dominated by cluster merging. Figure (f) shows the rate of increase in the number of clusters with respect to A c . Positive values indicate a regime where the growth of A c is dominated by cluster nucleation and negative values indicate that the merging of clusters dominates the contact-area growth. 79 4.3. Distribution of contact cluster areas. The main graph shows p(A/λ 2 s ,W ), the probability density function of normalized cluster areas for H = 0.8 and varying λ l /λ s , evaluated using twenty logarithmic bins. The inset shows the probability density function for values of H ∈ {0.6, 0.7, 0.8} with λ l /λ s = 128 . p(A,W ) can be approximated by a power-law within a certain cluster size interval, inside which it is independent of the applied load W . Increasing W or λ l /λ s increases the upper bound of the power-law interval. Normalization with the smallest wavelength λ s collapses all distributions to a single curve within the power-law interval. Varying the Hurst exponent has a limited effect on the resulting distribution. . . . . . . . . . . . . . . . . . . . . . . 4.4. Archard's wear coefficient in context of self-affine surface contact.

(a) shows the complementary cumulative probability distribution function (H = 0.8 ), (b) shows Archard's wear coefficient as a function of the applied load (for λ l /λ s = 128 and H = 0.8 ). Regardless of A * there exists a critical load at which the wear coefficient transitions from zero (i.e. no wear debris) to a constant value (i.e. steady-state mild wear regime). This critical load largely depends on the value of A * : the transition occurs at a higher critical load for contacts with lower interfacial shear strength and consequently larger A * (i.e. better lubrication condition). In (b) the y-axis is linear up to 10 -5 . . . . . .

4.5.

Wear coefficient and wear rate (H = 0.8 and λ l /λ s = 128). While the wear coefficient is non-constant with the load, its derivative decreases with the load. However, regardless of A * , the limit value of K is one, which does not correspond to experimental observations. Similarly, the wear rate is increasing non-linearly with the load, although its derivative stabilizes to a fixed value with increasing load. . . . . . . . . . . . . . . . . . . . . . . .

4.6.

Comparison between Archard's wear coefficient K PL and the proposed interpretation K PL for a power-law distribution of cluster areas (α = 1.5 ). As can be seen in the expressions of eqs. (4.9) and (4.10), the behavior of K is different from K when α < 2 : we have K → (a) shows that the local pressure can exceed values of 3σ y , the average pressure in contacts p = W /A c is closer to 2σ y , whereas the saturated model gives an average of 2.5σ y with large saturated portions of the micro-contacts. In this case, the normal load is W /W 0 is 6.5 • 10 -2 , and the saturated model predicts a contact area 20% smaller than the 2 prediction. As a result the connectivity of micro-contacts is different between the two models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Contact ratio and secant slope of the load/contact area relationship.

All models predict a sub-linear load/contact area law because the secant slope decreases with the load. On average the slope computed using the saturation pressure model is 20% lower than the slope computed with the 2 plasticity approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.3. Probability density function of surface pressures. Neither the elastic nor the saturated models qualitatively reproduce the pressure distribution of the elastic-plastic model. As expected, the pressure distribution of the saturated model tends to a Dirac distribution for p = p m , whereas the distribution for the elastic-plastic model tends to zero. Note that the results of [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF] have been renormalized (cf. companion notebook (Frérot, 2019b)). 95 5.4. Crack nucleation sites for the saturation pressure (a) and 2 plasticity (b) models. An "upwards" shear is applied on each contact. The true contact area is shown in light grey, the plastic zones in red and the crack nucleation sites in black. We can see that the 2 model has more crack nucleation sites than the saturation model. Since both models give approximately the same true contact area, this discrepancy must be due to plastic residual deformations which are not represented in the saturation approach. Figure (b) shows that there can be multiple crack nucleation sites per contact. . . . . . . . . 98 5.5. Crack nucleation likelihood (CNL) as a function of σ c and normal load. Figure (a) shows the results for an elastic contact, while (b) shows the saturation and 2 plasticity results. In (a) the CNL curves are uniformly shifted to the right when the load is increased, indicating an exponential increase in the CNL. The normalization of n crack by the true contact area makes explicit that this increase is due to stronger elastic interactions between contacts. The magnitude of the shift is given by the most frequent value of σ I called σ i/c . The two plasticity models in (b) have drastically different behavior: the crack nucleation is much more likely in the 2 approach because of plastic residual stresses, and the CNL experiences a faster decay for values of σ c > σ y . . The loading curves consist of two parts: an initial linear loading and a non-linear saturation of σ I . Although the initial slope is linear, the loading is not elastic as plastic deformations still evolve in the system. Their influence on the maximum tensile stress at the surface is however minimal. A sufficient shear stress may eventually create plastic strain at the surface, which causes a transition in the loading curve. Moreover, given a low enough yield stress, the initial indentation may cause surface yield, as seen for the curve where σ y /E = 2 • 10 -2 . . . . . . . . . . . 5.8. Failure regimes for a sheared spherical indentation. For the first row of graphs, the ratio σ y /E is varied for a perfectly plastic material. For the second row, the hardening ratio E h /E is varied for a yield ratio of σ y /E = 2 • 10 -2 .

On each graph the dashed line shows the yield stress. The white curve marks the transition between failure driven by slip rupture of the junction and failure by crack nucleation. The competition between the junction strength τ j and the critical stress σ c is influenced by σ y because of the saturation effect shown in fig. 5.7. While plasticity gives a failure mechanism independent of τ j , hardening allows the tensile stress to grow past the initial yield limit, giving a linear transition between failure mechanisms. . . . . . . . . . . . 5.9. Inter-contact stress as a function of contact ratio. Inset shows the probability density function of the largest stress eigenvalue at the surface. The stress value corresponding to the peak in the probability density is defined as the inter-contact stress σ i/c . We see that in the rough surface simulation and the Greenwood-Williamson model σ i/c evolves linearly with the contact ratio. The value of κ is estimated from the mean curvature of contacting summits in the rough surface. . . . . . . . . . . . . . . . . . . . . . . . . . 5.10. Crack nucleation likelihood as a function of re-normalized critical stress, normal load and junction shear strength. The material in contact behaves elastically. We normalize by h rms E * instead of κ because of the difficulty of defining asperity curvature on a rough surface. One can see on figure (a) that the curves corresponding to different normal loads are collapsed on a straight line, showing that the CNL does indeed follow the scaling established in eq. (5.12). When the strength of the junction is taken into account (figure (b)), or, equivalently, if a shear stress is applied, we observe qualitative deviations from eq. (5.12). There is a decrease in crack event density due to the interference of close contacts, which tends to unload the tensile stresses at the trailing edge of leading contacts. . . . . . . . . . . . . 5.11. Inter-contact stress as a function of contact ratio for the saturation and 2 plasticity models. While the curves do not match the analytical GW approach, their slopes are smaller than in the elastic case. The 2 model shows the smallest slope, indicating that the stresses due to plastic residual deformations have an influence on the inter-contact stress and actually reduce it compare to the underlying elastic stresses of the saturation pressure model.

5.12. Crack nucleation likelihood as a function of re-normalized critical stress and normal load for the saturation and 2 plasticity models.

The pressure saturation model reproduces a scaling similar to the elastic case in fig. 5.10a, while the 2 model shows a fundamentally different behavior.

The crack density is higher in the plastic case because of the additional tensile stresses caused at the edge of contacts by residual plastic deformations. . . . 107 C.1. Influence of L/λ l on the complementary cumulative distribution function of contact-cluster areas. Grey curves are results for each realization (20 in total), and red curves are ensemble averages. Surface sampled has λ l /λ s = 8 and λ s /∆l = 16 , so that only L is varied. As L/λ l increases, the dispersion of the individual realizations decreases, and the behavior of the ensemble average approximates accurately the behavior of each realization. We also note that as L/λ l increases, the distributions converge to a limit distribution. For the simulations carried out in the rest of this paper, a value of L/λ l = 8 , which is a good compromise between required number of realizations and computational cost, was selected. Note: we analyze the cumulative distribution instead of probability density to remove any bias due to binning. 125 C.2. Influence of λ s /∆l on the complementary cumulative distribution function of contact-cluster areas. Curves are ensemble averages of 20 realizations of a surface with L/λ l = 16 and λ l /λ s = 8 , so only ∆l is varied. As λ s /∆l is increased, the distributions converge to a limit distribution. Increasing λ s /∆l smoothens the distributions and reduces the systematic bias between the computed distribution and the limit of λ s /∆l → ∞. is often quoted for having said that "God made the bulk; surfaces were invented by the devil." Putting aside the grimness of this statement, an optimist would interpret that the challenging physics resides at interfaces. Although "the bulk" is not devoid of challenges, this quote holds some truth even today. While large scale, cutting edge experiments like CERN or LIGO are uncovering the most fundamental aspects of physics, the most trivial experiences, such as stopping a bike by pulling on the breaks or erasing a pencil mark are attributed to physical processes which are yet to be understood, namely friction and wear. This can also be said of large magnitude rare events such as earthquakes, which affect millions of lives each year, and are governed by … the very same interface phenomena.

The scientific investigation into friction and wear in the past sixty years has identified a quantity that plays a major role in the macroscopic friction and wear response: the true contact area. This was first propounded by [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF] in a different interface process: the conductivity of electric contacts. They show that nominally flat surfaces have an electric conductance four orders of magnitude smaller than the expected value based on the apparent contact area: in fact their results show that the conductance is independent of the apparent contact area and mainly depends on the normal load. They discuss the implications of this observation with Amontons' first law of friction 1 , i.e. the maximum friction force is proportional to the applied load [START_REF] Amontons | De La Résistance Causée Dans Les Machines[END_REF], or more classically F T ≤ µF N where µ is the proportionality coefficient. By comparing an elastic and a plastic model for the contact area, [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF] conclude that the electric conductivity magnitude cannot be explained by the elastic model, and plastic deformations must occur at the contacting asperities: since the contact load is supported on an area several orders of magnitude smaller than the apparent contact area, local contact pressures are expected to cause plastic flow of asperities. This was later corroborated using experimental methods allowing a direct imaging of the contact interface and a precise measurement of the true contact area (Dieterich and Kilgore, 1996;[START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF][START_REF] Zhang | Experimental and Finite Element Analyses of Contact Behaviors between Non-Transparent Rough Surfaces[END_REF]. Based on this observation, Bowden and Tabor (1942) propose an expression for the friction coefficient:

µ = shear strength flow pressure , (0.1)
that is based on the shear resistance of contact junctions and the plastic flow pressure of the material which governs the plastic contact behavior. Later, [START_REF] Dieterich | Modeling of Rock Friction: 1. Experimental Results and Constitutive Equations[END_REF] showed using rock sliding experiments that the kinematic friction force was dependent on the sliding velocity. His works were generalized in "rate-and-state" models [START_REF] Ruina | Slip Instability and State Variable Friction Laws[END_REF][START_REF] Gu | Slip Motion and Stability of a Single Degree of Freedom Elastic System with Rate and State Dependent Friction[END_REF] which include, via a state variable, a time-dependent effect which strengthens the contact. For example a possible formulation [START_REF] Sinai | Slow Rupture of Frictional Interfaces[END_REF] of the rate and state dependent friction coefficient is:

µ = µ 0 + a ln 1 + * + bµ 0 ln 1 + ϕ ϕ * ,
where is the sliding velocity, ϕ is the state variable and µ 0 , a, b, * , ϕ * are empirical quantities to be determined by experiments. This formulation is often coupled with an aging law [START_REF] Dieterich | Modeling of Rock Friction: 1. Experimental Results and Constitutive Equations[END_REF]:

dϕ dt = 1 - ϕ D ,
where D is also an empirical constant. This "aging", which was first observed by [START_REF] Coulomb | Théorie des machines simples en ayant égard au frottement de leurs parties et à la roideur des cordages[END_REF], is due to an increase over time of the real contact area [START_REF] Dieterich | Direct Observation of Frictional Contacts: New Insights for State-Dependent Properties[END_REF] and to the chemical aging of existing contacts [START_REF] Li | Frictional Ageing from Interfacial Bonding and the Origins of Rate and State Friction[END_REF][START_REF] Liu | Chemical Origins of Frictional Aging[END_REF].

Recently, the experiments of Jay Fineberg's group have given an unprecedented look at the contact interface and how it carries shear loading. They have in particular highlighted how fracture mechanics describes the propagation of slip pulses across the interface, even before the onset of macroscopic sliding [START_REF] Svetlizky | Classical Shear Cracks Drive the Onset of Dry Frictional Motion[END_REF][START_REF] Kammer | Linear Elastic Fracture Mechanics Predicts the Propagation Distance of Frictional Slip[END_REF]. The similarity between fracture and friction has been advanced further by [START_REF] Barras | The Emergence of Crack-like Behavior of Frictional Rupture: Edge Singularity and Energy Balance[END_REF] who have shown how fracture properties can emerge from rate-and-state friction models. This macroscopic view is necessarily the result of the collective behavior of individual contacts resulting from the unavoidable micro-scale roughness of surfaces2 .

Of course friction and wear are intimately linked: friction is necessary for wear particle formation, and the surface roughness evolution due to wear may in turn change the friction properties of the interface. The scientific investigation of wear of the mid-twentieth century [START_REF] Holm | Electric Contacts: Theory and Applications[END_REF][START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF][START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF][START_REF] Archard | The Wear of Metals under Unlubricated Conditions[END_REF] unveiled, similarly to friction, the basic laws of wear, which can be summarized in Archard's wear equation:

V = K W s H , (0.2)
which shows that the wear volume V is proportional to the normal load W times the sliding distance s and inversely proportional to the hardness H , with a proportionality constant K to be experimentally determined. Similarly to the expression of the friction coefficient by Bowden and Tabor (1942), it is apparent that this law was derived on the assumption of plastic deformations, as the term W /H represents the true contact area in a purely plastic contact model. Unlike friction however, where the µ is typically in the range of 0.1 -0.8 for unlubricated solids (with some exceptions, see e.g. [START_REF] Fontaine | Achieving Superlow Friction with Hydrogenated Amorphous Carbon: Some Key Requirements[END_REF], the proportionality factor K typically varies by orders of magnitude in 10 -7 -10 -1 (Rabinowicz, 1995).

Experiments being the only way to determine the proportionality factor K, also known as the wear coefficient, numerous experimental studies have proposed wear maps for various materials [START_REF] Hokkirigawa | Wear Mode Map of Ceramics[END_REF][START_REF] Riahi | Wear Map for Grey Cast Iron[END_REF][START_REF] Kong | An Interior-Point Algorithm for Elastoplasticity[END_REF][START_REF] Gahr | Sliding Wear of Ceramic-Ceramic, Ceramic-Steel and Steel-Steel Pairs in Lubricated and Unlubricated Contact[END_REF]. [START_REF] Meng | Wear Models and Predictive Equations: Their Form and Content[END_REF] have surveyed more than 180 wear equations, among which "no single predictive equation or group of limited equations could be found for general and practical use." They additionally point out as a contributing factor to this issue "the slow pace of translation of microscopic observations into macroscopic models", which is still the consensus today [START_REF] Vakis | Modeling and Simulation in Tribology across Scales: An Overview[END_REF]. Microscopic wear observations have flourished with the advent of nano-scale experimental techniques in the early 21 st century, which provided the means to develop an understanding of wear processes at the atomic scale [START_REF] Liu | Method for Characterizing Nanoscale Wear of Atomic Force Microscope Tips[END_REF]. Similarly, discrete modeling approaches have made recent breakthroughs in nano-scale wear. [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF] have uncovered a critical length-scale controlling wear particle formation at the asperity level. This length-scale stems from the Griffith energy balance between the deformation energy available in the contact of two asperities and the required energy to break off the wear particle from the two surfaces. Similarly, [START_REF] Candela | The Minimum Scale of Grooving on Faults[END_REF] have shown that the anisotropy in the roughness of slip surface in geological faults disappears below a given length scale (∼ 10 µm). This suggests an interplay between the roughness, the fracture processes that cause wear and the plastic deformation of the material due to the contact pressures (cf. fig. 0.2), which can have large-scale consequences on the mechanical properties of a fault. Indeed, plastic deformations dramatically increase the true contact area, as well as increase the local stress state of asperities in contact, which can modify how a crack develops.

These examples illustrate the need to understand the contact interface with plastic flow of the materials in contact. While the true contact area may control the magnitude of certain Because the contact load must be carried over an area typically two orders of magnitude smaller than the apparent area of contact, the local pressures may cause plastic activity. The principal shear stress is maximum below the surface [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], so plastic deformations first develop below contacting points, then grow in volume and eventually reach the surface. It is common that disjoint surface contacts "share" a plastic region below the surface. These plastic strains modify significantly the stress state compared to an elastic material, which drastically changes the macroscopic picture of stressdriven phenomena such as wear. Note that this figure is the result of an actual elastic-plastic contact simulation. The original method used to obtain it is detailed in Chapter 2 of this manuscript.

macroscopic quantities (such as thermal/electrical conductivity, sealing properties, etc.) and has been the subject of numerous studies in elastic contact [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF][START_REF] Bush | The Elastic Contact of a Rough Surface[END_REF]Persson, 2006;[START_REF] Hyun | Elastic Contact between Rough Surfaces: Effect of Roughness at Large and Small Wavelengths[END_REF]Yastrebov, Anciaux, and Molinari, 2017a;[START_REF] Pastewka | Contact between Rough Surfaces and a Criterion for Macroscopic Adhesion[END_REF], other aspects of the contact interface, e.g. the stress state in the vicinity of micro-contacts, may be of greater importance in the accurate description of tribological phenomena. Because the contact interface is where micro-scale mechanisms start interacting to give the macro-scale picture that we experience, its full understanding is necessary to bridge the gap between the micro-and macro-scales.

Objectives

While next-generation experimental techniques [START_REF] Garabedian | Quantifying, Locating, and Following Asperity-Scale Wear Processes Within Multiasperity Contacts[END_REF][START_REF] Sahli | Evolution of Real Contact Area under Shear and the Value of Static Friction of Soft Materials[END_REF][START_REF] Zhang | Experimental and Finite Element Analyses of Contact Behaviors between Non-Transparent Rough Surfaces[END_REF][START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF] allow unprecedented detail of the topography of the contact area and contact evolution during sliding, some important properties of the interface, such as contact pressures, underlying plastic deformations, etc., are difficult to gather experimentally. In addition, obtaining statistically significant results is experimentally involved and time consuming. For these reasons, modeling and simulations have proven efficient in the study of rough contact interfaces. While the earliest efforts in modeling contact of self-similar surfaces can be traced back to Archard (1957), simulations were not able to handle the multi-scale nature of rough surfaces until the late twentieth century [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF]Polonsky and Keer, 1999b). Since then, numerical investigations into the true contact area in elasticity have flourished (e.g. [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF][START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF][START_REF] Campañá | Practical Green's Function Approach to the Simulation of Elastic Semi-Infinite Solids[END_REF][START_REF] Putignano | A New Efficient Numerical Method for Contact Mechanics of Rough Surfaces[END_REF]. However, the previous works that form the basis for the undertaking of this thesis clearly show that some plastic activity must occur at the contact interface (e.g. [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF][START_REF] Archard | The Wear of Metals under Unlubricated Conditions[END_REF][START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF].

In this regard, the work of [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF] has pioneered the study of elastic-plastic rough surface contact. Although they were not the first to carry out such simulations (Mayeur, Sainsot, and Flamand, 1995, in 2D only), their analysis of the influence of plastic behavior on the true contact area in 3D as well as the contact topography and pressure distribution has shown that elastic contact models cannot realistically capture the relevant quantities for the understanding of tribological phenomena. The results presented unfortunately contain representativity and discretization inaccuracies, but are still unique in the rough contact community at large: a review by [START_REF] Vakis | Modeling and Simulation in Tribology across Scales: An Overview[END_REF] summarizes the need for realistic contact models incorporating plastic constitutive behavior in the comprehension of tribological processes.

The present thesis positions itself with the intent of developing a framework for the modeling of interface phenomena by focusing on the simulation of elastic-plastic rough contact systems, expanding on existing numerical methods to achieve this task. Precisely, the objectives of this thesis are: Contact modeling -to provide an efficient numerical method to model the realistic contact of elastoplastic solids with rough surfaces. Boundary integral methods (e.g. [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF] have been very successful in their application to elastic contact with rough surfaces.

While their volume counterparts have been applied to elastic-plastic contact [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF][START_REF] Wang | Numerical Simulation and Analysis for 3D Elastic-Plastic Rough Contacts[END_REF], it is unclear whether the existing formulations are competitive with a finite-element approach. The first concrete objective of this thesis is to develop a volume integral method that can handle the performance requirements of rough surface elastic-plastic contact.

collective behavior of individual contacts. How can one transfer knowledge of an asperity-level mechanism to the whole contact interface? How can one rationalize the wear coefficient of Archard's law (K in eq. (0.2)) in a multi-asperity contact setting? How is the local crack nucleation process in adhesive wear influenced by the global contact response?

Approach

As mentioned above, we wish to develop a volume integral method for the modeling of rough contact interfaces. Surface roughness in natural surface occurs on many length scales which cannot be separated [START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF][START_REF] Renard | Constant Dimensionality of Fault Roughness from the Scale of Micro-Fractures to the Scale of Continents[END_REF]: numerical methods must therefore enable the full-scale modeling of the contact of such surfaces. Boundary integral methods make this possible for elastic contact because of the translational invariance that arises in typical rough surface contact simulations. It is indeed common to model the rough surface as undeformable and the deformable solid as a flat elastic half-space [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]. This translational invariance turns the integral representation of the solution displacement into a convolution product which can be accelerated via a fast-Fourier Transform.

Such an approach has already been used for volume integral methods [START_REF] Sainsot | A Numerical Model for Elastoplastic Rough Contact[END_REF]Zhou, Jin, Wang, Wang, et al., 2016), but it relies either on the discrete Fourier transform of physical-space Green's functions and their full storage or 3D discrete Fourier transforms which have periodicity issues and require large discretized domains. We plan on developing a volume integral method that makes use of a Green's function that is directly available in the Fourier domain, thus reducing the computational burden.

The outcome of the developments on volume integral methods is an efficient simulation procedure for the contact of rough surfaces that we use in conjunction with recent findings on the adhesive wear process of pairs of contacting asperities. [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF] have shown the existence of a critical length-scale that governs the formation of a wear particle in the sliding contact of two hemispherical asperities. It gives a simple criterion on particle formation at the asperity scale: if the contact diameter of the two asperities is larger than the critical length-scale, the asperities will break an wear off. We will apply this knowledge to multi-asperity contact to rationalize the wear coefficient, which Archard interpreted as the probability of an asperity-pair forming a wear particle. Since the particle formation process has been shown to be deterministic, that probability can be epistemically understood within a stochastic rough surface contact.

We will also apply concepts of fracture mechanics to elastic-plastic rough surface contact to determine the location of crack nucleation sites and their density, and how the interasperity interactions and plastic deformations may influence them. Since fracture mechanics is dependent on the stress state, we expect that the true contact area may not be the major quantity driving the crack nucleation.

Outline

This dissertation is comprised of five chapters separated in two parts. Below is a brief description of each of them:

Chapter 1 -This chapter presents a state of the art on the modeling of rough surface elastic contact and its associated numerical approaches. I present the derivation of boundary integral methods applied to the equilibrium of an elastic half-space, and in particular their formulation in the Fourier domain, which allows for significant computational advantages.

We then show how the tools describing the half-space equilibrium are employed in the context of rough contact mechanics with the quadratic programming techniques typically used in this field. Finally, we provide a short review on the statistics of rough surfaces and the algorithms used to generate synthetic surfaces.

Chapter 2 -We present here the derivation of the Mindlin fundamental solution, an essential component of the volume integral method we wish to deploy, directly in the Fourier domain. We also provide a sound mathematical framework for the application of volume integral methods to periodic systems with the discrete Fourier transform. We validate our approach with analytical solutions and compare the methods performance to the finite-element method.

Chapter 3 -Two solutions strategies for the coupling of the contact and plasticity problems are presented here. We describe first a fixed point algorithm and validate it with finite elements simulations and published results from the literature. We also propose a new solving scheme based on the primal-dual interior point method that is often used in second-order conic programming. We present here the variational formulation of the elastoplastic contact problem which forms the starting point of the method. We then describe the method in its essential ideas. This is however ongoing work.

Chapter 4 -In this first application chapter, we use elastic rough surface mechanics to examine how the critical length-scale concept in adhesive wear proposed by [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF] can be up-scaled to obtain the macroscopic wear response of a contact interface. We focus our attention on the wear coefficient of Archard's model, how it can be interpreted at the macro-scale, and how it can be derived without relying on some of Archard's assumptions on the nature of adhesive wear at the micro-scale.

Chapter 5 -Finally, we apply our numerical methods for elastoplastic rough contact to study the number of cracks nucleating at the edge of contacts, a necessary part of adhesive wear, and how it depends on plastic behavior and shear resistance of the contacts. We propose analytical arguments to rationalize the observed behavior of the rough contact interfaces.

We finish this dissertation with a summary of the main results and discussion of the future work that could be carried out on the grounds of this thesis.

Part I.

Methods

Elastic rough contact: boundary integral methods

Abstract

A state of the art of continuum linear elastic rough surface contact is presented, with emphasis on boundary integral solution methods of the equilibrium of an elastic half-space. The Fourier approach taken in this thesis is introduced to represent the elastic response of a half-space to surface tractions. This lays the foundation for numerical methods used for the rough surface contact problem. The latter is formulated as a constrained quadratic program and associated solution strategies are described. Finally, rough surface generation algorithms are compared according to their statistical properties. This chapter provides the theoretical basis for the development of more complex numerical methods and result discussions.

Disclaimer

Parts of this chapter are reproduced from the article [START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF]. "A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact". In: Computer Methods in Applied Mechanics andEngineering 351, pp. 951-976. doi: 10.1016/ j.cma.2019.04.006, with permission of all authors. My personal contributions to this article include the development of the method, the code implementation, the running of simulations, the figure production and the writing of the article. T he field of contact mechanics can be said to have originated from the ground breaking works of Heinrich Hertz on the deformation of optical lenses in 1882. His solution of two elastic spheres in contact has seen countless uses in both research and industrial settings. In 1939, Westergaard provides the first multi-asperity contact solution by using a complex potential method for the contact of a sinusoidal surface with a deformable solid [START_REF] Westergaard | Bearing Pressures and Cracks[END_REF]. That same year, Bowden and Tabor show that friction is intimately linked to the true contact area of mating surfaces [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF], while Archard makes the same statement for adhesive wear more than a decade later [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF]. It is with the goal of quantitatively studying the physics at the true contact area that various methods for rough surface contact are then developed. The seminal works of Archard (1957) and [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF] paved the way for multi-asperity contact via fractal and statistical methods respectively, which evolved over time to include asperity curvature [START_REF] Bush | The Elastic Contact of a Rough Surface[END_REF], long-range elastic interaction [START_REF] Vergne | Elastic Behavior of Multiple Contacts: Asperity Interaction[END_REF], surface anisotropy [START_REF] Bush | Strongly Anisotropic Rough Surfaces[END_REF], surface power spectral density (Persson, 2006), etc. Numerical methods have also been developed and used to study rough surface contact. Although finite elements methods have been successfully employed [START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF][START_REF] Yastrebov | Rough Surface Contact Analysis by Means of the Finite Element Method and of a New Reduced Model[END_REF], boundary integral methods have proven very efficient in dealing with elastic rough contact. These include spectral methods [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF][START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF]Rey, Anciaux, and Molinari, 2017) and Green's functions methods (Polonsky and Keer, 1999b;[START_REF] Campañá | Contact Mechanics of Real vs. Randomly Rough Surfaces: A Green's Function Molecular Dynamics Study[END_REF][START_REF] Putignano | A New Efficient Numerical Method for Contact Mechanics of Rough Surfaces[END_REF]. They have been used to study the true contact area evolution in adhesionless [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF] and adhesive contact [START_REF] Carbone | Adhesive Contact of Rough Surfaces: Comparison between Numerical Calculations and Analytical Theories[END_REF][START_REF] Pastewka | Contact between Rough Surfaces and a Criterion for Macroscopic Adhesion[END_REF]Rey, Anciaux, and Molinari, 2017), as well as interfacial separation [START_REF] Almqvist | Interfacial Separation between Elastic Solids with Randomly Rough Surfaces: Comparison between Theory and Numerical Techniques[END_REF], the autocorrelation of the surface stresses and micro-contacts [START_REF] Campañá | Elastic Contact between Self-Affine Surfaces: Comparison of Numerical Stress and Contact Correlation Functions with Analytic Predictions[END_REF][START_REF] Ramisetti | The Autocorrelation Function for Island Areas on Self-Affine Surfaces[END_REF], the distribution of the areas of micro-contacts [START_REF] Frérot | A Mechanistic Understanding of the Wear Coefficient: From Single to Multiple Asperities Contact[END_REF][START_REF] Müser | Contact-Patch-Size Distribution and Limits of Self-Affinity in Contacts between Randomly Rough Surfaces[END_REF], the autocorrelation of sub-surface stresses [START_REF] Müser | Internal, Elastic Stresses below Randomly Rough Contacts[END_REF], etc. The contact mechanics challenge organized by Müser in 2016 references as many as 13 different techniques for adhesive elastic contact, from 12 different research groups, studying some of the quantities previously mentioned [START_REF] Müser | Meeting the Contact-Mechanics Challenge[END_REF].

In this introductory chapter, we present a state of the art on numerical methods used in rough surface elastic non-adhesive contact, with an emphasis on boundary integral methods and their formulation in the Fourier domain, as it will lay the theoretical grounds for the methods developed in Chapter 2. We discuss the quadratic programming techniques used for the contact problem as well as the characterization and generation of self-affine rough surfaces which are commonly employed to represent realistic contact interfaces.

Equilibrium of elastic bodies: problem introduction

Generally, solving a contact mechanics problem involves computation of unknown quantities located both within the compliant solids (e.g. the displacement field) and at the contact boundary (e.g. the contact tractions). The former are given as the solution of the problem of a deformable body in equilibrium, which is established in this section for the specific case of linear elasticity. This particular topic having a long history of both analytical and numerical solution methods, we will only focus on state-of-the-art procedures currently employed in the field of rough surface mechanics. 1.1.1. Variational form Let the open B ⊂ R 3 be a deformable elastic solid of boundary ∂B with outward normal n. Let u be the displacement vector field of B. The linearized strain tensor ε and the Cauchy stress tensor σ are respectively given by:

ε[u] := 1 2 ∇u + ∇u T , (1.1) σ [u] := C : ε[u], (1.2)
where C ∈ R 3×3×3×3 is the elasticity tensor, satisfying the usual ellipticity and (major and minor) symmetries. Here and thereafter, we follow the usual convention whereby the gradient operator ∇ acts "to the right", so that e.g. (∇u) i j = ∂ j u i = u i,j . Let n be the external normal of ∂B, we define the traction vector of the displacement field u as:

T [u] := σ [u] ∂ B
• n.

(1.3) Given a prescribed displacement u D applied on ∂B u and prescribed traction t D applied on ∂B t , the Dirichlet and Neumann boundary conditions are expressed as:

u ∂ B u = u D , and T [u] ∂ B t = t D .
(1.4)

The equilibrium of B subjected to the above boundary conditions and to the volumetric force density b can be expressed by the principle of minimum potential energy:

inf u ∈KA(B) I (u) = 1 2 ∫ B σ [u] : ε[u] dV - ∫ B b • u dV - ∫ ∂ B T [u] • u dS , (1.5)
which involves finding u ∈ KA(B), i.e. a kinematically admissible displacement field, that minimizes the total potential energy of the system. When complemented with the appropriate regularity hypotheses on B, b and boundary data (see e.g. Dacorogna, 2015, chap. 4), the above formulation is shown to be equivalent to the strong Euler-Lagrange equation: (1.6) which expresses the local balance of momentum. The reduced regularity requirements of eq. ( 1.5) make it more attractive than the above equation for use in numerical methods (see

div(σ [u]) + b = 0 a.e. in B,
1.
Elastic rough contact: boundary integral methods e.g [START_REF] Mason | Methods of Functional Analysis for Application in Solid Mechanics[END_REF]: indeed the function space KA is typically the Sobolev space H 1 := W 1,2 whose members are in L 2 as well as their first order partial derivatives (in the weak sense). The variational form (1.5) also has the added benefit of easily accommodating the contact boundary conditions that we will introduce in Section 1.3. Since C is symmetric positive definite, the elastostatic problem has one and only one solution u which minimizes the potential energy functional I . The stationarity condition of I yields the weak Euler-Lagrange equation (see e.g. Dacorogna, 2015, chap. 3) that expresses the global balance of momentum:

∫ B σ [u] : ε[v] dV = ∫ B b • v dV + ∫ ∂ B T [u] • v dS, ∀v ∈ KA(B). (1.7)
This variational equality is known as the "principle of virtual work" and expresses the fact that the potential energy of a system in equilibrium should stay constant under any admissible displacement perturbation v.

The variational equality (1.7) is the starting point of many families of numerical methods, including the finite element method and the boundary integral methods that we discuss in this chapter. The former classically applies the Galerkin approach which consists in solving eq. (1.7) for an unknown u h in a finite-dimensional subspace of functions that approximates KA(B), and using the same approximate space for the test functions. It is equivalent to solving the minimization problem (1.5) with u in a computationally tractable subspace of KA(B). The finite elements method has been lauded for its simplicity and flexibility, and is still today ubiquitous in all areas of computational solid mechanics.

Although it is widely used in contact mechanics (see e.g. Wriggers, 2006, for various contact methods), it has seen relatively limited application in the specific field of rough surface contact. This is in part due to the multi-scale roughness found in natural surfaces [START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF], which leads to discretization requirements too great for the finite elements method. In addition, contacts at the micro-scale are usually modeled with semi-infinite bodies, which in finite elements typically involves modeling a large domain to reduce the effect of the boundary on the zone of interest.

In this regard, boundary integral methods alleviate the above shortcomings of the finite element method for elastic rough contact. They are also based on the above variational equality, which is used as a starting point to derive the main principles of the method.

Boundary integral equation

The basic idea of integral methods is to take advantage of the linearities in the considered problem to express the sought solution as a combination of solutions known a priori (with appropriate incorporation of the problem data). To build a solution method on this principle, let

(u 1 , b 1 , t 1 D , u 1 D ) and (u 2 , b 2 , t 2 D , u 2
D ) be two solution states, each comprised of the resolved displacement field and the boundary data. One can apply eq. (1.7) to state 1 with u 2 as a test function and conversely. Taking the difference of the two resulting identities produces

∫ B σ [u 1 ] : ε[u 2 ] -σ [u 2 ] : ε[u 1 ] dV = ∫ B b 1 • u 2 -b 2 • u 1 dV + ∫ ∂ B T [u 1 ] • u 2 -T [u 2 ] • u 1 dS.
(1.8)

Due to the symmetry properties of C, the left-hand side cancels and we obtain the well-known Maxwell-Betti reciprocity identity:

∫ B b 1 • u 2 -b 2 • u 1 dV = ∫ ∂ B T [u 2 ] • u 1 -T [u 1 ] • u 2 dS. ( 1
div(C : ε[U k ]) + δ x e k = 0 ⇔ N[U k ](x, •) = δ x e k , (1.10)
where we have defined the Navier operator N[u] =div(C : ε[u]) (with differential operations acting on the variable of U k (x, ) unless otherwise noted).

Due to the convolution properties of the Dirac distribution, the Maxwell-Betti identity applied to a fundamental solution becomes the integral representation of the unknown displacement (dropping the exponents 1, 2 ):

u(x) • e k = ∫ ∂ B T [u]( ) • U k (x, ) -T [U k ](x, ) • u( ) dS + ∫ B U k (x, ) • b( ) dV , for x ∈ B.
(1.11)

The source point x cannot lie on ∂B because the function T [U k ] is singular at = x and not properly integrable [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF]. Nonetheless, this representation clearly shows the expression of the displacement as a superposition of known solutions. However, the boundary data (T [u], u ∂ B ) that is needed is not fully known on ∂B: the tractions are only known on ∂B t and the displacements only on ∂B u . Assuming some regularity on u and applying a limiting process to adequately account for the singularity in T [U k ], one can construct a regularized boundary integral equation suitable for determining the missing boundary data. This equation is at the heart of boundary element methods, which will not be developed in this manuscript.

The reader is nevertheless referred to the book of [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF] for a detailed picture of these methods, which extend far beyond the realm of contact mechanics.

The methods we wish to present now are also based on a boundary integral equation and therefore belong, like the boundary element methods, to the family of Boundary Integral Methods (BIM). We now introduce the hypothesis that ∂E ⊇ ∂B. We assume to have at our disposal a fundamental solution that additionally verifies the boundary condition T [U k ] = 0 , i.e. ∂E is a free surface. This cancels the non-integrable traction vector and gives an explicit expression for the unknown displacement at all x ∈ B ∪ ∂B: .12) As this integral representation underlies all the numerical methods for rough surface elastic contact that are based on continuum mechanics, we will now explore how it can be used.

u(x) • e k = ∫ ∂ B T [u]( ) • U k (x, ) dS + ∫ B U k (x, ) • b( ) dV . ( 1

Equilibrium of a half-space

It often happens in rough surface contact that the scales of interest permit the simplification of removing macroscopic effects and assuming the bodies in contact are two half-spaces with rough surfaces. In linear elastic friction-less contact one can apply to this simplified problem the following hypothesis [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]: the problem of two contacting isotropic elastic solids with rough surfaces is equivalent to the problem of a flat half-space with modified elastic properties in contact with a rigid indenter having an equivalent roughness profile. This simplified formulation of the contact problem has two major advantages when an integral equation solving strategy is employed:

1. the simplified half-space geometry means that fundamental solutions can be derived directly for the domain of interest (i.e. E = B),

2. the problem has a horizontal translational symmetry and the integral representation becomes a horizontal (partial) convolution.

We now provide a rigorous setting to the various mathematical tools that form the backbone of the solution method and fully exploit the above characteristics.

Integral operators and Fourier formulation

Let B := { ∈ R 3 :

• e 3 > 0}, according to the semi-infinite assumption previously described, with outward normal n := -e 3 . Points ∈ B will often be denoted as = ( ˜ , 3 ) with ˜ := ( 1 , 2 ). We rewrite eq. (1.12) in the form (assuming b = 0 ) u = M[p], (1.13) where p ( = T [u]) is the known traction distribution on ∂B2 and M is a boundary integral operator defined in terms of fundamental solutions. The constitutive uniformity and semiinfinite geometry of B together imply a translational invariance in any horizontal direction, i.e.

M[p(•

+ h)](x + h) = M[p](x), ∀x ∈ B, ∀ h = (h 1 , h 2 , 0) ∈ R 2 × {0}. (1.14)
Consequently, the integral operator M can be expressed as a convolution with respect to the horizontal coordinates. There is a tensor-valued function G(z, , 1900) such that:

x 3 , 3 ) = e k ⊗ G k (z, x 3 , 3 ) (Fredholm
M[p] = ∫ ∂ B G( ˜ -x, x 3 , 0) • p( ˜ ) d ˜ .
(1.15)

We accordingly define the partial convolution operation along the horizontal coordinates .16) where the operands f , may in addition depend on the vertical coordinates x 3 , 3 ; moreover, the convolution symbol will implicitly retain any tensor contraction involved in the integrals being recast in convolution form. Reformulating eq. ( 1.15) using definition (1.16) gives

x = (x 1 , x 2 ) by (f )( x) := ∫ R 2 f ( ˜ -x) ( ˜ ) d ˜ , ( 1
M[p] = G(•, x 3 , 0) p.
(1.17)

Remark 1. We emphasize that the partial convolution (1.16) makes sense only provided its operands f , are "convolvable", i.e. are in some sense mutually compatible. It is in particular well-defined in the classical sense if f , are either (a) integrable functions over R 2 (i.e. f , ∈ L 1 (R 2 )) or (b) locally integrable functions whose supports are convolutive3 , and only in the sense of distributions, a.k.a. generalized functions, or if (c) f is a tempered distribution and is a compactly supported distribution4 . We will shortly see that such considerations play an important role in this work, which involves convolutions of types (b) or (c) as kernels such as G derived from fundamental solutions may be locally-integrable but are not integrable.

Remark 2. Notice the slightly unconventional definition (1.16) of f , adopted here because of its convenience in the formulation of elastic potentials.

The (horizontal) convolutional form taken by the relevant integral operators, reflecting the aforementioned translational invariance, prompts the use of the two-dimensional partial Fourier transform, defined by

f (q, 3 ) = F [f ](q, 3 ) := ∫ R 2 f ( ˜ , 3 )e -iq • ˜ d ˜ (1.18)
for functions of the variable ˜ that are in L 1 (R 2 ). The coordinate space (q, 3 ) ∈ R 2 × R + underpinning the above Fourier transform will be referred to as the partial Fourier space, with q = (q 1 , q 2 ). Thanks to the celebrated Fourier convolution theorem, convolutions become mere multiplications upon application of the Fourier transform (1.18). In particular, for integral operators defined by horizontal convolution, we have Lemma 1. Under the partial Fourier transform defined by eq. (1.18) and its extension to distributions, the partial convolution product (1.16) obeys the identity 1.19) for any pair (f , ) that is convolvable in the sense of condition (a), (b) or (c) of Remark 1.

F f (q) = f (-q) (q), (
Consequently, subject to the same type of convolvability conditions, integral operators of the form of M have the expression

M[p](q, x 3 ) = G(-q, x 3 , 0) • p(q).
(1.20)

The integral in eq. (1.18) is well defined in the classical sense only in case (a), and the distributional extension of the Fourier transform must be used instead for cases (b), (c), see e.g. the appendix "Distributions" in [START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology: Volume 2 Functional and Variational Methods[END_REF].

This lemma makes systematic use of the Fourier transform the cornerstone of the methods presented in this manuscript: it has the dual advantage of leveraging the algorithmic complexity of the fast-Fourier transform (FFT) for the computation of convolutions and of providing a simplified mathematical framework for dealing with both fundamental solutions and discretized integral operators, as it will become apparent in Chapter 2 with the application of integral methods to the volume.

The Boussinesq-Cerruti fundamental solution

The fundamental solutions used in this work are, as mentioned in Definition 1, displacement fields created by singular point sources in the half-space B endowed with elastic properties. They are functionsU (x, ) of two spatial variables, namely the source variable x (i.e the location of an applied point force) and the field variable (i.e. where displacement is measured). Such fundamental displacements satisfy problems posed in terms of the Navier elastostatic operator N , defined for a generic displacement v by

N[v] := -div (C : ε[v]) (1.21)
(differential operators being understood in this work as acting on the field variable and defined in the weak sense of the theory of distributions). Fundamental solutions decay at infinity, i.e. verify U (x, ) → 0 asx → ∞ for a given x.

The integral operator M is defined in terms of the Mindlin fundamental tensor G = e k ⊗G k (created in B by a point force applied at x ∈ B, see [START_REF] Mindlin | Force at a Point in the Interior of a Semi-Infinite Solid[END_REF], which moreover satisfies the traction-free condition on ∂B necessary to provide the integral representation (1.12) of the displacement. However, M admits an alternative description based on the Boussinesq-Cerruti kernel (cf. Remark 3 for nomenclature details), which simplifies the treatment of the elastic equilibrium problem. The Boussinesq-Cerruti tensor B( x, ) = e k ⊗ B k ( x, ) satisfies for all ( x, 0) ∈ ∂B

N[B k ]( x, •) = 0 in B, T [B k ]( x, •) = δ x e k on ∂B, 1 ≤ k ≤ 3, (1.22)
where δ x is the Dirac distribution defined on ∂B and supported at ( x, 0) ∈ ∂B.

Lemma 2. The boundary integral operator M defined in terms of the Boussinesq-Cerruti kernel is written as:

M[p](q, x 3 ) = B T (q, x 3 ) • p(q).
(1.23)

Proof. [START_REF] Mindlin | Force at a Point in the Interior of a Semi-Infinite Solid[END_REF] has shown that the kernel G reduces to the Boussinesq-Cerruti fundamental tensor when x → ∂B, that is G(•, 0, 3 ) = B(•, 3 ). We make use of this property along with the property that U (x, ) = U T ( , x) for any fundamental tensor [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF]:

M[p] = G(•, x 3 , 0) p = ∫ ∂ B G( x -˜ , x 3 , 0) • p( ˜ ) d ˜ = ∫ ∂ B G T ( ˜ -x, 0, x 3 ) • p( ˜ ) d ˜ = ∫ ∂ B B T ( ˜ -x, x 3 ) • p( ˜ ) d ˜ .
Applying Lemma 1, we obtain the desired result (notice the change in sign in front of q due to the x ↔ ˜ swap).

Remark 3. The Boussinesq-Cerruti kernel cannot be considered a fundamental solution in the sense of Definition 1 as it solves a problem where the point force is introduced as a boundary condition and not as a second member in the local equilibrium equation, i.e. the point force has to lay on the surface an cannot be swapped for the field point in the domain.

As such, it does not satisfy the usual symmetries expected from fundamental solutions (e.g. we have B(q, 3 ) B T (-q, 3 ), whereas for the Mindlin tensor G(q, x 3 , 3 ) = G T (-q, 3 , x 3 )).

However, by abuse of terminology, we still refer to the Boussinesq-Cerruti displacement as a "fundamental solution" since it essentially fulfills the same role.

Since the Boussinesq-Cerruti tensor describes the deformation of a half-space subject to surface tractions, it is ubiquitous in elastic contact mechanics.

Partial Fourier solution

The first use of a spectral form of the operator M can probably be attributed to [START_REF] Westergaard | Bearing Pressures and Cracks[END_REF] who solved the line contact of wavy surfaces. Taking his contact solution to full contact gave the displacement due to a sinusoidal traction distribution on the surface. The generalization of this idea to surface contact is due to [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], and was first used for multi-scale rough surface contact by [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF].

Elastic rough contact: boundary integral methods

The use of the Westergaard solution in a numerical method is straightforward to conceive (and to execute), hence its popularity since the beginning of the century. Being defined for a periodic setting, the Westergaard solution is not exactly B in eq. (1.23), since B (and by extension B) does not describe a periodic medium. The first derivation of B was done by Amba-Rao (1969) 5 . We now show his derivation method, which we use as foundation for the volume integral operators introduced in Chapter 2 for the treatment of plastic flow in contact. The entire method is coded in a Jupyter Notebook [START_REF] Frérot | The Mindlin Fundamental Solution -A Fourier Approach[END_REF] and is free to (re-)use on binder.org.

In the partial-Fourier domain, the transformed isotropic Navier operator is expressed by Amba-Rao (1969) as:

N[u](q, 3 ) := N [ u](q, 3 ) = µ d 2 d 2 3 -q 2 I + c ∇ ⊗ ∇ u(q, 3 ), (1.24)
with c := (λ + µ)/µ = 1/(1 -2ν), ∇ := (iq 1 , iq 2 , d/d 3 ) and q 2 := q 2 = q 2 1 +q 2

2 . The process of deriving the Boussinesq-Cerruti fundamental tensor involves finding elements of ker( N ), the 6-dimensional space of functions u satisfying the ODE N [ u] = 0 , that fulfill the conditions given in eq. (1.22). [START_REF] Amba-Rao | Fourier Transform Methods in Elasticity Problems and an Application[END_REF] has derived a basis of ker( N ), defined as follows. Let the matrix-valued functions A + and A -be given by

A ± (q, 3 ) = e ∓q 3 I + c c + 2 q 3 ∆ ± ⊗ ∆ ± , (1.25)
where ∆ ± = e 3 ∓ ∆ and ∆ = (iq 1 /q, iq 2 /q, 0). Each column A ± j of A ± solves N [A ± j ] = 0 , so that we have ker( N ) = range A + + range (A -) .

(1.26) Moreover, the matrices A ± are invertible and linearly independent (see Appendix A.1): any element of ker( N ) therefore has six free coefficients, to be determined by additional conditions.

In the case of the Boussinesq-Cerruti fundamental solution, the coefficients are deduced from 1.27) which are the decaying and boundary conditions expressed in the partial-Fourier space. The former imposes that B k be solely expressed in terms of A + . Solving the latter for k = 1, 2, 3

lim 3 →∞ | B k (q, 3 )| = 0, T [ B k ] = e k , (
gives the solution for the tensor B: B(q, 3 ) = 1 q B 0,0 (q)e -q3 + B 1,0 (q)q 3 e -q 3 , (1.28) with

B 0,0 (q) = 1 2µ 2I + (1 -2ν)∆ + ⊗ ∆ -+ ∆ ⊗ ∆ -e 3 ⊗ e 3 , (1.29) B 1,0 (q) = 1 2µ ∆ + ⊗ ∆ + .
(1.30)

5 Actually B

When p( x) = A sin(ωx 1 )e 3 (a normal sine pressure), the above expression recovers the Westergaard solution in full contact. This can be verified with by taking the Fourier transform of the sine function (two Dirac distributions supported at ω and -ω), which when applied to B and subsequently transformed back to the physical space restores a sine wave of amplitude scaled by ω -1 . The application of the Dirac distributions to the continuous function B(q, •) yields discrete values of B that constitute a Fourier series representing a periodic function. This result is presented in the next section and discussed further in Chapter 2.

Horizontally-periodic setting

We have until now discussed the solution process for the equilibrium of an elastic half-space.

We have used the Fourier transform because of its ability to cast convolutions into function products. It also has the advantage of providing a straightforward discretization process by means of Fourier series. This requires the domain of interest to be periodic, and we therefore define its cell:

B p = - L 1 2 , L 1 2 × - L 2 2 , L 2 2 × R + , (1.31)
where L 1 (resp. L 2 ) is the horizontal dimension of the cell in the direction e 1 (resp. e 2 ). B pperiodic fields admit a representation as Fourier series. The surface traction is consequently expressed as

p( x) = k ∈Z 2 p(k) exp(2π i k • x), (1.32)
where k is the wave vector, ki = k i /L i and p(k) are the Fourier coefficients of the complex series. The resulting displacement also admits a representation as a Fourier series:

M[p]( x, x 3 ) = 1 4π 2 k ∈Z 2 B T (2π k, x 3 ) • p(k) exp(2π i k • x).
(1.33)

A more general form of this representation is discussed in Section 2.3. The machine computation of eq. ( 1.33) requires the truncation of the series in the set

Z 2 N = k ∈ Z 2 : - N 1 2 < k 1 < N 1 2 , - N 2 2 < k 2 < N 2 2 , (1.34)
where

N := (N 1 , N 2 )
is the number of surface points. The Fourier coefficients p(k) are approximated using a discrete Fourier transform values of p taken at finitely many points in the physical space. Values of p have to be represented in the physical space because contact constraints need to be enforced locally; a Fourier representation of p and cannot be used to enforce contact constraints because Fourier coefficients have a non-local nature.

Other integral equation methods

While the presented method has strong suits in both ease of implementation and computational efficiency, the periodicity requirements can be a drawback in some situations (e.g. rough ball on flat contact). Moreover, there are other fundamental solutions that one can use to represent the elastic response of a half-space. We present here some alternative methods.

Physical space methods

The literature on rough surface contact mechanics shows that the Boussinesq-Cerruti kernel can be employed directly in the physical space. Polonsky and Keer (1999b) discretize ∂B into a square grid where the pressure is represented by a function that is piecewise constant on each square element. Closed form solutions exist for the application of a normal pressure on a square patch [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], and the convolution product underlying M can be written as a matrix-vector product. They then employ a multi-level multi-summation technique to reduce the O(N 2 ) computational complexity to O(N ln 2 (N )), similar to the complexity of the fast-Fourier transform of O(N ln(N )). The advantage of this approach is the inherent nonperiodicity of the pressure (that now has a compact support) which makes the method ideal for non-nominally flat surface contact (e.g. rough ball on flat). Using a Fourier approach in such a case would require discretizing a domain larger than necessary. This cost can however be alleviated by a discrete FFT convolution technique [START_REF] Liu | A Versatile Method of Discrete Convolution and FFT (DC-FFT) for Contact Analyses[END_REF] that adapts the computational domain size to the physical contact size. This technique is often used with the pressure-patch solution defined in the physical space and transformed to the Fourier domain [START_REF] Zhou | A Mesh Differential Refinement Scheme for Solving Elastic Fields of Half-Space Inclusion Problems[END_REF][START_REF] Sahlin | A Mixed Lubrication Model Incorporating Measured Surface Topography. Part 1: Theory of Flow Factors[END_REF]). An adaptive non-Fourier method for normal contact has also been developed by [START_REF] Putignano | A New Efficient Numerical Method for Contact Mechanics of Rough Surfaces[END_REF]. The physical-space Boussinesq-Cerruti kernel was also applied to frictional contact [START_REF] Pohrt | Complete Boundary Element Formulation for Normal and Tangential Contact Problems[END_REF][START_REF] Li | A Semi-Analytical Approach to Three-Dimensional Normal Contact Problems with Friction[END_REF].

Outside the domain of rough surface contact, traditional boundary element methods can be applied to contact problems [START_REF] Man | Bem Frictional Contact Analysis: Load Incremental Technique[END_REF]. Classical formulations of BEM lead however to dense non-symmetric systems which are simply too expensive to solve for rough surface contact. However, the recent introduction of the multipole decomposition of the fundamental solutions (which allows a separation of the source and field variables in the integral equation) has led to the applicability of the multilevel fast multipole method (FMM) which drastically reduces the computational cost of BEM to O(N ) [START_REF] Liu | The Fast Multipole Boundary Element Method for Potential Problems: A Tutorial[END_REF]. The FMM has seen applications in wave propagation [START_REF] Chaillat | A Multi-Level Fast Multipole BEM for 3-D Elastodynamics in the Frequency Domain[END_REF][START_REF] Cao | A Fast Directional BEM for Large-Scale Acoustic Problems Based on the Burton-Miller Formulation[END_REF] but also in elastic contact mechanics [START_REF] Chen | The Fast Multipole Boundary Element Methods (FMBEM) and Its Applications in Rolling Engineering Analysis[END_REF], although it has to our knowledge never been applied to rough surface contact mechanics. The Boussinesq-Cerruti solution being derived from potential theory [START_REF] Love | A Treatise on the Mathematical Theory of Elasticity[END_REF], it possibly admits a multipole expansion. Another acceleration method for the BEM is the use of hierarchical matrices [START_REF] Chaillat | Theory and Implementation of H-Matrix Based Iterative and Direct Solvers for Helmholtz and Elastodynamic Oscillatory Kernels[END_REF][START_REF] Hodapp | Lattice Green Function Methods for Atomistic/Continuum Coupling: Theory and Data-Sparse Implementation[END_REF] which allow a representation of the dense matrix where the low-rank blocks are compressed, greatly decreasing storage and computation complexities. It also has the added advantage of being usable as a black box, whereas the FMM requires tedious problem-dependent developments.

The adaption of modern BEM acceleration techniques such as FMM and H -matrices is out of the scope of this thesis, but given the amount of literature on their applicability and efficiency in other areas of solid mechanics, they seem to give a promising and low-risk course of action in the development of rough surface contact numerical methods, especially in elastodynamics, where such an approach of contact has the potential to provide insight into, for example, earthquake nucleation [START_REF] Barras | Interplay between Process Zone and Material Heterogeneities for Dynamic Cracks[END_REF] and propagation [START_REF] Chaillat | A New Fast Multi-Domain BEM to Model Seismic Wave Propagation and Amplification in 3-D Geological Structures[END_REF]).

Green's function molecular dynamics

Although ubiquitous in contact mechanics, the Boussinesq-Cerruti kernel is not the only kernel that can be used to represent an elastic solid. More complex constitutive relations require different kernels [START_REF] Wang | Three-Dimensional Contact Analysis with Couple Stress Elasticity[END_REF], and other approaches for the derivation of kernels for semi-infinite solids can be taken. One such approach is the molecular dynamics based method of [START_REF] Campañá | Practical Green's Function Approach to the Simulation of Elastic Semi-Infinite Solids[END_REF]. They consider a harmonic lattice of n atoms, whose degrees of freedom U ∈ R 3×n are coupled with a stiffness matrix K ∈ R 3n×3n . Using statistical physics, one can show that at a given temperature T the stiffness can be obtained from the statistical correlation of the degrees of freedom:

G := U ⊗ U = k B T K -1 .
A Schur complement of K can be expressed to represent the full lattice response only as a function of the surface degrees of freedom. This is however not used in practice. The methodology to compute the Schur complement described by [START_REF] Campañá | Practical Green's Function Approach to the Simulation of Elastic Semi-Infinite Solids[END_REF] consists in evaluating the correlation matrix of the displacements in the Fourier domain G(q) := U(q) ⊗ U(q) on an atomic system of small size (i.e. of the order of the lattice spacing), use it to compute the block-diagonal form of the Schur complement in the Fourier domain, which is proportional to G -1 , and extrapolate the results to large system sizes. This procedure has the appeal of being able to formulate Green's functions for different lattice structures with inherent anisotropy and has the added advantage, given translational invariance, of computational attractiveness due to the use of fast-Fourier transforms. It is also well suited to represent infinite-body boundary conditions for atomistic simulation systems [START_REF] Pastewka | Seamless Elastic Boundaries for Atomistic Calculations[END_REF] and, like their continuous counterpart, to compute elastic contact solutions [START_REF] Müser | Meeting the Contact-Mechanics Challenge[END_REF].

Normal Contact

With the tools developed in the previous section, we are able to efficiently and accurately solve the equilibrium of an elastic half-space subject to Neumann boundary conditions. The contact conditions can now be introduced. We consider situations where the boundary ∂B is in contact with a surface S, which is defined as the graph of a scalar function h ∈ C 0 (∂B) and can be viewed as the boundary of an infinitely stiff solid. The resulting normal force bringing B and S together is denoted W := W e 3 . We ignore friction and adhesion between B and S and define the gap function as:

[u] := u • e 3 ∂ B -h, (1.35)
1. Elastic rough contact: boundary integral methods which is the separation between the deformed solid B and the surface S along n. The boundary conditions on ∂B are hence expressed with the Hertz-Signorini-Moreau conditions:

[u] ≥ 0,

(1.36a) p[u] := T [u] • e 3 ≥ 0, (1.36b) [u]p[u] = 0.
(1.36c)

Moreover we treat the contacting surface S as fixed, and enforce the applied resulting force W as an additional constraint that will allow to determine the mean vertical displacement:

∫ ∂ B p[u] dA = W .
(1.37)

The above contact conditions (1.36) are known as unilateral boundary conditions, as they do not behave symmetrically with respect to the sign of u6 . This makes the contact problem nonlinear even if the bodies in contact behave elastically. One can think of the contact problem as a solid mechanics problem with zones of the boundary where the displacement is known because the gap is null and the bodies touch (the pressure is unknown) and zones where the boundary is free of tractions because there is no contact (the gap is unknown). What condition (1.36c) expresses is that these two zones form a partition of ∂B: points where the gap is strictly positive have p = 0 and are disjoint from points that have a strictly positive contact pressure (those have = 0 ). The topography of this partition is however a priori unknown, hence the non-linearity of the problem.

Variational form

We have introduced earlier in this chapter the formulation of the general elastostatic problem as the minimization of the total potential energy. To introduce the contact constraints, one needs to define the space of kinematically admissible displacement fields KA(B). Apart from the usual regularity hypothesis (u ∈ H 1 ), one needs to carefully define KA because of:

• periodicity: analysis of the kernel B shows that for k = 0 the displacement is undefined. This is due to the periodic nature of the problem. A common workaround is to remove the constant displacement mode, as is discussed in Chapter 2;

• contact: an admissible displacement must satisfy the constraint on the gap.

The first consideration defines the function space

H 1 (B p ; R 3 ) = u + ū u ∈ H 1 (B p ; R 3 ), ū ∈ R 3 with u (0 , •) = 0 , (1.38)
which is a space of B p -periodic functions whose fundamental mode (i.e. the horizontal average) is constant with respect to x 3 . The space of admissible solutions to the contact problem is:

KA(B p ) := u ∈ H 1 (B p ; R 3 ) : [u] ≥ 0 .
(1.39)

The minimization principle is then written as:

inf u ∈KA(B p ) 1 2 ∫ B p σ [u] : ε[u] dV , (1.40)
which is a quadratic program with unilateral constraints. [START_REF] Kalker | Variational Principles of Contact Elastostatics[END_REF] has shown that the above problem can be expressed only in terms of boundary integrals, provided u additionally satisfies the local equilibrium eq. ( 1.6):

1 2 ∫ B p σ [u] : ε[u] dV = 1 2 ∫ B p σ [u] : ∇u dV = 1 2 ∫ B p (div(σ [u] • u) -div(σ [u]) • u) dV = 1 2 ∫ ∂ B p u • σ [u] • n dV , so that the contact problem is equivalent to inf u ∈KA(B p ) I u (u) = 1 2 ∫ ∂ B p T [u] • u dS . (1.41)
The above formulation is well suited to boundary integral methods because the unknown field u need only be resolved on the surface. Moreover, the integral representation (1.12) produces by construction displacement fields that satisfy the local equilibrium equation. This minimization problem can now be solved by quadratic programming (see Rey, Anciaux, and Molinari, 2017, for adhesive contact), but optimization techniques are usually applied to the dual of (1.41), which we derive by associating a Lagrange multiplier λ to the non-penetration constraint eq. (1.36a). The Lagrangian of the problem is:

L(u, λ) = 1 2 ∫ ∂ B p T [u] • u dS - ∫ ∂ B p λ(u • e 3 -h) dS. (1.42)
The Karush-Kuhn-Tucker optimality conditions imply stationarity of L with respect to u, so that we can express λ:

∂L ∂u = 0 ⇔ T [u] -λe 3 = 0 ⇔ λ = T [u] • e 3 ⇔ u = M[λe 3 ].
λ can be physically interpreted as the normal tractions acting on the surface. Replacing u = M[pe 3 ] in L (effectively computing the Legendre transform of I u ) leads to the quadratic program 1.43) where the primary unknown p is the normal surface traction and SA is the space of statically admissible tractions, i.e. tractions that satisfy p ≥ 0 (for non-adhesive contact) and

inf pe 3 ∈SA(B p ) I p (p) = 1 2 ∫ ∂ B p pe 3 • M[pe 3 ] dS - ∫ ∂ B p ph dS , (
∫ ∂ B p dS = W .

Solution methods

Most methods of (constrained) quadratic programming involve gradient descents [START_REF] Boyd | Convex Optimization[END_REF] with additional projection and/or constraint enforcement steps [START_REF] Bemporad | Optimization Algorithms for the Solution of the Frictionless Normal Contact between Rough Surfaces[END_REF]. The gradient computation of the functional I p is straightforward with the use of M as defined previously:

∇I p (p) = M[pe 3 ] • e 3 -h = u 3 -h = .
(1.44)

We will here provide an overview of two descent methods used in rough surface contact mechanics and discuss other methods found in the literature.

Steepest descent

Perhaps the simplest approach to solving (1.43) is the algorithm proposed by [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF]. It consists in taking steps in the -∇I p = -direction then shifting and truncating the normal traction so that both the unilateral and equilibrium constraints are satisfied. It is summarized in Algorithm 1. While only needing one application of M per loop, its convergence rate is sub-optimal. It however has the advantage of being easily adapted to pseudo-plastic methods such as the surface hardness phenomenological model proposed by [START_REF] Almqvist | On the Dry Elasto-Plastic Contact of Nominally Flat Surfaces[END_REF].

Algorithm [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF] steepest descent algorithm.

Input: normal total load W , surface profile H, tolerance ε, max. number of iterations N max P ← W /|∂B p | Average constant load initial guess h norm ← ||H|| k ← 0 repeat G ← M[Pe 3 ] • e 3 -H Gradient P ← P -G Find α 0 solution of (P + α 0 ) + -W = 0 (•) + is the ramp function P ← (P + α 0 ) + e ← P • (G -min(G))/(W h norm )
Error on complementarity

k ← k + 1 until e < ε or k > N max G ← G -min(G)
Ensure positive gap

Constrained conjugate gradient

The method proposed by Polonsky and Keer (1999b) is a variation of a conjugate gradient algorithm with an active set to accelerate the convergence. As expected from numerical linear algebra, it outperforms the steepest descent. Rey, Anciaux, and Molinari (2017) have shown that it can also be adapted to the primal problem where the unknown is the gap field in order to solve an adhesive contact problem7 . The primal and dual version of the constrained conjugate gradient are similar enough that one implementation can be done for both. In Algorithm 2, the main unknown is P which is the pressure, but can be swapped for the gap.

The active set strategy is relying on the definition of the set S = { x ∈ ∂B p : p > 0} which is the current set of points in contact, and the set R = { x ∈ ∂B p : p = 0 and < 0} which is the set of inadmissible points. The latter should be empty at the end of the algorithm.

Algorithm 2 Polonsky and Keer (1999b) constrained conjugate gradient algorithm.

Input: normal total load W , surface profile H, tolerance ε, max. number of iterations

N max P ← W /|∂B p | Average constant load initial guess T ← 0 Search direction h norm ← ||H|| G ← 0 , G old ← 1 δ ← 0 k ← 0 repeat S ← where(P > 0) Set of current points in contact G ← ∇I Gradient G ← G -mean(G S )
Centering on contact zone

G ← ||G S || 2 T S ← G S + δ G G old T S Update only in current contact zone G old ← G R ← M[T] R ← R -mean(R S ) τ ← G S • T S R S • T S Critical step size P ← (P -τ T) + R ← where(P = 0 and G < 0) Set of inadmissible points if R = ∅ then δ ← 1 else δ ← 0 end if P R ← P R -τ G R Positive pressure on inadmissible points P ← W P P Applied force constraint e ← P • (G -min(G))/(W h norm ) k ← k + 1 until e < ε ∨ k > N max G ← G -min(G)
Positive gap

Elastic rough contact: boundary integral methods

Other approaches [START_REF] Bemporad | Optimization Algorithms for the Solution of the Frictionless Normal Contact between Rough Surfaces[END_REF] have reviewed Algorithms 1 and 2 as well as adapted to algorithms typically used in quadratic programming: the non-negative least squares and the alternating direction method of multipliers. The latter uses an augmented Lagrangian formulation of problem (1.43), i.e. the unilateral constraint on p is incorporated in the Lagrangian with a Lagrange multiplier and a penalty term. This has the disadvantage of adding unknowns to the system, whereas the former (and the approaches we have previously presented) does not, as the Lagrange multiplier is known with the KKT conditions. The non-negative least squares algorithm proposed by [START_REF] Bemporad | Optimization Algorithms for the Solution of the Frictionless Normal Contact between Rough Surfaces[END_REF] seems to perform better than the constrained conjugate gradient of Polonsky and Keer (1999b), but it requires to solve a reduced linear system of the form M S [P S ] = U S whose number of unknowns P S is the number of points in the set S. This operation must be solved with an iterative solver, as the operator M S , which is a restriction of M on S cannot be readily inverted. The advantage of this algorithm is its performance when the number of points in contact is small. Like Algorithm 2, the active set strategy favors small contact situations when the dual formulation is employed.

Dynamic relaxation schemes have also been used to solve rough contact problems, both in a finite elements setting [START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF][START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF] and GFMD setting [START_REF] Campañá | Elastic Contact between Self-Affine Surfaces: Comparison of Numerical Stress and Contact Correlation Functions with Analytic Predictions[END_REF][START_REF] Zhou | Solution of Boundary-Element Problems Using the Fast-Inertial-Relaxation-Engine Method[END_REF]. While they are relatively simple, one needs to appropriately choose the fictional mass and damping parameters for optimal convergence, a process which is often problem dependent. However, like Algorithm 2, they are flexible enough to be adapted to non-convex problems like adhesive contact [START_REF] Müser | Meeting the Contact-Mechanics Challenge[END_REF].

Approaches based on a penalty or barrier term can also be employed. These generally tend to show a performance independent of the contact area, at the expense of approximately enforcing the optimality conditions [START_REF] Boyd | Convex Optimization[END_REF]). We will consider such an approach in Section 3.2 for elastic-plastic rough surface contact.

Rough surfaces and statistics

All solids have rough surfaces at a given range of scales. These vary dramatically depending on the system (e.g. geological faults scales compared to watch parts scales), but most natural (and some man-made) surfaces exhibit a self-affine scaling behavior over their characteristic range of scales8 [START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF]. The accurate representation of this scaling over as many length-scales as possible, as well as the associated contact response, is the main computational challenge of rough contact numerical methods. In this aspect, boundary integral methods outperform finite elements methods because of dimensionality reduction. Regardless of the method, using real surface profiles is important to obtain a realistic response of a contact interface. These are however difficult to obtain, and we turn to stochastic generation processes to obtain synthetic surface profiles with representative properties. The natural characteristic that artificial rough surfaces should reproduce is the self-affine scaling behavior. It is often expressed through the power-spectrum density ϕ of the surface (PSD):

ϕ(q) := |F [h](q)| 2 . (1.45)
To obtain a self-affine behavior, ϕ must scale as a power-law of its arguments (e.g. for an isotropic surface ϕ(q) ∼ q -τ ). The scaling is generally described with the Hurst exponent [START_REF] Meakin | Fractals, Scaling, and Growth Far from Equilibrium[END_REF], which is best defined with the height-height correlation function S(δr ) which characterizes changes in h over a distance δr smaller than the autocorrelation distance, .46) Intuitively, this shows that if one were to stand at a point x of the surface and move δr in any direction, the typical height variation would be of the order δr H (for an isotropic surface). This is equivalent to the fact that scaling space by a factor λ scales the surface heights by a factor λ H . Given a power-law behavior of the PSD ϕ(q) ∼ q -τ , τ can be expressed in terms of the Hurst exponent as τ = 2(H + 1) [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF]. Figure 1.1 shows an example of a self-affine rough surface.

S(δr

) := h(• + δx) |δ x |=δ r -h(•) 2 ∼ δr H . ( 1

Rough surface generation

Rough surfaces are often modeled as Gaussian random fields whose covariance encode the scaling properties previously mentioned [START_REF] Longuet-Higgins | Statistical Properties of an Isotropic Random Surface[END_REF][START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF][START_REF] Rey | Quantifying Uncertainties in Contact Mechanics of Rough Surfaces Using the Multilevel Monte Carlo Method[END_REF]. This is a simplifying assumption, as some classes of surfaces, such as those worn by a grinding process, do not follow a Gaussian distributions [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF]. We show here two common approaches for the generation of Gaussian fields with self-affine properties and finally discuss more complex approaches applicable for non-Gaussian surfaces.

Fourier filtering

The approach of [START_REF] Hu | Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis[END_REF] consists in using the fast-Fourier transform to construct a random field from a prescribed PSD ϕ(q) and a non-correlated Gaussian field. The procedure is summarized in Algorithm 3. The randomness introduced via the non-correlated field is transmitted to the PSD information by taking the discrete Fourier transform of the former and multiplying the square root of the latter. This has the effect that the PSD of the discrete surface H computed a posteriori includes spurious noise.

Algorithm 3 Filtering algorithm [START_REF] Hu | Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis[END_REF].

Input: surface sizes N := (N 1 , N 2 ), PSD ϕ : R 2 → R X ← N (0, I ) Gaussian uncorrelated vector S ← ϕ(k), with k ∈ Z 2 N Discrete PSD H ← FFT[X] H ← H • √ S
Point-wise product and square root

H ← | N | • FFT -1 [ H]
Output: H discrete surface heights Random phase [START_REF] Wu | Simulation of Rough Surfaces with FFT[END_REF] proposed a method that removes the randomness in the effective PSD of a rough surface. It is simply motivated by the fact that the PSD of an uncorrelated Gaussian field should be constant. The randomness is therefore introduced with a complex field of unit modulus and uniformly distributed phase. Since the phase disappears in the a posteriori computation of the PSD, the discrete PSD is recovered. The procedure is given in Algorithm 4.

Algorithm 4 Random phase algorithm [START_REF] Wu | Simulation of Rough Surfaces with FFT[END_REF].

Input: surface sizes N := (N 1 , N 2 ), PSD ϕ : R 2 → R X ← U(-π, π ) Uniformly distributed phases S ← ϕ(k), with k ∈ Z 2 N Discrete PSD H ← exp(Xi) H ← H • √ S
Point-wise product and square root

H ← | N | • FFT -1 [ H] Output: H discrete surface heights

Rough surfaces and statistics

Non-Gaussian surface generation

Rough surfaces can be auto-correlated non-Gaussian fields. This often happens in processed surfaces [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF], e.g. sanded surface that have worn out asperities. They may however conserve their scaling properties. In general, rough surfaces can be characterized by a probability density function and an auto-correlation function. In wear processes, both these functions may evolve as third body particles are created and work the surface [START_REF] Milanese | Emergence of Self-Affine Surfaces during Adhesive Wear[END_REF]. Modifications to the distribution of heights may induce changes to the frictional properties of the surface [START_REF] Chilamakuri | Contact Analysis of Non-Gaussian Random Surfaces[END_REF].

Generating non-Gaussian surfaces often involves the iterative solution of a linear system. Wu (2004) proposes a method based on Algorithm 4 where the distribution of phases is determined iteratively so that the generated surface matches the desired height distribution. The distribution of phases is controlled by a Johnson translator system, which is a way to construct an non-Gaussian random field with controlled skewness and kurtosis. [START_REF] Wang | A Simulation Method for Non-Gaussian Rough Surfaces Using Fast Fourier Transform and Translation Process Theory[END_REF] propose a different method that also uses a translation process to generate a non-Gaussian surface with a given distribution and auto-correlation from a gaussian field with a to-be-determined auto-correlation.

Although non-Gaussian rough surface may play a significant role in both friction and wear systems, they are not used in this thesis, and only discussed as reference for the reader.

Surface statistics

Some statistical properties of rough surfaces have measurable impact on contact-related quantities. For example, the density of surface summits governs the number of micro-contacts [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF], while the spectrum bandwidth and root-mean-square of slopes influence the load/contact area relationship (Yastrebov, Anciaux, and Molinari, 2017b;[START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF]. [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF] has shown that under the assumption of a Gaussian field, most surface measures can be derived from the PSD and its moments:

m st = ∫ R 2
q s 1 q t 2 ϕ(q 1 , q 2 ) dq 1 dq 2 .

(1.47)

We now compare the previously introduced generation algorithm on the root-mean-square of heights and root-mean-square of slopes:

h rms = |h| 2 = √ m 00 ,
(1.48)

h rms = |∇h| 2 = √ 2m 20 .
(1.49) Figure 1.2 shows the theoretical value for the two measures and histograms for 100 (1024×1024) realizations of each algorithm, with a reference spectrum given by: [START_REF] Wu | Simulation of Rough Surfaces with FFT[END_REF] gives consistent values around the analytical prediction, while the algorithm of [START_REF] Hu | Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis[END_REF] gives a much larger spread (approx. 10 times). Figure (b) shows the root-mean-square of slopes. Here the random phase algorithm does not reproduce the analytical value. This is due to the discrepancy between the continuous and discrete version of the PSD which is increased by the moment computation of eq. (1.47). This error decreases as the ratio L/λ l increases.

ϕ(q) =            C q l ≤ q < q r , C q q r -2(H +1) q r ≤ q < q s 0 otherwise, , (1.50
where C is typically set to one and (q l , q r , q s ) are respectively the long cutoff, roll-off and short cutoff spatial frequencies associated to the wavelengths (λ l , λ r , λ s ). In the case of fig. 1.2, q l = 4π , q r = 8π and q s = 256π , with a system size of L = 1 .

Statistics on surface can be computed in one of two ways: (a) geometric estimation, (b) spectral estimation. The former uses the generated height profile H directly to compute quantities of interest, while the latter computes them using the discrete PSD S. Depending on the quantity of interest, a discrepancy can arise between the two estimators. Figure 1.3 compares the estimated h rms using both techniques (with a finite differences approximation for the geometric estimator). One can see a significant bias of the geometric estimator for the lower surface discretization while the spectral estimator remains stable.

In linear elastic contact of rough surfaces, the slope of the load/contact area curve is proportional to h rms [START_REF] Bush | The Elastic Contact of a Rough Surface[END_REF][START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF]. It is therefore used in normalizing the applied load. In order to get an unbiased normalization and a correct estimate of the load/contact area slope, we use in the rest of this thesis the average value spectrally computed over a statistically significant number of surface realizations.

Periodicity and statistical significance

Algorithms 3 and 4 naturally construct periodic surfaces due to their use of the discrete Fourier transform. This conforms to the mathematical framework presented in Section 1.2, so the periodicity of the surface corresponds to the periodicity of the deformable half-space. One must however consider periodicity with respect to the physical problem at hand. The second part of this thesis is concerned with the up-scaling of microscopic wear mechanisms to the macro-scale. In this context, periodicity allows the reduction of finite boundary effects and approximates the response of an infinitely large system, provided the surface is representative. The latter issue has been studied by [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF], who have shown that the ratio of the periodic cell size L to the roll-off wavelength λ r of the surface spectrum controls the representativity in the sense of a representative volume element (commonly used in homogenization). An alternative criterion is that the autocorrelation length of the rough surface should be small compared to the period size. Representativity is however a non-objective notion, as it depends on the quantity that one wishes to evaluate with statistical accuracy (i.e. low variance). For example, Yastrebov, Anciaux, and Molinari (2012) report that L/λ r = 4 is acceptable for the evaluation of the true contact area, whereas in Chapter 4 we report a value of L/λ r = 8 to have a statistically converged distribution of micro-contact areas. This is discussed in further details in Chapter 4 and Appendix C.1. Some problems are sensitive to periodicity errors: for example, the contact of rough spheres may be polluted by periodicity [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF], or a contact moment can be applied in addition to a force. In those cases, the periodic solution can still be used with an appropriate padding of the computational domain. This entails larger computations, which may unfavorably compare to non-periodic simulation methods like the fast multi-pole approach.

Elastic rough contact: boundary integral methods

Case-by-case analysis should then be conducted to determine the most appropriate numerical method. However, because the physical problems considered in this thesis lend themselves to periodicity, we shall not undertake such analysis. A lthough the methods presented in last chapter for elastic rough surface contact are now mature, elastic theories have difficulties representing realistic contact behavior: local contact pressures can easily reach values higher than the Young's modulus [START_REF] Müser | Internal, Elastic Stresses below Randomly Rough Contacts[END_REF], and the true contact area vanishes as the rough surface spectrum grows wider (Persson, 2006). The works of [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF], [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF], and more recently [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF] have experimentally shown that elastic theories are not able to model realistic contact interfaces. It is the consensus of the tribology community [START_REF] Vakis | Modeling and Simulation in Tribology across Scales: An Overview[END_REF] that contact models need to evolve to include material non-linearities such as plasticity.

The work of [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF], while pioneering the study of elastic-plastic rough surface contact (using a finite element approach), suffers from discretization error (a single element is used to reproduce the smallest surface wavelength), inaccuracy of an explicit dynamic relaxation scheme to reduce the simulation cost of a static calculation, as well as statistical errors from the surface spectrum choice [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF] and number of realizations. Pei et al. nonetheless confirmed the key role of plasticity in quantifying the true contact area and the micro-contact distribution. Consequently, the objective of this chapter is to propose a high-performance, robust and quantitatively accurate method to study the contact behavior of elastic-plastic materials with statistically representative rough surfaces.

Our method is based on a volume integral approach [START_REF] Telles | On the Application of the Boundary Element Method to Plasticity[END_REF][START_REF] Telles | Implicit Procedures for the Solution of Elastoplastic Problems by the Boundary Element Method[END_REF][START_REF] Bonnet | Implicit BEM Formulations for Usual and Sensitivity Problems in Elasto-Plasticity Using the Consistent Tangent Operator Concept[END_REF]. Similarly to the boundary methods previously mentioned, the volume integral methods (VIMs) can exactly represent the elastic behavior of a semi-infinite solid. This limits the volume discretization to the potentially plastic regions, allowing better usage of computational resources. VIMs however require knowledge of singular fundamental solutions and the computation of their volume convolution with plastic deformations. This operation is costly with a "naive" implementation [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF], its algorithmic complexity being O(N 2 ) (N being the total number of plastic deformation unknowns). It has been accelerated using 2D fast-Fourier transform (FFT) of the fundamental solutions [START_REF] Sainsot | A Numerical Model for Elastoplastic Rough Contact[END_REF][START_REF] Wang | Numerical Simulation for Three Dimensional Elastic-Plastic Contact with Hardening Behavior[END_REF], at the cost of introducing a sampling error in the computed quantities of interest. Another approach for accelerating the convolution computation is to use a 3D-FFT [START_REF] Wang | A New Fast Method for Solving Contact Plasticity and Its Application in Analyzing Elasto-Plastic Partial Slip[END_REF], but this introduces a periodicity error and requires discretization of a domain more than twice the volume of the expected plastic zone to reduce these effects, thereby reducing the attractivity of VIMs for semi-infinite modeling.

Accordingly, we present a volume integral method well suited for periodic elastic-plastic contact problems that is based on a novel derivation of the half-space fundamental solutions directly in the 2D partial Fourier domain. The advantage is three-fold: it allows the use of the 2D-FFT to speed up the convolution computation in the plane parallel to the contact surface without introducing sampling error, it does not require storage of the discrete Fourier coefficients of the real-space fundamental solutions, and permits optimization of the convolution computation by exploiting the structure of the analytical solutions, while keeping the advantages described above for volume integral methods. Indeed, the algorithmic complexity of our treatment is only O(N log(N )) per convolution. Moreover, the use of the FFT renders the method matrix-free, making it attractive for use in conjunction with iterative solvers.

x 1 x 3 x 2

Problem statement and overview of solution methodology

Let B := { ∈ R 3 : • e 3 ≥ 0} be a deformable semi-infinite elastic-plastic solid of boundary ∂B, see fig. 2.1. Points ∈ B will often be denoted as = ( ˜ , 3 ) with ˜ := ( 1 , 2 ). Let u be the displacement vector field of B. The linearized strain tensor ε and the Cauchy stress tensor σ are respectively given by:

ε[u] := 1 2 ∇u + ∇u T , (2.1) σ [u, ε p ] := C : ε[u] -ε p , (2.2)
where C ∈ R 3×3×3×3 is the elasticity tensor, satisfying the usual ellipticity and (major and minor) symmetries, while ε p is the plastic strain. Here and thereafter, we follow the usual convention whereby the gradient operator ∇ acts "to the right", so that e.g. (∇u) i j = ∂ j u i = u i,j . The stress σ satisfies the conservation of momentum equation in the absence of body forces:

div σ = 0 a.e. in B.
(2.3)

Let n := -e 3 be the external normal of ∂B, we define the traction vector of the displacement field u as:

T [u, ε p ] := σ [u, ε p ] ∂ B
• n.

(2.4)

The evolution of the plastic strain is, for definiteness, assumed to obey a standard 2 plasticity model with an associated flow rule (see Section 2.4.1 for more details).

Horizontally-periodic setting

As we will employ Fourier methods to solve the elastic-plastic contact problem, we set the latter in a more natural setting for the application of the discrete Fourier transform. Let us define the periodic cell

B p = - L 1 2 , L 1 2 × - L 2 2 , L 2 2 × R + , (2.5)
where L 1 (resp. L 2 ) is the horizontal dimension of the cell in the direction e 1 (resp. e 2 ). It will become apparent in Section 2.3 that this helps defining the discretization procedure of the continuous operators presented below, since B p -periodic fields admit a representation as Fourier series.

Volume integral representation

Although integral equation methods are based on a linear superposition of states, they can still be applied to materially and geometrically non-linear problems (e.g. [START_REF] Telles | On the Application of the Boundary Element Method to Plasticity[END_REF] for plastic behavior, [START_REF] Jin | Boundary Element Formulation in Finite Deformation Plasticity Using Implicit Integration[END_REF] for large deformations). To this effect, eq. ( 2.3) can be rewritten in the form:

div(C : ε[u]) -div(C : ε p ) = 0 ⇔ N[u] = -div(C : ε p ), (2.6)
which can be solved with eq. ( 1.12) from Chapter 1 [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF]:

u(x) = ∫ ∂ B G(x, ) • T [u, 0 ]( ) dS - ∫ B G(x, ) • div(C : ε p )( ) dV = M[T [u, 0 ]] - ∫ ∂ B G(x, ) • (C : ε p ( )) • n dS - ∫ B ∇G(x, ) : C : ε p ( ) dV = M[T [u, ε p ]] + N [C : ε p ].
(2.7)

The volume integral operator N is defined for an eigenstress distribution w in B as:

N [w](x) := ∫ ∞ 0 ∫ R 2 ∇G( ˜ -x, x 3 , 3 ) : w( ) d ˜ d 3 = ∫ ∞ 0 (∇G w)( x, x 3 , 3 ) d 3 .
(2.8)

Much like the developments presented in Chapter 1, the translational invariance of the problem allows the application of the partial Fourier transform as described in Lemma 1, provided the functions involved satisfy the convolvability requirements of Remark 11 . However, unlike in Chapter 1, the Mindlin kernel gradient ∇G cannot be replaced by the Boussinesq-Cerruti kernel. Ergo any application of N in the partial Fourier domain necessitates the computation of ∇G.

Overview of solution methodology

In this thesis, we use a volume integral equation approach for solving the elastic-plastic contact problem. We make use of the fact established above that the elastostatic displacement u generated in B or on ∂B by any given traction distribution p on ∂B and eigenstress distribution w in B has the explicit representation u = M[p] + N [w], where the integral operators M, N are defined in terms of the Boussinesq-Cerruti and the Mindlin fundamental solutions respectively. This framework accounts for the semi-infinite geometry of the problem, and enforces automatic satisfaction of equilibrium and elastic constitutive relations. As a result, the solution of contact and elastic-plastic problems only requires the satisfaction of relevant relations on the contact surface (using evaluations of eq. ( 2.7) at surface points) and plastic regions (exploiting ∇u at internal points), respectively, with the latter task entailing computations of ∇M[p] and ∇N [w] for given p or w.

Our emphasis on convolvability restrictions underpinning eq. ( 2.8) (and eq. ( 1.20)) stems from the fact that the kernel G(z, x 3 , 3 ) of interest in this chapter, given by elastostatic fundamental solutions, is not in L 1 (R 2 ) as a function of z due to insufficiently fast decay at infinity; moreover, the kernel of ∇N also has a non-integrable singularity at z = 0 , and hence is not a locally integrable function. As a result, the convolution integral of eq. ( 2.8) applied to (G w)( x) is not always defined in the classical sense, even for "nice" (smooth, compactlysupported) densities w, while the properties of elastostatic fundamental solutions make it well defined as a convolution between distributions for any w having compact (i.e. spatially bounded) support, ensuring the validity of Lemma 1. The same features of G cause its partial Fourier transform to be well defined as a distribution but not in the classical sense of the integral in eq. (1.18). Finally, the very fact that fundamental solutions are fields created by singular loads (point forces) entails treating them, and their governing equations, in the sense of distributions. Summarizing, the framework of distribution theory is necessary for the main components of our proposed treatment to have clear meaning and validity.

Finally, our aim is to formulate and solve elastic-plastic contact problems where all fields have horizontal periodicity, whereas integral operators appearing in eq. ( 2.7) are a priori defined in a non-periodic setting. Indeed, the meaning of the relevant convolutions becomes a priori unclear for periodic (hence not compactly supported) densities. However, the evaluation at suitable discrete values of q of Fourier transforms of non-periodic convolutions will be found (see Theorem 3) to provide the required relations between a periodic displacement u and the periodic density w from which u emerges.

Displacement and displacement gradient

In this section, we discuss the different integral operators used and their partial Fourier representation.

Fundamental problems

The integral operators introduced in eq. ( 2.7) are expressed in terms of fundamental solutions as defined in Definition 1. We present here the fundamental problems of Kelvin and Mindlin, and their respective solutions needed for the integral operators. The convolutions using those kernels will be computed in the partial Fourier space using Lemma 1.

The Kelvin problem

The displacement caused by a point force in an infinite medium is at the heart of integral equation methods in solid mechanics. Here, the Kelvin tensor

U ∞ = e k ⊗ U k ∞ is used as a stepping stone to the Mindlin fundamental solution. The displacement U k ∞ satisfies N[U k ∞ ](x, •) = e k δ x in R 3 , for any x ∈ R 3 , 1 ≤ k ≤ 3 , (2.9)
where δ x is the three-dimensional Dirac distribution on R 3 supported at x, and decays at infinity. Due to its singular right-hand side, eq. ( 2.9) must be understood in the sense of distributions 2 , and the same will apply implicitly to the other fundamental solutions and their governing equations. The Kelvin fundamental solution possesses the full-space translational symmetry, that is U ∞ (x, ) = U ∞ (0 , -x), as well as the implied property for its gradients that

∇ x U ∞ (x, ) = -∇U ∞ (x, ). We now define the operator N ∞ N ∞ [w](x) := ∫ B ∇U ∞ (0 , -x) : w( ) d = ∫ ∞ 0 ∇U ∞ w ( x, x 3 , 3 ) d 3 .
(2.10)

The field N ∞ [w] decays at infinity and can readily be shown to satisfy

N N ∞ [w] = -divw in R 3 (2.11)
i.e. it is the elastostatic displacement created in an unbounded medium by an eigenstress distribution w.

The Mindlin problem

The problem of a point force in a semi-infinite isotropic elastic medium with a free surface was solved by [START_REF] Mindlin | Force at a Point in the Interior of a Semi-Infinite Solid[END_REF]. This fundamental solution allows to express the displacement created in B by an eigenstress distribution w in B and satisfying a traction-free condition on ∂B. The latter feature makes it very attractive for contact problems since it removes the need to solve an implicit boundary integral equation, which is a staple of conventional boundaryelement methods (Bonnet, 1995 and references therein). The Mindlin tensor x) after integrations by parts and recalling the self-adjointness of the Navier operator, for any test function ϕ ∈ C ∞ 0 (R 3 ; R 3 ) (with f , ϕ R 3 denoting the duality product). Similar interpretations implicitly apply for the other fundamental solutions.

G = e k ⊗ G k satisfies: N[G k ](x, •) = δ x in B, T [G k ](x, •) = 0 on ∂B. (2.12) 2 Equation (2.9) therefore means that U k ∞ must verify N[U k ∞ ](x, •), ϕ R 3 = e k δ x , ϕ R 3 , i.e. U k ∞ (x, •), N[ϕ] R 3 = ϕ k (
The operator N defined by

N [w](x) := ∫ ∞ 0 ∇G w ( x, x 3 , 3 ) d 3 .
(2.13) then yields the displacement created in B by an eigenstress distribution w. It is easy to see that, by virtue of linear superposition, G is given in terms of U ∞ and B (the Boussinesq-Cerruti tensor) by

G k (q, x 3 , 3 ) = U k ∞ (q, x 3 , 3 ) -B T (q, x 3 ) • T [ U k ∞ (q, x 3 , 3 )],
(2.14)

with the second term in eq. ( 2.14), expressed using Lemma 2, is canceling the traction vector on ∂B. Consequently, the operators M and N can be readily evaluated in the partial Fourier space once B and U ∞ are known.

Partial Fourier space solutions

The main novelty of this work is the derivation and use of fundamental solutions directly in the partial Fourier space. In addition to providing substantial memory and computational savings, knowing closed-form expressions of these fundamental solutions enables optimizations which were previously tedious (cf. Section 2.3.2). The fundamental solutions are found in this context by solving equations involving the transformed Navier operator for an isotropic medium, which has been expressed in eq. (1.24). We now present our general methodology for obtaining the desired partial-Fourier expressions of fundamental solutions. The details of this treatment, including the source code deriving the solutions, are available in the companion notebook [START_REF] Frérot | The Mindlin Fundamental Solution -A Fourier Approach[END_REF]. The process of deriving solutions for the Kelvin and Boussinesq-Cerruti fundamental problems involves finding elements of ker( N ), the 6-dimensional space of functions u satisfying the ODE N [ u] = 0 . This space has been studied in Section 1.2.2, as the Boussinesq-Cerruti tensor is constructed from members of ker( N ). Although the Kelvin fundamental solution is not in ker( N ), it can be constructed from it.

The Kelvin solution

The Kelvin problem in partial Fourier space consists in solving the distributional ODE:

N [ U k ∞ ](q, x 3 , •) = e k δ x 3 , for all (q, x 3 ) ∈ R 2 × R, (2.15)
where δ x 3 is the one-dimensional Dirac distribution supported at x 3 . To find U k ∞ , we follow the methodology of [START_REF] Chaillat | A New Fast Multipole Formulation for the Elastodynamic Half-Space Green's Tensor[END_REF] and seek the displacement vector separately in each semi-infinite interval extending from the source point: .16) We can, without loss of generality, set x 3 = 0 since U ∞ is invariant by translation along e 3 (we then have

U k ∞ (q, x 3 , 3 ) =        U k ,- ∞ (q, x 3 , 3 ) 3 ∈ ] -∞, x 3 [, U k ,+ ∞ (q, x 3 , 3 ) 3 ∈ [x 3 , +∞[. ( 2
U k ∞ (q, x 3 , 3 ) = U k ∞ (q, 0, 3 -x 3 )). Each contribution U k ,±
∞ satisfies the homogeneous Navier equation, and hence belongs to ker( N ). Using eq. ( 1.25) and eq. ( 1.26) with the requirement that U k ,± ∞ decay as 3 → ±∞, we obtain: (2.17) where C k ,± ∈ C 3 are the remaining free coefficients. The latter are determined by requiring that U k ∞ , expressed as a distribution by the single formula (employing the Heaviside function instead of eq. ( 2.16))

U k ,± ∞ (q, 0, 3 ) = A ± (q, 3 ) • C k ,± ,
U k ∞ (q, x 3 , 3 ) = H ( 3 -x 3 ) U k ,+ ∞ (q, x 3 , 3 ) + 1 -H ( 3 -x 3 ) U k ,- ∞ (q, x 3 , 3 ) (2.18)
should be continuous at 3 = x 3 and should satisfy the distributional ODE eq. ( 2.15) (recalling that H (• -x 3 ) = δ x 3 ). As can be seen in the companion notebook, this results in the following expression for

U ± ∞ = e k ⊗ U k ,± ∞ : U ± ∞ (q, x 3 , 3 ) = 1 q U ± 0,0 (q) ± 0 q( 3 -x 3 ) + U ± 1,0 (q) ± 1 q( 3 -x 3 ) .
(2.19a)

The symbols of eq. ( 2.19) are defined in Table 2.1. On observing that the functions ± 0 (z) and ± 1 (z) verify 2.19b) the regular parts of the distributional derivatives of U ∞ at any order are easily found to be given through the recurrence relations

∇ ± 0 (q 3 ) = ∓q∆ ± ± 0 (q 3 ), ∇ ± 1 (q 3 ) = ∓q∆ ± ± 1 (q 3 ) + qe 3 ± 0 (q 3 ), (
∇ n U ± ∞ (q, x 3 , 3 ) = q n-1 U ± 0,n (q) ± 0 q( 3 -x 3 ) + U ± 1,n (q) ± 1 q( 3 -x 3 ) , U ± 0,n (q) = ∓ U ± 0,n-1 ⊗ ∆ ± + U ± 1,n-1 ⊗ e 3 , U ± 1,n (q) = ∓ U ± 1,n-1 ⊗ ∆ ± .
(2.20)

The Kelvin tensor appears in the operators N ∞ and N only through its first-order gradient, see eqs. (2.10) and (2.14). Since U ∞ is continuous at 3 = x 3 , the representation in eq. (2.18) shows that ∇U ∞ can be identified with a discontinuous function (in the classical sense) and evaluated on the sole basis of formulas in eq. (2.20).

Then, use of Lemma 1 allows to express the operator N ∞ in Fourier space as

N ∞ [w](q, x 3 ) = ∫ ∞ 0 ∇U ∞ (-q, 3 -x 3 ) : w(q, 3 ) d 3 (2.21)
b 4(1 -ν) ∆, ∆ ± (iq 1 /q, iq 2 /q, 0), e 3 ∓ ∆ U - 0,0 1 2µb (bI + ∆ ⊗ ∆ -e 3 ⊗ e 3 ) U + 0,0 1 2µb (bI + ∆ ⊗ ∆ -e 3 ⊗ e 3 ) U - 1,0 - 1 2µb ∆ -⊗ ∆ - U + 1,0 1 2µb ∆ + ⊗ ∆ + ± 0 (z), ± 1 (z) e ∓z , ze ∓z

Boussinesq-Cerruti solution

The Boussinesq-Cerruti fundamental tensor given in Section 1.2.2 is recalled here:

B(q, 3 ) = 1 q B 0,0 (q) + 0 (q 3 ) + B 1,0 (q) + 1 (q 3 ) , with B 0,0 (q) = 1 2µ 2I + (1 -2ν)∆ + ⊗ ∆ -+ ∆ ⊗ ∆ -e 3 ⊗ e 3 , B 1,0 (q) = 1 2µ ∆ + ⊗ ∆ + .
Then, the recurrence relations (2.20) are also valid for the gradients of B. We can now readily construct G and its gradient ∇G using eq. ( 2.14) and the recurrence relations for the regular parts of U ∞ and B.

Displacement gradient computation

Due to the construction of the Mindlin fundamental solution in eq. (2.14), the evaluation of ∇u requires the computation of ∇N ∞ . This operator is singular (see e.g. [START_REF] Bui | Some Remarks about the Formulation of Three-Dimensional Thermoelastoplastic Problems by Integral Equations[END_REF]Bonnet, 2017;[START_REF] Gintides | Solvability of the Integrodifferential Equation of Eshelby's Equivalent Inclusion Method[END_REF], but the present distributional and partial-Fourier framework still allows for a very straightforward treatment. Indeed, we simply have

∇N ∞ [w](q, x 3 ) = -∇ 2 U ∞ (-q, •) w = - ∫ x3 0 q - 0 q( 3 -x 3 ) U - 0,1 (-q) : w(q, 3 ) ⊗ ∆ -+ U - 1,1 (-q) : w(q, 3 ) ⊗ e 3 + - 1 q( 3 -x 3 ) U - 1,1 (-q) : w(q, 3 ) ⊗ ∆ -d 3 - ∫ ∞ x3 q + 0 q( 3 -x 3 ) -U + 0,1 (-q) : w(q, 3 ) ⊗ ∆ + + U + 1,1 (-q) : w(q, 3 ) ⊗ e 3 -+ 1 q( 3 -x 3 ) U + 1,1 (-q) : w(q, 3 ) ⊗ ∆ + d 3 -∇U ∞ (q) : w(q, x 3 ) ⊗ e 3 , (2.23)
with the second equality resulting from the application of ∇ to ∇U ∞ as given by eqs. (2.18) and (2.20). The discontinuity jump in the gradient of U ∞ is computed, using eqs. (2.19) and (2.20), as:

∇U ∞ (q) = U + 0,1 (q) -U - 0,1 (q) = 1 µb (2e 3 ⊗ e 3 -bI ) ⊗ e 3 , (2.24)
with the vector e 3 post-multiplying the discontinuity term in eq. ( 2.23) coming from the fact that distributional terms can only arise from derivatives in the variable 3 using the present partial-Fourier framework.

Remark 4. It is interesting to compare the formulations of ∇N ∞ in the partial Fourier and physical spaces. In the latter case, the function ∇ 2 U ∞ (x, ) has a non-integrable | -x | -3 singularity. Consequently, the singular integral operator is the sum of a Cauchy principal value (CPV) integral and a free term, as pointed out e.g. in [START_REF] Bui | Some Remarks about the Formulation of Three-Dimensional Thermoelastoplastic Problems by Integral Equations[END_REF] (for example, a careful distributional interpretation of the application of ∇ to N ∞ [w] yields both contributions). Then, practical evaluations of ∇N ∞ in the physical space entail special methods for the integration of a CPV [START_REF] Guiggiani | A General Algorithm for Multidimensional Cauchy Principal Value Integrals in the Boundary Element Method[END_REF]. By contrast, the present method, which implicitly and indirectly accounts for both the CPV and the free term contributions (the partial-Fourier counterpart of the latter being the jump term (2.24)) is significantly easier to exploit numerically, as no specialized methods are required.

Discretized operators

In this section, we study the numerical evaluation of the integral operator N [w] for given discretized densities w (the evaluation of M[p], ∇M[p] and ∇N [w] can then be formulated similarly by adapting the considerations made for N [w]). This extends the developments of [START_REF] Zeman | A Finite Element Perspective on Nonlinear FFT-based Micromechanical Simulations[END_REF] to our partial Fourier representation. Although the use of discrete Fourier methods is widespread in simulating the contact of rough surfaces (e.g. Polonsky and Keer, 1999a;[START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF][START_REF] Wang | A New Fast Method for Solving Contact Plasticity and Its Application in Analyzing Elasto-Plastic Partial Slip[END_REF] [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], a theoretical basis for the discretization of continuous operators has, to the best of our knowledge, never been provided.

Spectral discretization and DFT

Let L = (L 1 , L 2 , L 3 ) ∈ R +3 be the three lengths defining the discretized domain and N = (N 1 , N 2 , N 3 ) ∈ N 3 the number of points considered in each direction. Let us also define the following sets:

Z 2 N = k ∈ Z 2 : - N 1 2 < k 1 < N 1 2 , - N 2 2 < k 2 < N 2 2 ,
(2.25a)

X 3 = x i 3 i=N 3 -1 i=0 ⊂ R + with L 3 = sup X 3 -inf X 3 , (2.25b) X = k 1 L 1 N 1 e 1 + k 2 L 2 N 2 e 2 + x 3 e 3 : k ∈ Z 2 N and x 3 ∈ X 3 ⊂ B p . (2.25c)
We recall from eq. ( 1.31) that B p is a semi-infinite cell of size L 1 × L 2 when projected on ∂B. X is a set of N 1 × N 2 × N 3 discrete points, which projected on the (Ox 1 x 2 ) plane forms a regular grid, while the projection on (Ox 3 ) gives the set of chosen positive values X 3 . Unlike full-space Fourier methods [START_REF] Moulinec | A Numerical Method for Computing the Overall Response of Nonlinear Composites with Complex Microstructure[END_REF][START_REF] Wang | A New Fast Method for Solving Contact Plasticity and Its Application in Analyzing Elasto-Plastic Partial Slip[END_REF][START_REF] Zeman | A Finite Element Perspective on Nonlinear FFT-based Micromechanical Simulations[END_REF], we are free to choose the spacing of points in X 3 . Any B p -periodic eigenstress w, as well as the B p -periodic displacement u = N [w], can be expressed as complex Fourier series:

w( x, x 3 ) = k ∈Z 2 w (k, x 3 ) exp(2π i k • x), (2.26a)
u( x, x 3 ) = k ∈Z 2 u (k, x 3 ) exp(2π i k • x), (2.26b)
where w (k, x 3 ) and u (k, x 3 ) are the Fourier coefficients of the series and ki = k i /L i . Because w is B p -periodic, it is no longer convolvable with ∇G, making it impossible to evaluate N [w] as a convolution (the same remark applies to N ∞ [w] and the other operators). The operator N [w] can nevertheless still be evaluated by means of the non-periodic partial Fourier representation of the fundamental solution obtained in Section 2.2.2, thanks to the following result:

Theorem 3. Let w be B p -periodic. Then N [w] is B p -periodic and N [w]( x, x 3 ) = 1 4π 2 k ∈Z 2 ∫ ∞ 0 ∇G(-2π k, x 3 , 3 ) : w (k, 3 ) d 3 exp(2π i k • x).
(2.27)

Proof. See Appendix A.2.
Remark 5. All the fundamental solutions presented in Section 2.2.2 have a O(q -1 ) weak singularity at q = 0 , which does not prevent normal use of their continuous inverse Fourier transforms. By contrast, discrete transform evaluation at k = 0 is not possible. For computing displacements, we can arbitrarily set e.g. G(0 , •) = 0 , following common practice [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF][START_REF] Zeman | A Finite Element Perspective on Nonlinear FFT-based Micromechanical Simulations[END_REF]. In this work, only the operator p → M[p] requires this adjustment, as all other operators involve gradients of fundamental solutions, which have no singularity in q, see eq. ( 2.20). The displacement resulting from the "regularized" operator is correct up to an additive constant and a linear displacement depending only on x 3 . The former will be determined by contact conditions (i.e. imposed average load or displacement), while the latter can be ignored: it diverges in a semi-infinite domain, and we are only interested in surface displacements for the contact problem. Note that it produces a gradient constant w.r.t x 3 , which is taken into account in ∇N. We suppose however that this gradient alone is not enough to trigger a plastic response of the material.

Unlike the method of [START_REF] Sainsot | A Numerical Model for Elastoplastic Rough Contact[END_REF] (and subsequent adaptations by Chen, Liu, and [START_REF] Chen | Fast Fourier Transform Based Numerical Methods for Elasto-Plastic Contacts of Nominally Flat Surfaces[END_REF][START_REF] Wang | A New Fast Method for Solving Contact Plasticity and Its Application in Analyzing Elasto-Plastic Partial Slip[END_REF], the direct use of the closed-form expression of ∇G implies that eq. ( 2.27) is exact. The discretization error inherent in the numerical calculation of ∇G based on the discrete Fourier transform of the real-space fundamental solution ∇G [START_REF] Firth | Discrete Transforms[END_REF][START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF] is therefore avoided, which is a definitive advantage over the previously-mentioned sampling methods. Moreover, storage of the discrete values ∇G(-2π k, x 3 , 3 ) is not necessary, yielding substantial memory gains, especially for higher order operators such as ∇N involving the fourth-order tensor ∇ 2 G.

Implementing the proposed method nevertheless entails unavoidable approximations. One stems from the necessary truncation of the Fourier series in eq. ( 2.27). Another appears when w results from plastic deformations, which require a local (i.e. physical space) representation in order to perform operations such as the return mapping procedure. Hence, the Fourier coefficients w are approximated using the discrete Fourier transform [START_REF] Firth | Discrete Transforms[END_REF][START_REF] Zeman | A Finite Element Perspective on Nonlinear FFT-based Micromechanical Simulations[END_REF]

w h := DFT[w X ] (2.28)
which is known to cause discretization errors [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]. Note that because operations like the computation of plastic deformations are intrinsically local in the physical space and the application of integral operators (e.g. N ) is local in the partial-Fourier space (for the x 1 and x 2 directions), the solving of the elastic-plastic problem will involve going back and forth between the physical and the partial-Fourier representations, using the discrete Fourier transform. Consequently, numerical evaluation of eq. (2.28) is done with the FFT algorithm [START_REF] Cooley | An Algorithm for the Machine Calculation of Complex Fourier Series[END_REF] because of its very attractive O(N 1 N 2 log(N 1 N 2 )) computational complexity, by means of the open-source library FFTW [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF] for the present implementation.

2.3.2. Discretization and integration in the x 3 direction Equation (2.27) involves integrals in the x 3 direction. The purpose of this section is to present the procedure developed to compute them. First, we introduce a generic interpolation of w in the x 3 direction:

w h (k, x 3 ) = N 3 -1 j=0 w j (k)ϕ j (x 3 ), (2.29)
where ϕ j (resp. w j ) are the interpolation function (resp. the Fourier coefficients of w) evaluated at x 3 ∈ X 3 . From eq. ( 2.27), the evaluation of N [w] takes the form:

N [w]( x, x 3 ) = 1 4π 2 k ∈Z 2 N 3 -1 j=0 ∫ ∞ 0 ∇G(-2π k, x 3 , 3 )ϕ j ( 3 ) d 3 : w j (k) exp(2π i k • x).
(2.30) Considering a specific node x i 3 ∈ X 3 , the associated Fourier coefficients in eq. ( 2.30) are computed as a weighted sum of N 3 convolution integrals. Equation (2.19) reveals, after rearrangements, that the integral in eq. ( 2.30) can be expressed in terms of integrals of the simpler form

∫ ∞ 0 ± k q( 3 -x i 3 ) ϕ j ( 3 ) d 3 .
(2.31)

with k = 0, 1 . Furthermore, N [w] has to be evaluated at every node x i 3 ∈ X 3 , so that in the worst case the total number of integrals eq. ( 2.31) to be computed is O(N 23 ).

This cost may however be mitigated, and we now propose a method whose efficiency is better than that of the naive approach consisting in evaluating eq. ( 2.31) for all 0 ≤ i, j < N 3 . We choose classical Lagrange polynomials as our basis of interpolation functions. Let

E i = [x i 3 , x i+1
3 ] (i ∈ {0, . . . , N -2}) be an element with Lagrange polynomials ϕ L j (j ∈ {0, . . . , n -1}) of degree n -1 ; the center x i c and half-length e i of E i are given by

x i c = 1 2 (x i + x i+1
) and

e i = x i c -x i .
Using the change of variables = x i c + ze i and the properties of exponential functions, we easily find

∫ E i ± 0 q( 3 -x 3 ) ϕ L j ( 3 ) d 3 = ± 0 q(x i c -x 3 ) G ± 0 (q, i, j), (2.32a
)

∫ E i ± 1 q( 3 -x 3 ) ϕ L j ( 3 ) d 3 = ± 0 q(x i c -x 3 ) G ± 1 (q, i, j) + q(x i c -x 3 )G ± 0 (q, i, j) , (2.32b)
with G ± k (q, i, j) given using the standard Lagrange polynomials φL j of degree n -1 defined on [-1, 1] by

G ± k (q, i, j) := e i ∫ 1 -1 ± k (qze i ) φL j (z) dz. (2.32c)
As mentioned, the brute-force computation of eq. ( 2.30) requires O(N 1 N 2 N 2 3 ) computations. However, since we have access to the analytical expression of ∇G, we can exploit it to devise complexity reduction methods. We present here two such methods that lend themselves to a straightforward implementation3 .

Cutoff method

It follows from eq. ( 2.32) that the inner integral in eq. (2.30) decays as exp(-q|x i c -x 3 |), allowing us to define a threshold criterion q|x i c -x 3 | < ε co for deciding which integrals need to be computed for the application of N . In practice, if |x i c -x 3 | > ε co /q, the integral over E i is not computed. This results in a drastic reduction of the number of operations needed to compute the integral of eq. ( 2.27), especially for large wave vectors k since the cutoff length

ε co /q is inversely proportional to q = 2π k . The time complexity is then O( N 2 1 + N 2 2 N 2 3 )
(see Appendix A.4). However it entails an approximation of integral (2.30) whose accuracy depends on the selected cutoff.

Linear integration method An alternative way4 to reduce the number of operations needed in the evaluation of eq. ( 2.30) is to rearrange the terms in eqs. (2.32a) and (2.32b) to separate the x 3 and 3 variables. Unlike the cutoff method, this approach does not introduce an error in the evaluation of eq. ( 2.30).

Theorem 4. For a non-periodic interpolation w h (q, x 3 ) defined by eq. (2.29) where ϕ j are linear shape functions, N ∞ [ w h ] evaluated on x l ∈ X 3 can be expressed as:

N ∞ [ w h ](q, x l ) = e -qx l U - 0,1 (-q) -qx l U - 1,1 (-q) : S ↑ l [ - 0 ](q) + U - 1,1 (-q) : S ↑ l [ - 1 ](q) + e qx l U + 0,1 (-q) -qx l U + 1,1 (-q) : S ↓ l [ + 0 ](q) + U + 1,1 (-q) : S ↓ l [ + 1 ](q) , (2.33a)
with the functionals S ↑ l and S ↓ l defined as:

S ↑ l [f ](q) = l -1 j=0 w j (q) ∫ E j f (q )ϕ j ( ) d + w j+1 (q) ∫ E j f (q )ϕ j+1 ( ) d , (2.33b)
S ↓ l [f ](q) = N -2 j=l w j (q) ∫ E j f (q )ϕ j ( ) d + w j+1 (q) ∫ E j f (q )ϕ j+1 ( ) d .
(2.33c)

Proof. See Appendix A.5.

Since the summed terms in S ↑ l and S ↓ l do not depend on x l , the application of N [ w h ] can be computed in O(N 1 N 2 N 3 ) computations if a two pass algorithm is used: the first term in eq. (2.33a) is computed in a forward pass for l = 0, . . . , N 3 -1 (involving only S ↑ l ), then a backward pass computes the second contribution for l = N -1, . . . , 0 . This however increases the memory cost of the evaluation of eq. ( 2.30) compared to the cutoff method, since the full Fourier representation of N [w] needs to be stored for the two-pass algorithm5 . In addition, because of the separation of variables, floating point precision overflow and underflow can occur in the exponential terms, while the cutoff method is stable in this regard.

Both integration the cutoff and linear integration bring the overall complexity of evaluating

N [w] to O(N ln N ) if N 1 = N 2 = N 3 ,
with the complexity of the operations in the Fourier 2.4. Elastic-plastic integral equation method domain being O(N ). However, since the scaling of the two integration methods differ if nonequal discretization is used, the choice of integration method depends on the desired scaling.

For rough surface contact, as larger surface spectra are often the goal, the O(

N 2 1 + N 2 2 N 2 3 )
scaling is more appropriate for "horizontal refinement." Finally, while the support of x 3 → w h (k, x 3 ) is a priori unknown in an elastic-plastic simulation, it is possible that just a few points in a given sub-surface layer show non-zero values of w h , in which case w i can be treated as a sparse vector, keeping track of the non-zero entries needed for evaluating integrals. This sparse treatment is applicable for both the cutoff and linear integral methods.

Elastic-plastic integral equation method

The use of domain integral equation methods for elastic-plastic analysis is now well established [START_REF] Telles | On the Application of the Boundary Element Method to Plasticity[END_REF][START_REF] Telles | Implicit Procedures for the Solution of Elastoplastic Problems by the Boundary Element Method[END_REF][START_REF] Bonnet | Implicit BEM Formulations for Usual and Sensitivity Problems in Elasto-Plasticity Using the Consistent Tangent Operator Concept[END_REF][START_REF] Gao | An Effective Boundary Element Algorithm for 2D and 3D Elastoplastic Problems[END_REF][START_REF] Yu | Development and Implementation of Some BEM Variants-A Critical Review[END_REF]. In this work, we use the implicit integral equation formulation proposed by [START_REF] Telles | Implicit Procedures for the Solution of Elastoplastic Problems by the Boundary Element Method[END_REF].

Von Mises plasticity

We limit, without loss of generality, the results of this chapter to von Mises plasticity, for which the yield function f y : R 3×3 sym → R is given by

f y (σ ) := 3 2 s
where s := σ -1

3 Tr(σ )I .

(2.34)

We define the cumulated equivalent plastic strain as:

e p := 2 3 ∫ t t 0 ε p dt, (2.35)
where ε p is the plastic strain rate. The plasticity conditions are then written as:

f y (σ ) -f h (e p ) ≤ 0, (2.36a
)

f y (σ ) -f h (e p ) e p = 0, (2.36b)
(2.36c)

where f h : R → R is the hardening function, with the associated flow rule

ε p = 3 e p 2f y (σ ) s(σ ).
(2.36d)

Although multiple choices are possible for f h , we only consider linear isotropic hardening, for which f h is given in terms of the initial yield stress σ y and the hardening modulus E h by

f h (e p ) = σ y + E h e p .
(2.37)

Implicit equilibrium equation

The elastic-plastic state S n of B at step t n is characterized by the cumulated plastic equivalent strain e p n and the total plastic strain ε p n ; we write S n = (e p n , ε p n ). At t n , the total strain in B can be expressed as a function of the applied boundary tractions t D n (known) and the plastic strain ε p n (unknown):

ε n = ∇ sym M[t D n ] + ∇ sym N [C : ε p n ].
(2.38)

Writing eq. ( 2.38) at step t n+1 and taking the difference of the two, we obtain the following implicit incremental equation:

∆ε n = ∇ sym M[∆t D n ] + ∇ sym N C : ∆ε p n (∆ε n ; S n ) .
(2.39)

We can see that this equation combines the balance of momentum (2.3), the kinematic compatibility (2.1) (through the use of the symmetrized gradient ∇ sym ) and the constitutive behavior through (2.2) and ∆ε p n (∆ε n ; S n ) (via a return mapping algorithm symbolized by ∆ε p n (∆ε n ; S n )).

Therefore, solving eq. ( 2.39) for ∆ε n gives the solution of the elastic-plastic problem. To achieve this, we seek the root of the residual function:

R[∆ε n ] := ∆ε n -∇ sym M[∆t D n ] -∇ sym N C : ∆ε p n (∆ε n ; S n ) .
(2.40) Equation ( 2.39) being non-linear in ∆ε n , we have to use an iterative method. Since the operators ∇ sym M and ∇ sym N applied to given arguments are very efficiently evaluated thanks to the use of the FFT, we want to avoid assembling them. Newton-Krylov solvers (see e.g. references in [START_REF] Knoll | Jacobian-Free Newton-Krylov Methods: A Survey of Approaches and Applications[END_REF])) are traditionally well suited for this type of approach, especially since the consistent tangent operator for eq. ( 2.39) is known [START_REF] Bonnet | Implicit BEM Formulations for Usual and Sensitivity Problems in Elasto-Plasticity Using the Consistent Tangent Operator Concept[END_REF]. However, we will present here a convenient Jacobian-less method (La [START_REF] Cruz | Spectral Residual Method without Gradient Information for Solving Large-Scale Nonlinear Systems of Equations[END_REF] to solve eq. ( 2.39), as it offers several advantages over Newton-Krylov solvers, and is readily available as part of the SciPy [START_REF] Jones | SciPy: Open Source Scientific Tools for Python[END_REF] library.

Jacobian-Free Spectral Residual Method

The DF-SANE algorithm developed by La [START_REF] Cruz | Spectral Residual Method without Gradient Information for Solving Large-Scale Nonlinear Systems of Equations[END_REF] is an attractive algorithm because of its low memory requirements compared to traditional Krylov solvers which require storage of span vectors for a subpart of Krylov space. It also simplifies the implementation, as it does not require knowledge6 nor evaluation of the consistent tangent operator. The iteration goes as:

δε i n = -α i σ i R[∆ε i n ], ∆ε i+1 n = ∆ε i n + δε i n (2.41) with σ i = δε i-1 n 2 δR i-1 , δε i-1 n L 2 , δR i = R[∆ε i+1 n ] -R[∆ε i n ]
(2.42) where •, • L 2 is the appropriate scalar product on L 2 (B; R 3×3 ) and α i is a step size determined by a line search on R[∆ε] 2 . More details can be found in (La [START_REF] Cruz | Spectral Residual Method without Gradient Information for Solving Large-Scale Nonlinear Systems of Equations[END_REF][START_REF] Birgin | Spectral Projected Gradient Methods: Review and Perspectives[END_REF].

Method validation

Method validation

Evaluations of integral operators N and ∇N by our approach are compared to the analytical solution of a hydrostatic eigenstrain in a spherical inclusion embedded in a half-space to validate the proposed methodology.

The displacement generated by a constant hydrostatic eigenstrain ε = αT I applied in a spherical region embedded in a half-space is given by u = N [w] with w = C : ε = 2µαT 1+ν 1-2ν I , and has been derived analytically in closed form by [START_REF] Mindlin | Thermoelastic Stress in the Semi-Infinite Solid[END_REF]. We use that reference solution to validate our implementation of operators N and ∇N . Let the inclusion support be the ball of radius a and center (0, 0, c) (with c > a). The surface displacements u r and u z are given, using cylindrical coordinates (r, θ, z), by:

u r (r ) = 4a 3 3R 3 1 β(1 -ν)r, u z (r ) = - 4a 3 3R 3 1 β(1 -ν)c, (2.43

Volume integrals

x /c -.

-.

-.

. One can observe a good agreement of the numerical approximation with the analytical solution outside the inclusion. Some Gibbs effect can be observed at the boundary of the inclusion, with spurious oscillations in the inclusion. This is due to the Fourier approximation of the discontinuous eigenstrain function.

σ z /( µβ) a Mindlin & Cheng ( ) Proposed Method x /c -. -. -. . . . σ θ /( µβ) b Mindlin & Cheng ( ) Proposed Method
The relevant analytical values are (for r = 0 ): .45) with β = αT 1+ν 1-ν . We observe on figures 2.3a and 2.3b a Gibbs effect at the inclusion boundary, caused by the Fourier approximation of the discontinuity of the eigenstrain, as mentioned in Section 2.3.1. This should however not have a significant effect on elastoplastic simulations, as plastic deformations should be continuous provided there is no shear band. Nonetheless, our method accurately represents the large discontinuity in the tangential stress, see fig. 2.3b.

σ θ =        2µ β a 3 3 4ν -3 (z+c) 3 + 6c (z+c) 4 + 1 |z-c | 3 z ∈ [0, a[∪]3a, +∞[, 2µ β a 3 3 4ν -3 (z+c) 3 + 6c (z+c) 4 -2 a 3 z ∈]a, 3a[, (2.44)
σ z =        2µ β a 3 3 6z+c (z+c) 4 -2 |z-c | 3 z ∈ [0, a[∪]3a, +∞[, 2µ β a 3 3 6z+c (z+c) 4 -2 a 3 z ∈]a, 3a[, (2

Algorithmic complexity

We now compare the computational cost of the application of the operator N to that of an elastic solve step of a finite elements simulation involving the same number N = N 1 N 2 N 3 of nodes. Irrespective of the numerical method used, the optimal complexity would be O(N ). In the proposed methodology, the evaluation of N [w] for given w is decomposed into two computational steps: the multiple 2D fast-Fourier transforms and the computation of eq. (2.31). Assuming

N 1 = N 2 = N 3 , the former has a complexity of O(N 1 N 2 N 3 log(N 1 N 2 )) ∼ O(N log(N )),
while the latter has a complexity of O(N 1 N 2 N 2

3 ) ∼ O(N 4/3 ) with a naive implementation. Using the integration methods proposed in Section 2.3.2 brings the algorithmic complexity of the evaluations of eq. ( 2.31) down to O(N ) (see Appendix A.4), so the asymptotic cost of evaluating N [w] is O(N log(N )). For the direct solve step of a finite element elastic simulation with N nodes, the algorithmic complexitiy of a sparse Cholesky factorization is O(N 3/2 )7 (Lipton, Rose, and Tarjan, 1979). We compare in fig. 2.4 the relative computation times of the application of N and the elastic solve step of Akantu, which uses the direct solver package MUMPS [START_REF] Amestoy | A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[END_REF] to perform the factorization. A regular mesh with N nodes was used8 . We can observe that for large problem sizes the computation times fit the theoretical asymptotic complexities, showing the clear advantage of the proposed VIM over FEM. One should also note that memory needed for the factorization of the stiffness matrix for 2 21 nodes was larger than 128 GB whereas the VIM only required 1.27 GB for this case. For large problems, both the memory imprint and (absolute) computing time are two orders of magnitude smaller with the proposed approach.

As a closing remark, evaluating N in physical space and without any acceleration method would entail a O(N 2 ) complexity, making its use unrealistic for 3D problems. This complexity can be brought down to O(N ) by means of a multi-level fast multipole (ML-FM) approach. However, implementing the latter is quite involved in general, and here would require intricate and expensive close-range numerical quadrature methods for dealing with the strong singularity and complex expressions of the kernel in the physical space. This may explain why ML-FM treatments of VIMs have received only limited attention to date. Hierarchical matrices can likewise be employed to reduce the complexity [START_REF] Zechner | A Fast Elasto-Plastic Formulation with Hierarchical Matrices and the Boundary Element Method[END_REF]. deformations, changing the residual displacements of the surface. The latter are used to compute a new contact solution, and so on until convergence. The second coupling method is an interior point algorithm. This class of algorithms is well suited for second-order conic programming. We show how the elastic-plastic contact problem can be recast into a secondorder conic program and devise the interior point algorithm.

N T/T

Fixed point strategy

The periodic elastic contact problem has been covered to great length in Section 1.3. The addition of plastic strain only changes the variational form:

inf u ∈Λ 1 2 ∫ B p σ [u, ε p ] : ε[u] -ε p dV . (3.1)
We can also write this principle of minimum energy on the boundary: one can show that the minimizer of eq. (3.1) can be written as

u = v + N [C : ε p ],
where v is a member of

Γ(ε p ) := v ∈ H 1 (B p ; R 3 ) : v + N [C : ε p ] ≥ 0 and div(σ [v]) = 0 , (3.2)
and minimises the potential energy written as surface integral

inf v ∈Γ(ε p ) 1 2 ∫ ∂ B T [v] • v dS .
(3.3) Remark 6. Note the modified contact condition v + N [C : ε p ] ≥ 0 , which explicitly accounts for the residual displacement at the surface caused by the plastic deformations. Since the contact problem is solved with fixed ε p , eq. (3.3) corresponds to an elastic contact problem with a modified contact surface

h mod := h -N [C : ε p ] • e 3 ∂ B
, a property that the fixed point coupling algorithm exploits.

Because problem (3.3) is the problem treated in Chapter 1, the solutions methods discussed there are directly applicable. Notably, the problem of contact with adhesion can be easily treated.

Plastic coupling

The full elastic-plastic contact coupling scheme, developed by [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF], is given in Algorithm 5. The algorithm leverages the modified contact condition (3.2) (see Remark 6). The central loop's purpose is to determine the increment of residual displacements ∆U p 3 at the surface of the elastic solid. These are incorporated into the surface profile H for the elastic contact solve step, during which plastic deformations are fixed. The result of the elastic contact is a traction distribution T on ∂B that acts as a Neumann boundary condition to the elastic-plastic problem. The solve step yields the total strain increment ∆E, which is used to compute the residual surface displacement. The convergence condition is established on the evolution of ∆U p 3 from one iteration to the next. Finally, when convergence is reached

Fixed point strategy

Algorithm 5 Elastic-plastic contact coupling algorithm [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF].

Data: W (normal load), H (surface profile), S (current state), ε tol (tolerance), N max (maximum iterations)

∆U p 3,prev ← 0 T prev ← T [S] Previous tractions k ← 1 H 0 ← H h norm ← H 0 repeat T ← elastic_contact(W , H)
Polonsky and Keer (1999b) with FFT

∆T ← T -T prev ∆E ← plasticity(∆T, S) Strain increment s.t. R[∆E] = 0 ∆U p 3 ← N [C : ∆ε p (∆E, S)] ∂ B • e 3 Surface residual disp. increment e ← ∆U p 3 -∆U p 3,prev / h norm Error ∆U p 3,prev ← ∆U p 3,prev + λ(∆U p 3 -∆U p 3,prev ) Relaxation H ← H 0 -∆U p 3,prev ∆U p 3,prev ← ∆U p 3 k ← k + 1 until e < ε tol or k > N max update(S, ∆E)
Increment e p and ε p within a specified tolerance, the state S is updated with the converged total strain increment ∆E and traction increment ∆T. The user of Algorithm 5 is free to use any elastic contact solver and non-linear plasticity solver as drop-ins for elastic_contact and plasticity. For the simulations presented in this chapter, the elastic contact solver uses the projected conjugate gradient proposed by Polonsky and Keer (1999b), with an FFT approach for the gradient computation [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF]Rey, Anciaux, and Molinari, 2017). For the non-linear plastic solver, we use the DF-SANE algorithm described in Section 2.4.3. The relaxation parameter λ can take values in ]0, 1] and helps the algorithm to converge on large loading steps.

Since Algorithm 5 is a fixed point, acceleration schemes such as those described by [START_REF] Ramière | Iterative Residual-Based Vector Methods to Accelerate Fixed Point Iterations[END_REF] can be applied to enhance the performance. However, these schemes require the unrelaxed fixed point (Algorithm 5 with λ = 1 ) to converge in order to guarantee convergence. This condition is however almost never satisfied, and there is no guarantee that either the relaxed or the accelerated fixed points will converge. In practice, convergence is affected by the magnitude of the influence of plasticity on the contact solution, which means that low yield stress or high contact stresses (due to a large normal load or large surface amplitude) make it more difficult for Algorithm 5 to find a solution.

Elastoplastic contact validation

We validate the complete method (including elastoplastic contact) against an axi-symmetric FEM analysis by [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] which provides the surface pressure distribution for a rigid spherical indenter on an elastic-perfectly plastic material. For Hertzian contact with ν = 0.3 , the maximum shear stress occurs at a depth x 3 ≈ 0.57a, with a the contact radius. Moreover, the von Mises stress reaches σ y for a maximum surface pressure of p y =1.6σ y [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]. From this, one can compute the total load W y and contact radius a y at the onset of yield [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF]):

(a) W y = π 3 R 2 6E 2 p 3 y , (b) a y = 3W y 2πp y , (3.4)
where R is the indenter radius and E := E/(1 -ν2 ) is the Hertz elastic modulus. We also define

τ y := σ y / √
3 as the shear yield stress. W y , a y and τ y are used to normalize loads, lengths and stresses respectively. Figure 3.1 compares the surface pressure p[u] computed using Algorithm 5 with the corresponding values obtained by [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] (dashed lines), for different load ratios W /W y . The dark continuous line is the circumferential average of the pressure, while the lighter zone shows the maximum and minimum pressure for a given radial coordinate r . The difference between the maximum and minimum pressures at the edge of contact is due to the discretization and the periodic boundary conditions, which make the contact region shape deviate from a disk. The simulation parameters are given in Appendix B.1.1. As mentioned in [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF], normalized results are independent of the E/τ y ratio.

Both solutions show a flattening of the initially-ellipsoidal pressure distribution near the axis of symmetry (also observed in [START_REF] Johnson | An Experimental Determination of the Contact Stresses between Plastically Deformed Cylinders and Sphere[END_REF][START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF][START_REF] Wang | A New Fast Method for Solving Contact Plasticity and Its Application in Analyzing Elasto-Plastic Partial Slip[END_REF], with a plateau that extends as the load is increased. There are however significant differences between the two sets of results: in the plastic range, the agreement on both the contact radius and the plateau value is poor. Also, the results of [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] feature oscillations in the pressure profile at the highest loading increment, which are likely due to the coarse mesh that was used, combined with large loading steps 2 .

To assess these effects and confirm the latter remark, we have effected additional comparisons, this time to an implementation of Algorithm 5 where the strain increment is computed using a first-order FEM approach with the open-source high-performance code Akantu [START_REF] Richart | Implementation of a Parallel Finite-Element Library: Test Case on a Non-Local Continuum Damage Model[END_REF], the rest of the algorithm using the same code as for the VIM. The geometry, material properties, loading and discretization are as given in Appendix B.1.2. Figure 3.2 shows the surface pressure and vertical displacement for a load ratio of W /W y = 5 , the corresponding Hertz elastic solution being also plotted for additional comparison. We observe an excellent agreement between the FEM and VIM solutions. In Increase of the applied load beyond the initial yield shows that the pressure profile deviates from the elliptic Hertzian profile by flattening of the curve at the axis of symmetry, with a plateau whose extent increases with the load. The results of [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] however show oscillations of the pressure profile at high plastic loads which is not reproduced by our simulation. As there is, to our knowledge, no physical reason to these oscillations, they are likely due to the coarse discretization of the finite-element mesh they used in their study

Interior point method

Disclaimer: the work presented here is an ongoing collaboration with Jérémy Bleyer (ENPC, Laboratoire Navier) and Marc Bonnet (ENSTA, Poems).

While the fixed point algorithm of Jacq et al. ( 2002) is a simple solving scheme, it has serious convergence issues which limit its applicability to realistic rough contact situations where a lot of plastic activity is expected. For this reason, we have developed an alternative approach to solve the elastic-plastic contact problem based on interior point methods (cf. Boyd and Vandenberghe, 2004, chap 11 for an introduction). This family of methods is well suited to second-order conic programs (SOCPs), a class of optimization problem that contains the elastic-plastic contact problem with an associated flow rule. Indeed, the nonpenetration contact condition of eq. (1.36a) and the yield condition of eq. (2.36a) are both conic constraints (cf. [START_REF] Yonekura | Second-Order Cone Programming with Warm Start for Elastoplastic Analysis with von Mises Yield Criterion[END_REF][START_REF] Kong | An Interior-Point Algorithm for Elastoplasticity[END_REF], who provide SOCP formulations for the plasticity problem). -.

-.

-.

u /a y The main idea of interior point methods is to introduce a repulsive barrier term for the conic constraints so that the computed iterates always remain within the admissible domain.

VIM FEM

Plastic strain norm

This allows to take large minimizing steps without having to project the solution on the space of constraints as is done by Algorithm 5. In an interior point method, the "strength" of the barrier is adjusted during the course of the algorithm so as to provide a balance between step size and accuracy of the solution, because a solution computed with a barrier only enforces an approximate version of the optimality conditions. We now describe the second-order conic program formulation of the elastic-plastic contact problem, followed by the application of the interior point procedure.

Variational form of elastic-plastic contact

To apply the interior point method to the elastic-plastic contact problem, one needs to write a minimization principle to a solid mechanics problem where the constitutive law has a plastic form. For simplicity, we will again only consider here classical von Mises associated plasticity. However, to derive a variational form of plasticity, we will start in a more general thermodynamic setting and follow the derivation of the properties of generic plastic constitutive laws of [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF]. Let us start by stating three fundamental equations relevant for our purposes:

∫ Ω σ : ε[v] dV = ∫ Ω b • v dV + ∫ ∂Ω t (σ • n) • v dS ∀v ∈ KA(Ω), (3.5)
d dt ∫ Ω e dV = ∫ Ω b • u dV + ∫ ∂Ω (σ • n) • u dS + ∫ Ω r dV - ∫ ∂Ω q h • n dS, (3.6) d dt ∫ Ω η dV ≥ ∫ Ω r θ dV - ∫ ∂Ω q h θ dS, (3.7)
which are respectively the balance of momentum equation, the first and the second law of thermodynamics for a subdomain Ω ⊂ B. The variables e, η, θ , r and q h are respectively the internal energy, the entropy per volume, the temperature, the inner heat rate per volume and the heat rate per area supplied to Ω. The latter equations can be written in local form because Ω is arbitrary:

e = σ : ε + r -div(q h ), (3.8) η ≥ r θ -div q h θ , (3.9)
where we have used eq. (3.5) with v = u in the conservation of energy to have the stress and strain rate appear in place of the power of external forces. We introduce the Helmholtz free energy f = eηθ which should be minimum when the system is in thermodynamic equilibrium, giving us a minimization principle. The conservation of energy allows to express the entropy rate as a function of f :

η = 1 θ σ : ε -f + r -div(q h ) -η θ .
(3.10)

Replacing the expression for η in the entropy inequality gives the dissipation inequality [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF]:

σ : ε ≥ f + η θ + q h θ ∇θ (3.11) ⇔ σ : ε ≥ f , (3.12)
whose second form comes from the assumption that θ is uniform in space and time. Next we assume that the free energy f is function of the kinematic variables ε (the linearized strain tensor), ξ and γ . The latter two are internal variables, later to be identified to the plastic strain tensor and cumulated plastic strain, which represent the material rearrangement due to crystal plasticity molecular mechanisms. The variables conjugate3 to (ε, ξ, γ ) are (σ, χ, k) derived from f :

σ = ∂ f ∂ε , χ = - ∂ f ∂ξ , k = - ∂ f ∂γ , (3.13)
and represent the force quantities related to the kinematic quantities (ε, ξ, γ ), with of course σ being the Cauchy stress tensor. The force variables χ and k, which will later be related to the deviatoric part of σ and the yield stress, represent the forces associated with the inner material rearrangement characterized by ξ and γ . Developing the time derivative of f = f (ε, ξ, γ ) in eq. (3.12) with the chain rule gives the following form of the dissipation equation:

χ : ξ + k γ ≥ 0.
(3.14) [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF] discuss in detail the nature of the postulates that are needed to derive a rate independent plasticity law, and in particular that the maximum dissipation postulate and the normality law (i.e. associated flow rule) are equivalent. We shall only give the main results that will allow us to formulate a principle of minimum free energy which integrates plasticity.

Let us first define the space S ⊂ R 7 of stress variables (χ, k) which are constitutively admissible. The existence of this set is physically motivated by experiments [START_REF] Bui | Étude de l'évolution de la frontière du domaine élastique avec l'écrouissage et relations de comportement élasto-plastique des métaux cubiques[END_REF] which show for metals that it is convex. When (χ, k) reach the boundary of S, dislocations start to nucleate and move in the atomic lattice. Since S is convex, it can be described by a support function ϕ [START_REF] Boyd | Convex Optimization[END_REF]:

ϕ( ξ, γ ) = sup (χ ,k)∈S χ : ξ + k γ . (3.15)
The function ϕ gives the maximum energy dissipated for a given couple ( ξ, γ ). An admissible stress should not dissipate more energy than ϕ( ξ, γ ) for all ( ξ, γ ). The Legendre transform of ϕ is the indicator function of S [START_REF] Boyd | Convex Optimization[END_REF].

Von Mises plasticity with isotropic linear hardening [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF] give as an application the example of a material with a von Mises elasticity criterion, linear isotropic and kinematic hardening. As in Chapter 2, we shall only consider isotropic hardening. For this purpose, the free energy is expressed as

f (ε, ξ, γ ) = 1 2 (ε -ξ ) : C : (ε -ξ ) + 1 3 E h γ 2 .
(3.16)

To link back with the variables that where introduced in the previous chapter, we have ε p = ξ and e p = 2/3γ . It can be seen from the definition of the internal static variables that χ = σ and k = -2γ E h /3 = -2/3E h e p . The von Mises criterion can be expressed with:

| dev(χ )| + k ≤ c 0 , (3.17)
where dev gives the deviatoric part of a tensor and c 0 = 2/3σ y is the initial radius of the von Mises yield cylinder in stress space. In addition, the maximum dissipation postulate is equivalent to the normal flow rule [START_REF] Simo | Computational Inelasticity[END_REF], so ξ is colinear to dev(χ ). The dissipation potential can be readily computed4 [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF][START_REF] Ladevèze | Constitutive Relation Error Estimators for (Visco)Plastic Finite Element Analysis with Softening[END_REF]:

ϕ( ξ, γ ) = sup (χ ,k )∈S χ : ξ + k γ = sup (χ ,k )∈∂S | dev(χ )|| ξ | + k γ = sup (χ ,k )∈∂S | dev(χ )|| ξ | -(| dev(χ )| -c 0 ) γ = sup (χ ,k )∈∂S | dev(χ )|(| ξ | -γ ) + c 0 γ = c 0 γ | ξ | ≤ γ +∞ | ξ | ≥ γ .
We now establish the principle of minimum free energy in a backward-Euler fashion [START_REF] Simo | Computational Inelasticity[END_REF]: the free energy at time t n is equal to the free energy at t n+1 plus the dissipated energy between times t n and t n+1 . The equilibrium is then sought by minimizing the total free energy in the system at t n with respect to

Ξ n+1 = (u n+1 , ε n+1 , ε p n+1 , γ n+1 , n+1 , U n+1 ): inf Ξ n+1 ∈KA(B) ∫ B f (ε n+1 , ε p n+1 , γ n+1 ) + ∫ t n+1 t n ϕ ε p , γ dt dV -W n+1 U n+1 , (3.18)
where U n+1 is the rigid body displacement of the surface corresponding to the applied normal load W n+1 . The function space KA(B) contains, in addition to the kinematic compatibility of ε and the gap condition, the incompressibility condition Tr(ε p n+1 ) = 0 .

Conic programming

The numerical solution of problem (3.18) requires approximation of the time derivative terms, which we replace with finite differences: ε p ≈ ∆ε p /∆t. This allows to express quantities at t n+1 as a function of known quantities at t n and to write the minimization principle in terms of increments. We develop the free energy5 :

f (ε n + ∆ε, ε p n + ∆ε p , γ n + ∆γ ) = f (ε n , ε p n , γ n ) + ∇f (ε n , ε p n , γ n ) • (∆ε, ∆ε p , ∆γ ) + f (∆ε, ∆ε p , ∆γ ), = f (ε n , ε p n , γ n ) + f (∆ε, ∆ε p , ∆γ ) + σ n : (∆ε -∆ε p ) + 2 3 E h γ n ∆γ
and approximate the dissipation integral:

∫ t n+1 t n ϕ ε p , γ dt ≈ ∆tϕ ∆ε p ∆t , ∆γ ∆t = ϕ ∆ε p , ∆γ .
We can now write the SOCP corresponding to the elastic-plastic contact problem:

inf ∆u ,∆ε ,∆ε p ,∆γ , n+1 ,U n+1 ∫ B f ∆ε, ∆ε p , ∆γ + σ n : (∆ε -∆ε p ) + 2 3 E h γ n ∆γ + 2 3 σ y ∆γ dV -W n+1 U n+1 , (3.19a)
under the constraints:

n+1 ≥ 0, (3.19b) ∆ε p ≤ ∆γ , (3.19c)
Tr(∆ε p ) = 0,

(3.19d) ∆ε = ε[∆u] (3.19e) n+1 = (∆u + u n ) • e 3 -h -U n+1 . (3.19f)
The last constraint is needed because we consider n+1 as a free variable. We introduce the Lagrange multipliers p n+1 and s = (s 0 , s) ∈ R7 for the first two (unilateral) constraints, such that p n+1 ≥ 0 and s ≤ s 06 . [START_REF] Yonekura | Second-Order Cone Programming with Warm Start for Elastoplastic Analysis with von Mises Yield Criterion[END_REF] have shown that the third constraint (3.19d) can be ignored if one replaces all occurrences of ∆ε p by dev(∆ε p ) in problem (3.19a). The last constraint is enforced with the Lagrange multiplier λ , so that the Lagrangian of the SOCP (3.19a) is:

L(∆ε, ∆ε p , ∆γ , n+1 , p n+1 , s, λ , U n+1 ) = F (∆ε, dev(∆ε p ), ∆γ , U n+1 ) - ∫ ∂ B p n+1 n+1 + λ ( n+1 -∆u 3 + U n+1 -[u n ]) dS - ∫ B (s 0 ∆γ + s : ∆ε p ) dV ,
with F being the objective function of (3.19a).

Optimality conditions

Next we assume eq. (3.19e) holds by construction and take variations of the Lagrangian with respect to the kinematic variables to obtain its stationarity conditions.

∂ ∆u L, v = ∫ B σ n + C : (∆ε -dev(∆ε p )) : ε[v] dV + ∫ ∂ B λ v • e 3 dS = 0, (3.20a) ∂ ∆ε p L = -dev σ n + C : (∆ε -dev(∆ε p )) -s = 0 , (3.20b) ∂ ∆γ L = 2 3 E h ∆γ + 2 3 E h γ n + 2 3 σ y -s 0 = 0, (3.20c) ∂ n+1 L = -p n+1 -λ = 0, (3.20d)
∂ U n+1 L = -W n+1 - ∫ ∂ B λ dS = 0. (3.20e)
These conditions respectively give the equilibrium equation, the interpretation of s as the negative deviatoric stress, the hardening rule and definition of s 0 as the hardened yield stress, the equivalence of λ (which enforces the relation betweenu and n+1 ) and the normal pressure, and finally the total applied load condition. To obtain the complete optimality conditions, we add the following complementarity conditions for the variables related to inequality constraints.

p n+1 n+1 = 0, (3.20f)
∆γs 0 + ∆ε p : s = 0.

(3.20g)

The last equation implies that (∆γ , ∆ε p ) and s lie on the boundary of their respective constraint cone. Hence, ∆γ = ∆ε p and s 0 = s , resulting in ∆ε p : s = -∆ε p s , which is the normality rule: the plastic strain increment is colinear to -s with a positive constant equal to ∆γ /s 0 . This is rewritten as:

s 0 ∆ε p + ∆γs = 0 . (3.20h)

Integral operators

The set of equations (3.20), complemented with the conic constraints and eq. (3.19f), describes completely the incremental elastic-plastic contact problem. Since the kinematic unknowns are considered independent, we can once again consider ∆ε p as initial eigenstrain and apply to

M 3 ∇ sym M[•e 3 ] ε from surface pressures N ∇ sym N [•]
ε from eigenstresses

M 33 e 3 • M[•e 3 ]| ∂ B
Surface normal displacements from pressures

N 3 e 3 • N [•]| ∂ B Surface normal displacements from eigenstresses R ∇ sym N [C : •] -• "Elastic" strain
eq. (3.20a) the treatment developed in Chapter 2. This actually allows us to exactly solve the equilibrium with:

∆ε = ∇ sym M[p n+1 e 3 ] + ∇ sym N [C : dev(ε p n + ∆ε p )] -ε n , (3.21a) ∆u = M[p n+1 e 3 ] + N [C : dev(ε p n + ∆ε p )] -u n . (3.21b)
The use of integral operators essentially removes ∆ε (and by extension ∆u) from the unknowns of the optimality conditions. Replacing ∆ε in eq. (3.20b) yields:

-dev C : ∇ sym M[p n+1 e 3 ] + ∇ sym N [C : dev(ε p n + ∆ε p )] -dev(ε p n + ∆ε p ) -s = 0. (3.22)
In order to simplify notations, we use the symbols defined in Table 3.1 in lieu of the linear operators used above. We also define the following operations for x, s ∈ R m+1 :

x • s = x 0 s 0 + x • s x 0 s + xs 0 = mat(x)s = mat(s)x, (3.23a) mat(x) = x 0 x T x x 0 I m . (3.23b)
We can then write the optimality conditions in compact form:

          -2 3 E h 0 0 DCRD I 0 0 DCM 0 0 -N 3 CD 0 I -M 33 a 0 0 0 a T 0                                  ∆γ ∆ε p s n+1 p n+1 U n+1                        =              2 3 E h γ n + 2 3 σ y -DCRDε p n N 3 CDε p n -h W n+1              , (3.24a) n+1 p n+1 = 0, (3.24b) (∆γ , ∆ε p ) • s = 0 , (3.24c)
where a is a vector of ones.

The basic idea of primal-dual interior point methods, which is the subcategory of IPMs that we employ here, is to operate a Newton-Raphson scheme on the non-linear system of equations (3.24) with the following perturbation: (3.25b) where x = (∆γ , ∆ε p ), µ > 0 and e = (1, 0 ). This effectively adds a logarithmic barrier to the conic constraints of the problem, whose parameter µ decreases during the algorithm so the perturbed problem approaches the original problem. Let

n+1 p n+1 = µ, (3.25a) x • s = µe,
z k = x k , s k , k n+1 , p k n+1 , U k n+1
T be the collection of unknowns at the Newton-Raphson iteration k and R(z k ) the residual of the perturbed optimality conditions. Algorithm 6 describes the primal-dual interior point method. We can see in this algorithm that the barrier parameter µ is updated at each step with the complementarity gap reduced by a parameter 0 < ω ≤ 1 . The value of this parameter determines the propensity of the algorithm to reduce the duality gap (low values of ω) or to follow the so-called central path, which is the set of optimal solutions for varying values of µ [START_REF] Boyd | Convex Optimization[END_REF]. The value of ω therefore affects the efficiency: heuristics can be used instead of a user-provided value, but [START_REF] Mehrotra | On the Implementation of a Primal-Dual Interior Point Method[END_REF] has proposed a predictorcorrector scheme which improves the algorithm convergence [START_REF] Bleyer | Advances in the Simulation of Viscoplastic Fluid Flows Using Interior-Point Methods[END_REF].

Algorithm 6 Primal-dual interior point method applied to elastic-plastic contact.

Data: W (normal load), H (surface profile), S (current state), ε tol > 0 (tolerance), N max (maximum iterations), ω ∈ (0, 1] (centering parameter)

P ← W Initial pressure G ← H Initial gap X ← (1, 0) S ← ( 2/3σ y , 0) U ← 0 µ ← P • G + X • S k ← 0 repeat Solve R (Z)δ Z = -R(Z) cf. Section 3.2.2 Line search on maximum α such that Z + αδ Z is admissible Z ← Z + αδ Z µ ← ω(P • G + X • S) k ← k + 1 until µ < ε tol and R(Z) < ε tol or k > N max

Linear system solving

The computationally intensive step of Algorithm 6 is the resolution of the Newton-Raphson direction. To reduce the computational burden, we actually solve a reduced version of the linear system where the only unknowns are (δε p , δp, δU ). This is possible because the matrices arising from the linearization of the complementarity conditions are invertible by block (cf. Appendix B.2). The solving of the resulting linear system requires particular attention: because we make use of integral operators, the assembled matrix corresponding to the discretized tangent operator is dense. The storing of this matrix becomes problematic rather quickly, so for large problem sizes we are restricted to iterative solvers, which are sensitive to conditioning (Quarteroni, Sacco, and Saleri, 2007). Since interior point approaches introduce a barrier, the linearized tangent contains terms of type P • G -1 and X • S -1 that result from the reduction of the complementarity conditions, cf. Appendix B.2. As µ → 0 and the unknowns get closer to the constraints, the components of these terms either diverge or go to zero, and the conditioning number of the tangent operator diverges, preventing convergence of whatever iterative solver is chosen if the system is not preconditioned. The issue of preconditioning the tangent in interior point methods is an active area of research [START_REF] Gondzio | Matrix-Free Interior Point Method[END_REF]. Since we wish to retain a matrix-free formulation of the algorithm, we cannot use incomplete LU factorization as is sometimes proposed in large sparse systems [START_REF] Bocanegra | Using a Hybrid Preconditioner for Solving Large-Scale Linear Systems Arising from Interior Point Methods[END_REF]. We have therefore developed a tentative preconditioner that inverses an approximation of the reduced tangent operator, cf. Appendix B.2. Finally, the predictor-corrector scheme proposed by [START_REF] Mehrotra | On the Implementation of a Primal-Dual Interior Point Method[END_REF] demands an additional linear solve with the tangent operator. With direct solver this can be done at a little expense since the tangent has already been factorized, but with iterative solvers this demands the full cost of an extra solve step.

The issues of performance of the primal-dual interior point method with integral operators are central questions of an ongoing research effort in collaboration with Dr. Jérémy Bleyer and Dr. Marc Bonnet.

Part II.

Applications

Wear of rough surfaces: understanding the wear coefficient

Abstract

Sliding contact between solids leads to material detaching from their surfaces in the form of debris particles, a process known as wear. According to the well-known Archard wear model, the wear volume (i.e. the volume of detached particles) is proportional to the load and the sliding distance, while being inversely proportional to the hardness. The influence of other parameters are empirically merged into a factor, referred to as wear coefficient, which does not stem from any theoretical development, thus limiting the predictive capacity of the model. Based on a recent understanding of a critical length-scale controlling wear particle formation, we present two novel derivations of the wear coefficient: one based on Archard's interpretation of the wear coefficient as the probability of wear particle detachment and one that follows naturally from the up-scaling of asperity-level physics into a generic multi-asperity wear model. As a result, the variation of wear rate and wear coefficient are discussed in terms of the properties of the interface, surface roughness parameters and applied load for various rough contact situations. Both new wear interpretations are evaluated analytically and numerically, and recover some key features of wear observed in experiments. This work shines new light on the understanding of wear, potentially opening a pathway for calculating the wear coefficient from first principles.

Disclaimer

This chapter is reproduced from the article Lucas Frérot, Ramin Aghababaei, and Jean-François Molinari (May 2018). "A Mechanistic Understanding of the Wear Coefficient: From Single to Multiple Asperities Contact". In: Journal of the Mechanics andPhysics of Solids 114, pp. 172-184. doi: 10.1016/j.jmps.2018.02.015, with permission of all authors. My personal contributions to this article include the development of the models, the code implementation, the running of simulations, the figure production and the writing of the article. [START_REF] Hatchett | Experiments and Observations on the Various Alloys, on the Specific Gravity, and on the Comparative Wear of Gold. Being the Substance of a Report Made to the Right Honourable the Lords of the Committee of Privy Council, Appointed to Take into Consideration the State of the Coins of This Kingdom, and the Present Establishment and Constitution of His Majesty's Mint[END_REF], but our current understanding was built upon research conducted in the middle of the last century [START_REF] Burwell | Survey of Possible Wear Mechanisms[END_REF]. Wear comes in various forms, with adhesive wear, the process of detachment of surface asperities tip by adhesive forces during the sliding contact of two solids, being one of the most prominent. Systematic wear experiments in the mid-20 th century [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF][START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF] suggested a general relation where the wear rate (i.e. wear volume per unit sliding distance) is linearly proportional to the applied normal load, in a certain range of the latter [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF][START_REF] Rabinowicz | Metallic Transfer between Sliding Metals: An Autoradiographic Study[END_REF], and related to the hardness of the material. Inspired by this experimental evidence, Archard [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF]) generalized Holm's concept of "atom removal" [START_REF] Holm | Electric Contacts: Theory and Applications[END_REF] to "debris removal" and pictured an adhesive wear model. He assumed that an asperity junction of radius a produces a debris volume proportional to a 3 over an effective sliding distance of 2a, giving a linear relationship between wear rate and real contact area at the asperity level. To extend this single-asperity relation to a multi-asperity contact, Archard argued that only a fraction of contacting asperities, a quantity referred to as the "wear coefficient", produces wear particles. This conception of the "wear coefficient" being key in understanding wear, [START_REF] Archard | The Wear of Metals under Unlubricated Conditions[END_REF] claimed that "…one of the most important problems in an understanding of wear is to explain the magnitude of the probability of the production of a wear particle at an asperity encounter." This long-standing problem has remained unresolved, and evaluation of the wear coefficient is still relying on empirical data, with no insight from a physical model.

Similarly to the friction coefficient [START_REF] Svetlizky | Classical Shear Cracks Drive the Onset of Dry Frictional Motion[END_REF], the wear coefficient is a system property that depends on many parameters including applied load [START_REF] Zhang | Transition between Mild and Severe Wear in Aluminium Alloys[END_REF][START_REF] Riahi | Wear Map for Grey Cast Iron[END_REF], material properties of sliding bodies (e.g. fracture toughness [START_REF] Fleming | The Relationship between Crack Propagation Rates and Wear Rates[END_REF][START_REF] Challen | Prediction of Archard's Wear Coefficient for Metallic Sliding Friction Assuming a Low Cycle Fatigue Wear Mechanism[END_REF]) and properties of interface (e.g. dry or lubricated contact [START_REF] Kato | Wear in Relation to Friction -a Review[END_REF], roughness parameters [START_REF] Kato | Wear in Relation to Friction -a Review[END_REF], chemical properties [START_REF] Mischler | Wear of CoCrMo Alloys Used in Metal-on-Metal Hip Joints: A Tribocorrosion Appraisal[END_REF]). However, all these effects are currently empirically merged into the wear coefficient, which limits the applicability of Archard's model. Therefore, our goal in this work is to further the understanding of adhesive wear in multiasperity setting, based on a physics-based understanding of the wear process. To this effect, we base our approach on the concept of critical length-scale governing the formation of wear particles [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF]. This concept stems from energy balance between the available deformation energy in an asperity encounter and the energy required to detach a wear particle [START_REF] Griffith | VI. The Phenomena of Rupture and Flow in Solids[END_REF][START_REF] Rabinowicz | The Effect of Size on the Looseness of Wear Fragments[END_REF]Rabinowicz, 1995). The balance states that contacts smaller than a critical length-scale d * plastically deform and contacts larger than d * break into a wear particle. This was recently shown with molecular dynamics simulations [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF]. The critical length-scale d * = λ • G∆w/σ 2 j is function of the shear modulus G, the fracture energy per unit area ∆w and the junction strength σ j with a shape factor λ, all of which can be determined by direct experiments or analytical predictions, with no fit parameter.

In this chapter, we present two new conceptions of the wear coefficient that is built upon the critical length-scale concept. The first concept incorporates Archard's interpretation of the wear coefficient as the "probability of production of a wear particle". The second is based on an up-scaling of single-asperity wear considerations to a multi-asperity contact setting. Both concepts are analytically and numerically studied in different contact situations. We compare them in the context of contact of self-affine surfaces and show that the concept based on Archard's interpretation leads to a constant coefficient wear within a certain range of load, while the second does not. We finally give possible explanations as to why this is the case, and potential improvements to the proposed models.

Review of Archard's wear model

Archard's wear model [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF] can be decomposed into two parts: the single-asperity wear model and the contact model. At the single asperity level, the amount of material removed in an asperity interaction is considered proportional to a 3 (a is the asperity contact radius), whereas the sliding distance required to break off the wear particle is proportional to a. This gives the general relationship for the wear rate (worn volume per sliding distance) of a single asperity (subscript 1, see Table 4.1):

R 1 = ωA, (4.1)
where A is the contact area of the asperity and ω is a generic shape factor, equal to 1/3 in the case of spherical asperities forming hemi-spherical wear particles. This hypothesis is discussed in (Rabinowicz, 1995) and has been verified for isolated debris with molecular dynamics [START_REF] Aghababaei | On the Debris-Level Origins of Adhesive Wear[END_REF]. At the multi-asperity level, Archard makes two hypotheses: I. the size and shape of the individual contact areas are given by a contact model considering a rough surface made of spheres with radius r uniformly distributed in depth, with density d (number of spheres per unit distance).

II. a probability factor K applies on each contact to account for the fact that not all asperity encounters result in a wear particle. Archard assumes K is independent of a.

Using these, the global wear can be related to the applied load W :

R(W ) = KωA c (W ) (4.2) = K ωbd 2 p + 1 cd W 2 p+1
, where A c is the true contact area, b, c and p are found in Table 4.2. In the case of a rigid-plastic material the wear equation is a linear relationship written as:

R(W ) = Kω W H .
(4.3) 2) shows that the total wear rate is proportional to the total contact area, with a proportionality factor K (modulo a shape factor), and eq. ( 4.3) recovers the experimentallyobserved linear relationship between wear rate and load [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF] plastic behavior. This limitation comes exclusively from Archard's contact model, as it governs theW → A c relationship. Note that other contact models [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF]Persson, 2006) would yield a linear W → A c relationship in the purely elastic case.

Wear coefficients based on critical length-scale

Archard interprets the proportionality factor as the probability that a given asperity encounter yields a wear particle, and introduces it at the asperity level by expressing the single asperity wear rate as R 1 = KωA. This suggests that the particle formation is a random process at the asperity level and independent of the micro-contact size. This is inconsistent with recent results [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF]) that exhibit a Griffith-like criterion governing the detachment of wear particles for homogeneous materials, the latter thus being deterministic.

We now fundamentally enrich Archard's interpretation of the wear coefficient by considering, within a multi-contact setting, a critical micro-contact area:

A * ∝ G∆w σ 2 j 2 , (4.4)
which is the square of the critical length-scale defined in [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF]. This length-scale is derived from the balance of available deformation energy and required energy to form a wear particle. In this expression, G is the shear modulus, ∆w is the fracture energy per unit area and σ j is the asperity junction's shear strength. We now complement the eq. ( 4.2) with two new statements:

A. the size and shape of individual contacts are the outcome of contact between random surfaces. The area of a single contact (also called contact cluster) is a random variable A characterized by a probability density function p(A,W ).

B. The process of debris formation is deterministic at the asperity level. It is governed by a critical area A * : if the area A of a cluster is larger than A * , a wear particle is formed (fig. 4.1).

With this, we are now in position to propose two alternative definitions of the wear coefficient. A contact cluster forms a wear particle upon sliding if its area is larger than A * , otherwise the asperity in contact deforms upon sliding without breaking (hypothesis B.). This hypothesis brings asperity-level physics in the wear particle formation process through the critical cluster size A * .

Wear coefficient based on Archard's interpretation

Archard's interpretation of the wear coefficient is the "probability of the production of a wear particle at an asperity encounter". With the critical area A * in mind, the probability of production of a wear particle is simply:

K ≡ P(A > A * ,W ) = ∫ ∞ A * p(A,W ) dA, (4.5)
which is the probability that a realization of the random variable A is larger than A * for a given load W . This definition maintains Archard's interpretation of the wear coefficient while defining it at the multi-asperity level.

Alternative formulation of wear coefficient

Although Archard's interpretation of the wear coefficient is widely accepted, and therefore of prime interest, it is not physically justified by Archard, and still relies on a direct sum of wear volume produced by every micro-contact (weighted by a constant probability coefficient), which is not compatible with the critical length-scale concept. Here, we propose a new interpretation of the wear coefficient: we directly use the single asperity wear relation of eq. ( 4.1) with no probability coefficient and sum the wear rate of all contacts forming a wear particle, thus truly up-scaling eq. ( 4.1) to multi-asperity contact.

Considering a system of finite size, there is a finite number N (W ) of clusters in the system. We define n(A,W ) ≡ N (W ) • p(A,W ). Provided A → An(A,W ) is integrable on [0, +∞[, we can compute the total wear rate, using Archard's single-asperity wear rate and the critical length-scale concept, as a weighted sum of all contacts larger than A * :

R(W ) = ∫ ∞ A * R 1 (A)n(A,W ) dA = ∫ ∞ A * ωAn(A,W ) dA = ωA * c (W ), (4.6)
where A * c is the cumulated area of all clusters forming a wear particle and ω is an average shape factor. We define the function K(W ) as:

K(W ) ≡ A * c (W ) A c (W ) = ∫ ∞ A * Ap(A,W ) dA ∫ ∞ 0 Ap(A,W ) dA , (4.7)
where A c (W ) is the total contact area. This allows us to write the wear rate as R(W ) = K(W )ωA c (W ), in the same form as eq. (4.2). We call K(W ) the "wear coefficient", but its definition is substantially different from Archard's, as we now have an area ratio instead of a probability of particle formation. The definition in eq. (4.7) results naturally from the upscaling of the single-asperity wear considerations to multi-asperity contact, unlike Archard's wear coefficient, which was introduced a posteriori into the wear rate equation.

Numerical results

Without making a priori assumptions on the distribution of cluster areas, numerical simulation allows the direct application of the wear models on realistic surfaces [START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF]. We evaluate the wear coefficient and the wear rate using a model for rough-surface contact consisting of a flat, semi-infinite elastic medium in contact with a rigid solid having a random rough surface. The contact problem is solved using a boundary-element approach [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF]Polonsky and Keer, 1999b) which fully accounts for elastic interactions. The wear coefficient is investigated through the analysis of the contact map across a representative sample of simulations. The area of each contact cluster is determined from the contact map using an 8-neighbors flood-fill algorithm1 . This approach allows a direct computation of the wear coefficient using eqs. (4.5) and (4.7).

The rough surfaces are random self-affine (fractal) isotropic surfaces [START_REF] Mandelbrot | Fractal Character of Fracture Surfaces of Metals[END_REF][START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF], with height function h(x, ) generated by a filtering algorithm [START_REF] Hu | Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis[END_REF]. They are defined through their power spectral density (PSD). The surfaces are isotropic, so the surface PSD depends only on the radial coordinate q = 2π /λ, where λ is a wavelength. The surface PSD is defined using λ s , λ l , respectively the short wavelength cut-off and the large wavelength cut-off [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF][START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF]. Between q l = 2π /λ l and q s = 2π /λ s , the PSD decays as q -2(H +1) . H is the Hurst exponent, governing the self-affine behavior of the fractal rough surface. We vary the Hurst exponent between 0.6 and 0.8 as is commonly occurring in natural surfaces [START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF]. There are three main properties governing the statistics of fractal surfaces:

• L/λ l where L is the size of the surface, governs the representativity of the surface [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF]. A large value yields a surface with many large asperities, allowing better statistics while also reducing the effect of the periodic boundary conditions.

• λ l /λ s controls the range of the PSD. This quantity influences the spectrum bandwidth of the surface [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF]).

• λ s /∆l, where ∆l is the discretization size, governs the discretization error which causes bias in evaluated mechanical and statistical quantities, but also governs the resolution of the details of the fractal contact clusters (Yastrebov, Anciaux, and Molinari, 2017a).

A large value reduces both of these sources of error. 4.2e. Black area is not in contact, yellow contact clusters are smaller than A * , and red clusters are larger. As the load is increased, clusters grow to span more of the available area, occasionally merging to form larger clusters, as shown in fig. 4.2e with the increase in contact area and maximum cluster size. The latter increases dramatically as the clusters merge, forecasting percolation. Figure 4.2f represents the increase in the number of clusters per increase of contact area ( dN / dA c ). Positive values indicate the regime where contact-area growth is dominated by nucleation of clusters, and negative values indicate the regime where cluster merging dominates contact-area growth. We focus our study to the former regime, in which the clusters are far enough to neglect interactions in the debris forming process. As postulated by [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF], this corresponds to the mild wear regime [START_REF] Wang | Wear and Wear Transition Modeling of Ceramics[END_REF][START_REF] Das | Ultra Mild Wear in Lubricated Tribology of an Aluminium Alloy[END_REF][START_REF] Hokkirigawa | Wear Mode Map of Ceramics[END_REF]. We let W /W 0 ∈ [0.001, 0.06], discretized in thirty load steps for the other simulations presented in this paper.

Statistical distribution of contact clusters

Figure 4.3 presents the probability density function of cluster areas at multiple applied loads for three different spectrum range parameters (λ l /λ s ) and a Hurst exponent of 0.8 . Results for different Hurst exponents and λ l /λ s = 128 are shown in the inset. In agreement with previous experimental observations (Majumdar and Bhushan, 1991;Dieterich and Kilgore, 1996) and numerical simulations [START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF][START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF][START_REF] Campañá | Elastic Contact between Self-Affine Surfaces: Comparison of Numerical Stress and Contact Correlation Functions with Analytic Predictions[END_REF], the probability density function of cluster areas follows a power-law in a given interval of A. The evaluated exponent value of 1.5 (using a maximum likelihood estimator [START_REF] Clauset | Power-Law Distributions in Empirical Data[END_REF]) is well within the measured range of 1.05 -2.69 from the experiments of Dieterich and Kilgore (1996), close to values of 1.6 predicted by the overlap model [START_REF] Mandelbrot | Stochastic Models for the Earth's Relief, the Shape and the Fractal Dimension of the Coastlines, and the Number-Area Rule for Islands[END_REF]) and 1.54 -1.56 measured in the experiments of [START_REF] Majumdar | Role of Fractal Geometry in Roughness Characterization and Contact Mechanics of Surfaces[END_REF]. Values predicted by various contact mechanics models (discrete and continuous) fully accounting for long-range asperity interactions are in the range of 1.45 -1.6 [START_REF] Hyun | Elastic Contact between Rough Surfaces: Effect of Roughness at Large and Small Wavelengths[END_REF][START_REF] Campañá | Using Green's Function Molecular Dynamics to Rationalize the Success of Asperity Models When Describing the Contact between Self-Affine Surfaces[END_REF][START_REF] Müser | Meeting the Contact-Mechanics Challenge[END_REF].

Figure 4.3 also shows that the behavior of p(A,W ) is independent of W within the powerlaw interval, and that the upper bound of this interval increases with the load, in agreement with experimental observations (Majumdar and Bhushan, 1991). The upper bound of the power-law interval is also increasing with the surface PSD range. Moreover, when cluster sizes are normalized with the shortest wavelength in the system, all probability density functions collapse to a single curve within the power-law interval. The inset of fig. 4.3 shows that the Hurst exponent has limited influence on the distribution of clusters and only affects the fall-off behavior at large cluster areas, making the power-law approximation less accurate in this range.

Wear coefficient

Figure 4.4a shows the wear coefficient, K num , as defined in eq. (4.5). It exhibits common features with fig. 4.3, namely a power-law behavior in a given interval of A * and the increase of the upper bound of that interval with W and with λ l /λ s . The load-invariant power law of fig. 4.4a signals a constant wear coefficient for a given load range. Figure 4.4b shows the evolution of the wear coefficient as a function of load (surface used has λ l /λ s = 128 ). Remarkably, regardless of A * , the wear coefficient transitions from zero (i.e. no observable wear) to a constant value given by the power law of fig. 4.4a, which corresponds to the proportionality constant observed in experiments [START_REF] Burwell | Survey of Possible Wear Mechanisms[END_REF][START_REF] Rabinowicz | Metallic Transfer between Sliding Metals: An Autoradiographic Study[END_REF][START_REF] Archard | The Wear of Metals under Unlubricated Conditions[END_REF]. This is the first time to our knowledge that a model derived from first principles predicts a constant wear coefficient within a given load range and a transition The main graph shows p(A/λ 2 s ,W ), the probability density function of normalized cluster areas for H = 0.8 and varying λ l /λ s , evaluated using twenty logarithmic bins. The inset shows the probability density function for values of H ∈ {0.6, 0.7, 0.8} with λ l /λ s = 128 . p(A,W ) can be approximated by a power-law within a certain cluster size interval, inside which it is independent of the applied load W . Increasing W or λ l /λ s increases the upper bound of the power-law interval. Normalization with the smallest wavelength λ s collapses all distributions to a single curve within the power-law interval. Varying the Hurst exponent has a limited effect on the resulting distribution. from no observable wear (i.e. wear coefficient is zero) to mild wear (i.e. wear coefficient is constant). This transition occurs at a critical load that depends on the value of A * , and is larger for systems with higher A * (as would be the case in lubricated contact). For systems with small A * values (e.g. poor lubrication conditions), the critical load may be lower than the lowest load we simulate. Note that in the presence of lubrication, the wear volume is affected by a change of A * as well as a reduction of the solid contact area.

Figure 4.5a and 4.5b show the wear coefficient K num and the wear rate as functions of the load, computed from eqs. (4.6) and (4.7) respectively. Figure 4.5a shows that the wear coefficient is zero up to a transition load that depends on A * . For A * in the low range of values we simulate, the transition load is smaller than 0.001 • W 0 . Similarly to Archard's wear coefficient, this new interpretation is able to exhibit the no-wear/wear transition that has been observed in experiments [START_REF] Colaço | An AFM Study of Single-Contact Abrasive Wear: The Rabinowicz Wear Equation Revisited[END_REF]. After the transition load, the wear coefficient increases monotonically up to one. In fig. 4.5b, the wear rate is quasi-linear after a transition region [START_REF] Kato | Wear Mechanisms[END_REF].

Although there is a qualitative agreement of eq. (4.5) with experimental observation, a quantitative agreement is difficult to obtain because measurements of hardness and RMS of slopes are difficult and not systematic in wear experiments. In the experiments of [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF], good care is taken in eliminating all possible sources of wear but adhesive wear, and they provide a good base for a quantitative comparison. However no precise estimation of the wear coefficient is given. They nonetheless give measurements of the transition load to severe wear, which is the limit of our model, and thus cannot be computed using our base Regardless of A * there exists a critical load at which the wear coefficient transitions from zero (i.e. no wear debris) to a constant value (i.e. steady-state mild wear regime). This critical load largely depends on the value of A * : the transition occurs at a higher critical load for contacts with lower interfacial shear strength and consequently larger A * (i.e. better lubrication condition). In (b) the y-axis is linear up to 10 -5 . While the wear coefficient is non-constant with the load, its derivative decreases with the load. However, regardless of A * , the limit value of K is one, which does not correspond to experimental observations. Similarly, the wear rate is increasing non-linearly with the load, although its derivative stabilizes to a fixed value with increasing load.

assumptions that neglect asperity interactions during the wear process.

Analytical results

In order to understand the properties of the proposed wear coefficients in the contact of self-affine surfaces, we develop analytical expressions of K and K for a power-law distribution of contact clusters, as suggested by fig. 4.3:

p PL (A, A m ) =              C A ∈ [0, A s ] C A A s -α A ∈ [A s , A m ] 0 A ∈ [A m , +∞) . (4.8)
C is a normalizing factor such that

∫ ∞ 0 p PL (A, A m ) dA = 1
and α is the power-law exponent, set to the value measured in fig. 4.3. The effect of increasing the contact load is taken into account via increasing A m . Applying the previously stated equations for the wear coefficient yields the following expressions (c.f. Appendix C.1 for details):

K PL =                1 - A -α s (1 -α)A * A 1-α m -αA 1-α s A * ∈ [0, A s ] 1 - (A * ) 1-α -αA 1-α s A 1-α m -αA 1-α s A * ∈ [A s , A m ] 0 A * ∈ [A m , +∞)
, (4.9)

K PL =                    1 - A -α s (1 -α 2 )(A * ) 2 A 2-α m -α 2 A 2-α s A * ∈ [0, A s ] 1 - (A * ) 2-α -α 2 A 2-α s A 2-α m -α 2 A 2-α s A * ∈ [A s , A m ] 0 A * ∈ [A m , +∞)
.

(4.10)

The graphs of these expressions are displayed in fig. 4.6. It is apparent that the behavior of K PL and K PL when A m increases is very different: while K PL tends to a limit whose value depends on A * , the limit value of K PL is one, regardless of A * . This discrepancy is caused by the value of α ≤ 2 . If α > 2 , then K PL has a finite limit dependent on A * , exhibiting a wear coefficient independent of the load.

Discussion

The values predicted for Archard's wear coefficient (fig. 4.4) are in the order of 10 -5 -10 -1 , which is in the range of values reported in experiments (10 -8 -10 -1 ; see e.g. Rabinowicz, 1995;[START_REF] Hsu | Wear Prediction of Ceramics[END_REF], especially values for dry sliding [START_REF] Archard | The Wear of Metals under Unlubricated Conditions[END_REF]. We have shown that the upper bound of the power-law interval in the wear coefficient increases with the range of the surface PSD (λ l /λ s ). Due to computational limits, the maximum value of λ l /λ s simulated was 128, much lower than values of 10 3 and above measured on real surfaces [START_REF] Power | Roughness and Wear during Brittle Faulting[END_REF][START_REF] Liang | Fractal Analysis of Engineering Ceramics Ground Surface[END_REF][START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF]. We postulate that the range of predicted values for the wear coefficient could reach the lower experimental values for larger PSD ranges: if λ s is decreased, both the power-law interval upper bound and A * /λ 2 s are increased. Moreover, it has been shown [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF]) that an elasto-plastic constitutive behavior also increases the upper bound of the power-law interval, and an elastoplastic behavior is expected to occur even at small loads [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF]Majumdar and Bhushan, 1991).

The newly proposed model for the wear coefficient (fig. 4.5) introduces an interpretation that does not rely on Archard's assumption that the probability of wear particle formation is the same for all contacts, but does not predict the behavior commonly observed in experiments, as the wear coefficient tends to a value of one regardless of the critical length-scale. The value of the power-law exponent α = 1.5 is the cause of this discrepancy: the integrals in eq. (4.7) are dominated by the value of A m , the largest cluster size. This can be seen in the expression for the self-affine model in eq. (4.10). There is no definite reasoning on the origins of α, as it depends on the spectrum of the surface, but also on constitutive behavior. A value of α > 2 would make the limit value of K dependent on the critical length-scale and strictly smaller than one. It may be the case that under elasto-plastic constitutive assumption this would be satisfied [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF], but it requires efficient numerical methods to be checked accurately. In addition, our contact model does not include transport of wear particles or relative movement of surfaces and possible reattachment. These aspects likely also change p(A,W ) and therefore change the properties of the wear coefficient.

The formation of a tribological layer [START_REF] Popov | Simulation of Wear in Combustion Engines[END_REF][START_REF] Scherge | Characterization of Wear Debris of Systems Operated under Low Wear-Rate Conditions[END_REF] may also change the debris formation process and geometrical as well as physical properties of the surfaces. This may induce lubrication and/or formation of a third body, which could be accounted for via the proper modifications of the contact model and the energies involved in the Griffith criterion for particle formation.

Conclusive experimental data for validation of our proposed models is difficult to obtain, as wear experiments often include several coupled physical phenomenon (e.g. adhesion, chemistry, temperature). The experiments of [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF] manage to focus on adhesive wear only, and show that the wear coefficient should remain constant for a specific load range before increasing sharply with the load. They relate the transition load to the hardness, arguing that when the average applied pressure reaches H /3 , individual contacts and detached wear particles start interacting to form larger debris. We have shown that our simulations remain in the regime where the contact area growth is dominated by cluster nucleation. Individual contacts are not close enough for interaction, and the wear coefficient computed could be compared with a measurement done with an average pressure less than H /3 . However, the experiments of [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF] lack precise measurements of the wear coefficient and surfaces' spectra, which make quantitative comparison impossible. Some wear experiments [START_REF] Power | Roughness and Wear during Brittle Faulting[END_REF][START_REF] Davidesko | Evolution of Slip Surface Roughness through Shear[END_REF] provide surface spectra, but do not observe the variations of the wear coefficient with the applied load, although the values reported by [START_REF] Power | Roughness and Wear during Brittle Faulting[END_REF] are close to 0.02 in a linear wear regime, which is in the range of values we observe.

Ideally, a validating wear experiment would provide spectrum measurement of the surface, study only adhesive wear of homogeneous materials (similarly to [START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF]) and give wear coefficient measurements with respect to the load. We hope that the ideas we put forth in our paper can encourage the experimental community to move towards this kind of experiments. In the meantime, we are optimistic that an elasto-plastic contact model may bring us closer to an experimental validation of those ideas.

Nonetheless, the new approach we describe both incorporates physical properties of the interface (e.g. surface energy, junction strength), geometrical and statistical properties of the surfaces in contact, while remaining generic with respect to the contact model. This leaves the possibility to tune the latter to advances in contact mechanics (e.g. elasto-plastic contact), but also real measurement of rough surfaces.

Conclusion

We have presented two physics-based interpretations for the adhesive wear coefficient. The use of the critical length-scale is the key ingredient in our approach, and the results obtained show common characteristics with experimental observations. In particular, the results based on Archard's assumption that the wear coefficient is the probability of wear particle formation exhibit a transition in the wear coefficient, after which it remains constant with load. We also propose a second interpretation based on a direct up-scaling of single asperity adhesive wear, and observe that the wear coefficient depends on load. We hypothesize that this is due to limitations in the contact model rather than the approach taken to derive the wear coefficient. We provide analytical results for both wear concepts in different simple contact situations that give a broad understanding of how the wear coefficient evolves with interface physical parameters and surface roughness.

Crack Nucleation in the

Adhesive Wear of an Elastic-Plastic Half-Space

Abstract

The detachment of material in an adhesive wear process is driven by a fracture mechanism which is controlled by a critical length-scale. Previous efforts in multi-asperity wear modeling have applied this microscopic process to rough elastic contact. However, experimental data shows that the assumption of purely elastic deformation at rough contact interfaces is unrealistic, and that asperities in contact must deform plastically to accommodate the large contact stresses. We therefore investigate the consequences of plastic deformation on the macro-scale wear response. The crack nucleation process in a rough elastic-plastic contact is analyzed in a comparative study with a classical 2 plasticity approach and a saturation plasticity model. We show that plastic residual deformations in the 2 model heighten the surface tensile stresses, leading to a higher crack nucleation likelihood for contacts. This effect is shown to be stronger when the material is more ductile. We also show that elastic interactions between contacts can increase the likelihood of individual contacts nucleating cracks, irrespective of the contact constitutive model. This is confirmed by a statistical approach we develop based on a Greenwood-Williamson model modified to take into account the elastic interactions between contacts and the shear strength of the contact junction.

Disclaimer

This chapter is reproduced from the article Lucas [START_REF] Frérot | Crack Nucleation in the Adhesive Wear of an Elastic-Plastic Half-Space[END_REF]. "Crack Nucleation in the Adhesive Wear of an Elastic-Plastic Half-Space". In: arXiv: 1910.05163 [cond-mat], with permission of all authors. My personal contributions to this article include the development of the models, the code implementation, the running of simulations, the figure production and the writing of the article. T he removal of a wear debris particle through adhesive forces is mainly driven by a fracture process, and as such obeys a balance between the energy release rate (i.e. the energy released by the crack front advancing) and the fracture toughness (i.e. the energy required to create new surfaces). This [START_REF] Griffith | VI. The Phenomena of Rupture and Flow in Solids[END_REF] energy balance has been verified in atomistic simulations [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF][START_REF] Aghababaei | Asperity-Level Origins of Transition from Mild to Severe Wear[END_REF][START_REF] Brink | Adhesive Wear Mechanisms in the Presence of Weak Interfaces: Insights from an Amorphous Model System[END_REF] and reduced to a critical length-scale d * ∝ G∆w/τ 2 j , with G being the shear modulus, ∆w the surface energy and τ j the junction shear strength. This gives a simple geometric criterion for the formation of hemispherical wear particles: if the contact diameter between two hemispherical asperities is larger than d * then a wear particle detaches from the surface upon shearing of the system. The issue of transposing asperity-scale wear mechanisms to multi-asperity contact is key in the goal of formulating predictive wear models [START_REF] Meng | Wear Models and Predictive Equations: Their Form and Content[END_REF][START_REF] Vakis | Modeling and Simulation in Tribology across Scales: An Overview[END_REF]. [START_REF] Popov | Adhesive Wear and Particle Emission: Numerical Approach Based on Asperity-Free Formulation of Rabinowicz Criterion[END_REF] and Pham-Ba, [START_REF] Brink | Adhesive Wear Mechanisms in the Presence of Weak Interfaces: Insights from an Amorphous Model System[END_REF] have recently proposed energy-based models for the formation of wear particles in multi-asperity settings. The former investigates the formation of hemispherical wear particles in an elastic rough surface contact by computing an energy-favored particle diameter based on the elastic deformation energy of the contact solution. The latter formulates the energetic competition between the formation of a single vs. multiple wear particles (for 2D line contacts), thus giving an energy approach to the crack shielding mechanism that leads to disjoint but sufficiently close contacts forming a single wear particle [START_REF] Aghababaei | Asperity-Level Origins of Transition from Mild to Severe Wear[END_REF].

In Chapter 4, we have applied the critical length-scale concept to rough elastic contact by defining a critical cluster area A * ∝ (d * ) 2 above which micro-contacts should form a wear particle. This work however suffers from two model inadequacies: (a) the contact solution is given by an elastic contact model, (b) it assumes that A * exists and is proportional to the square of d * . The latter is related to the topography and shape of contacts. Contacts resulting from interfaces with rough surfaces are not disk-shaped and the crack is not expected to produce hemispherical wear particles. Moreover, this does not account for disjoint contacts that may form a single particle [START_REF] Aghababaei | Asperity-Level Origins of Transition from Mild to Severe Wear[END_REF][START_REF] Pham-Ba | Adhesive Wear and Interaction of Tangentially Loaded Micro-Contacts[END_REF]. Thence, it is unclear if the Griffith balance can be characterized with a comparison as "naive" as A ? ≥ A * with A being the area of a single contact cluster. The former shortcoming (a) provides to the wear models developed in the previous chapter an unrealistic contact solution. Since both A * and the contact solution indirectly depend on σ y1 , the outcome of an elastic contact problem (which is independent of σ y ) leads to a paradox: more ductile materials (with lower σ y ) have a higher A * and thus wear less than more brittle materials (with higher σ y ). This is due to the contact solution being insensitive to changes in σ y , but also to the lack of surface roughness evolution in sliding. When the contact of two asperities does not create a wear particle, plastic smoothing of the asperities occurs, thus creating larger contacts. In any case, a contact model incorporating plastic effects is needed.

Experimental data clearly shows that some form of plasticity must occur at rough contact interfaces [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF][START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF][START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF][START_REF] Zhang | Experimental and Finite Element Analyses of Contact Behaviors between Non-Transparent Rough Surfaces[END_REF]. Modeling these interfaces with a non-linear constitutive behavior is however a challenge because of the multi-scale nature of rough surfaces. [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF] were the first to use the finite-element method to study elastic-plastic rough contact with a classical von Mises formulation [START_REF] Simo | Computational Inelasticity[END_REF]. [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF] have developed a volume integral method that we have refined with a Fourier approach to be able to handle the large discretization requirements of multi-scale rough surfaces, cf. Chapters 2 and 3. The majority of published works on elastic-plastic contact does not rely on classical formulations of plastic flow, but rather on the concept of surface flow pressure, which is associated to the hardness of a material [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF][START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF]Majumdar and Bhushan, 1991;Persson, 2001a). The surface flow pressure is usually taken as the maximum value of the mean pressure caused by an indenter of a given shape (it is therefore shape-dependent). [START_REF] Tabor | The Hardness of Metals[END_REF] has shown that for a spherical indenter, the mean pressure saturates at a value close to 3σ y (with σ y being the yield stress). The models previously mentioned are thereafter referred to as "saturation models", in the sense that they apply this concept of a maximum average pressure to a multi-asperity contact model see e.g Tabor, 1951, chap. 9 and assume that a given contact cannot have a pressure exceeding the saturation pressure noted p m2 . They have been used in conjunction with boundary integral approaches [START_REF] Almqvist | On the Dry Elasto-Plastic Contact of Nominally Flat Surfaces[END_REF] to study friction [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF], but have to our knowledge never been compared to classical plasticity formulations, and the relevance of the choice between the two plasticity models has never been studied. [START_REF] Akchurin | A Stress-Criterion-Based Model for the Prediction of the Size of Wear Particles in Boundary Lubricated Contacts[END_REF] and [START_REF] Li | Mesoscale Numerical Modeling for Predicting Wear Debris Generation[END_REF] have used a saturation plasticity model to compute the contact solution and applied a stress based criterion for the removal of debris particles: from the contact pressure profile, they computed the resulting von Mises stress caused in a purely elastic medium. Then the zones of the material where the von Mises stress exceeds the yield stress are removed, changing the surface profile. This has the advantage of foregoing any geometrical consideration, at the expense of providing an ad-hoc removal process that is not derived from the fracture energy balance, as well as using a stress distribution that does not account for plastic deformations.

In this work, we wish to investigate the multi-asperity wear process from a fracture mechanics perspective and understand the influence of plasticity in the contact model on the global wear response. To this end, we focus on the crack nucleation process in the contact of a rigid self-affine rough surface with an elastic-plastic flat half-space. One measure of particular importance is the crack spatial density. While it is not a measure of wear itself, crack nucleation is a necessary process of wear, and understanding what are the roles of the normal load, the critical nucleation stress, the junction resistance, and plastic behavior in crack nucleation is a fundamental step towards predictive wear models. We first highlight the importance of the choice of a plasticity model and the implications it may have on the contact response (Section 5.1). We then show how the crack nucleation density in a rough surface elastic-plastic contact depends on the fracture mechanics properties of the material, as well as the applied load and the junction shear strength (Section 5.2). To rationalize the differences between the elastic, the saturation and the von Mises plasticity approaches, we study the contact behavior of a single asperity to understand under which conditions a crack can nucleate and what is the influence of residual plastic deformations on this process (Section 5.3). These findings are applied to a simple multi-asperity contact model [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF] in order to obtain analytical predictions of the scaling of the crack density with respect to system parameters like the applied normal load (Section 5.4). These predictions are confronted to elastic rough contact simulations which are able to reproduce the contact shielding effect under shear loading, as seen in molecular dynamics simulations [START_REF] Aghababaei | Asperity-Level Origins of Transition from Mild to Severe Wear[END_REF][START_REF] Pham-Ba | Adhesive Wear and Interaction of Tangentially Loaded Micro-Contacts[END_REF]. Simulation results show that ductile materials in contact with rough surfaces produce more crack nucleation sites than brittle materials due to the residual stresses caused by plastic deformations. This effect is not captured by the elastic contact model nor the saturation plasticity model, indicating that the resolution of the aforementioned wear paradox should include the full plastic contact response. This further implies that the true contact area is not the only key quantity in wear modeling.

Elastic-plastic contact

At our disposal are (at least) two formulations of the elastic-plastic contact of solids, the choice of which may have an impact on the subsequent results we wish to obtain. The first formulation, which has been used in the finite-element studies of [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF], follows the classical modeling hypothesis of metal plasticity [START_REF] Simo | Computational Inelasticity[END_REF], which have both experimental [START_REF] Bui | Étude de l'évolution de la frontière du domaine élastique avec l'écrouissage et relations de comportement élasto-plastique des métaux cubiques[END_REF] and theoretical backgrounds [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF], and additionally are valid in other context than contact. The second, developed by [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF] and extended by [START_REF] Almqvist | On the Dry Elasto-Plastic Contact of Nominally Flat Surfaces[END_REF] in conjunction with a boundary integral approach, postulates that the surface contact pressure should nowhere exceed a maximum value p m . This is based on observations that for spherical indentation the mean contact pressure does not exceed a value around 3σ y [START_REF] Tabor | The Hardness of Metals[END_REF]. Recent finiteelement simulations [START_REF] Krithivasan | An Analysis of Three-Dimensional Elasto-Plastic Sinusoidal Contact[END_REF][START_REF] Song | Elastic-Plastic Spherical Indentation: Deformation Regimes, Evolution of Plasticity, and Hardening Effect[END_REF][START_REF] Ghaednia | A Review of Elastic-Plastic Contact Mechanics[END_REF] show that p m /σ y may depend on the ratio σ y /E * (with E * := E/(1 -ν 2 ) being the contact modulus) as well as the wavenumber in the case of sinusoidal contact surfaces. Despite these reports, saturation models are often used in computational tribology [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF][START_REF] Akchurin | A Stress-Criterion-Based Model for the Prediction of the Size of Wear Particles in Boundary Lubricated Contacts[END_REF][START_REF] Li | Mesoscale Numerical Modeling for Predicting Wear Debris Generation[END_REF] due to their simplicity and ease of implementation. Besides increasing the magnitude of the true contact area compared to elastic contact, plasticity influences other aspects of the contact interface (such as contact pressures). These additional aspects may be key ingredients in wear modeling. For this reason we wish to provide a comprehensive comparison between the von Mises associated plasticity and the saturation plasticity with p m = 3σ y in a rough contact situation, and determine the consequences of the choice of one model over the other. We start by giving the full mechanical formulation for both models, then proceed to the comparison.

Definitions In this work, we consider a deformable three-dimensional solid B spanning a half-space, with its (flat) boundary noted ∂B. Moreover, we suppose a horizontal periodicity in the cell B p = [0, L] 2 × R + . We note σ the Cauchy stress tensor, which is related to the small-strain tensor ε and the plastic strain tensor ε p by the relation σ = C : (ε -ε p ) where C is the usual isotropic linear elasticity tensor. The strain tensor is given by kinematic compatibility as a function of the displacement field: ε = ∇ sym u. Finally, σ is expected to be divergence-free to satisfy conservation of momentum without volume forces.

We additionally define some surface quantities: t and p := t • e 3 are respectively tractions and normal pressures applied on ∂B. Other surface quantities are noted with an over-bar • when not explicitly defined on ∂B, e.g. u is the surface displacement.

Saturation: perfect plasticity

The simplest form of saturation model, conceptually close to the notion of "perfect plasticity", is given as [START_REF] Almqvist | On the Dry Elasto-Plastic Contact of Nominally Flat Surfaces[END_REF]

min p 1 2 ∫ ∂ B p pM[p]dS - ∫ ∂ B p phdS , (5.1a)
which is a problem of finding the surface pressures p minimizing the complementary energy of the system under the constraints p ≥ 0,

(5.1b)

p ≤ p m , (5.1c)
∫ ∂ B p pdS = W .
(5.1d)

The linear operator M gives the normal surface displacement due to the applied pressure p if B is assumed elastic; h is a continuous function representing the rough surface brought in contact with ∂B and W is the total applied normal load in the periodic cell boundary ∂B p . The gap is defined as := M[p] -h and should satisfy weak Hertz-Signorini-Moreau conditions [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF]:

≥ 0 where p < p m , (5.2a) p ≥ 0,

(5.2b) p = 0 where p < p m .

(5.2c)

The solution to the above constrained optimization problem yields a negative gap where p = p m . The magnitude of the negative gap is often assumed to be the magnitude of the residual plastic displacements. Since the weak optimality conditions do not represent a physical system (the gap should be non-negative everywhere to avoid body interpenetration), it is necessary to replace h in eq. (5.1a) by h s mod := h + h pl , with h pl being in principle a correction due to residual plastic displacements, therefore h pl := -(M[p] -h) where p = p m . [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF] propose an iterative scheme to solve for h pl which we have implemented and made available in the open-source contact library Tamaas [START_REF] Frérot | Tamaas, a High-Performance Library for Periodic Rough Surface Contact[END_REF] (https://c4science.ch/tag/tamaas/).

The "perfect plasticity" aspect of the model comes from the fact that p m is homogeneous on ∂B and constant. [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF] have amended this hypothesis to include a form of hardening. The saturation pressure is simply expressed as a linear function of h pl (i.e. the initial saturation stress is zero, and rises in proportion with h pl ). We will however not discuss this particular model here.

von Mises plasticity

For the Cauchy stress tensor σ , the von Mises yield function f y is defined as

f y (σ ) = 3 2
||s ||, where s := σ -

1 3
Tr(σ )I .

(5.

3)

The equivalent cumulated plastic strain is expressed as the integral of the plastic strain rate ε p from some reference time t 0 :

e p := 2 3 ∫ t t 0 || ε p ||dt .
(5.4)

The admissibility and consistency conditions are written as:

f y (σ ) -f h (e p ) ≤ 0,
(5.5a)

f y (σ ) -f h (e p ) e p = 0, (5.5b)
where f h is the hardening function. In this work, we will only consider functions of the form f h (e p ) = σ y + E h e p , with σ y the initial yield stress and E h the hardening modulus3 . The associated flow rule that determines ε p is given by [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]:

ε p = 3 e p 2f y (σ )
s(σ ).

(5.6)

The numerical integration of the relations above is typically done with a backwards Euler scheme and is classical to the solid mechanics literature [START_REF] Simo | Computational Inelasticity[END_REF]. Its coupling with the equilibrium and contact conditions is however non-trivial.

Solution strategy [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF] established a numerical method for the solution of the elastic-plastic rough contact problem, which we summarize here. The method consists in solving the contact and the plasticity problems separately. The contact problem is solved for fixed plastic deformations: it is effectively an elastic contact problem with a rough surface

h mod := h -u p 3
, with u p 3 the vertical component of the actual4 plastic residual displacement. Various solution strategies for the elastic rough contact problem are available in the literature [START_REF] Bemporad | Optimization Algorithms for the Solution of the Frictionless Normal Contact between Rough Surfaces[END_REF], and we use here the modified conjugate gradient algorithm of Polonsky and Keer (1999b) coupled with the spectral approach of [START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF] for the gradient computation involving the operator M.

The plastic problem is solved with fixed boundary tractions, meaning that the contact area does not evolve during the resolution of the plastic strain increment. The procedure we employ, fully detailed in [START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF], relies on an implicit incremental volume integral equation formulation proposed by [START_REF] Telles | Implicit Procedures for the Solution of Elastoplastic Problems by the Boundary Element Method[END_REF]. The total strain increment is shown to be expressed as:

∆ε = ∇ sym M[∆t] + ∇ sym N [C : ∆ε p (∆ε; S)],
(5.7)

where S := (e p , ε p ) is the current plastic state, ∆t is the increment of surface tractions (in our case ∆t = ∆p e 3 as we are in a normal contact situation). The function ∆ε p (∆ε; S) represents the radial-return algorithm classically used in incremental plastic analysis [START_REF] Simo | Computational Inelasticity[END_REF]. Equation (5.7) is a non-linear equation that can be solved with the DF-SANE algorithm [START_REF] Cruz | Spectral Residual Method without Gradient Information for Solving Large-Scale Nonlinear Systems of Equations[END_REF] which has the advantage of being jacobian-free.

The operators M and N , which are at the heart of the method developed in [START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF], are linear integral operators which compute in B the displacement due to periodic distributions of surface traction and volume eigenstress respectively 5 . Their complete formulation and application in a discretized setting is extensively discussed in [START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF]. The coupling between the elastic contact problem and the plasticity problem is done with a relaxed fixed point strategy [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF][START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF]. The full implementation of the described solution method is also freely available in Tamaas.

Comparison: rough surface

While both plasticity models are phenomenological, associated plasticity is soundly grounded in experimental observations [START_REF] Bui | Étude de l'évolution de la frontière du domaine élastique avec l'écrouissage et relations de comportement élasto-plastique des métaux cubiques[END_REF] as well as thermodynamic principles [START_REF] Reddy | Internal Variable Formulations of Problems in Elastoplasticity: Constitutive and Algorithmic Aspects[END_REF][START_REF] Simo | Computational Inelasticity[END_REF], and expresses a macroscopic picture of dislocation systems at the micro-scale. This is not the case for the saturation models: they depend on the observation that the mean contact pressure saturates for spherical indentation [START_REF] Tabor | The Hardness of Metals[END_REF], which has been challenged by recent finite-element simulations [START_REF] Song | Elastic-Plastic Spherical Indentation: Deformation Regimes, Evolution of Plasticity, and Hardening Effect[END_REF][START_REF] Krithivasan | An Analysis of Three-Dimensional Elasto-Plastic Sinusoidal Contact[END_REF].

We aim here to provide a direct comparison for a rough surface between a perfectly plastic 2 model and the saturation model. The rough surfaces we use throughout this work are self-affine random surfaces. Their power-spectrum density is defined as where q l , q r , q s are the spatial frequencies associated to the long cutoff wavelength λ l , roll-off wavelength λ r and short cutoff wavelength λ s respectively, while H is the Hurst exponent.

ϕ(q) =              C q l ≤ |q| ≤ q r , C q r |q| -2(H +1) q r ≤ |q| ≤ q s , 0
For the purposes of comparison, we chose a rather narrow spectrum: L/λ l = 3 , λ l = λ r and λ r /λ s = 9 . The discretization size ∆l is such that λ s /∆l = 9 . The depth modeled in the 2 approach is L/5 . This gives a discretized system of 243 × 243 × 64 points. Naturally, the discretized saturated model has 243 × 243 points.

The normal loads applied in both models are adimensionalized by W 0 = E * L2 h rms where h rms := |∇h| 2 is the root-mean-square of the surface slopes. In elasticity, this normalization collapses load (W ) vs. true contact area (A c ) for different values of h rms (but the same spectrum parameters) [START_REF] Bush | The Elastic Contact of a Rough Surface[END_REF][START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF]. We do not intend here to modify the spectrum parameters but merely scale h rms , which therefore becomes a non-dimensional measure of surface summit amplitude. This is convenient to compare the results of the 2 and saturated models to an elastic reference, as we expect the contact behavior to depend on the surface peak amplitude because of plasticity. ), at the load W /W 0 = 6.5 • 10 -2 . The total contact ratio is 25% for 2 plasticity, 20% for saturation and 15% for elasticity, resulting in about 20% error in the contact area of the saturated model. Moreover, while the maximum pressure in the 2 model exceeds 3σ y (cf. fig. 5.1a), the average pressure on micro-contacts is closer to 2σ y , which the saturation model fails to capture with an average of 2.5σ . Local features of the contact patches also differ due to the three models being in different contact stages.

Figure 5.2 shows the contact area ratio and the secant slope of the load/contact area relationship, which is a quantity that has been extensively studied in elastic contact [START_REF] Bush | The Elastic Contact of a Rough Surface[END_REF][START_REF] Persson | Elastoplastic Contact between Randomly Rough Surfaces[END_REF][START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF][START_REF] Campañá | Contact Mechanics of Real vs. Randomly Rough Surfaces: A Green's Function Molecular Dynamics Study[END_REF][START_REF] Pastewka | Contact between Rough Surfaces and a Criterion for Macroscopic Adhesion[END_REF]Yastrebov, Anciaux, and Molinari, 2017b;[START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF]. We show here that all models have a sub-linear load/contact area law since the secant slope is decreasing. The saturation plasticity model fails to predict the 2 secant slope and is on average 20% smaller. Both models should give increased slopes when h rms is increased, but it is unclear if the relative error between them should stay constant with respect to h rms . One could nonetheless fit p m in the saturation model to match the elastic-plastic contact area ratio for a given value of h rms . Figure 5.3 shows the probability density function of the surface pressures for the three models, with the addition of reference data from [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF] (fig. 10a) which results from a 2 plasticity criterion used in a finite-element approach. The features of the curve corresponding to the 2 plasticity models are not qualitatively reproduced by the saturation model: the peak at p = 2.5σ y is non-existent and as expected the distribution of the saturated pressures tends to a Dirac at p = p m whereas the 2 distribution tends to zero. The difference between our results and those of [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF] can be explained by the coarseness of the mesh they used, as well as different spectrum parameters (e.g. H = 0.5 in their case).

To study macroscopic quantities (such as the true contact area), the saturation model is only suitable if one has reference data to fit p m . However, it fails to qualitatively reproduce local surface quantities (as can be seen in figs. 5.1 and 5.3), in addition to providing no information on the complete stress state at and below the contact surface. The saturation model may be useful in applications where quantitative errors in the contact area magnitude or its topography may be accepted, only to obtain qualitative relations (e.g. sealing, electric conductance). It can however give unequivocally wrong results when local stress based quantities drive the phenomenon one wishes to study. We shall see in the next section that this can be the case in wear modeling when we consider the crack nucleation process in an elastic-plastic rough contact.

Crack nucleation in rough surface contact

Most of the atomistic investigations of adhesive wear processes use geometries that contain stress concentrators [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF][START_REF] Milanese | Emergence of Self-Affine Surfaces during Adhesive Wear[END_REF][START_REF] Brink | Adhesive Wear Mechanisms in the Presence of Weak Interfaces: Insights from an Amorphous Model System[END_REF], such that in their model system the debris formation is only controlled by the Griffith energy balance. However, without a defect/stress concentration, a crack described by linear elastic fracture mechanics cannot nucleate. In a half-space geometry, with plasticity constitutive behavior, no such concentration can exist. For that reason, it is necessary to introduce a critical nucleation tensile stress σ c .

In linear elasticity the stresses are unbounded and depend linearly on the applied load. The picture however changes in plasticity, as one can expect the tensile stress to saturate, possibly below σ c , preventing crack nucleation altogether. [START_REF] Brink | Adhesive Wear Mechanisms in the Presence of Weak Interfaces: Insights from an Amorphous Model System[END_REF] have also shown that the resistance to shear of the contact junction plays a fundamental role in the wear particle formation. From a stress perspective, there is a competition between σ c and τ j (the junction shear resistance) for the formation of a crack: if the junction is strong enough, the maximum tensile stress may, under conditions depending on σ y , reach σ c and nucleate a crack. Conversely, if the junction is weaker, its slip may prevent the tensile stress from reaching σ c .

We investigate the interplay of these effects in the contact of a rigid self-affine rough surface with h rms = 0.1 and whose spectrum is defined with L/λ s = 16 , λ l = λ r = 8λ s , H = 0.8 , so that the surface is representative [START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF]. We have discretized the system with 729 × 729 points for the elastic and saturation models, and 729 × 729 × 32 for the 2 model. The applied normal mean pressure varies between 10 -2 h rms E and 8 • 10 -2 h rms E for the elastic case and between 10 -2 h rms E and 4 • 10 -2 h rms E for the saturation and 2 models. The yield stress is σ y = 10 -2 E and hardening modulus E h = 5 • 10 -2 E. Since we investigate the effect of the junction shear strength, we also apply a shear stress on the contacts at constant normal load. There is a simple correspondence between the applied shear stress and the junction strength τ j : any shear loading larger than τ j should not modify the stress state of the system since all contacts should be slipping and the stresses should not increase. We therefore interchangeably refer to the applied shear and the junction strength as τ j . In linear elasticity, the application of a constant shear stress on a patch of the surface creates a stress singularity at the edge of the patch [START_REF] Menga | The Surface Displacements of an Elastic Half-Space Subjected to Uniform Tangential Tractions Applied on a Circular Area[END_REF] because of the boundary condition discontinuity. In a physical system, a small amount of slip and rearrangement of the solids would occur at the edge of the contact junction so that the shear stress carried should be reduced on this zone. We therefore regularize6 the constant shear distribution over a transition zone of width ε τ small compared to the smallest surface wavelength (ε τ ≈ λ s /3 ), removing any numerical discrepancy due to the singularity.

We call a crack nucleation site a connected zone of the surface where the largest eigenvalue of the Cauchy stress tensor σ I is larger than σ c . Although the wear particle formation process is deterministic, the inherent randomness of the rough surfaces makes the process epistemically random [START_REF] Frérot | A Mechanistic Understanding of the Wear Coefficient: From Single to Multiple Asperities Contact[END_REF]. We apply a similar concept here and study the probability that a given contact nucleates a crack. However, because of the complex topography of the micro-contacts, there is no one-to-one correspondence between crack nucleation sites and contacts. This can be seen in fig. 5.4, where we show in grey the true contact area, in red the plastic zones and in black the crack nucleation sites (i.e. the zones where σ I > σ c ) for σ c = 0.1h rms E * and upwards applied shear stress τ j = 2 • 10 -3 E. Figure 5.4a shows the result for the saturation plasticity model and fig. 5.4b shows the 2 plasticity model. One can easily see on the latter that a single microcontact may have several crack sites. We also recognize the expected crescent shape crack nucleation on some contacts. Finally, the number of cracks is larger in the 2 model than in the saturation plasticity model. Since the contact areas predicted by both models are essentially the same, this difference can only be explained by the residual plastic deformations, which are not modeled in the saturation plasticity approach.

Because we still want to investigate a quantity akin to the probability that a contact nucleates a crack, we propose an adimensional measure called the crack nucleation likelihood (CNL) given by A 0 • n crack /A c , where A 0 = L 2 is the apparent contact area and n crack is the number of crack nucleation sites (i.e. the number of connected black zones in fig. 5.4). The CNL is conceptually a normalization of the number of crack nucleation sites by the density of contacts and relates to the probability of crack nucleation at a contact (this relationship will be detailed later on in this article). An "upwards" shear is applied on each contact. The true contact area is shown in light grey, the plastic zones in red and the crack nucleation sites in black. We can see that the 2 model has more crack nucleation sites than the saturation model. Since both models give approximately the same true contact area, this discrepancy must be due to plastic residual deformations which are not represented in the saturation approach. means that for a fixed value of σ c , the CNL increases exponentially when the load is increased.

The vertical lines indicate a quantity σ i/c , which is the stress for which the probability density function of σ I on the whole surface is maximum: in other words it is the most frequent stress value, and is typically found between contacts (hence the term "inter-contact stress"). The horizontal shift in the CNL curves corresponds to σ i/c . Since the latter depends on the spatial proximity of contacts, the CNL must depend on elastic interactions between contacts. Figure 5.5b shows that the two plasticity models have widely different behavior: the crack nucleation likelihood is much higher in the case of 2 plasticity, but also decays faster for values of σ c > σ y . It is surprising that despite being plastic with some hardening, the 2 model is more likely to lead to surface cracks than the saturation model. This may be a first step towards resolving the paradox highlighted in introduction.

Figure 5.6 shows the CNL when the junction strength increases (for the last normal load of fig. 5.5). Unlike previously, the elastic results are not shifted to the right when τ j increases but are instead scaled rightward. The same happens with the saturation plasticity model, whereas the 2 CNL seems relatively unaffected by τ j . This is due to plastic deformations having reached the contact surface (cf. fig. 5.4) and preventing an increase of σ I as fast as the elastic and saturation models, although hardening still allows some increase at a lower rate. In order to rationalize these results and provide evidence for the conclusions we have induced from our multi-asperity simulations, we now study the crack nucleation process for a single asperity. Although the initial slope is linear, the loading is not elastic as plastic deformations still evolve in the system. Their influence on the maximum tensile stress at the surface is however minimal. A sufficient shear stress may eventually create plastic strain at the surface, which causes a transition in the loading curve. Moreover, given a low enough yield stress, the initial indentation may cause surface yield, as seen for the curve where σ y /E = 2 • 10 -2 .

Single asperity crack nucleation

To investigate the effect of plasticity in the competition between σ c and τ j , we simulate a spherical indenter of radius R pushed onto an elastic-perfectly-plastic solid. The resulting contact junction is then subjected to a shear distribution (with ε τ = R/64 ), and the principal tensile stress σ I at the surface is recorded. Figure 5.7 shows the maximum tensile stress σ I as a function of the applied shear stress τ j across the contact. The different curves correspond to different yield stresses, with the dashed lines indicating the value of σ y for reference. Stresses here are normalized by the maximum tensile stress in Hertz contact σ Hertz I = (1 -2ν)p 0 /3 , with p 0 being the maximum hertzian contact pressure [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]. We can observe that the initial tensile stress (without applied shear) depends on the amount of plastic deformation: if σ y is decreased (or conversely the applied load increases), the initial tensile stress at the edge of the contact is higher. This is due to the residual stresses created by the plastic deformations that accommodate the indentation: the localized nature of the plastic strains causes the unloaded equilibrium position to not be stress-free. The additional stresses are tensile and add to the stress on the contact rim. In the case of σ y /E = 2 • 10 -2 , the plastic zone has reached the surface and the von Mises stress at the edge of contact has reached σ y by indentation alone (not shown here).

If the von Mises stress at the surface is below σ y , the application of a shear stress will cause an elastic loading phase, offset by the initial σ I value, as seen for the higher values of σ y /E. The loading continues until σ I reaches values close to σ y , as indicated by the dashed lines. Since the stress state at the edge of contact is triaxial, the maximum value σ I can reach is not σ y , as is seen for the most plastic case. After a certain point, eq. (5.7) becomes unsolvable because a plastic failure mechanism develops [START_REF] Drucker | Soil Mechanics and Plastic Analysis or Limit Design[END_REF]: we supposed that further loading will not increase the value of σ I . Of course, for a hardening material σ I should not saturate and instead increase further at a lower rate.

Figure 5.8 shows the competition between the junction strength τ j and the tensile strength σ c for different values of σ y (first row no hardening, i.e. E h /E = 0 ) and E h (second row with σ y /E = 2 • 10 -2 ) the hardening modulus. The zones below the white curve show for which values of (σ c , τ j ) a crack may nucleate, and the zones above show when the interface breaks (i.e. slips) before crack nucleation. The dashed lines indicate the values of σ y . We can see that for materials with high yield there is an affine boundary between the two mechanisms, which shows their competition. Of course materials with σ c < σ y are not plastic, so such transition can only happen for brittle materials. It does however exist for hardening materials, as opposed to perfectly plastic ones. With no hardening, the failure mechanism is purely determined by the value of σ c /σ Hertz I , which depends on the applied load. With this single-asperity analysis, we have explained why the crack nucleation likelihood is higher in the 2 model for no applied shear stress: the plastic residual deformations cause tensile stresses which combine with the contact stresses and increase σ I , thus increasing the CNL. This does not occur in the saturated model because it ignores plastic residual deformations in the stress computation. We have also explained why the 2 CNL is relatively insensitive to the applied shear/junction strength: plasticity has a saturation effect on σ I : when the system is sheared, the rim of contacts is in the plastic regime, and the increase in σ I is purely driven by hardening, which in the case of fig. 5.6 is only 5% of the Young's modulus. In order to rationalize the other aspects of the CNL highlighted by figs. 5.5 and 5.6, we resort to a statistical model for multi-asperity contact.

Multi-asperities

One can now apply a statistical approach to estimate the proportion of contacts that nucleate cracks in a multi-asperity setting. We thereafter use a Greenwood-Williamson (GW) model [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF] with an exponential distribution of asperity heights to obtain simple analytical results. Since the asperities are randomly distributed (all with the same radius R), σ Hertz I becomes in turn a random variable. There is however a significant difference between the single asperity case we have studied and the multi-asperity setting. Because of elastic interactions, the tensile stress at the edge of a contact depends on the proximity and magnitude of the neighboring contacts. In a traditional GW approach, contacts are independent of each other. We assume this is the case, but that the stress state is determined by the local contact with an additional contribution σ i/c , the inter-contact stress, determined from the neighboring contacts. The radial stress outside the area of a single contact of radius a is given by [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]:

σ r (r ) = 1 -2ν 3 • a 2 r 2 p 0 = κ sR r 2 (z * -h) 3 2 ,
(5.9)

with κ := 2(1 -2ν)E * 3π • s R , (5.10)
where r is the euclidean distance from the contact center, z * = z/s is the asperity height random variable normalized by the standard deviation of heights s and h is the normalized surface approach. We assume a spatial asperity density η and a contact density η c = ηe -h = A c /(A 0 πsR) [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF], with A c being the true contact area. To compute the inter-contact stress σ i/c , we assume the stress state outside each contact is given by the mean of eq. (5.9), averaged over contacting asperities, and we compute the largest stress eigenvalue at the center of a series of circles whose diameters are multiples of d c = 1/ √ η c which is the characteristic distance between contacts. The calculation process is detailed in Appendix D.1 and leads to the following expression: (5.11) with ξ ≤ 3ζ (3) ≈ 1.9 and ζ is the Riemann Zeta function. We use this upper bound for σ i/c in the rest of this work. We can see that σ i/c depends linearly on the contact ratio, and therefore is linear with the load. We suppose that σ i/c acts as a "background" stress, and that the maximum tensile stress is the sum of the local contact tensile stress σ Hertz I and the inter-contact stress. We can therefore quantify the probability that a contact nucleates a crack:

σ i/c = 3κ ξ √ π • A c A 0 ,
P crack = P σ c σ Hertz I + σ i/c < ω(τ j ) z * -h ≥ 0 = exp - σ c /ω(τ j ) -σ i/c κ 2 ,
(5.12)

where ω is the function describing the failure mechanism transition (white line in fig. 5.8).

The calculation details of eq. ( 5.12) can be found in Appendix D.2.

Remark 7. The crack nucleation likelihood A 0 • n crack /A c and P crack are related in our GW approach: indeed n crack = P crack η c A 0 = P crack A c /πsR, cf. (Greenwood and Williamson, 1966, p. 303).

Remark 8. The quantity σ c /ω(τ j ) -σ i/c strikingly explains the features of figs. 5.5 and 5.6 for the elastic model. In fig. 5.5, the curves are shifted to the right as the load increases, which is apparent in eqs. (5.11) and (5.12): σ i/c increases linearly with the load, and thus causes a rightward shift in the graph of the CNL. Similarly, as ω is linear in τ j , P crack is scaled horizontally, which can also be seen in the CNL on fig. 5.6.

Remark 9. When σ c /ω(τ j ) = σ i/c the probability is one, meaning that all contacts, regardless of size, nucleate cracks; in other words the whole surface should be cracking in a catastrophic breakdown. This does not happen in practice, as the normal loading process should nucleate and propagate cracks at single asperities before the breakdown is reached, thus relaxing the tensile stresses in the system.

Comparison to a rough surface

We wish to assess the validity of the above developments with simulations of self-affine rough surface contact. Because of the simplifying assumptions of a GW model, we do not hope to establish a quantitative agreement, especially since the asperity curvature is not unequivocally defined on a self-affine rough surface [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF]. Instead, we will focus on the qualitative relations between P crack , σ c , σ i/c and τ j highlighted above.

Elasticity Results

We first consider σ i/c for a rough surface. Recall that the inter-contact stress in a rough contact is the most frequent value of σ I , i.e. the value for which the probability density function p σ I of the surface tensile stress is maximum. This is illustrated in the inset The stress value corresponding to the peak in the probability density is defined as the inter-contact stress σ i/c . We see that in the rough surface simulation and the Greenwood-Williamson model σ i/c evolves linearly with the contact ratio. The value of κ is estimated from the mean curvature of contacting summits in the rough surface.

of fig. 5.9. In the latter, we plot the evolution of σ i/c for the elastic rough contact defined previously and for eq. ( 5.11). We can observe that both curves behave linearly with the contact area ratio, with different slopes. To compute the value of κ, which depends on s/R, cf. eq. (5.10), we have used Nayak's approach [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF] to estimate the mean curvature radius of the zones of the rough surfaces in contact, i.e. R = 3/m 4 /I (z * c ), where m 4 is the fourth moment of the surface spectrum, z * c is the normalized height of the surface in contact and I is a function defined in [START_REF] Nayak | Random Process Model of Rough Surfaces[END_REF]. As can be seen in fig. 5.9, the slope of the GW curve is approximately constant, showing the weak dependency on z * c . Note that although the values of R and A c /A 0 for the GW model are informed from the rough surface simulation, there is no fit parameter to the GW prediction. The lack of quantitative agreement between the theoretical approach and the rough contact simulation shows the prediction limit of GW. Using a more elaborate model like that of [START_REF] Bush | The Elastic Contact of a Rough Surface[END_REF] could improve the agreement, but as our aim is to establish a qualitative prediction we restrict the results presented to GW.

As previously mentioned, P crack is not directly measurable on a rough contact interface (cf. fig. 5.4). However, the crack nucleation likelihood acts as an alternative measure for P crack , cf. Remark 7. Figure 5.10a, confirms that this is indeed the case and we find the squared exponential dependency predicted by eq. (5.12), with all curves collapsed due to the shift caused by σ i/c (recall that ω(0) = 1 ). Figure 5.10b, on the other hand, shows that the CNL does not follow eq. (5.12) for non-zero τ j . While each curve remains close to a straight line, they do not overlap, but seem to converge to a master curve. More strikingly, the CNL decreases (5.12). When the strength of the junction is taken into account (figure (b)), or, equivalently, if a shear stress is applied, we observe qualitative deviations from eq. (5.12). There is a decrease in crack event density due to the interference of close contacts, which tends to unload the tensile stresses at the trailing edge of leading contacts.

as τ j increases, meaning that ω(τ j ), which was computed directly from the data of fig. 5.7, over-normalizes the data. This is again due to interactions between asperities. For our "singleasperity" analysis, because of periodicity, we in fact consider many interacting asperities on a square lattice, each separated by a distance L. When shear is applied, a positive σ I is created at the trailing edge of the contact and a negative σ I appears at the leading edge. Because the periodic images are equidistant and far apart, they weakly affect the stress distribution in the vicinity of the contact. However, when two contacts are close to each other, creating local anisotropy, the inter-contact stress distribution of each asperity is compensated by the other. The positive peak in σ I at the trailing edge of one contact is then reduced, thereby reducing P crack as seen in fig. 5.10b. This phenomena is akin to the crack shielding mechanism uncovered by [START_REF] Aghababaei | Asperity-Level Origins of Transition from Mild to Severe Wear[END_REF].

Plasticity Results

As for the elastic model, we compare the inter-contact stress computed from the plastic rough contact simulations to our GW approach. models still give a linear dependency of σ i/c on the contact ratio, although some more data may be required to draw an affirmative conclusion in this regard. Finally, fig. 5.12 shows the data of fig. 5.5b normalized to compare the results to eq. (5.12). While the saturation model seems to follow our GW prediction (which is based on elasticity assumption), it is clear that the 2 model does not conform to our scaling predictions for P crack . However, in light of fig. 5.11, it is interesting to note that although σ i/c is lowest for the 2 approach (indicating less interactions between contacts), the latter has the largest crack nucleation likelihood, because of the local effect of plastic residual deformations. As shown in fig. 5.7, this local effect of plastic deformations is stronger the more ductile a material is, as expected from experimental data which shows that softer materials wear more.

Conclusion

We have investigated in this work the nucleation of cracks at an elastic-plastic rough contact interface. This was motivated by the necessity for an accurate description of the process of crack nucleation for adhesive wear that includes plasticity. By comparing a classical 2 plasticity model and a saturation plasticity approach commonly used in tribology, we have concluded that the latter can only qualitatively reproduce the true contact area and fails to give satisfactory results on local quantities. This can be seen in the crack nucleation likelihood, which is much higher in the 2 approach. This is caused by plastic residual deformations which increase tensile stresses in the vicinity of contacts. In this regard, the saturation model is not applicable to study crack nucleation because it does not capture plastic deformations. We show with a single asperity analysis that more ductile materials can have larger surface tensile stresses and nucleate more cracks at the interface.

We have also showed that elastic interactions play a role in the crack nucleation likelihood of a single contact. They may increase the latter through proximity of contacts, or decrease it in shearing by elastic shielding. This was further supported by an analytical approach based on a Greenwood-Williamson model modified to take interactions into account.

Conclusion

This thesis has been concerned with the modeling of realistic rough surface contact and its applicability to up-scale asperity level observations on wear and friction. It has specifically focused on the mathematical aspects of efficiently characterizing rough surface contact between materials with elastic-plastic constitutive laws. Natural and man-made surfaces are seldom flat, and roughness plays a crucial role in phenomena such as friction [START_REF] Bowden | The Area of Contact between Stationary and between Moving Surfaces[END_REF], adhesion [START_REF] Pastewka | Contact between Rough Surfaces and a Criterion for Macroscopic Adhesion[END_REF], sealing [START_REF] Shvarts | Trapped Fluid in Contact Interface[END_REF], electric resistance of contacts [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF] and wear [START_REF] Archard | Contact and Rubbing of Flat Surfaces[END_REF]. Developing simulation tools for such problems is challenging because of the inherent non-linear material behavior and the multi-scale nature of rough surfaces [START_REF] Renard | Constant Dimensionality of Fault Roughness from the Scale of Micro-Fractures to the Scale of Continents[END_REF][START_REF] Persson | On the Nature of Surface Roughness with Application to Contact Mechanics, Sealing, Rubber Friction and Adhesion[END_REF].

To allow the efficient resolution of elastic-plastic rough contact problems, I devised a novel numerical approach for volume integral methods (cf. Chapter 2). These methods rely on fundamental solutions to construct, by superposition, stress and displacement fields that satisfy the required boundary conditions, kinematic compatibility and equilibrium equations [START_REF] Bonnet | Boundary Integral Equation Methods for Solids and Fluids[END_REF]. The breakthrough in computational efficiency presented in Chapter 2 was achieved with a novel formulation of the required fundamental solution directly in the Fourier domain, which drastically saves on memory resources as well as computing time. This novel method is able to handle the discretization requirements of rough surfaces while including a proper plasticity constitutive law [START_REF] Simo | Computational Inelasticity[END_REF]. As a byproduct of this research, I have implemented the method in a high-performance hybrid C++/Python open-source code called Tamaas, which has enabled me to conduct elastoplastic simulations of rough contact with upwards of 100 million degrees of freedom (729 × 729 × 32 × 6 components of strain) on a single compute node, which to my knowledge is an order of magnitude larger than simulations of the same kind reported in the literature [START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF]Zhou, Jin, Wang, Wang, et al., 2016).

To improve the performance of the fixed point algorithm developed by [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF], I have explored an interior-point approach to solve both the plasticity and the contact problems simultaneously (cf. Chapter 3). The proposed primal-dual interior point approach solves the conic program underlying the elastic-plastic contact formulation using the volume integral operators laid out in Chapter 2. This is however work in progress, as the efficient solving of the underlying linear system of equations is an active research topic in the optimization community.

To study the applicability of rough surface contact models in the rationalization of macroscopic wear models, I have adopted in Chapter 4 a simple approach to single-asperity wear: a contact of area A can only form a wear particle if A > A * ∝ (d * ) 2 in similar fashion to [START_REF] Aghababaei | Critical Length Scale Controls Adhesive Wear Mechanisms[END_REF]. I have shown that although the particle forma-tion process is deterministic at the single-asperity level, the wear coefficient, which Archard interprets as the probability of particle formation, can be understood at the multi-asperity level (i.e. macro-scale) as K = P(A > A * ), which is the probability that a contact A is greater than A * . With this interpretation, the wear coefficient is shown to remain constant with the applied load, as expected from experiments. However, this interpretation of the wear coefficient relies on assumptions that were not justified by Archard. I have therefore proposed a new interpretation of the wear coefficient which does not rely on Achard's assumption an directly up-scales asperity-level wear mechanisms.

Finally, I have used the elastic-plastic contact method developed in Chapters 2 and 3 to build and understanding of the crack nucleation process in the adhesive wear of a rough contact interface. In particular, I have showed that the density of contacts, and hence the density of elastic interactions, plays a key role in the number of crack nucleation sites. Using analytical arguments based on the statistical approach of [START_REF] Greenwood | Contact of Nominally Flat Surfaces[END_REF], I accurately predicted the dependency of the crack nucleation density on the normal load and the critical tensile stress in normal contact. With added shear stress, I have rationalized the departure from expected behavior with elastic shielding effects which prevent crack nucleation in-between contacts that are in close proximity. I have also showed that plastic residual deformations have a crucial influence of the stress state in the neighborhood of a micro-contact by comparing the results of the methods developed in this thesis with a saturation pressure approach, which qualitatively describes the relationship between the normal load and the contact area but fails to reproduce the stress state obtained with the more physically accurate von Mises plasticity.

Outlook

This thesis has developed two major axes of research which have been clearly outlined by the two parts of the present manuscript, and therefore two major topics for future research endeavors arise. The first topic is concerned with the improvement of the numerical methods that have been developed. As mentioned in Chapter 3, an ongoing collaboration aims to develop a more efficient solving scheme based on a primal-dual interior-point method and on the integral operators detailed in Chapter 2. While this method may indeed yield a better algorithm, more optimal approaches may be developed in the future. This is however not the only aspect that could be improved. We have proposed in Chapter 2 two algorithms for the integration of the remaining physical-space integral in the volume integral operator N . While the second algorithm has a linear complexity (as opposed to the quadratic complexity of the first), the numerical evaluation of the exponential terms may overflow or underflow the machine representation of floating point numbers because of the separation of variables that was operated (cf. Appendix A.5). Other approaches for the integration may be taken to alleviate this issue while retaining the complexity advantages (e.g. multi-level methods). On a more technical note, the present implementation of the integration method is sequential, and an investigation into the parallelizability of the different integration algorithms should also factor in the choice of integration approach.

As we have shown in this manuscript, the methods developed in the first half of this thesis are tools that can be used to efficiently solve contact problems with plastic constitutive behavior. However, some parts of these methods are not restricted to that specific class of problems: in particular, the volume integral operators are mathematical entities used to construct displacement and strain fields satisfying equilibrium and elastic behavior, which gives them a wider range of applications than what we show here. In fact volume integral methods have been used to model damage around stiff inclusions in contacting bodies [START_REF] Beyer | A Coupled Damage Model and a Semi-Analytical Contact Solver to Simulate Butterfly Wing Formation around Nonmetallic Inclusions[END_REF], and could be used in conjunction with gradient-based damage models, whose governing equation can also be solved with Fourier-based volume integral operators. Developing such a numerical method could be instrumental to the modeling of wear particle formation in rough surface (sliding) contact. This manuscript has brought in focus some key issues of wear modeling: the necessity of understanding wear laws as emerging from the interaction of micro-scale mechanisms and that the true contact area, while being a quantity of interest, is not sufficient to describe wear. These observations are also true of many other interface phenomena, including friction, sealing, etc. While it seems the major concern of the rough contact community at large has been the true contact area, I think future pursuits in tribology should consider the contact interface as a whole. For example, the experiments on friction conducted by Fineberg's group show that the fracture energy of a PMMA frictional interface is comparable to the bulk fracture energy [START_REF] Svetlizky | Classical Shear Cracks Drive the Onset of Dry Frictional Motion[END_REF][START_REF] Bayart | Slippery but Tough: The Rapid Fracture of Lubricated Frictional Interfaces[END_REF], indicating some level of extreme material behavior at the contact interface that cannot possibly be captured by the contact area alone. The simulation code that I have developed in my thesis is unique in its kind to help shed light on these issues, as well as provide insight into other applications of rough contact, such as vibrating systems, for which the dynamical response depends on the contact interface [START_REF] Gallego | A Comprehensive Method to Predict Wear and to Define the Optimum Geometry of Fretting Surfaces[END_REF][START_REF] Armand | On the Effects of Roughness on the Nonlinear Dynamics of a Bolted Joint: A Multiscale Analysis[END_REF] or sealing, where the percolation point of the contact interface is of prime importance [START_REF] Dapp | Self-Affine Elastic Contacts: Percolation and Leakage[END_REF][START_REF] Shvarts | Trapped Fluid in Contact Interface[END_REF] and plastic deformations are expected. New experimental methods are able to track the evolution of the rough contact interface during sliding [START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF][START_REF] Garabedian | Quantifying, Locating, and Following Asperity-Scale Wear Processes Within Multiasperity Contacts[END_REF], and the tools developed during this thesis should be instrumental in the rationalization of the refined experimental observations. Finally, bridging scales in friction can be attempted on the basis of this thesis. The issue in this domain is the choice of friction law at the interface: Coulomb friction is inherently ill-posed and should physically make sense only at the macro-scale, as like wear it is a phenomenon emerging from the micro-scale interactions. The tools developed in this thesis can be employed to study these interactions with great accuracy, but methods incorporating non-linear geometry may be required to gain sufficient understanding [START_REF] Caroli | Hysteresis and Elastic Interactions of Microasperities in Dry Friction[END_REF]. Micro-contact lifetimes when two rough surfaces are sliding are important to rationalize the state variable in modern friction models. Further investigation in this area may reinforce or challenge the current interpretation of the state variable in modern rate-dependent friction models as a "memory of contact", and perhaps give further leads as to how these friction models emerge from the contact interface. Finally, frictional slip events that are well described by linear-elastic fracture mechanics at the macro-scale are actually composed of micro-contacts breaking and interacting via elastic waves. Although we have provided a static picture of the contact interface in this thesis, it could be extended to include dynamic effects, for example by support of ϕ j on the interval [x j-1 , x j+1 ] (e.g. P1 shape functions). Then we have the following properties:

∫ x l 0 - k (q )ϕ j ( ) d =                    0 j > l ∫ x l x l -1 - k (q )ϕ j ( ) d j = l ∫ x j +1
x j -1 k (q )ϕ j ( ) d j < l for l > 0,

∫ ∞ x l + k (q )ϕ j ( ) d =                    0 j < l ∫ x l +1 x l + k (q )ϕ j ( ) d j = l ∫ x j +1
x j -1 + k (q )ϕ j ( ) d j > l for l < N 3 -1.

These allow the simplification of the sums, as some terms cancel: The same can be done for the second sum, and we have now eliminated the dependency of the sum terms on l, yielding the expression in Theorem 4.

N3-1 j=0 w j ∫ x l 0 - k (q )ϕ j ( ) d = w 1 ∫ x1 x0 - k (q )ϕ 1 ( ) d + l -1 j=1 w j ∫ xj+1 xj-1 - k (q )ϕ j ( ) d + w l ∫ x l x l -1 - k (q )ϕ l ( ) d = w 1 ∫ x1 x0 - k (q )ϕ 1 ( ) d + w l ∫ x l x l -1 - k (q )ϕ l ( ) d

B.2. Conic programming

with the block decomposition of E:

E = E 00 E T E E 11 .
One can now introduce δx 0 and δs in the second optimality equation: Finally, we can replace δ n+1 in the third optimality condition:

-N 3 CDδx + δ n+1 -M 33 δp + δU = -R 2 (z k ), ⇒ -N 3 CDδx -M 33 + GP -1 δp + δU = -R 2 (z k ) - µ p k n+1 + k n+1 .
We can now write the reduced linear system:

         DCRD -E 11 + E ⊗ E 2 3 E h + E 00 DCM 0 -N 3 CD -M 33 + GP -1 a 0 a T 0                   δx δp δU          = -R(z k ), (B.5)
where R is the reduced residual formed from the previous substitutions:

R(z k ) =                  R 1 (z k ) -s + µX -1 e - E 2 3 E h + E 00 R 0 (z k ) -s 0 + µX -1 e 0 R 2 (z k ) + µ p k n+1 -k n+1 R 3 (z k )                  .
To complete the formulation, we can explicit the expression X -1 e. We define the function det(x) := (x 0 ) 2 -x 2 , which helps expressing the inverse of X (Bleyer, 2018):

X -1 = 1 det(x)       x 0 -x T -x det(x)I m + x ⊗ x x 0       , X -1 e = x det(x)
with x = (x 0 , -x).

It remains to express E. [START_REF] Bleyer | Advances in the Simulation of Viscoplastic Fluid Flows Using Interior-Point Methods[END_REF] gives the expressions for F , which we develop here for E = F 2 :

E = θ 2 (-Q + 2w ⊗ w), θ = 4 det(s) det(x) , w = θ -1 s + θ x √ 2 • x • s + det(s) det(x) , Q = diag(1, -I m ).
The sub-parts of E can be expressed individually:

E 00 = θ 2 2w 2 0 -1 , E = 2θ 2 w 0 w, E ⊗ E = θ 2 w 0 2 w ⊗ w, E 11 = θ 2 (I m + 2w ⊗ w) .
Preconditioning Although we cannot proceed to an incomplete LU factorization on the reduced linear operator of eq. (B.5) because we wish to keep a matrix-free iteration scheme, we can explicitly invert the components of the reduced tangent that cause the condition number to diverge as µ → 0 , i.e.:

         E 11 - E ⊗ E 2 3 E h + E 00 0 0 0 -GP -1 a 0 a T 0         
.

Using the Sherman-Morrison formula, the first block can be inverted (β := (2/3 • E h + E 00 ) -1 ):

E 11 + βE ⊗ E -1 = E -1 11 -β E -1 11 E ⊗ E E -1 11 1 + βE • E -1 11 E , E -1 11 = θ -2 I m -2 w ⊗ w 1 + w • w .
The second block can also be inverted:

GP -1 -a T -a 0 -1 = A + 1 k Aaa T A 1 k Aa 1 k a T A 1 k A = GP -1 -1 = PG -1 k = -a T Aa
C.1.3. Wear coefficients based on up-scaling approach

D ≡ ∫ ∞ 0 Ap PL (A, A m ) dA = C ∫ A s 0 A dA + ∫ A m A s A A A s -α dA = C A s + A α s 2 -α (A 2-α m -A 2-α s )
(C.4) Grey curves are results for each realization (20 in total), and red curves are ensemble averages. Surface sampled has λ l /λ s = 8 and λ s /∆l = 16 , so that only L is varied. As L/λ l increases, the dispersion of the individual realizations decreases, and the behavior of the ensemble average approximates accurately the behavior of each realization. We also note that as L/λ l increases, the distributions converge to a limit distribution.

K PL = 1 D ∫ ∞ A * Ap PL (A, A m ) dA =                  C D ∫ A s A * A dA + ∫ A m A s A A A s -α dA A * ∈ [0, A s ] C D ∫ A m A * A A A s -α dA A * ∈ [A s , A m ] 0 A * ∈ [A m , ∞] =                    1 - (1 -α 2 )A -α s (A * ) 2 A 2-α m -α 2 A 2-α s A * ∈ [0, A s ] 1 - (A * ) 2-α -α 2 A 2-α s A 2-α m -α 2 A 2-α s A * ∈ [A s , A m ] 0 A * ∈ [A m ,
For the simulations carried out in the rest of this paper, a value of L/λ l = 8 , which is a good compromise between required number of realizations and computational cost, was selected. Note: we analyze the cumulative distribution instead of probability density to remove any bias due to binning. Curves are ensemble averages of 20 realizations of a surface with L/λ l = 16 and λ l /λ s = 8 , so only ∆l is varied. As λ s /∆l is increased, the distributions converge to a limit distribution. Increasing λ s /∆l smoothens the distributions and reduces the systematic bias between the computed distribution and the limit of λ s /∆l → ∞. For simulations carried out in the main paper, a value of λ s /∆l = 8 was selected, offering a reasonably low discretization bias, and making the simulations possible with our current code. We first compute σ r (r ), which is the average radial stress outside contacts:

D. Appendix to

σ r (r ) = 1 P(z * ≥ h) ∫ ∞ h κ sR r 2 (z * -h) 3 2 e -z * dz * = κ 3 √ π 4 • sR r 2 .
We then suppose for simplicity that all contacts have the same radial stress σ r (r ). Accordingly, their hoop stress is σ θ (r ) = -σ r (r ) [START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF]. We assume that all contacts are spatially uniformly distributed with density η c , so that the characteristic distance between contacts is d c = 1/ √ η c . We divide the infinite surface into concentric rings of width d c and diameters d i ∈ {d c , 2d c , 3d c , . . .}. Each ring can be reduced to a circle of contacts with linear density √ η c . We now wish to compute for the sum of all circles of diameter d 1,2,... the largest stress eigenvalue at the center. For a single contact positioned at an angle θ on a circle of diameter d i , the stress state in Cartesian coordinates is:

σ = σ r (d i /2)
cos(2θ )sin(2θ ) sin(2θ )cos(2θ ) .

For the circle number i, the expected number of contacts is πid c √ η c = i • π , but we simplify by assuming the expected number to be n i := 3i. The total stress state at the center, summing all contacts per circle and all circles:

σ = ∞ i=1 σ i r n i k =1 cos(2θ k ) -sin(2θ k ) -sin(2θ k ) -cos(2θ k ) ,
with σ i r := σ r (id c /2). The largest eigenvalue of σ is given by:

λ 2 = ∞ i=1 σ i r n i k=1 cos(2θ k ) 2 + ∞ i=1 σ i r n i k=1 sin(2θ k ) 2 = ∞ i=1 (σ i r ) 2 (c 2 i + s 2 i ) + ∞ i <j σ i r σ j r (c i c j + s i s j ),
where:

c i := n i k=1 cos(2θ k ), s i := n i k=1 sin(2θ k )
The angular position of each contact is assumed to be uniformly distributed in [0, 2π ]. The expected value of λ, which gives the inter-contact stress, is given by:

σ i/c = E[λ] ≤ E[λ 2 ] = ∞ i=1 (σ i r ) 2 E[c 2 i +s 2 i ] + ∞ i <j σ i r σ j r E[c i c j +s i s j ].
where we have used Jensen's inequality for a simple estimation. We note that:

c 2 i + s 2 i = n i k=1 cos(2θ k ) 2 + n i k =1
sin(2θ k )

σ i/c ≤ ∞ i=1 (σ i r ) 2 n i ≤ 3κ √ π sR d 2 c ∞ i=1 3 i 3 ≤ 3κ √ πsRη c 3ζ (3)
where ζ is the Riemann Zeta function. We can now use the GW contact model to replace sRη c = sRηe -h = A c /(πA 0 ):

σ i/c ≤ 3κ 3ζ (3) π • A c A 0 .
Note that only the 3ζ (3) term depends on the estimation from Jensen's inequality, so σ i/c is indeed linear with respect to the contact ratio.

D.2. Nucleation probability

P crack as defined in eq. (5.12) is a conditional probability. It expresses the question "knowing an asperity is in contact, what is the probability that a crack nucleates at the contact edge?" The final expression for this probability is obtained by manipulating the inequality. We evaluate σ r at the contact radius a = R • s √ z * -h:

σ Hertz I = σ r (R • s √ z * -h) = κ(z * -h) 1 2 ,
which is replaced in the inequality:

σ c σ Hertz I + σ i/c ≤ ω(τ j ) ⇔ σ Hertz I + σ i/c σ c ≥ 1 ω(τ j ) ⇔ σ Hertz I ≥ σ c ω(τ j ) -σ i/c ⇔ z * ≥ σ c /ω(τ j ) -σ i/c κ 2 + h.
Let us call X the event corresponding to the above inequality. We have:

1 .

 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Equilibrium of elastic bodies: problem introduction . . . . . . . . . . . . 1.1.

1 .

 1 1 whereas K plateaus at values in ]0, 1[. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.Pressure profiles for 2 plasticity (a), saturation (b) and elasticity (c) models. Although fig.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.6. Crack nucleation likelihood (CNL) as a function of σ c and junction strength. As for fig. 5.5, (a) shows the elastic model and (b) the plastic models. Unlike fig. 5.5a, the CNL for the elastic model is scaled to the right when the junction strength increases. The same can be said of the curves corresponding to the saturation plasticity in (b), but not of the 2 curves, which are relatively insensitive to changes in τ j . This is due to plasticity preventing increases in σ I . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 xv 5.7. Maximum tensile stress as a function of applied shear stress across single asperity contact with 2 plasticity.

  For simulations carried out in the main paper, a value of λ s /∆l = 8 was selected, offering a reasonably low discretization bias, and making the simulations possible with our current code. . . . . . . . . . . . . . . . . . . . . . .

Figure 0

 0 Figure 0.1.: Illustration of the contact area evolution. Contact area between elastic solids with self-affine rough surfaces. Darker shades correspond to larger applied normal loads. The complex contact topography and the associated stress state governs virtually all tribological phenomena.

Figure 0

 0 Figure 0.2.: Illustration (3D) of the sub-surface plastic zones. Lighter shades are zones deeper below the contact interface.Because the contact load must be carried over an area typically two orders of magnitude smaller than the apparent area of contact, the local pressures may cause plastic activity. The principal shear stress is maximum below the surface[START_REF] Johnson | The Contact of Elastic Regular Wavy Surfaces[END_REF], so plastic deformations first develop below contacting points, then grow in volume and eventually reach the surface. It is common that disjoint surface contacts "share" a plastic region below the surface. These plastic strains modify significantly the stress state compared to an elastic material, which drastically changes the macroscopic picture of stressdriven phenomena such as wear. Note that this figure is the result of an actual elastic-plastic contact simulation. The original method used to obtain it is detailed in Chapter 2 of this manuscript.

Figure 1 .

 1 Figure 1.1.: Generated self-affine rough surface. Color code indicates the asperity height. The "fractal" aspect of the roughness is represented by small asperities sitting on top of larger ones.

) 1 .

 1 Figure 1.2.: Histogram comparison of rough surface generation algorithms. Figure (a) showsthe root-mean-square of heights, with the average value in dashed line and the colored zone containing ±σ around the average. The random phase algorithm by[START_REF] Wu | Simulation of Rough Surfaces with FFT[END_REF] gives consistent values around the analytical prediction, while the algorithm of[START_REF] Hu | Simulation of 3-D Random Rough Surface by 2-D Digital Filter and Fourier Analysis[END_REF] gives a much larger spread (approx. 10 times).Figure (b) shows the root-mean-square of slopes. Here the random phase algorithm does not reproduce the analytical value. This is due to the discrepancy between the continuous and discrete version of the PSD which is increased by the moment computation of eq.(1.47). This error decreases as the ratio L/λ l increases.

Figure 1

 1 Figure 1.3.: Measures of the root-mean-square of slopes. Due to the finite difference approximationof the surface gradient, the geometric estimator is dependent on discretization and is therefore not a good measure of the root-mean-square of slopes. The spectral estimator computes the discrete version of eq. (1.47) and does not depend on discretization if the shortest wavelength is larger than twice the sampling size (i.e. if the PSD can be represented in the Fourier domain).

  2.6. Algorithmic complexity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

Figure 2 .

 2 Figure 2.1.: Schematic representation of a periodic elastic-plastic contact problem, with geometrical contact quantities represented by blue arrows and plastic deformation zones with dashed red contours. The elastic-plastic body B is represented in deformed shape. Note that this schematic is taken from a real two-dimensional simulation.

-Figure 2 . 2 .:

 22 Figure 2.2.: Normalized surface displacements due to a hydrostatic spherical inclusion in a half-space. Comparison between the Mindlin and Cheng (1950) solution and the proposed VIM.The displacements are shown along the x 1 = x 3 = 0 line. The agreement between the analytical and computed solutions is very good in the central part of the periodicity cell, while distortion induced by periodic conditions becomes apparent near its boundary.

2 +

 2 c 2 and β = αT 1+ν 1-ν . For this example, we use c = 2a and ν = 0.3 , while the periodicity cell B p and its discretization N are defined by B p =] -15a, 15a[ 2 ×]0, 10a[ and N = (128, 128, 126). Figures 2.2a and 2.2b show the horizontal displacement u r and the vertical displacement u z , respectively, evaluated along the x 1 = x 3 = 0 line. The computed and reference displacements are in good agreement in the central zone, the expected distortion induced by the periodic boundary conditions becoming apparent only close to the boundary of B p . Recalling now that σ = C : ∇N [w] -w, fig. 2.3 shows the validation of ∇N by way of a comparison of stresses σ θ , and σ z produced by the inclusion, whose values along the vertical line r = 0 going through the inclusion center are shown in figures 2.3a and 2.3b, respectively.

Figure 2 .

 2 Figure 2.3.: Normalized stresses due to a hydrostatic spherical inclusion in a half-space. Comparison between the Mindlin and Cheng (1950) solution and the proposed volume integral method. The stresses are shown along the x 3 axis. The red region is where the eigenstrain ε = αT I is imposed.One can observe a good agreement of the numerical approximation with the analytical solution outside the inclusion. Some Gibbs effect can be observed at the boundary of the inclusion, with spurious oscillations in the inclusion. This is due to the Fourier approximation of the discontinuous eigenstrain function.

Figure 2

 2 Figure 2.4.: Relative computation times for VIM and FEM. We compare the application of the operator N to an elastostatic FEM solve step (Cholesky factorization). The reference time T 0 is T 0 = 1.41 • 10 -3 s for the VIM and T 0 = 1.77 • 10 -2 s for the FEM. The scaling for large problem sizes agrees with the theoretical algorithmic complexities. For N = 2 21 the stiffness matrix factorization needed over 128 GB of memory, an amount two orders of magnitude larger than what is required for the VIM.

  Figure 3.1.: Elastic perfectly-plastic Hertzian contact, comparison with Hardy, Baronet, and Tordion(1971). Increase of the applied load beyond the initial yield shows that the pressure profile deviates from the elliptic Hertzian profile by flattening of the curve at the axis of symmetry, with a plateau whose extent increases with the load. The results of[START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] however show oscillations of the pressure profile at high plastic loads which is not reproduced by our simulation. As there is, to our knowledge, no physical reason to these oscillations, they are likely due to the coarse discretization of the finite-element mesh they used in their study 1 . The simulation parameters are given in Appendix B.1.1

1 .

 1 Figure 3.1.: Elastic perfectly-plastic Hertzian contact, comparison with Hardy, Baronet, and Tordion(1971). Increase of the applied load beyond the initial yield shows that the pressure profile deviates from the elliptic Hertzian profile by flattening of the curve at the axis of symmetry, with a plateau whose extent increases with the load. The results of[START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] however show oscillations of the pressure profile at high plastic loads which is not reproduced by our simulation. As there is, to our knowledge, no physical reason to these oscillations, they are likely due to the coarse discretization of the finite-element mesh they used in their study 1 . The simulation parameters are given in Appendix B.1.1

Figure 3

 3 Figure 3.2.: Elastoplastic Hertzian contact. Comparison between the elastic Hertz solution (Johnson, 1985), a simulation using Algorithm 5 where the surface residual displacement is calculated by FEM, and a full VIM simulation. Both simulations have identical surface discretization. As expected (Johnson, 1968), the plastic pressure distribution deviates from the elliptical shape of the Hertzian distribution. The contact radius is larger in the plastic case. The simulation parameters are given in Appendix B.1.2.

Figure 3

 3 Figure 3.3.: Plastic strain norm |ε p | in elastoplastic Hertzian contact. Results obtained with Algorithm 5, with the plastic problem solved using our VIM approach (left) or a first order FEM (right). Note that in the FEM case, plastic deformations are piece-wise constant, whereas in the volume integral result they are interpolated between nodal values (cf. Section 2.3). Nonetheless, the solutions give similar plastic zone size and maximum plastic deformation norm. The simulation parameters are given in Appendix B.1.2.

Figure 4 . 1 .:

 41 Figure 4.1.: Schematic representation of rough contact and corresponding wear mechanisms.A contact cluster forms a wear particle upon sliding if its area is larger than A * , otherwise the asperity in contact deforms upon sliding without breaking (hypothesis B.). This hypothesis brings asperity-level physics in the wear particle formation process through the critical cluster size A * .

Figures C. 1

 1 Figures C.1 and C.2 show a sensitivity analysis that justifies the use of L/λ l = λ s /∆l = 8 in the work presented. All data obtained from simulation is normalized by λ s for lengths and W 0 = E |∇h| 2 L 2 for loads, where E is the effective Young's modulus and |∇h| 2 is the standard deviation of surface slopes.Figure 4.2a shows one realization of this surface. Figures 4.2b-d show the contact map for different load steps in the range W /W 0 ∈ [0.001, 0.2], indicated by dashed lines in fig. 4.2e. Black area is not in contact, yellow contact clusters are smaller than A * , and red clusters are larger. As the load is increased, clusters grow to span more of the available area, occasionally merging to form larger clusters, as shown in fig. 4.2e with the increase in contact area and maximum cluster size. The latter increases dramatically as the clusters merge, forecasting percolation. Figure 4.2f represents the increase in the number of clusters per increase of contact area ( dN / dA c ).Positive values indicate the regime where contact-area growth is dominated by nucleation of clusters, and negative values indicate the regime where cluster merging dominates contact-area growth. We focus our study to the former regime, in which the clusters are far enough to neglect interactions in the debris forming process. As postulated by[START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF], this corresponds to the mild wear regime[START_REF] Wang | Wear and Wear Transition Modeling of Ceramics[END_REF][START_REF] Das | Ultra Mild Wear in Lubricated Tribology of an Aluminium Alloy[END_REF][START_REF] Hokkirigawa | Wear Mode Map of Ceramics[END_REF]. We let W /W 0 ∈ [0.001, 0.06], discretized in thirty load steps for the other simulations presented in this paper.

  Figures C.1 and C.2 show a sensitivity analysis that justifies the use of L/λ l = λ s /∆l = 8 in the work presented. All data obtained from simulation is normalized by λ s for lengths and W 0 = E |∇h| 2 L 2 for loads, where E is the effective Young's modulus and |∇h| 2 is the standard deviation of surface slopes.Figure 4.2a shows one realization of this surface. Figures 4.2b-d show the contact map for different load steps in the range W /W 0 ∈ [0.001, 0.2], indicated by dashed lines in fig. 4.2e. Black area is not in contact, yellow contact clusters are smaller than A * , and red clusters are larger. As the load is increased, clusters grow to span more of the available area, occasionally merging to form larger clusters, as shown in fig. 4.2e with the increase in contact area and maximum cluster size. The latter increases dramatically as the clusters merge, forecasting percolation. Figure 4.2f represents the increase in the number of clusters per increase of contact area ( dN / dA c ).Positive values indicate the regime where contact-area growth is dominated by nucleation of clusters, and negative values indicate the regime where cluster merging dominates contact-area growth. We focus our study to the former regime, in which the clusters are far enough to neglect interactions in the debris forming process. As postulated by[START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF], this corresponds to the mild wear regime[START_REF] Wang | Wear and Wear Transition Modeling of Ceramics[END_REF][START_REF] Das | Ultra Mild Wear in Lubricated Tribology of an Aluminium Alloy[END_REF][START_REF] Hokkirigawa | Wear Mode Map of Ceramics[END_REF]. We let W /W 0 ∈ [0.001, 0.06], discretized in thirty load steps for the other simulations presented in this paper.

  Figure 4.2.: Evolution of contact area with increasing normal load. Simulation of forty realizations with λ l /λ s = 8 and fifty load steps in [0.001, 0.20] • W 0 . A * is taken as 2λ 2 s . (a) shows a rough surface sample. (b), (c) and (d) show the state of one realization at the loads indicated by dashed lines in (e). The video provided with the supplementary material shows the contact evolution for this realization. It can be seen from (b), (c) and (d) that the number and the size of clusters increase with W . (e) shows the combined effect of those two contact-area growth mechanisms on the total contact area A c . It also shows the size of the largest cluster, which increases dramatically when the growth of A c is dominated by cluster merging. Figure (f) shows the rate of increase in the number of clusters with respect to A c . Positive values indicate a regime where the growth of A c is dominated by cluster nucleation and negative values indicate that the merging of clusters dominates the contact-area growth.

Figure 4

 4 Figure 4.3.: Distribution of contact cluster areas. The main graph shows p(A/λ 2s ,W ), the probability density function of normalized cluster areas for H = 0.8 and varying λ l /λ s , evaluated using twenty logarithmic bins. The inset shows the probability density function for values of H ∈ {0.6, 0.7, 0.8} with λ l /λ s = 128 . p(A,W ) can be approximated by a power-law within a certain cluster size interval, inside which it is independent of the applied load W . Increasing W or λ l /λ s increases the upper bound of the power-law interval. Normalization with the smallest wavelength λ s collapses all distributions to a single curve within the power-law interval. Varying the Hurst exponent has a limited effect on the resulting distribution.

  Figure 4.4.: Archard's wear coefficient in context of self-affine surface contact. (a) shows the complementary cumulative probability distribution function (H = 0.8 ), (b) shows Archard's wear coefficient as a function of the applied load (for λ l /λ s = 128 and H = 0.8 ).Regardless of A * there exists a critical load at which the wear coefficient transitions from zero (i.e. no wear debris) to a constant value (i.e. steady-state mild wear regime). This critical load largely depends on the value of A * : the transition occurs at a higher critical load for contacts with lower interfacial shear strength and consequently larger A * (i.e. better lubrication condition). In (b) the y-axis is linear up to 10 -5 .

  Figure 4.5.: Wear coefficient and wear rate (H = 0.8 and λ l /λ s = 128). While the wear coefficient isnon-constant with the load, its derivative decreases with the load. However, regardless of A * , the limit value of K is one, which does not correspond to experimental observations. Similarly, the wear rate is increasing non-linearly with the load, although its derivative stabilizes to a fixed value with increasing load.

2 s

 2 Figure 4.6.: Comparison between Archard's wear coefficient K PL and the proposed interpretation K PL for a power-law distribution of cluster areas (α = 1.5 ). As can be seen in the expressions of eqs. (4.9) and (4.10), the behavior of K is different from K when α < 2 : we have K → 1 whereas K plateaus at values in ]0, 1[.

Contents 5 . 1 .

 51 Elastic-plastic contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.2. Crack nucleation in rough surface contact . . . . . . . . . . . . . . . . . . 96 5.3. Single asperity crack nucleation . . . . . . . . . . . . . . . . . . . . . . . . 100 5.4. Multi-asperities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5

  For reference, we can express with M and N both the surface vertical displacement due to an applied pressure M[p] = M[p • e 3 ] ∂ B • e 3 and the residual vertical displacement ūp 3 = N [C : ε p ] ∂ B • e 3 .

  Figure 5.1.: Pressure profiles for 2 plasticity (a), saturation (b) and elasticity (c) models. Although fig. (a) shows that the local pressure can exceed values of 3σ y , the average pressure in contacts p = W /A c is closer to 2σ y , whereas the saturated model gives an average of 2.5σ y with large saturated portions of the micro-contacts. In this case, the normal load is W /W 0 is 6.5 • 10 -2 , and the saturated model predicts a contact area 20% smaller than the

Figure 5

 5 Figure 5.2.: Contact ratio and secant slope of the load/contact area relationship. All models predict a sub-linear load/contact area law because the secant slope decreases with the load.On average the slope computed using the saturation pressure model is 20% lower than the slope computed with the 2 plasticity approach.

  Figure 5.4.: Crack nucleation sites for the saturation pressure (a) and 2 plasticity (b) models.An "upwards" shear is applied on each contact. The true contact area is shown in light grey, the plastic zones in red and the crack nucleation sites in black. We can see that the 2 model has more crack nucleation sites than the saturation model. Since both models give approximately the same true contact area, this discrepancy must be due to plastic residual deformations which are not represented in the saturation approach. Figure (b) shows that there can be multiple crack nucleation sites per contact.

Figure 5

 5 Figure 5.4.: Crack nucleation sites for the saturation pressure (a) and 2 plasticity (b) models.An "upwards" shear is applied on each contact. The true contact area is shown in light grey, the plastic zones in red and the crack nucleation sites in black. We can see that the 2 model has more crack nucleation sites than the saturation model. Since both models give approximately the same true contact area, this discrepancy must be due to plastic residual deformations which are not represented in the saturation approach. Figure (b) shows that there can be multiple crack nucleation sites per contact.

  Figure 5.7.: Maximum tensile stress as a function of applied shear stress across single asperitycontact with 2 plasticity. The loading curves consist of two parts: an initial linear loading and a non-linear saturation of σ I . Although the initial slope is linear, the loading is not elastic as plastic deformations still evolve in the system. Their influence on the maximum tensile stress at the surface is however minimal. A sufficient shear stress may eventually create plastic strain at the surface, which causes a transition in the loading curve. Moreover, given a low enough yield stress, the initial indentation may cause surface yield, as seen for the curve where σ y /E = 2 • 10 -2 .

Figure 5

 5 Figure 5.8.: Failure regimes for a sheared spherical indentation. For the first row of graphs, the ratio σ y /E is varied for a perfectly plastic material. For the second row, the hardening ratio E h /E is varied for a yield ratio of σ y /E = 2 • 10 -2 . On each graph the dashed line shows the yield stress. The white curve marks the transition between failure driven by slip rupture of the junction and failure by crack nucleation. The competition between the junction strength τ j and the critical stress σ c is influenced by σ y because of the saturation effect shown in fig. 5.7. While plasticity gives a failure mechanism independent of τ j , hardening allows the tensile stress to grow past the initial yield limit, giving a linear transition between failure mechanisms.

  Figure5.9.: Inter-contact stress as a function of contact ratio. Inset shows the probability density function of the largest stress eigenvalue at the surface. The stress value corresponding to the peak in the probability density is defined as the inter-contact stress σ i/c . We see that in the rough surface simulation and the Greenwood-Williamson model σ i/c evolves linearly with the contact ratio. The value of κ is estimated from the mean curvature of contacting summits in the rough surface.

Figure 5

 5 Figure 5.10.: Crack nucleation likelihood as a function of re-normalized critical stress, normal load and junction shear strength. The material in contact behaves elastically. We normalize by h rms E * instead of κ because of the difficulty of defining asperity curvature on a rough surface. One can see on figure (a) that the curves corresponding to different normal loads are collapsed on a straight line, showing that the CNL does indeed follow the scaling established in eq. (5.12). When the strength of the junction is taken into account (figure (b)), or, equivalently, if a shear stress is applied, we observe qualitative deviations from eq. (5.12). There is a decrease in crack event density due to the interference of close contacts, which tends to unload the tensile stresses at the trailing edge of leading contacts.

Figure 5

 5 Figure5.11.: Inter-contact stress as a function of contact ratio for the saturation and 2 plasticity models. While the curves do not match the analytical GW approach, their slopes are smaller than in the elastic case. The 2 model shows the smallest slope, indicating that the stresses due to plastic residual deformations have an influence on the inter-contact stress and actually reduce it compare to the underlying elastic stresses of the saturation pressure model.

Figure 5

 5 Figure 5.12.: Crack nucleation likelihood as a function of re-normalized critical stress and normal load for the saturation and 2 plasticity models. The pressure saturation model reproduces a scaling similar to the elastic case in fig. 5.10a, while the 2 model shows a fundamentally different behavior. The crack density is higher in the plastic case because of the additional tensile stresses caused at the edge of contacts by residual plastic deformations.

  ϕ j+1 ( ) d .

2 3E 2 3 2 3

 222 δs = -Eδx 0 -E 11 δxs + µX -1 e, ⇒DCRDδx -Eδx 0 -E 11 δx + DCM δp = s -µX -1 e -R 1 (z k ), ⇔ (DCRD -E 11 ) δx -E R 0 (z k ) -E • δxs 0 + µX -1 e 0 E h + E 00 + DCM δp = s -µX -1 e -R 1 (z k ), ⇔ DCRD -E 11 + E ⊗ E h + E 00 δx + DCM δp = s -µX -1 e -R 1 (z k ) + E E h + E 00 R 0 (z k ) -s 0 + µX -1 e 0 .

  Figure C.1.: Influence of L/λ l on the complementary cumulative distribution function of contact-cluster areas.Grey curves are results for each realization (20 in total), and red curves are ensemble averages. Surface sampled has λ l /λ s = 8 and λ s /∆l = 16 , so that only L is varied. As L/λ l increases, the dispersion of the individual realizations decreases, and the behavior of the ensemble average approximates accurately the behavior of each realization. We also note that as L/λ l increases, the distributions converge to a limit distribution.For the simulations carried out in the rest of this paper, a value of L/λ l = 8 , which is a good compromise between required number of realizations and computational cost, was selected. Note: we analyze the cumulative distribution instead of probability density to remove any bias due to binning.

  Figure C.2.: Influence of λ s /∆l on the complementary cumulative distribution function of contact-cluster areas.Curves are ensemble averages of 20 realizations of a surface with L/λ l = 16 and λ l /λ s = 8 , so only ∆l is varied. As λ s /∆l is increased, the distributions converge to a limit distribution. Increasing λ s /∆l smoothens the distributions and reduces the systematic bias between the computed distribution and the limit of λ s /∆l → ∞. For simulations carried out in the main paper, a value of λ s /∆l = 8 was selected, offering a reasonably low discretization bias, and making the simulations possible with our current code.

2 i + s 2 i

 22 k ) cos(2θ l ) + n i k=1 sin 2 (2θ k ) + 2 n i k <l sin(2θ k ) sin(2θ l ) = n i + 2 n i k <l cos(2(θ k -θ l )) c i c j + s i s j = θ k -θ l ))Computing the expected value of the above expressions gives integrals of the form: θ -γ )) dθ dγ = 0, and we simply obtain E[c ] = n i = 3 • i and E[c i c j + s i s j ] = 0 . Therefore:

P

  crack = P(X z * ≥ h) = P(X and z * ≥ h) P(z * ≥ h) = P(X ) P(z * ≥ h) = exp -σ c /ω(τ j ) -σ i/c κ

  Equilibrium of a half-space . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.2.1. Integral operators and Fourier formulation . . . . . . . . . . . . . 16 1.2.2. The Boussinesq-Cerruti fundamental solution . . . . . . . . . . 1.2.3. Other integral equation methods . . . . . . . . . . . . . . . . . . 1.3. Normal Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.1. Variational form . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.3.2. Solution methods . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.4. Rough surfaces and statistics . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.1. Rough surface generation . . . . . . . . . . . . . . . . . . . . . . 29 1.4.2. Surface statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3. Periodicity and statistical significance . . . . . . . . . . . . . . . . Illustration of the contact area evolution. Contact area between elastic solids with self-affine rough surfaces. Darker shades correspond to larger applied normal loads. The complex contact topography and the associated stress state governs virtually all tribological phenomena. . . . . . . . . . .

	List of Figures
	0.1.
	2. Volume integrals
	2.1. Problem statement and overview of solution methodology . . . . . . . . .	37
	2.1.1.	Horizontally-periodic setting . . . . . . . . . . . . . . . . . . . .	38
			vii

1. Variational form . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2. Boundary integral equation . . . . . . . . . . . . . . . . . . . . . 14 1.2. 3 0.2. Illustration (3D) of the sub-surface plastic zones. Lighter shades are zones deeper below the contact interface.

representation of a periodic elastic-plastic contact prob- lem, with

  .47) and does not depend on discretization if the shortest wavelength is larger than twice the sampling size (i.e. if the PSD can be represented in the Fourier domain). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . geometrical contact quantities represented by blue arrows and plastic deformation zones with dashed red contours. The elastic-plastic body B is represented in deformed shape. Note that this schematic is taken from a real two-dimensional simulation. . . . . . . . . . . . . . . . . . . . . . . .

2.1. Schematic 2.2. Normalized surface displacements due to a hydrostatic spherical inclusion in a half-space. Comparison between the Mindlin and Cheng

  Equilibrium of a half-space . . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.2.1. Integral operators and Fourier formulation . . . . . . . . . . . . . 16 1.2.2. The Boussinesq-Cerruti fundamental solution . . . . . . . . . . 18 1.2.3. Other integral equation methods . . . . . . . . . . . . . . . . . . . 21 1.3. Normal Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 1.3.1. Variational form . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.3.2. Solution methods . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.4. Rough surfaces and statistics . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.4.1. Rough surface generation . . . . . . . . . . . . . . . . . . . . . . 29 1.4.2. Surface statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 1. Elastic rough contact: boundary integral methods 1.4.3. Periodicity and statistical significance . . . . . . . . . . . . . . . . 33

Contents 1.1. Equilibrium of elastic bodies: problem introduction . . . . . . . . . . . . 13 1.1.1. Variational form . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.1.2. Boundary integral equation . . . . . . . . . . . . . . . . . . . . . . 14 1.2.

  .9) This identity allows the expression of an unknown state as a function of a known one. Let us now define the (known) state number 1 : b 1 represents a unit point force in the direction e k applied on point x ∈ E ⊇ B, while u 2 = U k (x, ) is the displacement vector at point caused by the point force. Note that E, the domain of state 1 , should be chosen so that U k can be obtained analytically. The simplest choice is E = R 3 , but other considerations may enter the choice of E, as we shall see next. Due to the nature of the volume force density b 1 , it is necessary to enlarge the mathematical framework defined until now, as b 1 is not a function in the traditional sense. Indeed, it is best expressed with the Dirac distribution δ x supported on {x }. The theory of distributions is a more appropriate setting: it generalizes classical functions while retaining properties of differentiation, integration, convolution and integral transformation, which allow for a rigorous development of the methods we present in the first part of this manuscript 1 .

	Definition 1. We call a fundamental solution a displacement field U k defined on E that
	satisfies

Table 2 .

 2 1.: Symbols for the full-space fundamental solution

	Symbol Expression

  which use DFT of real-space fundamental solutions and[START_REF] Yastrebov | From Infinitesimal to Full Contact between Rough Surfaces: Evolution of the Contact Area[END_REF] Rey, Anciaux, and Molinari, 2017;[START_REF] Weber | Molecular Probes Reveal Deviations from Amontons' Law in Multi-Asperity Frictional Contacts[END_REF] which use Fourier-space fundamental solutions) and dates back to[START_REF] Stanley | An FFT-Based Method for Rough Surface Contact[END_REF] (which implicitly use a Fourier space fundamental solution via Johnson,

Table 3

 3

		.1.: Simplified notation for linear operators
	Notation	Definition	Description
	D	dev(•)	Deviatoric operator

  Contents 4.1. Review of Archard's wear model . . . . . . . . . . . . . . . . . . . . . . . 73 4.2. Wear coefficients based on critical length-scale . . . . . . . . . . . . . . . . 75 4.2.1. Wear coefficient based on Archard's interpretation . . . . . . . . . 76 4.2.2. Alternative formulation of wear coefficient . . . . . . . . . . . . . 76 4.3. Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.3.1. Statistical distribution of contact clusters . . . . . . . . . . . . . . 78 4.3.2. Wear coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.3.3. Analytical results . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.4. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.5. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 T he scientific study of wear dates back to the early 19 th century

Table 4

 4

		.1.: Symbols and notations
	Symbol	Description	Physical dimension
	A	single asperity contact area	area
	R	wear rate	volume/distance
	W , δ	normal applied load, indentation depth	force, distance
	d * , A *	critical length-scale, critical area	distance, area
	(•) 1	quantity relative to a single asperity con-	-
		tact
	K	wear coefficient based on Archard's inter-	dimensionless
		pretation
	K	wear coefficient based on up-scaling ap-	dimensionless
		proach
	(•) PL	quantity computed with power-law con-	-
		tact model
	(•) num	quantity numerically computed	-
	A c , A * c	real contact area, cumulated area of con-	area
		tacts larger than A *
	p(X, )	probability density function of random	dimension of 1/X
		variable X with parameter
	P(X > x, ) probability of the event X > x with pa-	dimensionless
		rameter
	N	number of contacts	dimensionless
	λ l , λ s	largest and smallest wavelengths	distance
	H	Hurst exponent	dimensionless
	L, ∆l	system size, discretization	distance
	d	height density	1/distance
	A s , A m	smallest and largest micro-contact areas	area
	α, C	power-law exponent and plateau value	dimensionless,
			1/area
	E, ν, H	Young's modulus, Poisson's ratio, indenta-	pressure, dimension-
		tion hardness	less, pressure
	E	effective Young's modulus (E ≡ E/(1 -	pressure
		ν 2 ))
	G, ∆w, σ j	shear modulus, adhesive energy per unit	pressure,	en-
		area, junction shear strength	ergy/area, pressure
	Equation (4.

  (within a certain load interval). The wear equation is however non-linear in any other case than rigid-4.2. Wear coefficients based on critical length-scaleTable 4.2.: Parameters for Archard's multi-asperity contact model. E is the Young's modulus, H is the hardness of the material and r is the radius of the spherical asperities.

	Behavior	b	c		p
	Elastic	πr	4.25E	√ r 3 2
	Plastic	2πr	2πr H	1

To expand on biblical imagery, "the devil is in the detail."

The reader is referred to the books ofDavid and Gosselet (2015, chap. 4) andDautray and Lions (2000, Appendix "Distributions") for an introduction to the theory of distribution and summary of key properties.

We only focus here on Neumann problems.

Meaning that f and have supports such that ˜x remains bounded, which is in particular the case when one of the operands is compactly supported.

i.e. f ∈ S and ∈ E , using standard notation for spaces defined by the theory of distributions, see e.g. the appendix "Distributions" in[START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology: Volume 2 Functional and Variational Methods[END_REF]

in plane-strain, but Amba-Rao's derivation method is valid for arbitrary distributions of surface tractions (including the point force).

Cables are an intuitive examples of unilateral conditions: unlike bars, they behave differently if their ends are brought together or pulled apart.

As adhesive potentials are formulated in terms of distance between surfaces, a primal approach is more convenient. However, the problem looses convexity: non-unique and unstable solutions can arise.

The lower bound of the scaling behavior is unfortunately often arbitrarily determined by the resolution of the apparatus used to measure the roughness profile.

I.e. if both distributions are (a) in L 1 (R

), (b) locally integrable with convolutive supports, (c) one is a tempered distribution and the other is a compactly supported distribution[START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology: Volume 2 Functional and Variational Methods[END_REF].

In contrast, the application of these methods to volume integral approaches that compute the DFT of ∇G (e.g.[START_REF] Sainsot | A Numerical Model for Elastoplastic Rough Contact[END_REF] are intricate if not completely impossible.

Not present in[START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF].

The memory complexity remains unchanged since the full Fourier representation of w needs to be stored anyway.

Even approximation of the Jacobian via finite differences.

This bound is established for 2D finite elements, but we expect it to remain representative in the present 3D case.

Note that actually applying the FEM to the present mechanical problem would in fact require discretization of a much larger domain to model the fields away from the potentially plastic zone.

Their stiffness matrix fits in the 512K RAM of the IBM/360 used for simulations, which is impressive for 1971. The smoothness of the results (extracted from figure5) is likely due to the figure being drawn by hand.

Our own experiments with finite-element simulations exhibited the same behavior, as well as disagreement in the plateau value, for too-large element sizes.

In the sense of the Legendre-Fenchel transform. They are alternatively called "dual" variables.

Note that Tr(a) = 0 ⇒ a : b = a : dev (b).

Since the free energy in eq. (3.16) is quadratic, the second order term of the Taylor expansion is expressed directly with f and the expansion is exact.

Note that s and (∆γ , ∆ε p ) both belong to the same Lorentz cone K = {(x 0 , x) ∈ R

: x 0 ≥ x }.

Neighborhood definition has no influence on distribution of cluster areas when discretization is fine enough.

We have A * ∝ σ y -4 if one assumes τ j ∝ σ y and the total contact area A c ∝ σ -1 y according to saturation models

More often than not, the saturation pressure is referred to as "hardness". As[START_REF] Burwell | On the Empirical Law of Adhesive Wear[END_REF] discuss, the saturation (or flow) pressure cannot be absolutely known but is of the same order of magnitude as the value given by usual hardness tests. We therefore keep separate notations for clarity.

This corresponds to linear isotropic hardening.

In this approach the residual displacement is directly computed from ε p , whereas in the saturation plasticity model it is merely assumed equal to the negative gap.

prediction. As a result the connectivity of micro-contacts is different between the two models.

Regularization is done by convolution with a function of the form ϕ ε (r ) = exp(-1/(1 -(r /ε)))/ε 2[START_REF] David | Équations aux dérivées partielles: cours et exercices corrigés[END_REF].

, since z * follows the canonical exponential distribution.

Acknowledgments

Acknowledgments All authors acknowledge the insightful discussions with T. Brink and E. Milanese, as well as the financial support of the Swiss National Science Foundation (grant #162569, "Contact mechanics of rough surfaces").

Supplementary data All codes used in this work are available on Zenodo (Frérot, Anciaux, Rey, et al., 2019; Frérot, 2019b).

Volume integrals

Abstract

A Fourier-based volume integral method for the solution of elastic-plastic problems is presented. The major novelty of this work consists in the derivation of required fundamental solutions directly in the partial Fourier domain introduced in Chapter 1, which allows for drastic memory savings and optimization of the remaining integral. The approach is validated against analytical solutions from the literature and the computational costs are compared with a standard finite element method.

Disclaimer

Parts of this chapter are reproduced from the article [START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF]. "A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact". In: Computer Methods in Applied Mechanics andEngineering 351, pp. 951-976. doi: 10.1016/ j.cma.2019.04.006, with permission of all authors. My personal contributions to this article include the development of the method, the code implementation, the running of simulations, the figure production and the writing of the article.

Contents

Abstract Two solution strategies for the elastic-plastic contact problem are presented. The first is a fixed point approach dealing with the plasticity problem (cf. previous chapter) separately from the contact problem, which is treated as elastic (cf. Chapter 1). The implementation of this algorithm is validated against results from the literature and reference finite-element simulations. The second solution strategy deals with the contact and plasticity problems in a monolithic fashion. The elastic-plastic contact problem is formulated as a second-order conic program which can be solved with the interior point method. Disclaimer Parts of this chapter are reproduced from the article [START_REF] Frérot | A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact[END_REF]. "A Fourier-Accelerated Volume Integral Method for Elastoplastic Contact". In: Computer Methods in Applied Mechanics andEngineering 351, pp. 951-976. doi: 10.1016/ j.cma.2019.04.006, with permission of all authors. My personal contributions to this article include the development of the method, the code implementation, the running of simulations, the figure production and the writing of the article. I n order to solve the full elastic-plastic contact problem, one needs to solve the unknown boundary tractions and the plastic deformations. The displacement field satisfying equilibrium once these quantities have been resolved is obtained with eq. (2.7). The difficulty of the elastic-plastic contact problem resides in the coupling between the contact and the plastic problems: a change in plastic deformation will displace the surface, changing the gap function and therefore the contact solution, while a change in surface traction will influence the plastic deformations.

We present in this work two different coupling methods: the first is an alternative coupling strategy [START_REF] Jacq | Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code[END_REF] where the contact problem is solved with a fixed distribution of plastic deformations, then the resulting contact tractions are used to update the plastic coupling with another boundary integral method [START_REF] Barras | When Dynamic Cracks Meet Disorder: A Journey along the Fracture Process Zone[END_REF].

Appendices

A. Appendix to Volume Integrals

A.1. Invertibility and linear independence of A + and A - Due to its structure, A ± has two unit eigenvalues. The determinant of A ± is therefore:

so does not vanish ∀(q, 3 ) ∈ R 2 × R. Therefore A ± is invertible and rank(A ± ) = 3 . The linear independence of A + and A -then stems from their being proportional to different exponential functions.

A.2. Proof of Theorem 3

The displacement u = N [w] (N being the Mindlin integral operator) solves the problem

In the non-periodic case, the partial Fourier version of the above problem is the ODE system

where q ∈ R 2 acts as a parameter. Solving (A.4) for arbitrary sources w(q, 3 ) yields u(q, •) = N [w](q, •) for any q ∈ R 2 , N [w] being given by (1.20) with H = ∇G. Consider now the case where w is B p -periodic, so has the Fourier series form (2.26a). Since

Using this in the right-hand side of (A.4a), we deduce that u must also assume the form of a series of weighted Dirac distributions with the same supports. Hence u is also B p -periodic and can be expressed as a Fourier series (2.26b) with its coefficients u (k, x 3 ) to be determined:

Applying the partial Fourier transform to problem (A.3) and inserting w(q, 3 ), u(q, 3 ) as given above thus yields

Using (for example) the distributional equality N (q)δ (q -a) = N (a)δ (q -a), we deduce the relations for

which are formally identical to (A.4) with the replacements q → 2π k, u(q, 3 ) → u (k, x 3 ) and w(q, 3 ) → w (k, x 3 ). This shows that u (k,

where N [w] is still given by (1.20) with H = ∇G. This completes the proof of the theorem.

A.3. Simulation data

In this appendix, we give the detailed geometry, discretization and loading of each simulation.

We note A.3.1. Comparison with [START_REF] Mindlin | Thermoelastic Stress in the Semi-Infinite Solid[END_REF] With a the radius of the inclusion, the free parameters are fixed to:

• Discretization N = (128, 128, 126).

• ν = 0.3

Values of E, α or T do not influence the normalized results.

A.3.2. Scaling simulations

With a cubic domain of side L = 1 , the number of points in the discretized system direction (with N 1 = N 2 = N 3 ) in both FEM and VIM simulations vary in {4 3 , 8 3 , 16 3 , 32 3 , 64 3 , 128 3 }.

The source code for the finite elements simulation can be found here: 10.5281/zenodo.2613614 (Frérot, 2019a).

A.4. Complexity of integration with cutoff

The cutoff condition on the integration of an element of center x i c for a point of interest x is that q|x

2 . So for a given value of x 3 , the number of operations needed for the computation of integral (2.31) for all discrete values of q is of the order of:

(A.10)

Since the cutoff is smaller at high wavenumbers the number of terms to be accounted for decreases. We can approximate this series with an integral, which gives the value

A.5. Proof of Theorem 4

Let the function w h : R 2 × R + → R 3×3 be of the form:

w j (q)ϕ j (x 3 ).

(A.11)

For conciseness, we will drop the q dependency on w j as well as

for k = 0, 1 , which are the "bits" constituting the kernel of N ∞ , cf. equation (2.20). The correct q-dependency for the latter can be deduced from Lemma 1. The application of N ∞ to w h yields:

Here we have made use of the properties of ± 0 and ± 1 to separate the x l and 3 variables. There still remain a dependency to x l in the sums, which can be removed if we assume a compact [START_REF] Hardy | The Elasto-Plastic Indentation of a Half-Space by a Rigid Sphere[END_REF] With R the radius of the indenter, the free parameters are fixed to:

• Discretization N = (81, 81, 32)

• ν = 0.3

• LoadingW /W y ∈ [0.9, 15.5] 20 ∪{1, 2.1, 6.4, 15.5} (sorted for monotonic loading history)

Values of E or σ y do not influence the normalized results. The tolerance of Algorithm 5 is set to 10 -9 .

B.1.2. Comparison with Akantu [START_REF] Richart | Implementation of a Parallel Finite-Element Library: Test Case on a Non-Local Continuum Damage Model[END_REF] With R the radius of the indenter, the free parameters are fixed to:

• Discretization N = (81, 81, 32)

• FEM discretization: 81×81 surface nodes, 54 nodes in the x 3 direction, 85975 total nodes and 460246 linear tetrahedron elements. The mesh is refined at the contact interface.

• ν = 0.3

The tolerance of Algorithm 5 is set to 10 -9 . The source code of Akantu is available on https://akantu.ch. The finite-element mesh is generated with GMSH [START_REF] Geuzaine | Gmsh: A 3-D Finite Element Mesh Generator with Built-in Pre-and Post-Processing Facilities[END_REF] and results are visualized with Paraview [START_REF] Ayachit | The ParaView Guide: Updated for ParaView Version 4[END_REF].

B.2. Conic programming

We construct the sequence

T with a Newton-Raphson iteration such that R(z k ) → 0 where R(z k ) is the residual of the system (3.24) with the perturbed complementarity (3.25).

The tangent operator R is comprised of the linear operator of eq. (3.24a) augmented with the linearized complementarity conditions:

where

), S := mat(s k ) and X := mat(x k ). In order to simplify the resolution of the Newton-Raphson iteration R (z k )δz = -R(z k ), we operate two changes to the linear system: (1) introduce a rescaling, proposed by [START_REF] Nesterov | Primal-Dual Interior-Point Methods for Self-Scaled Cones[END_REF], which provides a symmetric version of the complementarity equation (3.25b) [START_REF] Bleyer | Advances in the Simulation of Viscoplastic Fluid Flows Using Interior-Point Methods[END_REF]; (2) reduce the linear system size by inverting R by block. [START_REF] Nesterov | Primal-Dual Interior-Point Methods for Self-Scaled Cones[END_REF] have shown that for all admissible (x, s) there exists a matrix F such that Fx = F -1 s =: v and eq. (3.25b) is equivalent to v • vµe = 0 . As consequence, we have X = V F -1 and S = V F with V = mat(v). The equations associated with the linearized complementarity can be manipulated to explicitly obtain the unknowns δ and δs (cf. Appendix B.2):

The above unknowns can be replaced in the Newton-Raphson linear system. We can additionally isolate δx 0 by replacing δs 0 in the first optimality equation: Recent contact simulations [START_REF] Hyun | Finite-Element Analysis of Contact between Elastic Self-Affine Surfaces[END_REF][START_REF] Pei | Finite Element Modeling of Elasto-Plastic Contact between Rough Surfaces[END_REF][START_REF] Hyun | Elastic Contact between Rough Surfaces: Effect of Roughness at Large and Small Wavelengths[END_REF][START_REF] Campañá | Using Green's Function Molecular Dynamics to Rationalize the Success of Asperity Models When Describing the Contact between Self-Affine Surfaces[END_REF] and experiments [START_REF] Majumdar | Role of Fractal Geometry in Roughness Characterization and Contact Mechanics of Surfaces[END_REF]Dieterich and Kilgore, 1996) have shown that, in contact with self-affine surfaces, the probability density of cluster areas is constant up to an area A s , then follows a power-law of exponent -α in the interval

where A m is the size of the largest cluster in the system and is the only parameter depending on the load (cf. fig. 4.3 for additional numerical evidence). Consider the following probability density function for cluster areas:

where C is chosen to satisfy

C.1.2. Wear coefficients based on Archard's interpretation

Lucas Frérot

Curriculum Vitae

Supervision of student projects