Heterogeneity of material structure determines the stationary surface topography and friction

The character of surface roughness and the force of friction in the stationary state after a sufficiently long run-in process are of key importance for numerous applications, e.g. for friction between road and tire. In the present paper, we study theoretically and experimentally the asymptotic worn state of a bi-phasic material that is arbitrarily heterogeneous in the contact plane, but homogeneous in the direction of the surface normal. Under the assumption of Archard’s wear law in its local formulation, the asymptotic shape is found in the closed integral form. Given the surface profile, the coefficient of friction can be estimated, since the coefficient of friction is known to be strongly correlated with the mean square root value of the surface slope. The limiting surface profiles and the corresponding coefficient of friction are determined as functions of size, relative concentration and wear ratio of the phases. The results of numerical calculations are compared to and validated by experiments carried out on simplified model systems. The main conclusion is that the rms value of the surface slope is not influenced by the characteristic linear size of inclusions and depends solely on the relative concentration of phases, as well as the ratio of their wear coefficients.

Friction and wear occur everywhere in our daily life. One prominent and economically very important example is the contact of a vehicle's tire with the road. Many important functional properties are determined by the exact nature of this contact, e.g. skid resistance, which is closely related to the driving safety 1,2 . In the last few decades, a lot of experimental and numerical modeling efforts have been undertaken to study the tire-pavement interaction, which led to the development of new measurement devices 3 , improved characterization of surface texture 4 , development of mixed aggregates and asphalt binder 5 , and assessment of environmental factors 6 . A number of tests have shown that the properties of the aggregate, including hardness and wear resistance, as well as the mixture of coarse and fine aggregates have a strong influence on the frictional properties of the pavement surface [7][8][9] . The tire-road contact is only one example and a typical representative of a broad class of tribological contacts where a relatively stiff rough surface (which, however, is rarely replaced and able to be worn away over time) is in sliding contact with an elastomer. The viscoelastic contact partner also undergoes wear, but this is of lesser importance, either because it is replaced more often or because its surface profile has lesser influence on the frictional properties of the contact (both of which is the case with tires). Due to wear, the surface topography of the stiff contact partner will evolve over time, which, even more importantly, may lead to changes in the coefficient of friction. The coefficient of friction between a rough rigid substrate and an elastomer is known to depend on a number of system and loading parameters including sliding velocity, temperature, normal load, etc 10 . However, under conditions relevant for most practical applications (corresponding to the region of the long plateau of the coefficient of friction as a function of velocity) it can be loosely estimated as the root mean square (RMS) of the surface gradient 11,12 . In the present paper, we will therefore consider the root mean square value of the surface gradient (surface slope) as a parameter roughly characterizing friction coefficient.
The system we study in the present paper is schematically shown in Fig. 1. It represents a cylindrical punch with multiphasic composition sliding horizontally with velocity v on an elastic half space in the presence of a normal load F N . Similar problems have been studied in the past using the finite element method 13,14 , the discrete element method 15 and molecular dynamics 16 . However, if we are only interested in the limiting shape, solution can be obtained in closed integral form, which enables extensive parameter studies.

Results
We consider a cylindrical punch composed of a multiphasic material, whose phases i have the wear coefficients k i (defined by eq. (2)) and the cross-sectional areas A i . Let us denote the initial shape of the indenter face as f 0 (x, y). To achieve a qualitative understanding of the wear process, we make the following simplifying assumptions: (a) only the indenter is subject to wear, (b) the elastic properties of all phases are equal, while their wear coefficients may have arbitrary ratios, (c) wear is governed by Archard's law in its local form, which means that the depth change df at the coordinate (x, y) due to wear is proportional to the local pressure p(x, y) and the sliding distance ds, and is inversely proportional to the material hardness σ 0 wear 0 k wear is a constant that depends on the material couple and is called the wear coefficient. Introducing the notation k = k wear /σ 0 , we can rewrite (1) in the form df (x, y) = kp (x, y) ds. For a multiphasic indenter, this equation is valid for each particular phase: df i (x, y) = k i p i (x, y) ds. Dividing by the time increment dt, the linear wear rate (wear depth per unit time) is obtained: The steady state is achieved when the linear wear rate at all points of the indenter is the same. This means that the pressure inside the area of each particular phase reaches a constant value, which we just continue to denote p i = p i (x, y) (despite the fact that the pressure no longer depends on the coordinates inside one phase; it is still a function of coordinates in the sample as a whole), while pressures in different phases fulfil the condition: The constant C can be determined from the equation for the total normal force With this pressure distribution, the normal deformation of the elastic half space can be calculated as 11 is the effective modulus of elasticity of the elastic half-space, E is its Young's modulus, ν is its Poisson ratio, and = − ′ + − ′ r x x y y ( ) ( ) 2 2 is the in-plane distance between the points (x, y) and (x′, y′). As we consider the indenter to be much stiffer than the counter-body, the deformation of the latter is equal to the worn shape of the indenter. Thus, the limiting worn profile is given by This equation provides the final worn profile in explicit integral form, so that both the profile and all parameters depending on the profile (as e.g. the rms value of the slope) can be calculated very effectively.
Limiting worn profile of a homogeneous cylinder. Let us start consider some specific examples. We begin with the simplest case of a monophasic cylinder with radius a and wear coefficient k. In the final state, the pressure over the whole contact will be constant and equal to π = p F a /( ) N 2 . The integral (8) has, in this case, the known analytical solution 11,17 Here r is the polar coordinate and E(k) denotes the complete elliptical integral of the second kind The limiting profile (10) Fig. 2. The rms value of the surface gradient of this profile is equal to A 0 2 Note that neither the limiting shape nor the rms slope depend on the wear coefficient, and are completely determined by the average pressure in the contact and the effective elastic modulus.
Limiting worn profile of a biphasic cylinder. In the following, we consider bi-phasic indenters. We continue to consider a cylindrical sample with radius a, which, however, is now composed of two different materials: phase 1 having the wear coefficient k 1 and area A 1 , and phase 2 with k 2 and A 2 (an example of a typical phase distribution considered in this paper is shown in Fig. 3a; the phases are marked with black and grey colors). Let us introduce the "fill factors" of phases With this notation, the final surface topography according to (9) is where p is the average pressure in contact. Taking into account (12), one can introduce a normalized shape which does not depend on the normal force and elastic properties of the contacting bodies and is completely determined by the geometry of the phase distribution and the ratio of wear coefficients. The generation of phase distributions and numerical evaluation of the integrals (15) is described in the "Methods" section. An example phase distribution and the resulting limiting worn profile are shown in Fig. 3a,b correspondingly. The main target parameter of our study is the rms value of the surface slope, ∇f. Let us start by elucidating the role of scaling in the formation of the limiting value of ∇f. For this sake, we generated four samples with the same statistical properties of phase distribution but different characteristic scale, as shown in Fig. 4(a) for the case of ρ 1 = ρ 2 = 0.5. Similar structures have been produced for fill factors ρ 2 varying linearly from  0.05 to 0.95 and for a number of ratios of wear coefficients, k 2 /k 1 , varying logarithmically from 10 −3 to 10 3 . The value of the rms slope for each set of parameters was averaged over 10 random realizations of the phase distribution. The results of these simulations are shown in Fig. 4b, where the normalized rms slope ∇ ∇ ∞ f f / 0 is plotted. Figure 5a, shows detailed dependencies of the rms slope on the fill factor ρ 2 for five different values of the ratio of wear coefficients and for all four samples with different coarseness of phase distribution (which are marked with different symbols). One can see that the points corresponding to the same values of wear ratio but different coarseness collapse onto a common curve, thus clearly showing that the scaling has no influence on the rms slope. Therefore, in Fig. 5b we only show the results for the sample 2 as a representative of structures with arbitrarily coarse grain.
Thus, the asymptotic rms slope is a function of only two parameters: the fill factor (e.g. ρ 2 ) and the ratio of wear coefficients k 2 /k 1 . It is clear that this dependency remains invariant by transformation ρ 1 →ρ 2 , k 1 /k 2 → k 2 /k 1 , which just changes notation without changing the configuration of the contact: (1 , / ) 2 2 1 2 1 2 . From the overview in Fig. 4(b), one can see that the largest values of the rms surface gradient (and thus coefficient of friction) correspond to the case when both ρ 2 and k 2 /k 1 are large or when both are small. In other words, large final rms slopes correspond to structures with a small fill factor of the phase with the smaller wear coefficient. In this case, the low-wear phase represents a dilute solution of particles in an easily worn matrix. In this extreme case the resulting surface topography resembles a set of "pillars". It is intuitively clear that such a brush-like structure will have a large coefficient of friction. In the other extreme case of large concentration of the low-wear phase, the worn surface will consist almost entirely of this phase with only small number of shallow "dimples", where the softer inclusions have been worn away. Two surface topographies with the same fill factor ρ 3 = 0.3 but different ratios of wear coefficients (k 2 /k 1 < 1 and k 2 /k 1 > 1 correspondingly) are shown in the insert of Fig. 5b.
Let us take a closer look at the extreme "pillars" case, assuming that A A 1 2 2 1 2 1 1 2 In this case, the profile (14) takes the form where we used relation ρ ≈ p p 2 2 which, under the above assumptions, follows from (7). Assuming for simplicity that the hard inclusions have a circular cross-section, we can find the rms slope by immediate use of eq. (12) In this limiting case, the rms slope does not depend any more on the ratio of wear coefficients and is inversely proportional to the fill factor of the hard phase. This asymptote can be easily seen in Fig. 5(a) for k 2 /k 1 = 10 −3 . Note that the rms gradient depends only on the overall fill factor, but not on the radius of the inclusions, which is consistent with the lack of scale-dependence found above in the study of random configurations. We thus come to the conclusion that the best way to maintain a high coefficient of friction in the stationary worn state is to use a matrix with a small fraction of hard (low-wear) inclusions. As an example, consider the problem mentioned in the introduction -wear and friction in the tire-road contact. Assuming that the deformation of tire in the contact area is of the order of magnitude of 8%, and thus ≈ .
⁎ p E / 008, a stationary rms slope of the order of unity will be Experimental results. To experimentally validate the results of our numerical simulations, we investigated the worn shape of homogeneous and heterogeneous cylinders. We did not aim to carry out comprehensive parameter studies as in the foregoing theoretical analysis, and confined ourselves to three model systems, which we consider to be representative generic cases: (i) a homogeneous sample, (ii) high-wear matrix with low-wear inclusions, and (iii) low-wear matrix with high-wear inclusions. We first studied wear of a homogeneous cylinder made of epoxy with a radius r = 4 mm. The cylinder was pressed against a rubber band with the normal load F N = 10 N and pulled back and forth with a constant velocity v = 0.01 m/s. The worn surface, measured with a scanning laser interference microscope (see Section "Methods"), is shown in Fig. 6a, and the profiles of four cross sections, marked in Fig. 6(a) with numbers 1 to 4, are shown in Fig. 6b. For comparison, the theoretical solution given by equation (10) is shown in the same plot (green line). It is seen that the theoretical curve is quite close to the experimentally obtained profiles. The small asymmetry is likely the result of a slightly inclined sample.
In the next step we prepared bi-phasic brass-steel cylinders. In Figs 7 and 8, one can see the final worn profiles of bi-phasic cylinders consisting of five rods with radius r = 2.5 mm embedded in a cylinder with radius r = 10 mm: Fig. 7 shows the results for brass rods in a steel matrix and Fig. 8, correspondingly, steel rods in a brass matrix. In Fig. 7a, the steel rods are less worn and produce visible protrusions. In the case of brass in steel, the rods are worn more intensively and clearly result in the formation of dimples (Fig. 8a). The numerical results with parameters F N = 35 N and k steel /k brass = 1.7 for both cases are shown in the corresponding subplots. The elastic modulus of the elastic counterpart was assumed to be E * = 1.3 MPa in all cases, as determined from experiments with homogeneous cylinders. To allow for a more detailed comparison of the effects of heterogeneity on the worn  profile, in Figs 7b and 8b differential profiles of four cross sections are shown. Differential profiles were determined as the difference between the profile of some cross section and the "averaged profile" obtained by averaging the surface height over ring segments at any given radius r. In this way, the macroscopic curvature of the sample is subtracted and the resulting profile shows more clearly the effect of the phase heterogeneity (for a homogeneous material the differential profile would be approximately constant.) It can be observed that the experimental findings match theoretical predictions reasonably well in both cases of steel rods in brass and brass rods in steel. A possible cause for the observed small discrepancies between experimental findings and theoretical results could be the use of the glued abrasive paper, while in simulation we assumed the counter body to be a homogeneous elastomer.

Discussion
Under the assumption of Archard's wear law, we proposed a method for calculating the surface topography of a multiphasic punch in the stationary worn state. The limiting shape of the body after a long run-in process has been obtained in explicit integral form. The evaluation of the integrals was carried out using an in-house BEM code described in 18 . The obtained solution is applicable to systems with an arbitrary number of different phases and arbitrary geometrical in-plane configuration of the phases. A detailed numerical study was carried out only for the case of bi-phasic indenters.
Using the obtained limiting profiles, the rms slope has been determined as a topographic property most directly related to the coefficient of friction. Both simulations and analytical estimations show that the limiting shape does not depend on the coarseness of the phase structure. Thus, compositions that only differ in coarseness of the grain will have the same limiting rms slope. The complete set of parameters determining the final state and the limiting rms slope was determined to be: average (apparent) pressure in the contact area, fill factors (area fractions) of the phases, and the ratio of the wear coefficients. The largest limiting rms slope is obtained in the case of a small fraction of hard (low-wear) fibers, which leads in the final state to a "pillar-like" structure. In the case of a composite with a small concentration of very hard inclusions, the final rms slope is roughly inversely proportional to the area fill factor of the hard phase. The obtained result provides a rule for the design of composite structures with specified frictional properties after run-in.
For validation of the underlying contact mechanics, we carried out experiments on homogeneous and heterogeneous samples and compared the results with theoretical predictions. The experimental results show very good qualitative and acceptable quantitative agreement with theoretical predictions. They confirm that the use of Archard's law of wear in its local formulation is a sensible assumption, allowing correct understanding of wear of composite structures.
It should be stressed that the above results are not universal, and rely on a number of assumptions: (1) The wearable composite contact partner (cylindrical punch in Fig. 1) was assumed to be heterogeneous in the contact plane but homogeneous in the z-direction (depth), effectively consisting of a bundle of fibers or rods oriented normal to the surface; (2) the normal load is kept constant during the wear process; (3) the sliding velocity is assumed to be constant over the entire contact area. This means that that gross sliding is assumed. However, the actual velocity is not relevant and may even vary in time; (4) the counter-body is considered to be elastic and either non-wearable or regularly replaced; (5) finally we assume the validity of Archard's wear law in a local form. However, the results will be similar for any local wear law that is not necessarily linear in normal load.
The assumption (4) may look strange at first glance, but it was chosen with the motivation of pavement wear due to contact with tires described in the introduction. With regard to elasticity, the road can be considered rigid, while tires have significant elasticity. With regard to wear, the situation is the opposite: We assume that only the road is worn, and tires not. With the described conceptual background, it is not important, whether the tires wear or not, as they are changed out regularly, while the road remains the same over long intervals of time. In the experiments, we sought to give our model system the same conceptual properties: a rigid and wearable sample and soft and replaceable counter body.
In real systems one or many of these assumptions may be violated. Thus, the assumption (1) (cylindrical geometry) is of course very artificial. In real composites such as asphalt pavement the composite components have a general three-dimensional geometry. For such geometry, a stationary state will never be reached, since some inclusions are worn down completely, while others become exposed. However, if the ratio of the wear coefficients is not too large, we expect that the worn state will be very similar to that described in the present paper at any given time, although the exact statistical realization of the topography will change slowly with progressive wear. However, if the ratio of wear coefficients is very large, then the system with three dimensional inclusions will be far from the stationary state described here. For such systems, additional analysis, probably using different numerical approaches, would be needed. With respect to the wear law, the results are more robust. As a matter of fact, the only required assumption used in the paper is that the wear rate is local and depends only on pressure. The pressure-dependence can be nonlinear, but if the wear law is non-local (e.g. we consider the finite size of wear particles), then the results will be violated and additional analysis will again be needed. Since the wear law tends to be strongly system-dependent, one cannot make a more general statement of the range of validity of the presented results. The present study has thus to be considered as a mathematical model under some strict conditions, which have to be checked for particular applications.

Methods
Boundary Element Method. The final surface topography (9) was calculated using the Boundary Element Method (BEM), with the given pressure distribution (7). This Method is based on the evaluation of integrals of the type of (8) using a fast convolution based on the Fast Fourier Transform [18][19][20] . In the present contact problem, the evaluation of integrals (14) was implemented as a module for our in-house BEM code.
Generation of surface patterns. The biphasic surface patterns have been generated using the method of Voronoi diagrams, which has many practical applications in various scientific fields, such as geometry, biology, meteorology, architecture, etc. 21 . In this algorithm, seed points are distributed randomly in a plane, after which cells are formed around each seed by assigning every point in the plane to its closest seed point. In the four samples of Fig. 4a, cells with a total number of 20 2 , 40 2 , 80 2 and 160 2 points have been generated. Afterwards, seeds were randomly labelled as phase 1 or phase according to the defined fill factor. The "average diameter" of the cells in the four samples is approximately 1, 1/2, 1/4 and 1/8 relative to each other. The produced structures can be considered self-similar under magnification.

Experimental investigation.
The experiment was carried out with the experimental set-up shown in Fig. 9. The rubber block with a size of 380 × 40 × 40 mm was glued to a linear stage that moved horizontally with controlled velocity. The material of the band is polyurethane with elastic modulus of about 1.3 MPa. The indenter was fixed on another linear stage that could move vertically. The normal load was measured by a force sensor placed between the actuator and the sample. To increase and precisely control the wear rate, abrasive fine sandpaper was glued to the rubber band.
The indenter was pressed into the rubber band with a fixed load; subsequently, the linear stage moved horizontally with a constant velocity 0.01 m/s. After 50 times of back and forth linear movement, the sandpaper was replaced to provide constant wear conditions. The orientation of the indenter remained unchanged during the whole experiment.
Measurement of surface topography. The topography of the worn samples was measured with the confocal 3D laser scanning microscope VKX150/160 of Keyence. This microscope uses a laser light receiving element with 16 bit PMT, a color CCD image sensor with a recording resolution of 3072 × 2304 pixels and confocal optics with pinhole. A single recorded image has a size of 168 × 126 μm with a height resolution of 5 nm. To measure the entire sample 340 images were taken and merged them to a single image 20 × 20 mm in size with the help of the "Auto Stage for intelligent xy control-program". Photos of the two worn samples that were measured in Figs 7 and 8 are shown in Fig. 9.