Elastic-Perfectly Plastic Contact of Rough Surfaces: An Incremental Equivalent Circular Model

In this paper, an incremental eqivalent contact model is developed for elastic-perfectly plastic solids with rough surfaces. The contact of rough surface is modeled by the accumulation of circular contacts with varying radius, which is estimated from the geometrical contact area and the number of contact patches. For three typical rough surfaces with various mechanical properties, the present model gives accurate predictions of the load-area relation, which are verified by direct finite element simulations. An approximately linear load-area relation is observed for elastic-plastic contact up to a large contact fraction of 15%, and the influence of yield stress is addressed.

contact pressure, Pullen and Williamson [12] adopted an energy balance method to determine the load-area relation for plastic contact of rough surfaces. The original GW model has been widely extended by applying proper elastic-plastic contact model to each isolated contacting asperity. For example, Chang et al. [13] furthered this model by considering the volume conservation of plastically deformed asperities. In view of the fractal characteristic of rough surfaces, Majumdar and Bhushan [14] presented a fractal contact model, in which contact patch with area larger than a critical value obeys the elastic Hertz theory, while that with smaller area is assumed deforming fully plastically. Inspired by Archard's model [1], Jackson and Streator [15] proposed an iterative multi-scale approach to calculate the contact area of elastic-plastic solids. Adopting the concept of scale magnification, Persson [16] accounted for the plasticity through introducing the yield stress as the upper bound of local contact pressure. Using an explicit dynamic Lagrangian finite element calculation, Pei et al. [17] modelled the contact between a rigid plane and an elasto-plastic solid with fractal surface. Scale dependent strain gradient plasticity has also been incorporated into the contact of rough surfaces [18]. Even with the presence of appreciable plastic deformation, these works displayed a linear rise of contact area with load. And plastic deformation always takes place throughout the contact process, which was observed in a recent experiment [19].
In addition to multi-asperity models, the profilometric model was introduced by Abbott and Firestone in 1933 [20]. In this model, the contact area was taken as the geometrical intersection of the rough surface by a virtual plane. For fully plastic deformation, the load was given by the contact area times the flow pressure or the indentation hardness [21]. When the contact area is an appreciable fraction of the apparent area, the whole concept of multi-asperity model fails, while the profilometric theory can still be employed to determine the contact area for realistic separation [22].
Recently, Wang et al. [23] adopted the profilometric model to calculate the contact area, and then used an incremental equivalent model to analyze the elastic contact of rough surfaces, which showed good agreement with direct finite element simulation.
In the present work, this approach will be used for the contact of elastic-perfectly plastic solids with rough surfaces.
The paper is organized in the following order. In section 2, the incremental equivalent contact model for rough surfaces is summarized. The indentation of elastic-perfectly plastic solids by a circular punch is simulated in section 3 to obtain the general expression of the contact stiffness. In section 4, the elastic-plastic contacts for three typical rough surfaces are modelled by finite element method. Finally, the comparisons between the predictions of new model and the finite element simulations are discussed in section 5.

The incremental equivalent contact model
Recently, Wang et al. [23] developed an incremental elastic contact model for rough surfaces, in which the contact of rough surfaces was modelled by a group of equivalent circular contacts with varying radius. Here the basic concept and formulations of this model are presented briefly, and are adopted to deal with the elastic-plastic contact of rough surfaces.
Consider an elastic-perfectly plastic substrate with rough surface being compressed by a rigid plane, as schematically shown in Fig. 1(a). The Young's modulus of the substrate is E, the Poisson's ratio ν and the yield stress σ y . Choose the mean plane of the rough surface as the reference plane, and use the separation between the rigid plane and the reference plane indicating the compressive process.
Given the height information of a rough surface, the total section area A c and the number of section patches N generated by a virtual plane at a surface separation z can be determined through geometrical analysis. The size and shape of these sections may be irregular as shown in Fig. 1(b). They are further simplified by N identical circular contacts having the same total contact area A c , as shown in Fig. 1 (1) In this model, the purely geometrical generated section area will be treated as real contact area, and then it is to determine the load required to achieve such contact area reversely by an incremental form.
At a specific surface separation z, denote the current load by P as shown in Fig. 1.
Then for a decrement of separation dz, we try to calculate the required increment of load dP using the contact stiffness defined by dP/dz. For a circular contact of radius R, the contact stiffness in elastic deformation is given by 2E * R with E *  E(1ν 2 ) being the plain-strain elastic modulus of the substrate [24]. While for elastic-perfectly plastic deformation, the contact stiffness is a little complicated, but can be expressed in terms of the current load, the yield stress σ y in addition to E * and R to be given in section 3.
For N circular contacts with identical radius R, the total contact stiffness is given in which g is a dimensionless function given in section 3. Substitution of Eq. (1) into ( The initial load is taken as P  0 for the surface separation z at infinity. To achieve a contact area A c (z), the total load is obtained by solving the differential equation (3).
The present model demands both the contact area A c (z) and the number of contact patches N(z). From the Abbot-Firestone curve [20], A c (z) can be immediately extracted. However, the number of contact patches N(z) is too complicated to be given in an explicit expression [21]. As an option, numerical techniques are adopted to obtain both A c (z) and N(z) from rough topographies in the present study.
In our operations, a rough surface is discretized by nodes with identical lateral spacing, then each node contributes equally to the projected nominal contact area A 0 .
For a given surface separation z, any node with a height higher than z are marked as one contact node. Then the contact fraction A c (z)/A 0 can be calculated as the ratio of the number of contact nodes to the amount of surface nodes. Meanwhile, any two surface nodes with their projections on the reference plane being adjacent are recognized as within the same contact patch, then the amount of the contact patches N(z) can be obtained.
Eq. (3) is a differential equation, and the contact stiffness dP/dz is a function of z and P itself. At a sufficiently large surface separation z, the normal load P vanishes.
Such an initial value problem could be solved numerically by explicit iteration methods as the fourth-order Runge-Kutta method. Firstly, the surface separation z is discretized as z n with identical interval h in a descending order, and z 0 is large enough such that no contact happens and P(z 0 )  0. Then, by using Eq. (3), the contact stiffness dP/dz at points (z n , P n ), (z n h/2, P n k 1 h/2), (z n h/2, P n k 2 h/2) and (z n h, P n k 3 h) are calculated successionally as k 1 , k 2 , k 3 and k 4 , respectively. As surface separation decreasing from z n to z n , the load increment P n1 P n is then calculated as (k 1 k 2 k 3 k 4 )/6 times the interval h. Finally, the load sequence P n and the subsequent load-area relation are obtained. The interval of surface separation h is chosen as 0.002 times the root-mean-square (rms) roughness σ of the rough surface.
Since the numerical method adopted here is a fourth-order method, the total accumulated error of load P is on the order of o[(0.002σ) 4 ], which is negligible.

Contact stiffness of elastic-perfectly plastic substrate
Riccardi and Montanari [25]  process. The accuracy of numerical results has been guaranteed through convergence tests.
According to Riccardi and Montanari's analysis, the overall contact response can be described by the variation of the dimensionless load L/(σ y R 2 ) with respect to the normalized indentation depth 2E * δ/(σ y R). For the indentations on various substrates, it is found that the load-depth relation does follow this scaling law, as shown in Fig. 3.
When the mean pressure p  L/(R 2 ) is below the yield stress, the linear relation holds between the load and indentation depth. As the depth increases further, the curve will deviate from the linear elastic relation and the slope will drop. At very large depth, the mean pressure maintains almost a constant (i.e. 3σ y ).  For the sake of testing the efficiency of our model, the contact process of a rigid plane and an elastic-perfectly plastic substrate with various rough surfaces is simulated through ABAQUS, as shown in Fig. 4. The rms roughness σ is calculated as 0.500 m, 1.00 m and 0.103 m for surface A, B and C, respectively.

Results and discussions
For those three rough surfaces, we use the numerical approach in section 2 to calculate the contact fraction A c /A 0 and the amount of contact patches N, as displayed in Fig. 7. As the surface separation z/σ decreases, the contact area continuously mean pressure over the rough surfaces. Fig. 11 displays the variation of the mean pressure P/A c normalized by the plain-strain elastic modulus E * with respect to the contact fraction. For small contact fraction less than 1%, statistical fluctuations appear, in which the total number of contact nodes is quite small and the mean normal pressure is quite sensitive to the contact area. When the contact fraction is larger than 1%, the mean pressure gets stable throughout the contact process. Since the beginning of contact, the mean pressure is beyond the yield stress. For examples, the mean pressure is about 2.5σ y to 3.5σ y for material with σ y /E = 0.001, and is about 2σ y to 2.5σ y for substrate with σ y /E = 0.007. Therefore, plastic deformation takes place throughout the contact process of rough surfaces, which is clearly distinct from the indentation by a single punch with elastic deformation firstly and plastic deformation subsequently in section 3. Recently, experiments on metals have verified this characteristic [19]. Fig. 11 demonstrates that the magnitude of mean pressure depends on the mechanical properties and topography of rough surface. With the yield stress enhanced, the mean pressure will rise.
For a given contact fraction (i.e. 5%), the dependence of mean pressure on yield stress is plotted in Fig. 12. For a specific rough surface, the mean pressure is bounded by an upper limit corresponding to purely elastic contact. While for fully plastic contact at small σ y /E * , the mean pressure is often given by the hardness of substrate, which is about 3 times the yield stress [11,21]. For a high yield stress, the mean pressure moves to the elastic limit and the topography plays an evident effect. While for a low yield stress, the influence of topography becomes weak. For σ y /E * under

Conclusions
In this paper, we develop a new method to determine the load-area relation for contact of elastic-perfectly plastic solids with rough surfaces. For various material properties and rough surfaces, the predictions of our model agree reasonable well with direct finite element simulations. Even for elastic-plastic contact, the contact area is linearly proportional to the load till for a large contact fraction as 15%. Since the beginning of contact, the plastic deformation plays its role and tends to increase contact area. The mean pressure is beyond the yield stress at the initial contact and keeps almost constant throughout the contact process. The enhancement of yield stress could improve the mean pressure, following a power law relation. The new proposed model demonstrates its accuracy, conciseness and general applicability, and is a powerful tool to deal with contact of rough surfaces.