Effect of surface roughness on normal contact compression response

Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology http://pij.sagepub.com/ Effect of surface roughness on normal contact compression response A Kumar, F Sadeghi and C. M. Krousgrill Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology 2006 220: 65 DOI: 10.1243/13506501JET71 The online version of this article can be found at: http://pij.sagepub.com/content/220/2/65 Published by: http://www.sagepublications.com On behalf of: Institution of Mechanical Engineers Additional services and information for Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology can be found at: Email Alerts: http://pij.sagepub.com/cgi/alerts Subscriptions: http://pij.sagepub.com/subscriptions Reprints: http://www.sagepub.com/journalsReprints.nav Permissions: http://www.sagepub.com/journalsPermissions.nav Citations: http://pij.sagepub.com/content/220/2/65.refs.html >> Version of Record - Feb 1, 2006 What is This? Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 65 Effect of surface roughness on normal contact compression response A Kumar, F Sadeghi , and C M Krousgrill School of Mechanical Engineering, Purdue University, West Lafayette, Indiana, USA The manuscript was received on 24 January 2005 and was accepted after revision for publication on 3 February 2006. DOI: 10.1243/13506501JET71 Abstract: This work presents a numerical model for normal engagement of two rough surfaces in contact. In this study, the Johnson translator system with a linear filter is used to transform a Gaussian white-noise input to an output surface with prescribed moments and autocorrelation function. The rough surface contact model employs influence coefficients obtained from finite element analysis of the contacting bodies. The contact solution accounts for the effects of macroscopic geometry and boundary conditions, and can be used to simulate engagement at a wide range of loads, including loads at which bulk effects dominate the response. A description of bulk deflection in terms of the displacement of the surface mean plane is also presented. The effects of surface topography on normal engagement stiffness are discussed. Keywords: normal contact, rough surfaces, roughness deformation, bulk response 1 INTRODUCTION Even the smoothest engineering surfaces when viewed at a microscopic scale exhibit significant roughness. Rough surfaces can be considered as two-dimensional random processes [1]. They are described completely, in a statistical sense, by their height distribution and autocorrelation function (ACF). When two surfaces are brought in contact, the contact occurs at discrete junctions located near the surface peaks (asperities). The sum of the areas of the contact spots is the true, or real, area of contact. This is much smaller than the nominal area, which is the macroscopic geometry of the contacting bodies. Modelling rough surface contact requires an understanding of material properties and surface topography. The topographical characteristics that are most important are the heights and curvatures of asperities which come in contact with the opposing surface. The asperity heights influence the contact status of the regions, whereas the curvatures affect the local contact stiffness.  Corresponding author: Department of Mechanical Engineering, Purdue University, 1288 Mechanical Engineering Building, West Lafayette, IN 47907-1288, USA. email: sadeghi@ecn.purdue.edu JET71 # IMechE 2006 Statistical models have been developed [2 – 5] to study the effects of surface topography on contact behaviour. Such models usually employ the halfspace theory on which the classical Hertz contact equations are based. Although the simplicity of these statistical models has many advantages especially in obtaining a general understanding of the effects of roughness variation, they are not used in this work. Instead, this investigation describes a deterministic approach to studying normal contact. In recent years, models for rough surfaces have been developed that consider interactions among surface asperities. Majumdar and Tien [6] show that the power spectra of a variety of surfaces follow a power law within the length scales considered and suggest characterizing surface roughness by fractal geometry. In a different approach, Patir [7] describes a method to generate any threedimensional surface roughness with prescribed statistics. In this method, the statistics are described by the height distribution and surface ACF. The method transforms a Gaussian white-noise sequence to the output surface using a linear filter. The filter depends on the desired ACF of the output surface. In simulating Gaussian surfaces, this method can be used for arbitrary ACF. This is because linear transformations of Gaussian random surfaces result in Gaussian random surfaces. However, Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 66 A Kumar, F Sadeghi, and C M Krousgrill non-Gaussian random surfaces do not have this property and their linear transformation usually results in surfaces with different distributions. Therefore in simulating non-Gaussian surfaces, the relationship between the statistical characteristics of the input and output sequences must be established. Patir illustrates this procedure for the special case of a bilinear ACF, in which both the input random numbers and the output roughness heights follow a Gamma distribution. A similar approach for rough surface generation has been used by Hu and Tonder [8] and Chilamakuri and Bhushan [9], who used the Johnson translator system of statistical distributions to transform the Gaussian white-noise sequence. By providing a relationship between the input and output height distributions for surfaces with arbitrary ACF, the use of the Johnson translator system extends Patir’s approach such that surfaces with arbitrary ACFs can now be generated. This approach is also discussed by Bakolas [10]. Using the above-described approaches to generate rough surfaces with realistic topography, the requirements on the number of mesh elements make the finite element approach for solving the contact problem impractical. However, specialized methods to simulate discrete surface contact have been developed. Such simulations study the effects of bringing two bodies with rough surfaces in contact. In addition to surface topography, these methods require knowledge of the relationship between surface tractions and displacement. In these methods, only the surfaces are discretized. The roughness profile is defined as a sequence of heights defined over the entire nominal area. The contact problem is usually formulated in terms of the distributions of surface traction, and an iterative approach is required to obtain the solution. Once the pressure distribution is known, the displacement field can also be obtained. This approach has been used by Karpenko and Akay [11] to develop a model to study contact with friction between two arbitrary rough surfaces. Bhushan [12] reviews various discrete contact models and discusses their implications in friction and wear studies. Borri-Brunetto et al. [13] describe an algorithm to study the evolution of the contact domain at different scales, using a fractal description of surface roughness. All these studies rely on stress – displacement relations that are derived from the elastic half-space solution. The discrete formulation of the contact problem has been discussed previously by Kalker [14], Bhushan [12], Peng and Bhushan [15], and Stanley and Kato [16]. In previous literatures, the models typically use pressure –displacement relations (influence coefficients) derived from the elastic half-space solution. Therefore, the contacting surface is considered as a linearly elastic half-space that is loaded over its real contact area. In this approach, it is assumed that the contact stresses can be treated separately from the more general stress distribution arising in the body from its geometry and boundary conditions. The use of the elastic half-space solution requires the following restrictions [17]. 1. The dimensions of the local contact areas must be small compared with the relative radii of curvature at the individual contact regions. This is necessary so that the regions just outside the contact region can be approximately modelled by the plane surface of a half-space. For a given height distribution, the accuracy of the half-space solution is higher for roughness profiles with long spatial wavelengths, because such profiles are better approximated by a plane surface. 2. The strains within the contact region must be small enough such that the behaviour can be considered linear elastic. This justification is implicitly contained within the requirement of small contact regions. 3. The dimensions of the potential contact zone must be small compared with the dimensions of the contacting bodies. This is necessary to ensure that the stress field, based on the theory of a semi-infinite half-space in extent, is not affected by proximity to actual boundaries. This restriction justifies the consideration of contact stresses independent of the macroscopic stress distribution that is affected by geometry and boundary conditions. The above restrictions also apply to statistical models, which also employ the half-space theory. The third restriction, which requires the independence of contact stresses from bulk geometry and boundary effects, is prohibitive for the purposes of studying contact response at higher loads when bulk effects become important. Furthermore, in flat surfaces, although the real contact area is generally much smaller than the nominal area, the potential contact region is equivalent to the nominal area of the friction material. This precludes the use of halfspace solutions, which assume bodies with infinite size. In this work therefore, influence coefficients are derived using a finite element model of the contacting bodies. 2 GENERATION OF ROUGH SURFACES In this investigation, surfaces are defined by specifying the probability distribution of surface heights and the ACF. This fully defines the surface statistics, implying that all joint distributions of surface geometric variables (such as the slopes, peak heights, Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 JET71 # IMechE 2006 Effects of surface roughness and curvatures) can be determined. The approach is based on the work of Patir [7], Chilamakuri and Bhushan [9], and Bakolas [10]. For a variety of surfaces, the ACF is found to exhibit exponential decay with distance from the origin. For such surfaces, the ACF is defined by its autocorrelation lengths bx and by along the orthogonal x – y-coordinate directions. The ACF can then be written as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 3 u 2 u x y 5 R(x, y) ¼ s2 exp4 2:3t þ bx by 2 (1) where s is the standard deviation of surface heights. It follows that bx (or by ) is the distance from the origin along the x (or y)-coordinate direction at which the ACF decays to 10 per cent of its value at the origin. In the procedure adopted below, the height distribution of the surface is defined by the first four central moments of the corresponding probability density function (PDF). These moments are given by m1 ¼ 1 X iP½z(x, y) ¼ i Š ¼ i¼ 1 m2 ¼ 1 X i2 P½z(x, y) N X M 1 X zij ; mz (2) NM i¼1 j¼1 mz ¼ i Š 67 may depend on a number of factors, such as the characteristics of the class itself, or the process used to generate members of this class. The choice in this work was made because of convenience and also because of the precedents set by previous literature. In recent works [9, 10] on contact of non-Gaussian surfaces, the Johnson translator system was used to generate surfaces with prescribed first four moments. This system uses one of the three functions to transform a Gaussian random sequence to a sequence with prescribed non-Gaussian moments m1 , m2 , m3 , and m4 . The choice of which function to use depends on the moments m3 and m4 of the desired surface. 1. In this approach, a Gaussian white-noise sequence hi,j is transformed using the Johnson system to a non-Gaussian white-noise sequence xi,j . The translation is performed using one of the following functions h g i,j SU : xi,j ¼ j þ l sinh d h g i,j SL : xi,j ¼ j þ l exp d l SB : xi,j ¼ j þ l 1 þ exp ((hi,j (6) (7) g)=d) (8) i¼ 1 ¼ m3 ¼ N X M 1 X (zij NM i¼1 j¼1 1 X i3 P½z(x, y) mz )2 ; s2z (3) mz ¼ i Š i¼ 1 ¼ m4 ¼ N X M 1 X (zij NM i¼1 j¼1 1 X i4 P½z(x, y) mz )3 ; Skz s3z (4) mz ¼ i Š i¼ 1 ¼ N X M 1 X (zij NM i¼1 j¼1 mz )4 ; Kuz s4z (5) Here P[ ] is the probability operator, mz the mean, sz the standard deviation, Skz the skewness, and Kuz the kurtosis of the sequence z. For a Gaussian PDF, Skz ¼ 0 and Kuz ¼ 3. In general, the first four moments alone do not uniquely define a PDF. Two surfaces could have the same first four moments, but quite different distributions in which the higher moments differ. In the procedure used, uniqueness of PDF is guaranteed by restricting the set of possible PDFs to a class in which the members are completely defined by the first four moments. There are infinitely many such classes of PDFs defined by a finite set of moments. The choice of a particular class JET71 # IMechE 2006 Here, j, l, g, and d are the parameters of the system, which depend on the desired moments of xi,j . Solving this system then reduces to first finding the appropriate curve and then solving for the parameters, based on the iterative approach described by Hill et al. [18]. A MATLAB program was developed by Levy [19] to select the appropriate function and to calculate the parameters of the Johnson system. 2. In this approach, the non-Gaussian white-noise sequence xi,j is transformed to the desired output sequence zi,j using a linear filter hi,j . The input –output relationship of the filter is given by zi,j ¼ XX k hk, l xi k,j l (9) l In the spatial frequency domain, this relationship can be written as Z(vx , vy ) ¼ H(vx , vy )I(vx , vy ) (10) where I and Z are the Fourier transforms of xi,j and zi,j , respectively, and H is the transfer function of the filter hi,j . The power spectral density of the Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 68 A Kumar, F Sadeghi, and C M Krousgrill input and output are related as Sz (vx ,vy ) ¼ jH(vx , vy )j2 Sx (vx , vy ) (11) The basis for this procedure is as follows. 1. The output of the transformation that occurs in any of the equations (6), (7), or (8) is a whitenoise process. Therefore, the magnitude of the system transfer function jH(vx , vy )j can be determined directly from the desired spectral characteristics of the output sequence zi,j , using equation (11). 2. If the filter hi,j is known, there is a relationship between the moments of xi,j and zi,j , which is valid for any moving average process [10]. The input and output skewness and kurtosis are related as P P 4 i j hi,j Kuz ¼ P P (12) 2 (Kux 3) þ 3 2 h i j i,j P P 3 i j hi,j (13) Skz ¼ P P 3=2 Skx 2 h i j i,j Therefore, knowing hi,j the desired moments of xi,j can be calculated. The moments Skx and Kux are then used to identify the Johnson translator that is to be used among equations (6) to (8) and to calculate its parameters. Next, the translator can be used to generate xi,j from a Gaussian white-noise sequence. Surfaces are defined only in a statistical sense by their moments and ACF. Repeated generation of surfaces with the same desired statistics results in random surfaces. In this approach, this randomness can occur in two ways. 1. The initial Gaussian white-noise input hi,j is randomly generated each time. 2. When the filter is calculated using equation (10), phase information is lost because only the magnitude of the system transfer function jH(vx , vy )j is known. Phase can, therefore, be arbitrarily assigned in calculating H(vx , vy ), resulting in different output surfaces zi,j for the same starting input hi,j . In this work, a constant phase of zero is assigned to H(vx , vy ) and all randomness in the output zi,j arises from randomness in the input hi,j . The choice of output skewness and kurtosis values is not arbitrary. They must satisfy the condition Kuz 5 Sk2z þ 1 (14) Fig. 1 Regions in the (Skz , Kuz ) space covered by different frequency curves. In the above illustration, bx ¼ by ¼ 2 mm. The nominal dimensions of the surface are Lx ¼ Ly ¼ 1 cm A similar relationship must also be satisfied for the intermediate non-Gaussian white-noise sequence xi,j . This is written as Kux 5 Skx2 þ 1 (15) These relationships are illustrated in Fig. 1. Also shown in this figure are the regions in the skewness – kurtosis space covered by the different Johnson curves. The ST curve is the special case of the SB curve and is defined by the equality in equation (15). The SL curve marks the boundary between the SB region to one side, and the SU region to the other. Some theoretically attainable output (Skz , Kuz ) combinations that satisfy equation (14) cannot, in practice, be obtained because of the additional constraint (15) on the skewness and kurtosis of the intermediate sequence. Therefore, the ST curve defines the allowable (Skz , Kuz ) range. The allowable region is not constant but depends also on Rz (x, y). This is illustrated in Fig. 2, for values of b ranging from 1 to 3 mm. For fixed surface dimensions and grid spacing, it is found that increasing b decreases the range of allowable (Skz , Kuz ) combinations. Although it is not illustrated here, the location of the SL curve also changes with autocorrelation length. Figure 3 shows generated Gaussian surfaces with the same roughness sz , but having different autocorrelation lengths. Figure 4 depicts non-Gaussian surfaces with the same spectral characteristics, but with different values of skewness and kurtosis. The negatively skewed surface shown in the figure has a height distribution close to that of a measured surface whose properties Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 JET71 # IMechE 2006 Effects of surface roughness Fig. 4 Fig. 2 Variation of Kuz,min with Skz , corresponding to the ST curve, for different values of bx ¼ by ¼ b. The nominal dimensions are Lx ¼ Ly ¼ 1 cm were of interest. Height distributions of these generated surfaces are illustrated in Fig. 5. 3 FRICTIONLESS NORMAL CONTACT MODEL In the previous section, a method was described to generate rough surfaces that are defined by their first four central moments and ACF. In this section, Fig. 3 Generated surfaces with different autocorrelation lengths. All surfaces have Gaussian height distributions, and sz ¼ 40 mm. All dimensions are in metres JET71 # IMechE 2006 69 Generated surfaces with different values of skewness and kurtosis. All surfaces have sz ¼ 40 mm and are isotropic with bx ¼ by ¼ 1 mm. All dimensions are in metres the numerical approach to study normal contact between rough bodies is discussed. In this method, only the surfaces of the contacting bodies are discretized. The profiles of the surfaces are generated using the method described in the previous section. In addition to the surface profiles, the contact model requires a relationship between surface normal pressures and the normal displacements at all locations on the surface grid. It is assumed that the contacting surfaces are frictionless, so that only normal pressure acts between them. The primary objective of this contact model is to obtain the relationship between normal load and displacement. It is expected that the effects of friction on this relationship is small. The method used to simulate contact takes into account the Fig. 5 Probability density of heights for generated surfaces shown in Fig. 4 Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 70 A Kumar, F Sadeghi, and C M Krousgrill potential interactions of all surface points and predicts the contact regions and pressures under prescribed normal load. The contacting region is discretized as a mesh of small elements over the entire domain of potential contact. The problem then reduces to finding the displacements ui,j and pressures pi,j at each element (i, j) for a surface whose undeformed topography is given by zi,j . The relation between pressure and displacement is given by the influence coefficients, C(i,j),(k,l) , as ui,j ¼ XX k (16) C(i,j),(k,l) pk,l l 1 V ¼ 2 ðð where zm is the initial level of the smooth surface at the onset of contact and is equal to the maximum height of the rough surface. The prescribed displacement dz (x, y, d) is derived from the geometric interference condition of equation (19), and is written as dz (x, y, d) ¼ d f (x, y) (22) ~ i,j ¼ p pi,j Enc (23) where Enc is equal to the effective Young’s modulus of the contacting bodies. The displacement is nondimensionalized as u~ i,j ¼ ui,j sz (24) where sz is the standard deviation of rough surface heights. Similarly p(x, y)(uz1 (x, y) þ uz2 (x, y) dx dy) V di,j d~ i,j ¼ sz p(x, y)(dz1 (x, y) þ dz2 (x, y) dx dy) ðð p(x, y)uz (x, y) dx dy ð ðV (25) (17) where uz1 and uz2 are the displacements of the two bodies and dz1 and dz2 are the prescribed displacements. As only normal frictionless contact is considered in which the contact pressure acts unidirectionally, the results depend only on relative separation between the contacting surfaces. Therefore, the contact is modelled as occurring between a rigid smooth surface and an equivalent deforming rough surface. Equation (17) can then be rewritten as 1 V ¼ 2 (21) z(x, y) V ðð  f (x, y) ¼ zm To non-dimensionalize the problem, let If the pressures pi,j on each element are known, equation (16) can be used to find the deflections and solve the contact problem. However, the individual pressures are usually unknown, and an iterative procedure has to be adopted for the solution. The above contact problem can also be formulated in its variational form as one of the minimization of complementary potential energy [20]. The total complementary potential energy for two contacting bodies which behave linearly is  of the smooth and rough surfaces given by p(x, V y)dz (x, y) dx dy (18) where uz ¼ uz1 þ uz2 is the total displacement and dz ¼ dz1 þ dz2 is the total prescribed displacement. The geometrical constraints on uz are uz (x, y) þ f (x, y) uz (x, y) þ f (x, y) d ¼ 0 (within contact area) (19) d . 0 (outside contact area) (20) Here d is the rigid body displacement under the applied normal load and f (x, y) the initial separation The non-dimensional influence coefficients are written as Enc C~ (i,j),(k,l) ¼ C(i,j),(k,l) sz (26) The variational formulation of the problem in nondimensional form is then XX 1XX ~ i,j ~ k,l Min V~ ¼ p C~ (i,j),(k,l) p 2 i j k l XX  ~ i,j 5 0 ~ i,j d~ i,j s.t. p p  i ! (27) j The method of steepest descent is used to solve the above equation, using the procedure illustrated in Fig. 6. This yields displacement u~ i,j and pressure p~ i,j . Solving the problem in non-dimensional form has the following advantages. 1. The numerical solution is easier to implement. 2. The effects of varying standard deviation sz and Young’s modulus Enc across a wide range of values can be obtained without having to perform additional contact simulations. Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 JET71 # IMechE 2006 Effects of surface roughness 71 roughness heights and Young’s modulus. This can be written as (30) pi,j / sz Enc If the dimensional pressure field p0i,j corresponding to a particular set of parameters s0 and E0 is known, then the solutions corresponding to arbitrary parameters sz and Enc are related to p0i,j as 0 pi,j (d=sz ) pi,j (d=s0 ) ¼ sz Enc s0 E 0 (31) The above relationship also applies to the engagement pressure P required to obtain a particular value of non-dimensional displacement. Similarly, the fractional area in contact Ar for arbitrary parameters is given by Ar Fig. 6 Illustration of steps used to solve normal contact problem C(i,j),(k,l) ¼ C^ (i,j),(k,l) Enc (28) where C^ (i,j),(k,l) are fixed and reflect the fact that all properties affecting influence coefficients (bulk geometry, boundary conditions, and Poisson’s ratio) except Enc are fixed. The above equation therefore reflects the effect of Young’s modulus on the pressure – deflection relationship. This relationship can be made more explicit by rewriting equation (16) in terms of the non-dimensional displacement field as sz u~ i,j ¼ X X C^ (i,j),(k,l) k l Enc pk,l (29) It follows from this equation that the dimensional pressure corresponding to a given pressure and displacement solution in non-dimensional form is proportional to both the standard deviation of JET71 # IMechE 2006    d d ¼ A0r sz s0 (32) which is independent of elastic modulus. It is to be emphasized that the above discussion and the relationships in equations (31) and (32) apply only to surfaces with fixed skewness and kurtosis. 4 The effects of varying sz and Enc , keeping other properties fixed, is seen by writing the influence coefficients as if they depend only on these two parameters. For a linear material, the influence coefficients C(i,j),(k,l) can then be written as  CONTACT SIMULATION RESULTS The influence coefficients C(i,j),(k,l) of equation (7) are obtained using a finite element model which has the nominal geometry and material properties of the contacting bodies, but excludes the roughness effects. In this work, the contact model was applied to the type of contact that occurs during automobile disc brake engagement, in which a relatively lowmodulus friction material comes in contact with a high-modulus brake rotor. Because the elastic moduli of the brake rotor and friction material differ by two orders of magnitude, only the more compliant body is modelled. The features of this model are shown in Table 1. In this application, one is mainly interested in the variation of normal stiffness with engagement load. Table 1 Characteristics of ANSYS model used to calculate influence coefficients Dimensions Grid spacing Material properties Length: 9 mm Number of elements in length: 60 Number of elements in width: 60 Number of elements in thickness: 5 Linear, isotropic Width: 9 mm Thickness: 1 cm Young’s modulus: 5.11  108 Pa Poisson’s ratio: 0 Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 72 A Kumar, F Sadeghi, and C M Krousgrill This stiffness is normalized as kP ¼ L dF An dd (33) where F is the normal force, d the displacement of the rigid surface, L the thickness of the compliant body, and An the nominal area. The superscript P on k denotes that it has units of pressure. The variation of normalized stiffness with pressure is shown in Fig. 7 for surfaces with different height distributions. The results are based on the averaged load – deflection curves for 60 surfaces with each set of statistics. The variance in stiffness for each case is much smaller than the average difference between the cases. All surfaces have the same bulk properties. At high loads, the normalized stiffness approaches the elastic modulus of the bulk material. The differences in stiffness results are then because of the differences in topography. Negative skewness generally results in stiffer surfaces, whereas positive skewness results in surfaces that are more compliant compared with Gaussian surfaces. It might appear surprising that for Kuz ¼ 4:5, the zero skewness surface is more compliant than the surface with high positive skewness. This is because for Kuz ¼ 4:5, the zero skewness surface has sharper, and therefore more compliant, peaks than the surface with high positive skewness. This can be seen from the height distributions in Fig. 5. The total deflection d is what is measured during a force – deflection test. Although the underlying displacement field is a result of asperity and bulk Fig. 8 Lumped spring model of normal contact effects, it is physically meaningless to separate the total deflection into an asperity and bulk deflection. However, in this work it was necessary to separate the total measured deflection into an asperity and a bulk component. This was required to estimate the material properties from the experimental force – deflection curve. This is illustrated in Fig. 8. The total deflection is written as the sum of the asperity deflection da and the bulk deflection db as (34) d ¼ da þ db In this model, the bulk depends only on the material properties and geometry, whereas the asperity deflection depends is also affected by surface roughness. To make this separation into asperity and bulk components, it is necessary to assume a model for both the asperity and bulk springs. The deflections are then obtained from least-squares estimation using the experimental results. 4.1 Mean plane deflection and bulk behaviour Although the bulk deflection defined previously does not correspond to a physically observed deflection, a suitable candidate exists that is independent of surface roughness. This is the deflection of the surface mean plane. This is calculated as db ¼ Fig. 7 Variation of normalized stiffness with pressure for different surfaces. Each curve corresponds to results averaged over 60 surfaces with similar statistics PN PM i¼1 j¼1 NM ui,j (35) The variation of db with nominal pressure P is shown in Fig. 9 for different surfaces. Each colour corresponds to a surface with different statistics (skewness and kurtosis) and there are five surfaces for each pair of skewness and kurtosis values. It was determined that the relationship between P and db for every surface lies on the straight line Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 JET71 # IMechE 2006 Effects of surface roughness 73 Interchanging the order of summation, equation (40) can be rewritten as db ¼ PN k¼1 PM l¼1 pk,l PN PM i¼1 j¼1 C(k,l),(i,j) NM (41) Using equation (38), equation (41) becomes db ¼ L Enc PN k¼1 PM l¼1 pk,l NM (42) Because the areas of all elements on the contact grid are equal to a constant DA N X M X pk,l DA ¼ F ¼ NMDAP (43) k¼1 l¼1 Variation of mean plane and asperity deflections with nominal pressure for different surfaces. All results are for surfaces with sz ¼ 30 mm and bx ¼ by ¼ 1 mm Fig. 9 shown in the figure. Assuming that the relationship between P and db is linear, the modulus corresponding to the linear relationship can be written as Eu ¼ PL db (36) This is calculated to be 5.11  108 Pa, the same value as the Young’s modulus Enc. This behaviour can be explained as follows: for the case of n ¼ 0 used in calculating influence coefficients, a uniform pressure P applied over the entire nominal area results in a constant displacement field ui,j ¼ PL Enc (37) Comparing with the pressure – displacement relation in equation (16), it is found that M N X X C(i,j),(k,l) ¼ k¼1 l¼1 L Enc (38) Because the influence coefficients are for a linear, conservative system and the contact grid is uniform (39) C(i,j),(k,l) ¼ C(k,l),(i,j) Therefore, the bulk displacement in equation (35) can be written as db ¼ PN PM PN i¼1 j¼1 JET71 # IMechE 2006 k¼1 PM NM l¼1 C(k,l),(i,j) pk,l (40) Using equation (43) in equation (42) yields db ¼ PL Enc (44) which is the relationship observed in Fig. 9. Although this derivation is for a uniform contact grid with constant element area DA, it applies also to a general grid. For a more general grid, equation (16) must be reformulated in terms of the normal force on each element. The changes, starting from equation (38), result in the same conclusions. Equation (44) then refers to the properties of a linear, conservative material with n ¼ 0, and not any particular choice of contact grid. For the case of n . 0, equations (37) and (38) are not exactly satisfied and the deviation from the equality depends on the geometry and the value of n. To study the difference for n . 0, simulations were performed with the features described in Table 1 but with n ¼ 0:3. It was found that the relationship between db and P was still linear, but the modulus calculated from equation (36) differed slightly between surfaces. The deviation of this modulus from the material modulus Enc was approximately 5 per cent. 4.2 Asperity behaviour As described previously, the separation of the total deflection into asperity and bulk deflections according to equation (34) is primarily used as a tool to estimate the Young modulus and does not have any physical significance. However, in the previous section it was also shown that all surfaces with given influence coefficients, independent of surface topography, exhibit the surface mean plane deflection behaviour of equation (44). In this section, the effects of surface topography total deflection d are illustrated. Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 74 A Kumar, F Sadeghi, and C M Krousgrill This is done, however, by studying the behaviour of d that is in excess of mean plane deflection. This difference is defined as the asperity deflection da ¼ d (45) db where db is defined in equation (35). The variation of da for different surfaces is also shown in Fig. 9. Asperity deflection as defined in equation (45) dominates the response at low loads, but as the load is increased the asperities become stiffer and the constant-stiffness bulk accounts for a larger fraction of the deflection. It is seen that da not only varies with skewness and kurtosis as expected, but also varies for surfaces with very similar skewness and kurtosis. The differences in the behaviour of da for different surfaces can be best explained by investigating the differences in their undeformed height distributions. Figure 10 shows the height distributions of the same surfaces in the undeformed state. In this figure, the level at which the cumulative distribution is studied is defined as the depth below the highest location on the surface zmax . It is seen by comparing Figs 9 and 10 that surfaces with negative skewness are stiffer than surfaces with positive skewness. For zero skewness, increasing kurtosis decreases the stiffness of the surface. It is also seen in Fig. 10 that the height distributions of surfaces with approximately the same statistics (e.g. with Sk ¼ 0 and Ku ¼ 3:0) have the same slope near the middle of the distribution, but differ at the ends. This accounts for the large differences in asperity force –deflection behaviour shown for surfaces with similar statistics. Although the statistics calculated over the entire range of heights for two surfaces might be very similar, the statistics of points near the tallest peaks most affect stiffness. The variation in this region is likely much higher than the averaged statistics of the surface reflect. It is also possible that slight variations in the spectral characteristics of the output surfaces account for differences in the asperity deflections shown in Fig. 9. However, for the contact model used in this work it is not expected that spectral characteristics have much effect on force – deflection response. This point is made by Fig. 11, which shows the normalized stiffness for surfaces with the same height distribution but with two different autocorrelation lengths. Even though results such as real contact area, maximum contact pressure, and average size of contact spots are different (these results are not shown), the stiffness of the entire surface is independent of autocorrelation length. 4.3 Real contact area and contact pressure In the previous section, it was shown that differences in total (or asperity) force –deflection behaviours are best explained by differences in the cumulative height distribution of surfaces. This raises the question of whether quantitative predictions of the force – deflection (or stiffness) behaviour can be made directly from knowledge of the cumulative height distribution (or height PDF). Fig. 11 Fig. 10 Height distributions of different generated surfaces in the undeformed state. The surfaces used are the same as in Fig. 9, with sz ¼ 30 mm and bx ¼ by ¼ 1 mm Normalized stiffness for contact model with different autocorrelation lengths. Each curve shows average results for 60 Gaussian surfaces with sz ¼ 24:5 mm and Lx ¼ Ly ¼ 9 mm Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 JET71 # IMechE 2006 Effects of surface roughness Thus far, no such relationship has been found. One reason for doubting the existence of such a relationship is that the range and distribution of local contact pressures varies with surface characteristics. For fixed autocorrelation length, the variation of maximum pressure with applied nominal pressure for different surfaces is shown in Fig. 12. Negatively skewed surfaces have lower contact pressures as compared with positively skewed surfaces at the same engagement load. The zero skewness surface has the highest contact pressures because it has the tallest peaks, as can be seen from the height distribution in Fig. 10. The limited explanatory ability of the moments skewness and kurtosis is illustrated by the behaviour of the surfaces with zero skewness, emphasizing further the importance of studying the complete distributions shown in Fig. 10. Figure 13 illustrates the variation in real contact areas for different surface statistics. Negatively skewed surfaces, in general, have a higher real contact area than positively skewed surfaces, at the same nominal engagement pressure. Furthermore at low pressures, the order of increasing contact area exactly corresponds to the order of decreasing maximum pressure shown in Fig. 12. However, at a pressure of approximately 3 MPa, the two real area curves cross while the corresponding pressure curves remain distinct. This indicates that the PDFs of contact pressure distributions are quite different for different surfaces. Figure 14 depicts the effect of changing nominal area on the stiffness calculated using the contact model. Both models have Gaussian surfaces with Fig. 12 Variation of maximum pressure with nominal pressure for surfaces with different skewness and kurtosis. Each curve shows the average results for 60 surfaces with similar statistics JET71 # IMechE 2006 Fig. 13 75 Variation of real contact area with nominal pressure for surfaces with different skewness and kurtosis. Each curve shows the average results for 60 surfaces with similar statistics the same roughness statistics. For the case of n ¼ 0 shown, normalized stiffness is the same for both areas. This means that both asperity and bulk stiffness scale linearly with area. For n . 0, it is expected that the bulk behaviour will not exactly satisfy equation (44) and the magnitude of the variation would depend on the geometry and material properties. However, asperity stiffness will scale linearly with area. This makes it unlikely that bodies with non-zero Poisson’s ratio have stiffness that scale as in Fig. 14. Fig. 14 Normalized stiffness results for contact models with different nominal areas. Each curve shows average results for 60 Gaussian surfaces with sz ¼ 24:5 mm and bx ¼ by ¼ 1 mm Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 76 5 A Kumar, F Sadeghi, and C M Krousgrill SUMMARY In this work, a normal contact model was described that can be used for a variety of applications. The contact model requires discretized profiles of the contacting surfaces and relations between the normal pressure and surface normal displacements at all locations on the surface as inputs. The rough surfaces were generated by transforming a Gaussian white-noise sequence using one of the four Johnson functions and then transforming the resulting output through a linear filter. The spectral characteristics of the output surface are controlled by the choice of a linear filter, whereas the height distribution of the surface is controlled by the choice of the Johnson transformation. The resulting height distribution is defined by the first four central moments of the surface. As the height distributions generated in this work are defined uniquely by its first four moments, the choice of height distribution for the output surface is not arbitrary. Rather, the possible height distributions that can be generated using this approach are restricted to a class of distributions whose members are defined by the first four moments. Replacing the Johnson translator system in this work by a different set of functions will result in a different class of height distributions. In general, two surfaces with the same first four moments, one generated using the Johnson system and the other using an alternate system, will not have the same height distribution. In general, negatively skewed surfaces are stiffer than Gaussian surfaces, whereas positively skewed surfaces are more compliant. However, to associate a particular combination of standard deviation, skewness, and kurtosis with its characteristic force – deflection curve, it is necessary to restrict the study to a particular set of height distributions such as the one generated by the Johnson system. More generally, when comparing surfaces which might not belong to such a restricted set of height distributions, it is more appropriate to compare the height distributions of the surfaces than to compare only the standard deviation, skewness, and kurtosis. Force– deflection results for normal contact are more closely related to the height distribution than to any single combination of standard deviation, skewness, and kurtosis. Even for a fixed class of distributions such as the Johnson system, the lower moments have limited explanatory ability as illustrated by the zero skewness surfaces, which had sharper peaks than surfaces with positive skewness. Furthermore, the contact stiffness is most strongly affected by the statistics of surface heights near the tallest peaks. The normal contact model can be used to obtain detailed results such as the distribution of contact pressure, the normal displacement field of the surface, and the real contact area for a given normal load. These detailed results also depend, in addition to the height distribution of the undeformed surface, on the spectral characteristics of the surfaces in contact. However, the macroscopic load –deflection behaviour is independent of small changes in autocorrelation length. This result, however, does not mean that spectral characteristics do not affect the macroscopic force – deflection behaviour. In this work, when the autocorrelation length was varied the sampling length used in the contact model was kept fixed. Furthermore, owing to limitations in calculating the influence coefficients the contact grid was coarse at 60  60 points. It is conceivable that a smaller sampling length that varies to account for effects at different length scales will illuminate differences in force – deflection response for surfaces with different spectral characteristics. The approach used for normal contact in this and the related work referenced in the survey of literature has the advantage over statistical models that it accounts for elastic interactions between points on the surface. The influence coefficients in this work are calculated using a finite element model of the bodies in contact. In this work, Poisson’s ratio was set equal to zero. As a consequence, the load – deflection results obtained from contact simulations can be normalized so that the results are independent of nominal area of the contacting rough surface. The bulk deflection is defined as the deflection of the mean plane. This means that when a Poisson’s ratio of zero is used to calculate influence coefficients in the contact model, the mean surface displacement at any applied pressure depends only on the pressure and not on the surface topography. For non-zero Poisson’s ratio, this relationship is not exactly satisfied and the difference depends on the nominal geometry and Poisson’s ratio. For such a case, the mean surface displacement also varies slightly from surface to surface. The bulk deflection is the total deflection that would be observed for perfectly smooth contacting bodies. ACKNOWLEDGEMENT The authors express their deepest appreciations to Bosch for their support and technical expertise during the course of this project. REFERENCES 1 Whitehouse, D. J. and Archard, J. F. The properties of random surfaces of significance in their contact. Proc. R. Soc. London, Series A, 1970, 316(1524), 97 – 121. Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014 JET71 # IMechE 2006 Effects of surface roughness 2 Greenwood, J. A. and Williamson, J. B. P. Contact of nominally flat surfaces. Proc. R. Soc. London, Series A, 1966, 295, 300 – 319. 3 Greenwood, J. A. and Tripp, J. H. The contact of two nominally flat rough surfaces. Proc. Instn Mech. Engrs, 1971, 185, 625 – 633. 4 Bush, A. W., Gibson, R. D., and Thomas, T. R. The elastic contact of a rough surface. Wear, 1975, 35, 87 – 111. 5 Chang, W. R., Etsion, I., and Bogy, D. B. An elastic – plastic model for the contact of rough surfaces. ASME J. Tribol., 1987, 109, 257 – 263. 6 Majumdar, A. and Tien, C. L. Fractal characterization and simulation of rough surfaces. Wear, 1990, 136, 313 – 327. 7 Patir, N. A numerical procedure for random generation of rough surfaces. Wear, 1978, 47, 263 – 277. 8 Hu, Y. Z. and Tonder, K. Simulation of 3-dimensional random rough surface by 2-dimensional digital filter and Fourier analysis. Int. J. Mach. Tools Manuf., 1992, 32, 83– 90. 9 Chilamakuri, S. K. and Bhushan, B. Contact analysis of non-Gaussian random surfaces. Proc. Instn Mech. Engrs, Part J: J. Engineering Tribology, 1998, 217, 19 – 32. 10 Bakolas, V. Numerical generation of arbitrarily oriented non-Gaussian three-dimensional rough surfaces. Wear, 2003, 254(5 – 6), 546 – 554. 11 Karpenko, Y. A. and Akay, A. A numerical model of friction between rough surfaces. Tribol. Int., 2001, 34, 531 – 545. 12 Bhushan, B. Contact mechanics of rough surfaces in tribology: multiple asperity contact. Tribol. Lett., 1998, 4, 1– 35. 13 Borri-Brunetto, M., Chiaia, B., and Ciavarella, M. Incipient sliding of rough surfaces in contact: a multiscale numerical analysis. Comput. Meth. Appl. Mech. Eng., 2001, 190, 6053 – 6073. 14 Kalker, J. J. Variational principles of contact elastostatics. J. Inst. Math. Appl., 1977, 20, 199 – 219. 15 Peng, W. and Bhushan, B. A numerical threedimensional model for the contact of layered elastic/ plastic solids with rough surfaces by a variational principle. ASME J. Tribol., 2001, 123, 330– 342. 16 Stanley, H. M. and Kato, T. An FFT-based method for rough surface contact. ASME J. Tribol., 1997, 119, 481 – 485. 17 Johnson, K. L. Contact mechanics, 1985 (Cambridge University Press, Cambridge). 18 Hill, I. D., Hill, R., and Holder, R. L. Fitting Johnson curves by moments. Appl. Stat., 1976, 25, 180 – 189. 19 Levy, B. Johnson curve fitting. March 2004. MATLAB Central File Exchange on the World Wide Web: http:// www.mathworks.com/matlabcentral/fileexchange/ loadCategory.do. 20 Bhushan, B. Contact mechanics of rough surfaces in tribology: multiple asperity contact. Tribol. Lett., 1998, 4, 1– 35. JET71 # IMechE 2006 77 APPENDIX Notation An C(i,j),(k,l)  di,j Enc Eu fi,j F hi,j H(vx, vy) I(vx, vy) ka, l kb kP Kuz L mi pi,j p~ i,j P Rz (x, y) Sz (vx, vy) Sx (vx , vy ) Skz ui,j u~ i,j u s V zi,j zm Z(vx, vy) bx , by d da , db hi,j mz, sz n j, l, g, d xi,j nominal area of bulk material influence coefficients prescribed displacement field Young’s modulus used in calculating influence coefficients Young’s modulus calculated from mean plane deflection initial separation between rough surfaces normal engagement load linear filter used in surface generation Fourier transform of hi,j Fourier transform of xi,j asperity spring constants bulk spring constant normalized stiffness Kurtosis of zi,j thickness of bulk material ith central moment of un-deformed surface height distribution contact pressures non-dimensional contact pressures nominal pressure autocorrelation function power spectral density of zi,j power spectral density of xi,j skewness of zi,j normal displacements of surface non-dimensional displacements of surface normal deflection of the surface mean plane complementary potential energy of contacting bodies rough surface heights maximum height of rough surface Fourier transform of zi,j autocorrelation lengths deflection of the rigid surface asperity and bulk deflections in lumped model Gaussian white-noise input sequence mean and standard deviation of zi,j Poisson’s ratio Johnson system parameters non-Gaussian white-noise sequence Proc. IMechE Vol. 220 Part J: J. Engineering Tribology Downloaded from pij.sagepub.com at INDIAN INST SCI on July 27, 2014