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