A Finite Element-Based Elastic-Plastic Model for the Contact of Rough Surfaces

Three-dimensional elastic-plastic contact of two nominally flat rough surfaces is considered. Equations governing the shouldershoulder contact of asperities are derived based on the asperity constitutive relations from a finite element model of the elasticplastic interaction proposed by Kogut and Etsion (2002), in which asperity scale constitutive relations are derived using piecewise approximate functions. An analytical fusion technique is developed to combine the piecewise asperity level constitutive relations. Shoulder-shoulder asperity contact yields a slanted contact force consisting of two components, one in the normal direction and a half-plane tangential component. Statistical summation of the asperity level contact force components and asperity level contact area results in the total contact force and total contact area formulae between two rough surfaces. Approximate equations are developed in closed form for contact force components and contact area.


Introduction
The GW theory [1] of contact between nominally flat rough surfaces has been preferred by numerous researchers as it benefits from relatively simple representation of a rough surface.It is based on a statistical account of a rough surface in which three parameters are identified.These include (1) standard deviation of asperity height distribution, σ; (2) average asperity summit radius of curvature, β; (3) area asperity density, η.The GW model treats both elastic and plastic contacts and it presumes that asperity contacts occur independent of each other, that is, no influence from adjacent local contacts on a given asperity contact.In the treatment of elastic interaction, GW model relies on the presumption of the Hertz contact.The GW model has been followed by numerous other studies, as summarized in the review paper by Adams and Nosonovsky [2], which take into account various aspects of surface topography such as contact between two rough surfaces, nonuniform radii of the asperities, non-Gaussian distributions of the asperity summit heights, anisotropy, and plasticity.The work proposed by Greenwood and Tripp [3] extended the GW model to contact between two rough surfaces.Greenwood and Tripp (GT) demonstrated that the contact between two rough surfaces could be treated as that between a flat and a rough surface if the composite statistics of the two surfaces are employed.Namely, Gaussian distribution of heights is in terms of the height sum distribution of the surfaces and the standard deviation of asperity height sum distribution is employed in the formulation of contact.This simplification required a modified function related to the interference of asperities involving the integration of interference function over the range of asperity tangential offset.McCool [4] extended GW microcontact model to include skewness in the distribution of surface summit heights and the presence of a surface coating of prescribed thickness and compliance.Recently, Sepehri and Farhang [5] developed an elastic model for two nominally flat rough surfaces in which asperity shoulder-shoulder contact was permitted to derive formulae for elastic contact of two rough surfaces.
A major contribution to the modeling of nominally flat rough surfaces is the work in 1987 by Chang et al. [6], who proposed a method for treating elastic-plastic contact of rough surfaces.This model, widely known as the CEB model, is based on volume conservation of an asperity during its plastic flow.The CEB model enjoys the simplicity of the GW model while providing a predictive tool for contact problems not amenable to an elastic contact assumption.Many publications have appeared since the CEB model that are based on the CEB or are inspired by the method employed by the CEB model .Many researchers have employed statistical models for the elastic-plastic contact of rough surfaces [7][8][9][10][11][12][13][14][15][16][17].Others have advocated the use of deterministic methods based on fractal characterization of roughness [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35].
Another approach is to use the finite element method (FEM) to study the elastic-plastic contact of a single asperity contact.Kogut and Etsion (KE) [36] performed such an FEM analysis of an elastic-perfectly plastic spherical asperity in contact with a rigid flat.The KE model was then used to give empirical expressions for the contact area, the contact force and the average contact pressure as functions of the interference.Jackson and Green [37] also studied an elasticperfectly plastic hemisphere in frictionless contact with a rigid flat using the FEM and with material yielding based on the Von Mises criterion.This model went farther into the elastic-plastic regime and also examined a wider range of conditions.The finer meshes provided more accurate results over the entire range of deformation.Etsion et al. [38] and then Jackson et al. [39] analyzed different aspects of single unloading of an elastic-plastically loaded sphere in contact with a rigid flat for a wide range of sphere material properties and radii.Jackson et al. [40] used a semianalytical model and finite element model to generate empirical equations describing the tangential and normal contact forces between sliding elastic-plastic spheres.
The FEM based models can be used as building blocks to study multi-asperity contacts with mixed elastic-plastic deformation.Kogut and Etsion [41] and Jackson and Green [42] used the FE models in [36,37] in conjunction with the GW methodology [1] to present an elastic-plastic model for the contact of rough surfaces.Similarly, Kucharski et al. [43] investigated elastic-plastic contact between a hemisphere and a rigid plane using the FEM and combined the resulting relations with a statistical description of rough surfaces.
In this paper, we consider elastic-plastic contact of nominally flat rough surfaces.Equations governing the shouldershoulder contact of asperities are derived based on the asperity-asperity constitutive relations from a finite element model of the elastic-plastic interaction proposed by Kogut and Etsion [36].Shoulder-shoulder asperity contact yields a slanted contact force consisting of both tangential (parallel to mean plane) and normal components.An analytical fusion technique is developed to combine the piecewise asperity level constitutive relations for contact force and real contact area.Statistical summation of tangential contact force component along an arbitrary tangential direction yields the half-plane tangential contact force.Similarly, statistical summation of contact force along the normal direction obtains the elastic-plastic normal contact force formulae for two rough surfaces.Approximate equations are developed in closed form for contact force components and contact area as a function of mean plane separation, sum of curvature radii of asperity summits, and plasticity index.

Elastic-Plastic Contact
Consider the elastic-plastic contact of two nominally flat rough surfaces.As shown in Figure 1, let w c1 and w c2 be the critical interferences of the Surface 1 (S 1 ) and the Surface 2 (S 2 ), respectively.Note that the critical interference on a surface defines the plastic asperities on that surface as illustrated by the dashed curves in Figure 1.
Let w c2 < w c1 , then elastic-plastic behavior would be primarily by the asperities on S 2 .Hence, we consider w c2 as the critical interference for the inception of plastic deformation for the contact of the two rough surfaces.For simplicity we denote the lower critical interference by w c .It should be noted that prior to interference of S 1 with plastic asperities of S 2 there is only elastic contribution.Any elasticplastic contribution would be due to the interference of the asperities on S 1 and the plastic asperities of S 2 .
Since in general asperities meet in a shoulder-to-shoulder contact, a contact force between two asperities would be slanted, giving rise to both normal and tangential force.This is illustrated in Figure 2 wherein the interference between shoulders of two asperities and the resulting contact force are depicted.It can be shown by considering the geometry of interference between surface asperities (Figure 2) that the interference is [5]  where β s is the sum of curvature radii of asperity summits and r the tangential offset of the mating asperities so that when r = 0 the asperities interfere along the normal to the mean planes.In ( 1) and (2) the parameters have been normalized with respect to the standard deviation of asperity height sum σ, so that s is (z 1 + z 2 )/σ, h is d/σ and r and β s are the normalized values using σ as the normalization parameter.
Kogut and Etsion [36], using an FEA model, obtained the following piecewise fits for contact load and area of contact between a deformable sphere and a rigid flat: Elastic Range: w cr = 0-1 Elastic-Plastic Range (1): w cr = 1-6 P cr (w cr ) = 1.03(w cr ) 1.425 , Elastic-Plastic Range (2): w cr = 6-110 P cr (w cr ) = 1.40(w cr ) 1.263 , A cr (w cr ) = 0.94(w cr ) 1.146 , (5) where w cr is the ratio of interference to the critical interference and the critical interference is that corresponding to the onset of plastic flow proposed by Greenwood and Williamson [1] where β = β 1 β 2 /β s is the equivalent radius of curvature of asperity summit and H is hardness of softer material, that is, Surface 2. The hardness coefficient, K, is related to the Poisson ratio by K = 0.454 + 0.41ν and the hardness is assumed H = 2.8S y .Alternatively, from the Jackson and Green model [37] where C is related to the Poisson ratio by C = 1.295 exp(0.736ν).E in (7a) and (7b) is the combined Young's modulus for the two surfaces where E 1 , E 2 and ν 1 , ν 2 are Young's Moduli and Poisson ratios of two contacting materials, respectively.P cr in (3)-( 5) is the ratio of contact load to the load at critical interference, P/P c .Likewise, A cr is the ratio of contact area to the contact area at critical interference, A/A c ; where P c and A c are, respectively, Here we propose, the following continuous form of the (3)- (5).
The asperity contact force in (11) is directed along the normal to the contact patch.It yields two components as shown in Figure 2. Denoting f N and f t the components of the asperity contact force along the normal and tangential (parallel to the mean plane) direction, respectively, we find, with the help of ( 2) and ( 9), The asperity contact area with the help of (10) can be found as

Normal Force
The normal components of various contact forces are parallel and can be algebraically summed by statistical means to obtain the total normal force of one surface on another.Statistical summation of asperity normal force components yields the total normal contact force between the two rough surfaces as follows: where η 1 and η 2 are the number of asperity per unit nominal area on S 1 and S 2 , respectively, and A n is the nominal area.I N is the statistical integral.For a Gaussian distribution of asperity height sum it is e −s 2 /2 r dr ds. ( It is noteworthy to mention that in ( 15) and ( 16) the parameters are normalized with respect to the standard deviation of asperity height sum, σ, so that s is z/σ, h is d/σ, and r, β, β s , and w c are the normalized values using σ as the normalization parameter.

Tangential Force
The tangential components due to various interactions cannot be algebraically added as they are projections of contact force onto the mean plane and depend on circumferential position of asperities on surface S 2 (Figure 4).In considering the tangential component of contact force, we seek the components of the tangential contact force along an axis of interest, for instance tangential force component along the x-axis, depicted in Figure 4. We are interested in formulating the cumulative effect of x-component of tangential force along the positive x direction (as shown in Figure 4).
Hereafter, as we generate result for the x-component of the tangential force due to positive contact slope, we will refer to this as the "tangential force" and denote by the force component F x .The goal here is to account for the tangential force of an asperity that would be experienced on each side, and therefore accumulation or summation of such forces would establish the tangential load on a surface from each side, that is, due to all contacts at positive slope.Tangential force due to all asperities at height z 2 confined in area dA and at radial distance r can be found as where φ 2 (z 2 ) is the density function associated with asperity heights on the surface S 2 and η 2 , the number of asperities per unit nominal area on S 2 .The component of this force along +x is By considering Figure 4, the force due to all asperities in the +x half plane at height z 2 and distance r would be obtained by integration of ( 18) over θ = −π/2 to π/2, resulting in Using (19), accounting for the contribution of all asperities and considering a Gaussian distribution of asperity height sum, it can be shown that the component of the tangential force between surfaces S 1 and S 2 , along +x, may be found using the following equation: where e −s 2 /2 r 2 dr ds. (21)

Contact Area
All the asperity contact areas can be algebraically summed by statistical means to obtain the total contact area of one surface on another.Statistical summation of asperity contact area yields the total contact area between the two rough surfaces as follows: where, η 1 and η 2 are the number of asperity per unit nominal area on S 1 and S 2 , respectively, and A n is the nominal area.For a Gaussian distribution of asperity height sum so that A cr (w cr )e −s 2 /2 r dr ds. (23)

Approximate Equations
In this section we introduce approximate equations for the integral functions of normal and tangential forces as well as contact area.Based on the dominant physical interaction, we define three ranges for critical interference or corresponding plasticity index (ψ = 1/ √ w c ) to be able to find the most accurate fitting functions.
The approximate function for each integral is denoted using an additional letter "a" in the subscript to signify approximation.For instance, the approximations to dimensionless normal contact force component, I N , is denoted I Na and is given as follows: Elastic Range: w c = 2.75-200 or ψ = 0.07-0.6Elastic-Plastic Range ( 1): w c = 0.16-2.75 or ψ = 0.6-2.5 Elastic-Plastic Range ( 2  The approximate function in (24) yields accuracy to within 7 percent (7%) over the entire domain of h, β s , and w c .The approximate equation for the dimensionless tangential contact force component, I x , is where for the elastic range for the elastic-plastic range ( 1) and for the elastic-plastic range ( 2) Figure 8 illustrates I x over h = 1 to 4 and β s = 100 to 2000 for the values of plasticity index of 0.1.A similar observation applies to the results relevant to the elastic-plastic ranges.Assess the accuracy of the approximation in ( 29) by defining the following error between the dimensionless load component and its approximation: Similar accuracy (7%) is obtained by ( 29) for the half-plane tangential force component.
In the same way, we find the approximate equation for the dimensionless contact area, I A , as follows where for the elastic range for the elastic-plastic range ( 1) and for the elastic-plastic range ( 2) Figure 9 depicts I A over h = 1 to 4 and β s = 100 to 2000 for the values of plasticity index of 0.1.Define the percent error as follows: Using the above, we find that the approximate function in (34) yields accuracy to within 8 percent (8%) over the entire domain of h, β s , and w c .

Comparison with CEB-Based Model
The model based on CEB [17] extended the CEB model to handle the oblique contact of asperities on two rough surfaces in contact.From extension to CEB model [17] we have where e −s 2 /2 r dr ds, where e −s 2 /2 r 2 dr ds, e −s 2 /2 r 2 dr ds, e −s 2 /2 r 2 dr ds, where e −s 2 /2 rdr ds.
(44) (45) As shown in Figures 10,11,and 12, for elastic contact (w c > 3 or ψ < 0.6), both the present and the extension to CEB models yield identical results as would be expected.However, large differences (of up to 45% in the contact load and contact area for a given separation) are found for w c < 1 or ψ > 1.It is interesting to note from Figures 9 to 11 that error between the two models does not depend on asperity summit radius of curvature sum, β s .

Concluding Remarks
The asperity level constitutive equations were presented based on the work by Kogut and Etsion [36].Kogut and Etsion [36] developed a finite element model of an elastic plastic sphere in contact with a rigid flat.Based on the FEA results, they established the relation between contact force and interference and contact area and interference for different ranges of interference ratio.This paper developed continuous constitutive asperity equations relating (1) the asperity contact force to interference and (2) asperity area of contact to the interference.This was accomplished by devising an analytical fusion   half-plane tangential contact force component resulted in the formulation of the tangential force impinged upon one surface by the other due to the cumulative effect of interactions in a half plane.In the absence of an applied tangential force the net tangential force transferred between the two surfaces is zero due to symmetry of interactions about an asperity.Similarly, statistical summation of the asperity contact force along the normal direction and asperity contact area, respectively, yielded the total normal contact force and contact area formula for two rough surfaces.Approximate equations were forwarded for the integral functions of contact force components and contact area.These equations were shown to provide accuracy within seven and eight percent, respectively, for contact force components and contact area over ranges of mean plane separation, asperity summit radius of curvature sum, and plasticity index.The approximate equations greatly simplify solution of problems involving elastic-plastic contact of rough surfaces.
A comparison with the approximate elastic-plastic CEB model showed identical results for elastic contacts having plasticity index values below 0.6 but substantial differences for elastic-plastic contacts with plasticity index values above 1.

Nomenclature
w: Dimensionless interference w 1 : Dimensionless interference defined in GT [3] α: Contact angle between two asperities r: Dimensionless asperity tangential offset A n : Nominal area E : Combined Young's modulus H: Hardness of the softer material (= H 2 ) K: Maximum contact pressure factor σ: Standard deviation of asperity height sum h: Dimensionless mean separation z 1 , z 2 : Heights of asperities on the surfaces 1 and 2 measured from the mean asperity heights s: Dimensionless asperity heights sum β 1 , β 2 : Dimensionless average summit radius of asperities on the surfaces 1 and 2 β: Combined asperity summit radius of curvature β s : Dimensionless asperity summit radius of curvature sum d: Mean separation w c : Smaller dimensionless critical interference ψ: Plasticity index w cr : Ratio of interference to the interference for onset of plastic flow S y : Yield strength η 1 , η 2 : Asperity areal density for the surfaces 1 and 2 P cr : Dimensionless contact load in KE [36] A cr : Dimensionless contact area in KE [36]

Figure 1 :
Figure 1: Elastic-plastic contact of two rough surfaces-for w c2 < w c1 , elastic-plastic behavior would be dominated by the Surface 2.

Figure 2 :
Figure 2: Asperity contact: Overlap region showing normal and oblique interferences-Elastic-plastic force and its components.

Asperities on S 2 AsperityFigure 4 :
Figure 4: Schematic showing the tangential components of contact force exerted by asperities of S 2 on an asperity on S 1 .

Figures 5 , 6 ,
Figures 5, 6, and 7 illustrate I N over h = 1 to 4 and β s = 100 to 2000 for the values of plasticity index of 0.1 in the elastic range, 1 in the elastic-plastic range (1), and 5 in the elasticplastic range (2), respectively.To assess the accuracy of the approximation in (24), we define the following error between the dimensionless normal contact force component and its approximation in percent error form E N h, β s , w c = I N h, β s , w c − I Na h, β s , w c I N h, β s , w c × 100.(28)
E 1,2 : Young's moduli ν 1,2 : Poisson ratios P c : Contact load at critical interference A c : Contact area at critical interference f N : Component of the asperity contact force along the normal direction f t : Component of the asperity contact force along the tangential direction A s : Asperity contact area F N : Total normal contact force I N : Dimensionless total normal contact force F x : Total half-plane tangential contact force I x : Dimensionless total half-plane tangential contact force A: T otalcontactarea I A : Dimensionless total contact area I Na : Approximate function for I N I xa : Approximate function for I x I Aa : Approximate function for I A E () : Percent error between I () and I ( )a .