Finite Element Analysis of Three-Dimensional Elasto-Plastic Sinusoidal
Contact and Inclusion in a Multi-Scale Rough Surface Contact Model
Except where reference is made to the work of others, the work described in this thesis is
my own or was done in collaboration with my advisory committee. This thesis does not
include proprietary or classi ed information.
Vijaykumar Krithivasan
Certi cate of Approval:
Pradeep Lall
Thomas Walter Professor
Mechanical Engineering
Robert L. Jackson, Chair
Assistant Professor
Mechanical Engineering
Hareesh V. Tippur
Alumni Professor
Mechanical Engineering
George T. Flowers
Interim Dean, Graduate School
Finite Element Analysis of Three-Dimensional Elasto-Plastic Sinusoidal
Contact and Inclusion in a Multi-Scale Rough Surface Contact Model
Vijaykumar Krithivasan
A Thesis
Submitted to
the Graduate Faculty of
Auburn University
in Partial Ful llment of the
Requirements for the
Degree of
Master of Science
Auburn, Alabama
May 10, 2008
Finite Element Analysis of Three-Dimensional Elasto-Plastic Sinusoidal
Contact and Inclusion in a Multi-Scale Rough Surface Contact Model
Vijaykumar Krithivasan
Permission is granted to Auburn University to make copies of this thesis at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Vijaykumar Krithivasan, son of Mr. B. Krithivasan and K. Sathyabama was born in
Bangalore, India on Macrh 26, 1981. He graduated high school from Sri Sankara Senior
Secondary School, Chennai in April 1999. He graduated from University of Madras, Madras,
India with a Bachelor of Engineering degree in July 2003. He began his graduate studies in
Mechanical Engineering at Auburn University in January 2005.
iv
Thesis Abstract
Finite Element Analysis of Three-Dimensional Elasto-Plastic Sinusoidal
Contact and Inclusion in a Multi-Scale Rough Surface Contact Model
Vijaykumar Krithivasan
Master of Science, May 10, 2008
(B.E., University of Madras, 2003)
104 Typed Pages
Directed by Robert L. Jackson
Researchers have developed many models to simulate the elasto-plastic contact of
spheres. However, there does not appear to exist a closed-form analytical model for elasto-
plastic three dimensional sinusoidal contact. This work uses a  nite element model (FEM)
to characterize elasto-plastic sinusoidal contact. Although at initial contact the spherical
and sinusoidal cases are very similar and can both be described by the classic elastic Hertz
contact case, once the contact is pressed past a certain range of inelastic deformation the
two cases are very di erent. The FEM model is used to produce equations which can be em-
ployed to approximately relate the area of contact to the contact pressure for elasto-plastic
sinusoidal contact. The equations are obtained by  tting to the FEM results and existing
elastic solutions to sinusoidal contact. An empirical expression for the average pressure
which causes complete contact between elasto-plastic sinusoidal contacts was also devel-
oped. The results showed that the required pressure for complete contact is signi cantly
less in the elasto-plastic regime than the elastic regime. In addition, this pressure is shown
to be greater than the traditional hardness, of 3?Sy.
v
One of the major motivations for this work was to generate a model that could be
used in sinusoidal and frequency based rough surface contact models. A multiscale model
is a non-statistical model and non-fractal that is used to describe normal contact between
rough surfaces featuring multiple scales. The empirical equations developed in the sinusoidal
contact model are used to characterize asperity contact in the multiscale contact model.
Based on this, predictions are made for contact area as a function of applied load. It
was interesting to note that the real area of contact versus the applied load exhibits a
linear relationship for both elastic and elasto-plastic cases up until the surface is completely
 attened out. As expected, the real area of contact undergoing elasto-plastic deformation
predicted by the multiscale model was higher than when the surface is undergoing elastic
deformation. For a given applied load it was also found that the lower frequency ranges, as
opposed to higher frequency ranges, dictated the predicted level of real contact area.
vi
Acknowledgments
I am deeply indebted to Dr. Robert.L.Jackson for his technical, educational, and
moral support throughout my master?s program. I am deeply in uenced by his research
and teaching philosophy which made my graduate work an educational and professionally
enriching experience. I am grateful to Dr. Hareesh.V.Tippur for providing a great deal of
encouragement and technical support. I would also like to thank Dr. Pradeep Lall for his
encouragement and for agreeing to be on my master?s thesis committee.
I would like to thank my colleagues Mr. Everett Wilson, Mr. Jeremey Dawkins and
Mr. Santosh Angadi at the Multiscale Tribology Laboratory here at Auburn university.
I am grateful to my friends Ms. Anjeli Singh, Dr. Vivek Krishnan, Mr. Kashyap
Yellai, Mr. Anand Sankarraj, Mr. Karthik Narayanan, Ms. Kavita Arumugam, Mr. Sankar
Balasubramanian, Mr. Raghu Viswanathan, Mr. Naren Parinar, Ms. Rebecca Ibrahim,
and Mr. Sumit Sen for there support and encouragement.
Finally, I would like to thank my family for their unwavering support and encourgament
throughout my master?s program. Many sacri ces made by my parents Krithivasan and
Satyabama cannot be described in words. My brother, Dr. Ramkumar Krithivasan, has
been a pillar of support throughout my life. His achievements have always been a motivation
for me to move forward. I am forever indebted to my family.
vii
Style manual or journal used LATEX: A Document Preparation System by Leslie
Lamport (together with the style known as  aums ) and Bibliography as per Tribology
Transactions
Computer software used TEX (speci cally LATEX), VossPlot, ANSYS, MATLAB
7.0.4, MS O?ce Excel 2003 and the departmental style- le aums.sty
viii
Table of Contents
List of Figures xi
List of Tables xiii
nomenclature xiv
1 Introduction 1
1.1 Organization of thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2 Literature review 3
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Spherical contact models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Elastic sinusoidal contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Elasto-plastic sinusoidal contact . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 Statistical contact model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.6 Fractal contact model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3 Finite element model 20
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.2 Building the solid model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3 Mesh convergence test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4 Simulation methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.5 Veri cation of model accuracy . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4 elasto-plastic fem results 29
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Parametric study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.3 Critical elasto-plastic pressure estimation . . . . . . . . . . . . . . . . . . . 34
4.4 Real contact area estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.5 The empirical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5 Multiscale Rough Surface Contact Model 50
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.2 Multiscale model foundation . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3 Elasto-plastic multiscale model . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.4 Frequency spectrum for di erent surfaces . . . . . . . . . . . . . . . . . . . 54
5.5 Frequency level iteration versus contact area . . . . . . . . . . . . . . . . . . 58
5.6 Elasto-plastic and elastic comparison . . . . . . . . . . . . . . . . . . . . . . 62
ix
6 Conclusions 67
Bibliography 69
Appendices 72
A Derivation of critical values for sinusoidal contact 73
B matlab code for contact area versus frequency iteration 76
C matlab code for frequency spectrum and load versus contact area 84
x
List of Figures
2.1 Topographical contour plot of the sinusoidal surface geometry. . . . . . . . . 4
2.2 Cross-section of sinusoidal type contact. . . . . . . . . . . . . . . . . . . . . 6
2.3 Diagram of the hypothetical progression of change in hardness with geometry
for spherical and sinusoidal contact. . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Compressed cavity of material. . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5 Islands de ned by a horizontal section. . . . . . . . . . . . . . . . . . . . . . 18
3.1 Uniform mesh on the rigid  at formed by 21x21 array of elements. . . . . . 22
3.2 Schematic of degree of freedom restraints used for the one quarter sinusoidal
FEM model (actual modeled geometry is much longer in the z direction). . 24
3.3 Various stages of contact. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 Deformation history. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5 Comparison of elastic FEM results with JGH model. . . . . . . . . . . . . . 28
4.1 pSy versus contact area ratio (A?F2) for di erent yield strength. . . . . . . 31
4.2 Contact area ratio (A?F2) versus ?p=p? for di erent yield strength. . . . . . 32
4.3 The e ect of Sy=E0 on pressure to cause complete contact. . . . . . . . . . . 35
4.4 The e ect of  = on pressure to cause complete contact. . . . . . . . . . . . 37
4.5 Contact area ratio (A?f2) versus ?p=p?ep for di erent yield strength values. . 39
4.6 Contact area ratio (A?f2) versus ?p=p?ep for di erent values of Young?s Modulus
E. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.7 Contactarearatio(A?f2)versus ?p=p?ep fordi erentvaluesofthedimensionless
geometric quantity  = . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
xi
4.8 Comparison of the FEM results and the elasto-plastic sinusoidal contact
model (marked by lines) given by Eq. 4.8 for di erent yield strength val-
ues. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.9 Comparison of the FEM results and the elasto-plastic sinusoidal contact
model (marked by lines) given by Eq. 4.8 for di erent Young?s modulus values. 48
4.10 Comparison of the FEM results and the elasto-plastic sinusoidal contact
model (marked by lines) given by Eq. 4.8 for di erent sinusoidal geometries. 49
5.1 Flow chart of iterative asperity contact model. . . . . . . . . . . . . . . . . 53
5.2 Frequency spectrum of surface 1. . . . . . . . . . . . . . . . . . . . . . . . . 55
5.3 Frequency spectrum of surface 2. . . . . . . . . . . . . . . . . . . . . . . . . 56
5.4 Frequency spectrum of surface 3. . . . . . . . . . . . . . . . . . . . . . . . . 57
5.5 Contact area ratio versus frequency iteration for surface 1. . . . . . . . . . . 59
5.6 Contact area ratio versus frequency iteration for surface 2. . . . . . . . . . . 60
5.7 Contact area ratio versus frequency iteration for surface 3. . . . . . . . . . . 61
5.8 Elastic and Elasto-Plastic contact area versus applied load for surface 1. . . 64
5.9 Elastic and Elasto-Plastic contact area versus applied load for surface 2. . . 65
5.10 Elastic and Elasto-Plastic contact area versus applied load for surface 3. . . 66
xii
List of Tables
3.1 Material properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.1 Parametric study based on the Yield Strength Sy . . . . . . . . . . . . . . . 30
4.2 Parametric study based on the Young?s Modulus E0 . . . . . . . . . . . . . 33
4.3 Parametric study based on the Geometric ratio  = . . . . . . . . . . . . . 33
4.4 Variation of p?epp? (FEM) and p?epp? (Fit Eq. 4.1) with respect to Sy=E0 . . . . 36
4.5 Variation of p?epp? (FEM) and p?epp? (Fit Eq. 4.1) with respect to  = . . . . . 38
xiii
Nomenclature
?A individual asperity contact area
?Ai single asperity contact area
?Fi single asperity contact load
?P contact load
?p average pressure over entire surface
?i amplitude at the ith level
?i areal asperity density at the ith level
? scaling parameter
 wavelenth (1f)
!c critical interference between hemisphere and surface
! interference between hemisphere and surface
 (z) gaussian distribution
 m;n random phase
?lateral lateral stress
? standard deviation of surface heights
?y Shear strength
xiv
4 amplitude of sinusoidal surface
? Poisson?s ratio
Ac critical area of contact value at onset of plastic deformation
AE area of contact
AJGH area of contact from model by Johnson et al. [1]
AKE area of contact from model by Kogut and Etsion [2]
Ar real area of contact
Cs non-dimensional load constant
D fractal dimension 2<D<3
E elastic modulus
E0 reduced elastic modulus
F contact load
f spatial frequency (reciprocal of wavelength)
F? Friction force
FE contact force
fi reciprocal of wavelength at the ith level
G fractal roughness
h height of sinusoidal surface from base
xv
L sample length
M number of superimposed ridges
n frequency index
nmax upper limit of n
p average contact pressure
p? average pressure for complete contact
p?ep critical pressure for elasto-plastic complete contact
R radius of curvature of tip of sinusoidal surface or hemispherical asperity
z height of asperity density measured from the mean of asperity heights
Sy yield strength
xvi
Chapter 1
Introduction
Computational modeling of contact between surfaces, or contact mechanics, has been
realized as an important  eld for many years, probably starting sometime around the classic
solution of elastic spherical contact developed by Hertz in 1888 (originally developed to
model optical contacts) [3]. Since then the  eld has evolved to consider rough surface
contact between mechanical devices and more recently in MEMS. Many earlier models [4] 
[10] employed mostly spherical contact models to simulate individual asperity peaks on a
surface. Several used alternative approaches to model periodic surface roughness as sine
waves [1, 11]. These  rst models considered only elastic contact. However, due to the high
loads often seen in contact situations, the stresses can cause plastic deformation as well.
The FFT based contact methodology outlined by Stanley and Kato [12] made use of these
elastic contact models to model rough surface contact by  rst transforming the surface data
into the frequency domain.
In keeping with the need for a new deformation model which could be easily incorpo-
rated into the frequency based rough surface contact models, the current work develops an
elasto-plastic sinusoidal contact model. Empirical equations were developed that fully char-
acterize a sinusoidal surface undergoing elasto-plastic deformation. The developed equations
were used inside the multiscale model foundation to predict the real area of contact as a
function of applied load. Since the multiscale model uses Fourier series, and it consists series
of sine waves, it would be logical to model the asperities as an array of stacked sine waves.
1
Therefore sinusoidal contact model should be used instead of spherical contact models in
the multiscale model.
1.1 Organization of thesis
Chapter 2 of this thesis presents a thorough background on di erent types of contact
models ranging from spherical type contact to sinusoidal type contact.
Chapter 3 (also published in [13]) presents the modeling methodology employed in
formulating the Finite Element Model (FEM). This chapter is divided into three sections.
The  rst section of this chapter describes how the solid model is built for for the  nite
element analysis. The second section deals with the mesh convergence that was achieved.
The third section titled Simulation methodology deals with how the  nite element method
employed the displacement method for the analysis.
Chapter 4 (also published in [13]) is divided into four sections. The  rst describes
in detail the elasto-plastic FEM results obtained from the model. The remaining sections
present the various mathematical  ts that were generated in estimation of real contact area
and average contact pressure.
In Chapter 5 the multi-scale model for elasto-plastic contact of rough surfaces is pre-
sented. Here the mathematical  t is incorporated into the multi-scale model.
Chapter 6 concludes the thesis with a summary of the  nite element analysis results.
2
Chapter 2
Literature review
2.1 Introduction
Elasto-plastic sinusoidal contact has recently become more important with the devel-
opment of several multiscale rough surface contact models [14] [20]. This is because some of
these multiscale rough surface contact models use a Fourier series to consider the di erent
scales. Since a Fourier series consists of a series of sine waves with di ernt amplitudes and
periods, it is logical to use sinusoidal contact models to consider individual asperity contact.
Sinusoidal contact in some cases may be a more realistic depiction of surface contact than
spherical contact. Intutively, surface asperities should over all be shaped more like sine
waves than hemispheres. A model of elasto-plastic sinusoidal contact would be very useful
to these multiscale contact models. Gao et al. [21] modeled two-dimensional sinusoidal
contact using the  nite element method, and made note of some interesting results which
will also be considered in the current analysis. Gao et al. [21] reported a H=Sy = 5.8 for the
elasto-plastic case. The current work  nds trends similar to [21] resulting from the three
dimensional elastic case. It was interesting to note that the H=Sy for the 3D elasto-plastic
case was not restricted to 5.8 as reported by Gao et al. [21].
This chapter elucidates various rough surface contact models such as the statistical
contact model, the fractal contact model and the multiscale contact model. A justi cation
for inclusion of sinusoidal deformation model into the multiscale model foundation is also
provided.
3
Figure 2.1: Topographical contour plot of the sinusoidal surface geometry.
4
2.2 Spherical contact models
Although the following section may at  rst seem o subject, in reality spherical contact
is very useful for modeling sinusoidal contact, especially when the deformation is restricted
to the tips of the peaks (see Fig. 2.2). For this reason, Johnson et al. [1] provides two limiting
solutions to three dimensional sinusoidal contact, and the  rst being based upon the Hertz
elastic spherical contact solution. The Hertz solution provides closed-form expressions to
the deformations and stresses of two spheres in a purely elastic contact. The two spheres
may have di erent radii and di erent elastic properties. However, the closed-form solutions
render an equivalent case where a single elastic sphere, having an equivalent elastic modulus,
E0 , and an equivalent radius, R, is in contact with a rigid  at (see Eqs. 2.1 2.4 that follow).
5
Figure 2.2: Cross-section of sinusoidal type contact.
6
The interference, !, can be described as the distance the sphere is displaced normally
into the rigid  at. The Hertz solution assumes that the interference is small enough such
that the geometry does not change signi cantly. The solution also approximates the sphere
surface as a parabolic curve with an equivalent radius of curvature at its tip. The resulting
equations for contact radius and load from the Hertz solution are:
AE = ?R! (2.1)
FE = 43E0pR(w)32 (2.2)
where
1
E0 =
1 ?21
E1 +
1 ?22
E2 (2.3)
1
R0 =
1
R1 +
1
R2 (2.4)
and E1, ?1, R1, E2, ?2, R2, are the elastic properties and radii of sphere 1 and 2, respectively.
Again, since the sinusoidal surface and the tip of a sphere are very similar at initial
contact, the initiation of plastic deformation (de ned as the critical interference) can be
derived from Hertz contact theory, but can also be used for sinusoidal contact for cases of
relatively small loads and interferences. The critical values of interference, area, and load
7
are derived for the case of initial sinusoidal contact (when the contact geometry is still
similar to a parabola or sphere) in Appendix A.
Kogut and Etsion [2] performed a  nite element analysis of the case of an elastic-
perfectly plastic sphere in contact with a rigid  at. Their work gives a very detailed analysis
of the stress distribution in the contact region, and piecewise empirical expressions are
provided for the contact area, the contact force, and the average contact pressure. At
values !=!c < 1 the Hertz contact solution is assumed. Two sets of equations are provided
for the ranges 1 < !=!c < 6, and 6 < !=!c < 110. They describe the deformation only
up to !=!c = 110, at which point full plasticity is assumed. In their analysis, the value
of hardness, H, is set to be  xed at 2:8?Sy. Jackson and Green [22] also provide a more
complicated model which may describe large deformations better than the Kogut and Etsion
(KE) model. The most notable di erence being that H is often much less than 2:8?Sy for
high loads. However, Quicksall et al. [23] found that the Jackson and Green and KE
model agree fairly well for small deformations and the KE model is also much easier to
algebraically manipulate. For these reasons, the KE model is used in the current work to
model the initial elasto-plastic deformation of the sinusoidal surfaces.
The KE model is used in the current work to obtain an empirical relationship between
contact pressure and area. Since the current work is mostly concerned with spherical contact
for small load and interferences, the pressure-interference and the area-interference relation?s
for 1 < !=!c < 6 from Kogut and Etsion [2] are used:
8
: p
Sy
9
; = 1:19
8
: !
!c
9
;0:289 (2.5)
8
8
:AKE
Ac
9
; = 0:93
8
: !
!c
9
;1:136 (2.6)
With some manipulation Eq. 2.5 becomes
!
!c =
h 1
1:19
8
: p
Sy
9
;i
1
0:289 (2.7)
The average contact pressure, p, is related to the average pressure, ?p, over the sinusoidal
surface by:
p =
8
: ?p
2?AKE 
2
9
; (2.8)
The ?p denotes average pressure over the entire surface, including areas that are not in
contact while p denotes average contact pressure, just in the contact area. When complete
contact occurs, ?p = p . Substituting Eq. 2.7 and 2.8 into Eq. 2.6 the following is obtained:
AKE
Ac =
8
: ?p
2?AKE ?Sy ?1:19 
2
9
;3:93 (2.9)
Further solving for AKE results in the equation below:
AKE = (Ac) 14:93
8
: ?p
2:38?Sy 
2
9
;
3:93
4:93 (2.10)
9
An equation similar to Eq. 2.10 is used in the current work to model initial elasto-plastic
contact of the sinusoidal surfaces. This is explained in more detail later.
2.3 Elastic sinusoidal contact
The analysis for the case of 3-D waviness developed by authors Johnson, Greenwood
and Higginson [1] (hereafter referred to as "JGH") provides a relation between pressure and
contact area. First ?p is de ned as the average pressure on the surface (considering both
contacting and non-contacting regions) and p? is the amplitude of a sinusoidal pressure
superimposed on the mean pressure that yields complete contact. Complete contact is
de ned as when the entire sinusoidal surface is  attened and in contact with the opposing
surface (i.e. there is no gap between the surfaces). The amplitude of the pressure pro le,
p? is given as:
p? = p2?E0 f (2.11)
where E0 is the reduced elastic modulus,  is the amplitude of the sinusoidal surface, and
f is the frequency or reciprocal of wavelength,  . Thus when ?p  p?, the pressure loads the
surfaces so that there is no gap between them. Alternatively, when ?p < p? the contact is not
complete, and a closed form solution for the three-dimensional waviness contact problem is
not available. However, Johnson et al. [1] provides two asymptotic solutions to the problem.
For ?p ? p? the following equation derived from Hertz contact theory applies:
(AJGH)1 = 2?f2
h 3
8?
?p
p?
i2
3 (2.12)
10
and when ?p approaches p? (i.e., contact is nearly complete) the following equation applies:
(AJGH)2 = 1f2
8
>:1 3
2?
h
1 ?pp?
i9>
; (2.13)
Since no general analytical solution is available, an equation linking Eqs. 2.12 and 2.13 is  t
by Jackson and Streator [14] to the experimental and numerical data provided by Johnson
et al. [1]:
For ?pp? < 0:8:
A = (AJGH)1
8
>:1 h ?p
p?
i1:519>
;+(AJGH)2
8
>: ?p
p?
9
>;1:04 (2.14)
For ?pp?  0:8:
A = (AJGH)2 (2.15)
2.4 Elasto-plastic sinusoidal contact
The previous cases of elasto-plastic spherical contact and elastic sinusoidal contact
are now su?ciently expanded using  nite element results so that a model of elasto-plastic
sinusoidal contact can be formulated. Since it has been established that at low deformations
the sinusoidal contact will behave similarly to the spherical contact, this work will focus
mostly on the case of heavy deformations in sinusoidal contact.
Intuitively, the hardness of a sinusoidal shaped contact will follow a much di erent
trend than the spherical case. As shown in the Fig. 2.4, the slope of the surface at the
edge of contact for the sinusoidal case will be much di erent than the spherical case as the
11
amount of deformation is increased. While the spherical case approaches a rod or cone type
problem [22, 24], the sinusoidal case appears to reduce to a  at against a  at problem as
interference increases and the sinusoidal surface is  attened out. Thus, it would appear that
the hardness of the sinusoidal surface will initially stay constant or slightly decrease with
interference (depending on the value of  = ) and then eventually start increasing again as
 attening occurs.
For sinusoidal surfaces with large values of  = , the hardness to yield strength ratio
H=Sy will probably decrease initially because the geometry is very similar to that of a post
or rod. However, for smaller values of  = , H=Sy will probably never decrease because
the surface is more  at and the post or rod shape is never approached. Then, when the
contact is more  at, the e ective H=Sy becomes larger (as shown in the following sections).
It should be noted that the surfaces simulated in this work are for fairly low values of  = ,
and so examples of the  rst case are never seen.
When the sinusoidal surface is almost completely  attened, one may initially conclude
that the contact condition is similar to the spherical case at low a=R values and thus H=Sy
returns to approximately 2.84 or higher values (as was concluded by Gao et al.[21]). How-
ever, if this case is examined closely, it becomes unclear what the actual limiting pressure
is during the complete contact case. For example, once complete contact has occurred, the
case becomes very similar to that of a compressed cavity of material (see Fig. 2.4). In this
case the material cannot deform in the lateral direction parallel to the surface because there
are neighboring identical sinusoidal contacts. During elastic contact and using Hooke?s Law
12
Figure 2.3: Diagram of the hypothetical progression of change in hardness with geometry
for spherical and sinusoidal contact.
13
Figure 2.4: Compressed cavity of material.
14
the normal pressure is related to the lateral stresses by
?lateral = 1 3?1 ? P (2.16)
However, the current case is not elastic, but it is unclear what the material does when large
plastic deformations occur, or if it can reach the plastic regime once complete contact has
occurred. One way to model a plastically deforming material is to assume that it conserves
volume (? = 0:5) then ?lateral =  P. For the cavity case, plastic deformation then cannot
occur since the shape of the material is restricted, and so the only way for the material
to deform is for it to compress volumetrically (hydrostatic stress). Thus the volume of
material can only deform elastically and H=Sy approaches in nite values (see Fig. 2.4). In
the work by Gao and Bower [18] and Gao et al. [21] it also appears that H=Sy increases
dramatically when a loaded 2-D surface approaches the case of complete contact. However,
in the current 3-D model the authors?  nd that p? during fully plastic deformation is not
limited by a value of 5:8?Sy that appears to limit the 2-D case [21].
Greenwood and Rowe [25] cover the topic of plastic crushing of serrated surfaces and
also note that the force to cause complete contact becomes drastically larger as the surfaces
come closer together. The case of a conical indenter correlates to this case as well, since
as the angle of the cone tip approaches 180o the stress beneath the indenter will then
become hydrostatic (see Marsh [26] and Johnson [27]). This is also a similar concept to that
employed by the fractal contact models of Majumdar and Bhushan [28], in that they predict
as the contact area increases, the contacting asperities will e ectively become more elastic.
However, it should be noted that in reality the material in the current case did initially
undergo severe plastic deformation, so that even though the contact appears to behave
15
elastically near complete contact, it has undergone a large amount of plastic deformation
to arrive at that condition.
2.5 Statistical contact model
Greenwood and Williamson [29]  rst introduced the statistical approaches for surface
topography characterization in contact mechanics. This method has long been viewed as a
simple and easy to use rough surface contact model. However the simplicity of this model
has a few short comings such as the dependence of spectral moments on the resolution of
the surface measuring apparatus and the sample length as shown by [30] and [31]. This
dependence can skew the results for the contact parameters, such as the contact area and
contact load, which depend on the surface parameters.
Greenwood and Williamson in their work [29] showed that rough surfaces can be mod-
eled as a set of mutually exclusive asperities with constant radii and a variable height based
on a particular height distribution function. The parameters that characterize the surfaces
are the standard deviation of asperity heights ?s, the areal asperity density ?, and the asper-
ity radius R. These parameters are a function of the spectral moments which are extracted
from the surface using [32]. The value of ?s is related to ?, ?, R given by (Eq. 2.17).
?2s = ?2  3:717e
 4
?2 ?R2 (2.17)
16
AGaussiandistribution (z)isusuallyassumedfortheheightdistribution. Thecontact
parameters of rough surfaces are obtained as given by
A(d) = ?An
Z 1
d
?A(z  d) (z)dz (2.18)
P(d) = ?An
Z 1
d
?P(z  d) (z)dz (2.19)
where An is the nominal or apparent area of contact, which is de ned by the overlap of
surfaces in contact. The individual asperity contact area, ?A, and the corresponding contact
load ?P, are functions of !, which is the individual asperity interference. This original GW
model assumes elastically deforming hemispherical asperites de ned by Hertz solution [3].
Later models also expand the statistical methodology using elasto-plastic asperity contact
models [22] and [33].
2.6 Fractal contact model
In order to overcome the drawback of spectral moments dependence on surface param-
eters in the statistical model, the fractal contact model proves bene cial in characterization
of surface topography. The fractal geometry proposed by [30] has been utilized to char-
acterize the surface topography in contact mechanics. Although the fractal contact model
captures the surface topography in a multi-scale nature, it is important to note that not all
engineering surfaces have pro les which exhibit fractal behavior.
17
The feature of  ner and  ner detail becoming apparent as the measurement length scale
is reduced is a characteristic of the fractal model. Many authors [34], [35], have proposed
to replace the asperity model concept by a description of a surface as a fractal. The trouble
with the fractal model, apart from what as mentioned earlier is that there is no obvious way
to solve the contact problem for two fractal surfaces, or for one fractal surface contacting
a plane. The most widely accepted fractal theory is that due to Majumdar and Bhushan
[36] in which the distribution of contact areas is tied to the distribution of ?islands? cut o 
from the surface by a given horizontal plane as shown in Fig. 2.5. Therefore, the contact
area is simply a truncation of the fractal surface geometry and the plane. This is a classical
quantity in fractal theory. Majumdar and Bhushan [36] replace the material in each island
by a smooth parabolic asperity, this e ectively de nes a fractal distribution of asperities
that can then be treated as in a conventional asperity model.
Figure 2.5: Islands de ned by a horizontal section.
A typical fractal surface pro le can be generated using the Weierstrass-Mandelbrot
function [34]. A truncated two-variable fractal pro le can be described as [34]
18
z(x;y) = L(GL)d 2(ln?M )1=2
MX
m=1
nmaxX
n=0
?(cos m;n cos[2??
n(x2 +y2)1=2
L cos(arctan(y=x) 
?m
M )+ m;n])
(2.20)
where L is the sample length, D is the fractal dimension, G is the fractal roughness, ?
(=1.5) is a scaling factor, M is the number of superimposed ridges, n is a frequency index,
with nmax = R[log(L=Ls=log?)] represents the upper limit of n, where Ls is the cut-o 
length and  m;n is a random phase. Majumdar and Bhushan [36] have extensively used
fractal methods to characterize surface topography and also in scale independent rough
surface contact models. Kogut and Jackson [31] in the paper showed a series of results
comparing short falls of the statistical and fractal models for rough surface contacts. To
address these issues, several new multiscale models have been developed and are outlined
in Chapter 5.
19
Chapter 3
Finite element model
3.1 Introduction
This chapter describes in detail the methodology employed in building the solid  nite
element model. The solid model is divided into two parts. The  rst being the rigid  at plate
that acts as the contact surface and the second being the sinusoidal surface on which target
elements are created. Initially, the elastic case will be compared to the known solution to
con rm the model accuracy. once veri ed the model will be given elastic perfectly plastic
material properties.
3.2 Building the solid model
The current analysis will examine the case of three-dimensional elastic perfectly plastic
sinusoidal contact by building on these previous works, using fundamental solid mechanics
theory, and conducting a parametric study using the  nite elements method. The sinusoidal
surface considered by the current work is described by
h = 4?
?
1 cos
?2?? ?x
 
?
cos
?2?? ?y
 
??
(3.1)
and is shown (Fig. 2.1), where h is the height of the sinusoidal surface from its base. This
is very similar to the surface used by Johnson et al. [1] and results in the same analytical
equations for elastic contact.
20
A three-dimensional model was developed and the commercial ANSYSTM 8.1 package
was used to further analyze the elasto-plastic sinusoidal contact problem. Owing to the
symmetry of the sinusoidal surface (shown in Fig. 2.1 and 3.2), only a quarter section
of the whole problem is modeled. The sinusoidal pro le was generated based on equation
(Eq. 3.1). Keypoints were created for various values of x, y, and h. These keypoints
were then connected through lines. A conformal area was  t to generate the sinusoidal
surface. In all, over 42,000 elements where used in the analysis. Solid45, which is an 8-node
brick element, was used within the entire volume of the model. Contact174 and Target170
elements collectively formed the contact pair to model interaction between the surfaces.
3.3 Mesh convergence test
The number of elements was increased iteratively by a factor of 2 until mesh conver-
gence was obtained. The  nal value corresponded to concordant values over two successive
iterations. Upon  nal convergence the rigid  at surface, comprising of Contact174 elements,
was a uniform mesh formed by a 21x21 array (Fig. 3.1) of elements. In all, there were ap-
proximately 1200 contact elements in the model. The contact sti ness value was also altered
until the model agreed with the elastic solution (see section 3.5).
The uniform mesh on the rigid surface is used to predict the real contact area with the
sinusoidal surface. By determining the contact status of each node during post-processing
the total number of nodes in contact and the corresponding contact pressures are obtained
from the nodal solution for incremented values of displacement. The ratio of the number of
nodes in contact to the total number of nodes over the surface gives the real area of contact
normalized by the apparent or nominal area of contact.
21
Figure 3.1: Uniform mesh on the rigid  at formed by 21x21 array of elements.
22
3.4 Simulation methodology
The displacement method was employed to simulate the contact problem. This method
applies a  nite displacement to the rigid  at surface in the z-direction toward the sinusoidal
surface, and then solves the contact problem. This non-linear analysis accounts for large
displacements, hence additional measures were taken to enable the non-linear geometry
option during the course of the simulation. Displacement boundary conditions were enforced
on the surface areas of the solid model, and not in the contact area (see Fig. 3.2). The base
surface along the xy plane was  xed in all directions. The  at rigid surface was constrained
to move only along the z axis. Constraints were applied to surfaces located on the xy, xz
and yz planes in the direction perpendicular to the plane (thus enforcing a cyclic boundary
condition). For example, a surface area along the yz planes was constrained in the x-
direction, as shown in Fig. 3.2.
The simulation methodology can be understood better from Fig. 3.3. A contour of
the displacement in the z-direction is plotted on the deformed sinusoidal structure for the
elastic sinusoidal case. The material properties that were used to model this case is given in
Table 3.1. The surface on which the zero displacement boundary conditions were imposed
can be clearly seen during the course of simulation.
23
Figure 3.2: Schematic of degree of freedom restraints used for the one quarter sinusoidal
FEM model (actual modeled geometry is much longer in the z direction).
24
Figure 3.3: Various stages of contact.
25
The deformation history of the sinusoidal contact for the elasto-plastic benchmark
case is shown in Fig. 3.4. The material properties used for the benchmark case are E =
200GPa, Sy = 1GPa, and ? = 0:3. The dimensionless geometric ratio  = was set as
0.02. The contour plot (Fig. 3.4) represents the nodal displacements in the z-direction.
Interpenetrations that were observed during the course of simulation runtime were well
within accepted values of interference. Hence, these interpenetrations did not a ect the
 nite element analysis results.
Figure 3.4: Deformation history.
26
3.5 Veri cation of model accuracy
Model accuracy was achieved by comparing the elastic results of the  nite element
model (FEM) to that of Johnson et al.?s data [1] (see Fig. 3.5). From Fig. 3.5, FEM data
from the elastic sinusoidal case lie well within the two asymptotic solutions provided by
Johnson et al. [1]. It is also clear that the curves for both the data sets followed almost
identical paths. Table. 3.1 shows the parameters describing the JGH [1] data and the FEM
results from the elastic sinusoidal model [13]. An average error of 6% was found to exist
between the FEM and JGH data over the entire range of pressure.
Table 3.1: Material properties
Treatment JGH data FEM data
Young?s Modulus E 200 GPa 200 GPa
 = 0.02 0.02
Average pressure for complete contact (p?) 19.5 GPa 19.5 GPa
27
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
A.f
2
JGH data
FEM data
(A
JGH
)
1
(A
JGH
)
2
A(0? 
?
p/p
*
?1)
Figure 3.5: Comparison of elastic FEM results with JGH model.
28
Chapter 4
elasto-plastic fem results
4.1 Introduction
The material properties of the solid 3D model used to run the elastic case, were modi ed
for the elasto-plastic case. To enhance convergence, as outlined by Brizmer et al. [37], the
material of the sinusoidal surface was assumed elasto-plastic bi-linear isotropic hardening
with a yield stress, Sy, of 1 GPa and a tangent modulus, ET, 2% of the Young?s modulus E.
This isotropic hardening signi cantly improves the convergence without causing a signi cant
change in the results.
4.2 Parametric study
In order to formulate a  t for the FEM contact pressure a parametric analysis of the
elasto-plastic sinusoidal surface contact problem was conducted. A benchmark case was set
to analyze the contact problem. The material properties used for the benchmark case are
E = 200GPa, Sy = 1GPa, and ? = 0:3. The material properties along with the  = ratio
were then individually varied to perform a parametric study. First, a range of yield stresses
were considered in the model (see Figs. 4.1 and 4.2). The yield strength was varied from
Sy = 0:75GPa to Sy = 2:25GPa.
First the results are presented as a plot of the average pressure divided by the yield
strength, ?p=Sy, as a function of the normalized area, A?f2 (see Fig. 4.1). As the area in-
creases the contact is becoming more complete and the amount of deformation is increasing.
29
Table 4.1: Parametric study based on the Yield Strength Sy
Sy Sy  = E0 Sy=E0 p? p?ep
7.50?108 7.5?108 0.02 2.2?1011 3.41?10 3 1.95?1010 4.35?109
1.00?109 1.0?109 0.02 2.2?1011 4.55?10 3 1.95?1010 5.00?109
1.25?109 7.5?108 0.02 2.2?1011 5.69?10 3 1.95?1010 5.66?109
1.50?109 7.5?108 0.02 2.2?1011 6.82?10 3 1.95?1010 6.31?109
1.75?109 7.5?108 0.02 2.2?1011 7.96?10 3 1.95?1010 6.97?109
2.00?109 7.5?108 0.02 2.2?1011 9.10?10 3 1.95?1010 7.63?109
2.25?109 7.5?108 0.02 2.2?1011 1.02?10 2 1.95?1010 8.28?109
Traditionally this would be considered to be the fully plastic regime and the average pres-
sure, ?p, would be the hardness, H. As shown, the ratio of ?p=Sy increases past the typical
H=Sy value of 3 and even past the value of 5.8 found in [21] for 2-D sinusoidal surfaces in
contact. This agrees with the earlier theory that the contact becomes more elastic as the
contact area becomes more complete (A?f2 ! 1).
Next, the average contact pressure, ?p, resulting from the elasto-plastic model was nor-
malized using p?. Then the normalized contact area A?f2 was plotted versus ?p=p?. It can
be seen that the normalized contact area increases steadily as the yield strength decreases.
It can be seen from Fig. 4.2 that for a ?p=p? = 0:42, the normalized contact area ratio is 1
for all the cases (the contact is complete). Thus for the elasto-plastic cases, complete con-
tact occurs much earlier than when it occurs in elastic contact, as is intuitively expected.
Although the trend captured partially resembles the perfectly elastic behavior, in order to
 nd empirical expressions for the elasto-plastic case, the  ndings suggest that a new average
pressure to cause complete contact during elasto-plastic contact, p?ep, can be found.
30
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
A.f
2
p/S
y
S
y
=0.75 GPa
S
y
=1.00 GPa
S
y
=1.25 GPa
S
y
=1.50 GPa
S
y
=1.75 GPa
S
y
=2.00 GPa
S
y
=2.25 GPa
Figure 4.1: pSy versus contact area ratio (A?F2) for di erent yield strength.
31
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
A.f
2
?
p/p
*
S
y
=0.75 GPa
S
y
=1.00 GPa
S
y
=1.25 GPa
S
y
=1.50 GPa
S
y
=1.75 GPa
S
y
=2.00 GPa
S
y
=2.25 GPa
Figure 4.2: Contact area ratio (A?F2) versus ?p=p? for di erent yield strength.
32
Table 4.2: Parametric study based on the Young?s Modulus E0
E0 Sy  = E0 Sy=E0 p? p?ep
1.50?1011 1.00?109 0.02 1.65?1011 6.07?10 3 1.46?1010 4.41?109
2.00?1011 1.00?109 0.02 2.20?1011 4.55?10 3 1.95?1010 5.00?109
2.50?1011 1.00?109 0.02 2.75?1011 3.64?10 3 2.44?1010 5.60?109
3.00?1011 1.00?109 0.02 3.30?1011 3.03?10 3 2.93?1010 6.20?109
3.50?1011 1.00?109 0.02 3.85?1011 2.60?10 3 3.42?1010 6.79?109
Table 4.3: Parametric study based on the Geometric ratio  = 
 = Sy  = E0 Sy=E0 p? p?ep
0.010 1.00?109 0.010 2.20?1011 4.55?10 3 9.76?109 3.54?109
0.015 1.00?109 0.010 2.20?1011 4.55?10 3 1.35?1010 4.30?109
0.020 1.00?109 0.020 2.20?1011 4.55?10 3 1.95?1010 5.00?109
0.032 1.00?109 0.032 2.20?1011 4.55?10 3 3.12?1010 6.33?109
0.040 1.00?109 0.040 2.20?1011 4.55?10 3 3.91?1010 7.08?109
33
4.3 Critical elasto-plastic pressure estimation
The average pressure, p?ep, that causes complete contact is extracted from the  nite
element model data for each modeled case. This value corresponds to the average pressure
when the area ratio, A=f2 = 1. At this stage complete contact occurs between the rigid
 at and the sinusoidal surface. An equation is  t to the values of p?ep from the  nite
element data. The goal of the empirical  t is to obtain a single expression that takes
into consideration both the material and geometric properties of the modeled elasto-plastic
sinusoidal contact problem. Since the material properties and geometric properties are each
varied independently from the benchmark case, an equation can easily be  t to each trend.
It was also found that the e ects Sy and E0 appear to be almost exactly inverse and can be
combined into one normalized variable Sy=E0. The resulting equation  t to the FEM data
is given as
p?ep
p? =
8
>:4:172? Sy
E0 +0:0173
9
>;?
r
 
 (4.1)
The e ect of Sy=E0 on p?ep=p? is shown by the plot in Fig. 4.3. The relationship of
p?ep=p? to Sy=E0 appears to be linear in nature. As shown in Fig. 4.3, Eq. 4.1 di ers from
the FEM data by an average error of 3.82% with the maximum and minimum errors being
7.89% and 0.17% respectively. The empirical  t generated thus seems to be e ective.
34
0 0.002 0.004 0.006 0.008 0.01 0.012
0
0.05
0.1
0.15
0.20
0.25
0.30
0.35
0.40
0.45
S
y
/E
p
*
ep
/p
*
Fit Eq
FEM data
Figure 4.3: The e ect of Sy=E0 on pressure to cause complete contact.
35
Table 4.4: Variation of p?epp? (FEM) and p?epp? (Fit Eq. 4.1) with respect to Sy=E0
Sy=E0  = p?epp? (FEM) p?epp? (Fit Eq. 4.1) Error
3.41?10 3 0.02 0.2344 0.2226 4.9936
4.55?10 3 0.02 0.2429 0.2562 5.5666
5.69?10 3 0.02 0.2832 0.2897 2.3355
6.82?10 3 0.02 0.3305 0.3233 2.1570
7.96?10 3 0.02 0.3712 0.3568 3.8621
9.10?10 3 0.02 0.3784 0.3904 3.1821
1.02?10 2 0.02 0.4189 0.4240 1.2284
6.07?10 3 0.02 0.3268 0.3009 7.8912
4.55?10 3 0.02 0.2429 0.2562 5.4853
3.64?10 3 0.02 0.2382 0.2293 3.7062
3.03?10 3 0.02 0.2009 0.2114 5.2895
2.60?10 3 0.02 0.1983 0.1987 0.1769
Average Error 3.8228
36
0 0.005 0.015 0.025 0.035 0.045
0
0.05
0.15
0.25
0.35
0.45
0.55
?/?
p
*
ep
/p
*
Fit Eq
FEM data
Figure 4.4: The e ect of  = on pressure to cause complete contact.
37
Table 4.5: Variation of p?epp? (FEM) and p?epp? (Fit Eq. 4.1) with respect to  = 
Sy=E0  = p?epp? (FEM) p?epp? (Fit Eq. 4.1) Error
4.55?10 3 0.010 0.3615 0.3623 0.2340
4.55?10 3 0.015 0.2995 0.2958 1.2202
4.55?10 3 0.020 0.2427 0.2562 5.5666
4.55?10 3 0.032 0.2027 0.2025 2.3742
4.55?10 3 0.040 0.1892 0.1811 4.2418
Average Error 2.7274
The e ect  = has on p?ep=p? appears to be non-linear, as shown in Fig. 4.4. The p?ep=p?
value also decreases with an increase in the  = ratio. As apparent in Eq. 4.1, a square
root function appears to produce a reasonable  t. The average error resulting from the  t
is 2.72%, with maximum and minimum error values of 5.56% and 0.23% respectively.
Next the p?ep values obtained from Eq. 4.1 are used as a normalization factor for the
contact pressure, ?p, predicted by the elasto-plastic FEM model (see Figs. 4.5- 4.7). This
normalization is useful because it is successful at collapsing the curves onto almost the same
curve (see Fig. 4.5 in comparison to Fig. 4.2). Since similar plots are not shown for the
results when E and  = are varied, the values are varied by as much as 100.5% between
the E cases and 91.1% between the  = cases, which is large in comparison to the error
values given above. The results for the cases of various yield strengths, Sy, are shown in the
plot of A?f2 versus in Fig. 11. In this case the equivalent elastic modulus, E0 = 220GPa
and the geometric property  = = 0:02 are held constant for all the cases modeled. From
a lower Sy value of 0.75 GPa to a higher value of 2.25 GPa the values of curves seem to
increase slightly. As the yield strength, Sy, increases the curves appear to converge to a
single curve. Towards complete contact it can be seen that although the curves are for
di erent yield strengths, they converge to the same point due to the normalization by p?ep.
38
Since the trend of the normalized contact area A?f2 as a function of ?p=p?ep should always
begin at A?f2 = 0, ?p=p?ep = 0 (just before initial contact) and end at A?f2 = 1, ?p=p?ep = 1
(when complete contact occurs), the normalization is successful just at collapsing these end
points together. However, the curvature of the normalized trend can still depend on the
material and geometric properties. Therefore, additional measures are taken to  t equations
to the FEM data.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
ep
A.f
2
S
y
=0.75 GPa
S
y
=1.00 GPa
S
y
=1.25 GPa
S
y
=1.50 GPa
S
y
=1.75 GPa
S
y
=2.00 GPa
S
y
=2.25 GPa
Figure 4.5: Contact area ratio (A?f2) versus ?p=p?ep for di erent yield strength values.
From Fig. 4.6 it is seen that the elastic modulus E0 and the yield strength Sy both
display similar trends. This case is modeled by varying only the equivalent elastic modulus,
E0, while the yield strength, Sy =1 GPa and  = =0.02 are held constant for all the di erent
39
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
ep
A.f
2
E=350 GPa
E=300 GPa
E=250 GPa
E=200 GPa
E=150 GPa
Figure 4.6: Contact area ratio (A?f2) versus ?p=p?ep for di erent values of Young?s Modulus
E.
40
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
ep
A.f
2
?/?=0.01
?/?=0.02
?/?=0.032
?/?=0.04
Figure 4.7: Contact area ratio (A?f2) versus ?p=p?ep for di erent values of the dimensionless
geometric quantity  = .
41
cases. As the elastic modulus increases, the curve also increases. At lower values of ?p=p?ep
the curves seem to fall onto a single curve. When the ?p=p?ep value reaches 0.1 the curves
seem to diverge. Near complete contact the curves converge again.
The geometric property  = is now varied from 0.01 to 0.04, as shown in Fig. 4.7.
The equivalent modulus E0=220 GPa and the yield strength, Sy=1 GPa, are held constant
in these various cases that are modeled. The contact pressure, ?p , from the FEM data
is again normalized by the average pressure for complete elasto-plastic contact, p?ep. This
ratio is then plotted versus the contact area ratio. It can be seen that normalized area
(Af2) steadily increases when the  = ratio increases. Another point should be noted,
that the lower values of the ratio  = allow for the solution to converge faster, while some
higher values of  = do not converge at all. The reason for this dependence on  = for
convergence could be from the fact that with higher  = ratios larger deformations are
required to  atten out the sinusoidal surface.
4.4 Real contact area estimation
The area described by the KE model is next modi ed to better suit the current elasto-
plastic sinusoidal contact problem. The premise for the argument is based on the fact that
when A=Ac=1, the pressure to critical pressure ratio p=pc should equal 1. Notice that p is
now used to denote the average contact pressure over the contact area, while ?p denotes the
average pressure over the entire surface area. Based on the KE model, a logical function to
model initial elasto-plastic contact is described below
A
Ac =
8
: p
pc
9
;d (4.2)
42
where pc from Eq. (A5) in Appendix 1 is
pc =
8
:2?C ?Sy
3
9
;
substituting this into Eq. 4.2 we get,
A
Ac =
8
: p
2?C ?Sy ?3
9
;d (4.3)
where average pressure over the contact area, p, in Eq. 4.3 is converted to the average
pressure over the entire surface, ?p , by
p =
8
:?p? 2
2?A
9
; (4.4)
By substituting in Eq. 4.4 into Eq. 4.3 the following expression results
A
Ac =
8
>: ?p? 2 ?3
4?Sy ?C ?A
9
>;d (4.5)
Solving for A in Eq. 4.5, one obtains the modi ed expression for contact area, which
is now termed as Aep. Also noting that the two sphere-like peaks occur in the  2 area of
sinusoidal contact than a factor of two must be included. The result is given by
Aep = 2(Ac) 11+d
8
>: 3? ?p
4?C ?Sy 
2
9
>;
d
1+d (4.6)
Although this expression (Eq. 4.6) for area is derived di erently, it is interesting to note
that it is a similar result to the area to pressure relation from the KE model (Eq. 2.10).
43
The value of d is obtained empirically by  tting a curve to the FEM results. To allow a
good  t with the FEM results, the value of d found in Eq. 4.6 for Aep is allowed to vary
with E0=Sy and  = . The expression for d is given by
d = C1 ?
8
>:E0
Sy ?
 
 
9
>;C2 (4.7)
where, C1 = 3:8 and C2 = 0:11 and are constants which are obtained empirically. E0 and
Sy are the e ective Young?s modulus and yield strength respectively. For the benchmark
case of E0=220 GPa, Sy=1 GPa, and  = =0.02, then d=4.47. As expected, this is close to
the value of 4.93 that is predicted by the KE model.
To model the contact area as a function of load, the results of Jackson and Streator
[14] are now modi ed with the elasto plastic results now given by Eqs. ( 4.1and 4.6). The
modi ed version of the model to consider elasto-plastic deformation is
A = (Aep)
8
>:1 h ?p
p?ep
i1:519>
;+(AJGH)2
8
>: ?p
p?ep
9
>;1:04 (4.8)
where Aep is the spherical elasto-plastic contact area as predicted by the model given by
Eqs. ( 4.6 and 4.7) and based on the KE model. Then as the load increases, the contact
will diverge from the spherical case and asymptotically approach the sinusoidal case. Since
the H=Sy will theoretically increase as the contact becomes more complete (see previous
discussion concerning Fig. 2.4), the contact area may theoretically approach the elastic
limit (AJGH)2 as ?p=p?ep approaches 1. However, when the solution approaches (AJGH)2 the
average pressure is p?ep instead of p? as in the elastic case. Also, if A, as predicted by Eq. 4.8,
increases past  2 then A= 2.
44
4.5 The empirical model
The FEM results are presented alongside the model given by Eq. 4.8 in Figs. 4.8, 4.9,
and 4.10. In Fig. 4.8 we see the  t between the two extreme cases that were modeled based
on the yield strength of the material. The lowest Sy value considered being 0.75 GPa and
the higher value being 2.25 GPa. We can see clearly that the modi ed KE model for contact
area agrees reasonably well with the numerical results. The average error between the new
model given by Eq. 4.8 and the FEM results when the yield strength is varied independently
is less than 5%.
45
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
ep
A.f
2
Model S
y
=0.75 GPa,E=200 GPa
FEM S
y
=2.25 GPa,E=200 GPa
Model S
y
=2.25 GPa,E=200 GPa
FEM S
y
=2.25 GPa,E=200 GPa
Figure 4.8: Comparison of the FEM results and the elasto-plastic sinusoidal contact model
(marked by lines) given by Eq. 4.8 for di erent yield strength values.
46
Likewise, the plots in Fig. 4.9 and 4.10 show that the model also compares well for the
cases where the elastic modulus, E, and the geometry,  = , are varied (the error in these
cases is also less than 5%). Fig. 4.9 a comparison of the FEM results and the elasto-plastic
sinusoidal contact model (marked by lines) given by Eq. 4.8 for di erent Young?s modulus
values is shown. Similarliy in Fig. 4.10 a comparison of the FEM results and the elasto-
plastic sinusoidal contact model (marked by lines) given by Eq. 4.8 for di erent sinusoidal
geometries is shown.
47
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
ep
A.f
2
Model E=150 GPa,S
y
=1.00 GPa
FEM E=150 GPa,S
y
=1.00 GPa
Model E=350 GPa,S
y
=1.00 GPa
FEM E=350 GPa,S
y
=1.00 GPa
Figure 4.9: Comparison of the FEM results and the elasto-plastic sinusoidal contact model
(marked by lines) given by Eq. 4.8 for di erent Young?s modulus values.
48
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
?
p/p
*
ep
A.f
2
Model ?/?=0.01,S
y
=1.00 GPa
FEM ?/?=0.01,S
y
=1.00 GPa
Model ?/?=0.04,S
y
=1.00 GPa
FEM ?/?=0.04,S
y
=1.00 GPa
Figure 4.10: Comparison of the FEM results and the elasto-plastic sinusoidal contact model
(marked by lines) given by Eq. 4.8 for di erent sinusoidal geometries.
49
Chapter 5
Multiscale Rough Surface Contact Model
5.1 Introduction
In this chapter a non-statistical multiscale model using the elastic and elasto-plastic
sinusoidal contact, based on the work by Jackson and Streator [14], is presented. The equa-
tions generated emperically using FEM in earlier chapter?s are included in the formulation
of the multiscale model. This model incorporates the e ect of asperity deformations at
multiple scales into a simple foundation for modeling rough surfaces. This model considers
the e ect of having smaller asperities located on top of larger asperities in repeated fashion,
very similar to the "protuberance upon protuberance" theory proposed by Archard [40].
5.2 Multiscale model foundation
The current multiscale model uses the same direction as Archard [40], but di ers in
the ease of use when applied to real rough surfaces. The model assumptions derived from
Jackson and Streator [14], which are di erent from the statistical and fractal models, are
as follows:
1. Asperities are stacked so that smaller asperities (higher frequency) are on top of larger
asperities (lower frequency).
2. Each iteration or frequency level of asperities carries the same total load.
3. The load is shared equally among all the asperities at a given frequency level.
50
4. At a given frequency level, each asperity deforms according to Hertz theory or to
a chosen elasto-plastic asperity contact model, irrespective of the presence of higher
frequency asperities upon it. In the current work the sinusoidal contact models will
be used.
5. At a given frequency level the contact area cannot increase beyond what is experienced
by the frequency level below it.
These assumptions set up the following foundation of equations for the contact model:
Ar =
?imaxY
i=1
?Ai?i
?
An (5.1)
F = ?Fi?iAi 1 (5.2)
where Ar is the real area of contact, F the contact load, An the nominal contact area,
and the subscript i denotes a frequency level, with imax denoting the highest frequency
level considered. A  ow chart of the method is given by Fig. 5.1. This illustrates how
the contact model foundation is used to model the contact between rough surfaces. First
FFT of the acquired surface data is performed. If a spherical contact model is used to
consider individual asperity, the areal density and radius of curvature are computed for
each frequency level from the resulting Fourier series. Although in the current work these
quantities are not used, they can be computed as follows:
?i = 2f2i (5.3)
51
Ri = 14?2?
if2i
(5.4)
where ?i is the amplitude at the given frequency level and fi is the frequency. The source of
factor two is from [1]. This factor only means that there are two density of asperity peaks per
reference area. However, the abpve methodology di ers when a sinusoidal contact model
is used. because the amplitude and wavelength can be directly used from the FFT. By
using the sinusoidal contact model assumption 5 will also never be violated. The nominal
contact area An is de ned at the zeroth frequency level or at i = 0. For a chosen level of
frequencies, the number of asperities is calculated for every frequency level. At a given level
the total load is then equally divided among all the asperities of that level. From the given
asperity load, dimension, the material property and the chosen deformation model, the
single asperity area is determined. A temporary total contact area for the given frequency
level is computed by multiplying the number of asperities at that level. Assumption 5
mentioned above is then veri ed at every frequency level. If the contact area predicted at
that level is greater than the contact area at frequency level below it, the smaller contact
area serves as the nominal contact area for the given frequency level. This procedure is
repeated iteratively until all the frequency levels are considered, resulting in a prediction of
the real contact area as a function of applied load.
52
Figure 5.1: Flow chart of iterative asperity contact model.
53
5.3 Elasto-plastic multiscale model
The sinusoidal asperity model developed earlier is used within the multiscale model
foundation described above to relate the contact area to the applied load. The modi ed
KE contact area is used to de ne the asperity deformation model. Aep is given by
Aep = 2(Ac) 11+d
8
>: 3? ?p
4?C ?Sy 
2
9
>;
d
1+d (5.5)
To model the contact area as a function of load inside the multi-scale model founda-
tion, the JGH contact area is modi ed with the elasto plastic results now given by Eqs.
( 4.1and 5.5). The modi ed version of the model to consider elasto-plastic deformation is
A = (Aep)
8
>:1 h ?p
p?ep
i1:519>
;+(AJGH)2
8
>: ?p
p?ep
9
>;1:04 (5.6)
where Aep is the spherical elasto-plastic contact area as predicted by the model given by
Eqs. ( 5.5 and 4.7) and based on the KE model. Then as the load increases, the contact
will diverge from the spherical case and asymptotically approach the sinusoidal case. A
limiting condition for elastic sinusoidal case was provided by Johnson et al. [1].
5.4 Frequency spectrum for di erent surfaces
Three surfaces surface 1, surface 2 and surface 3 were selected for the multiscale model.
A stylus pro lometer was used to measure the pro les of arbitrarily machined metal samples
is shown in Figs. ( 5.2, 5.3, and 5.4). The displayed pro le each comprised of 496 to 3900
data points. A frequency spectrum of the surface pro le is also shown in Figs.( 5.2, 5.3,
and 5.4). The frequency spectrum (Fig. 5.2) of the surface 1 shows asperities of little
54
variation of amplitude with scale. The frequency spectrum plots for surface 2 and surface
3 show decreasing amplitudes with higher frequencies, as would normally be expected.
Since the motivation behind this work was to develop a multiscale model incorporating
the sinusoidal deformation model, the same material (Aluminum) with di erent levels of
roughness was considered. The typical material properties that were de ned were Young?s
modulus E = 200GPa, yield strength Sy = 0:1GPa, and Poisson?s ratio ? = 0:33. The
sample length was 400?.
10
0
10
1
10
2
10
3
10
?9
10
?8
10
?7
10
?6
i
Log Amplitude (m)
Figure 5.2: Frequency spectrum of surface 1.
55
10
0
10
1
10
2
10
3
10
4
10
?12
10
?11
10
?10
10
?9
10
?8
10
?7
10
?6
i
Log Amplitude (m)
Figure 5.3: Frequency spectrum of surface 2.
56
10
0
10
1
10
2
10
3
10
4
10
?14
10
?13
10
?12
10
?11
10
?10
10
?9
10
?8
10
?7
10
?6
i
Log Amplitude (m)
Figure 5.4: Frequency spectrum of surface 3.
57
5.5 Frequency level iteration versus contact area
To further emphasize how the multi-scale model works, the contact area as predicted by
the model as a function of frequency level iteration, i, for two load cases FAnE0 = 9:774x10 5
and FAnE0 = 9:774x10 9 is shown in Figs. ( 5.5, 5.6, and 5.7). To see the trends from
the elastic and the elasto-plastic multi-scale sinusoidal contact model, results from both the
models are plotted alongside in the same  gure. From Fig. ( 5.5) for the random surface
sample, we can see that at higher loads the contact area for elasto-plastic case is close to 1
for the frequency range 10-100. The elastic case varies from .1 to 1 for the same load and
frequency range. The drops in area occur due to the assumption 5 stated in the multiscale
model foundation. For any given load level the contact area predicted by the immediate
higher level cannot be more than the contact area predicted by the previous level. Based
o of a visual inspection it is clear that both lower and higher frequency level ranges dictate
the contact area prediction. This is a stark di erence from the Jackson and Streator [14]
model which considered surface 1 for the multiscale model. The contact area predicted by
the Jackson and Streator model [14] is dictated by the lower frequency level ranges and
the higher frequency level ranges have little to no e ect on the contact area predicted. A
possible reason for this di erence could be the use of spherical asperity deformation model
in their multiscale frame work. The surface pro les surfaces 2 and 3 that were considered
in the multiscale model exhibit similar trends to the Jackson and Streator [14] multiscale
model. This similarity is shown in Figs. ( ?? and ??) for these surface pro les.
58
10
0
10
1
10
2
10
3
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
Frequency iteration,i
A
r
/A
n
Elastic Load case F/(E*A
n
)=100
Elastic?plastic Load case F/(E*A
n
)=100
Elastic Load case F/(E*A
n
)=10
Elastic?plastic Load case F/(E*A
n
)=10
Figure 5.5: Contact area ratio versus frequency iteration for surface 1.
59
10
0
10
1
10
2
10
3
10
4
10
?6
10
?5
10
?4
10
?3
10
?2
Frequency iteration,i
A
r
/A
n
Elastic Load case F/(E*A
n
)=100
Elastic?plastic Load case F/(E*A
n
)=100
Elastic Load case F/(E*A
n
)=10
Elastic?plastic Load case F/(E*A
n
)=10
Figure 5.6: Contact area ratio versus frequency iteration for surface 2.
60
10
0
10
1
10
2
10
3
10
4
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
Frequency iteration,i
A
r
/A
n
Elastic Load case F/(E*A
n
)=100
Elastic?plastic Load case F/(E*A
n
)=100
Elastic Load case F/(E*A
n
)=10
Elastic?plastic Load case F/(E*A
n
)=10
Figure 5.7: Contact area ratio versus frequency iteration for surface 3.
61
5.6 Elasto-plastic and elastic comparison
Acomparisonofelasto-plasticandelasticroughsurfacecontactareasshowsanexpected
trend. In Figs. ( 5.8, 5.9, and 5.9) we can see that model predicts a larger non-dimensional
contact area as a function of non-dimensional surface load for the elasto-plastic case. This is
in agreement with the fact that large contact areas are produced when a surface is deforming
elasto-plastically. We can also see that the predicted real area of contact for both elastic and
elasto- plastic as function of applied is linear in nature. The Greenwood et al. [29] statistical
model and the Archard?s [40] protuberence theory also show similar trends. Furthermore
Amonton?s law of friction of can also be explained from the linear relationship between
contact area and applied load. The real contact area predicted by the multiscale model
is a function of the applied load. This relationship can be mathematically equated by
introducing a non-dimensional contact factor.
Ar = Cs ?Fy (5.7)
The friction force F? can be represented as a product of real area of contact and the Shear
strength ?y.
F? = Ar ??y (5.8)
Now, susbstituing for Ar from Eq. 5.7 into Eq. 5.8 yields the following:
F? = Cs ??y (5.9)
62
Eq. 5.9 can be viewed as the Amonton?s law of friction for a body experiencing normal
reaction on a rough surface. In other words, the friction force, F?, appears to be related
to the normal force, Fy, by a constant friction coe?cient. This is therefore in accordance
with one of the empirical ?law of friction?. This model can also be used to predict electrical
contact resistance and thermal contact resistance, adhesion and wear.
63
10
?9
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
F/EA
n
A
r
/A
n
Elastic
Elasto?plastic
Figure 5.8: Elastic and Elasto-Plastic contact area versus applied load for surface 1.
64
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
F/EA
n
A
r
/A
n
Elastic
Elasto?plastic
Figure 5.9: Elastic and Elasto-Plastic contact area versus applied load for surface 2.
65
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?8
10
?7
10
?6
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
F/EA
n
A
r
/A
n
Elastic
Elasto?plastic
Figure 5.10: Elastic and Elasto-Plastic contact area versus applied load for surface 3.
66
Chapter 6
Conclusions
The contact problem comprised of an elasto-plastic deformable sinusoidal surface and a
rigid  at was analyzed and modeled using the  nite element method. The numerical analysis
of the current  nite element model provided an approximate solution for the elasto-plastic
regime of sinusoidal contact. An empirical equation based on previous sinusoidal elastic
contact models and elasto-plastic spherical contact models was formulated to  t the FEM
data. Dimensionless expressions were empirically derived for the mean contact pressure
which causes complete contact during large elasto-plastic deformations. By normalizing
the contact area and the mean contact pressure, the current model provides analytical
expressions of contact area as a function of pressure for the elasto-plastic regime. The
expression for the contact area for the sinusoidal contact problem was obtained by modifying
the Kogut and Etsion model [2] and the elastic sinusoidal models [1]. The errors in the
numerical  ts are fairly low (less than 5%) and suggest that the  ts are reasonable. These
equations should prove useful for modeling contact between rough periodically structured
surfaces and in multiscale rough surface contact models.
The multiscale contact model utilizing the sinusoidal contact deformation model was
developed for the elastic and the elasto-plastic cases. The contact area as a function of
applied load was predicted using the model. The model predicted a linear trend as expected
both for the elastic and the elasto-plastic sinusoidal asperity deformation models. The
predicted area of contact was compared between the elastic and the elasto-plastic cases. As
expected the elasto-plastic case predicted more area of contact. The frequency distribution
67
of surface pro les in a multiscale model is a Fourier series which is a series of sine and
cosine wave forms. In keeping with inherent property of the model, the sinusoidal contact
model seems to be a good choice for a asperity deformation model. The multi-scale model
by Jackson and Streator [14] uses a spherical asperity deformation model for both elastic
and elasto-plastic cases, which is di erent from the current multi-scale contact model. The
most signi cant weakness of the multiscale model is the assumption that all asperites at a
given frequency level exhibit deformation characteristics under indentical load conditions.
This means that a asperities of given wavelength may not reside at di erent heights.
68
Bibliography
[1] Johnson,K.L., Greenwood,J.A, and Higginson, J.G., The contact of elastic regular
wavy surfaces., Int.J.Mech.Sci., 1985. 27(6): p. 383-396.
[2] Kogut,L., and Etsion,I., Elastic-plastic contact analysis of a sphere and a rigid
 at., J. of Applied Mechanics, Trans. ASME,2002. 69(5): p. 657-662.
[3] Johnson,K.L., Contact mechanics., Cambridge University Press., 1985.
[4] Dyson, A., Approximate calculations of Hertzian compressive stresses and contact
dimensions, J. of Mech. Engg. Sci., 1965. 98(7): p. 214-224.
[5] Galin, L.A., Contact problems in the theory of elasticity, Moscow., 1963.
[6] Garg, V.K., Anand, S.C., and Hodge, P.G., Elastic-plastic analysis of a wheel rolling
on a rigid track, Int.J.Solids and Structures, 1974. 10(945): p. 290-295.
[7] Goodman, L.E., Contact stress analysis of normally loaded rough surfaces, J. of
Applied Mechanics, Trans. ASME,1962. 29(515): p. 119-121.
[8] Goodman, L.E., and Keer, L.M., Elastic-plastic analysis of a wheel rolling on a
rigid track, Int.J.Solids and Structures, 1965. 1(407): p. 110-117.
[9] Johnson, K.L., Surface interaction between elastically loaded bodies under tan-
gential forces, Proceedings, Royal Society, 1955. A320(531): p. 75-220.
[10] Greenwood, J.A. and Tripp, J.H The elastic contact of rough spheres, J. of Applied
Mechanics, Trans. ASME,1967. 34(300): p. 417-420.
[11] Westergaard,H.M., Bearing pressure and cracks., ASME J. of Appl. Mech., 1939. 6:
p. 49-53.
[12] Stanley,H.M. and Kato,T., FFT-based method for rough surface contact., Journal
of Tribology, Transactions of the ASME, 1997. 119(3): p. 481-485.
[13] Krithivasan.V., and Jackson,R.L., An analysis of three-dimensional elasto-plastic
sinusoidal contact., Trib.Letters., 2007. 27(1): p. 31-43.
[14] Jackson,R.L. and Streator,J.L., A multiscale model for contact between rough sur-
faces., Wear, 2006. 261(11-12): p. 1337-1347.
69
[15] Bora,C.K, Flater,E.E, Street,M.D., Redmond,J.M., Starr,M.J., Carpick,R.W., and Ple-
sha,M.E., Multiscale roughness and modeling of MEMS interfaces., Tribology Let-
ters, 2005. 19(1): p. 37-48.
[16] Ciavarella, M., Demelio, G., Barber, J.R., and Jang, Y.H., Linear elastic contact of
the weierstrass pro le., Proc. R. Soc. Lond. A, 2000(456): p. 387-405.
[17] Ciavarella, M., Murolo, G., Demelio, G., and Barber, J.R., Elastic contact sti ness
and contact resistance for the weierstrass pro le., J. Mech. Phys. Solids, 2004.
52(6): p. 1247-1265.
[18] Gao, Y.F. and Bower, A.F., Elastic-plastic contact of a rough surface with weier-
strass pro le., Proc. R. Soc. A, 2006. 462: p. 319-348.
[19] Persson, B.N.J., Elastoplastic contact between randomly rough surfaces., Physical
Review Letters, 2001. 87(11): p. 116101.
[20] Persson, B.N.J., Bucher, F., and Chiaia, B., Elastic contact between randomly rough
surfaces: Comparison of theory with numerical results., Physical Review B, 2002.
65: p. 184106-1.
[21] Gao, Y.F., Bower, A.F., Kim, K.S., Lev, L., and Cheng, Y.T., The behavior of an
elastic-perfectly plastic sinusoidal surface under contact loading., Wear, 2006.
261(2): p. 145-154.
[22] Jackson, R.L. and Green, I., A  nite element study of elasto-plastic hemispherical
contact., ASME J. Tribol., 2005. 127(2): p. 343-354.
[23] Quicksall, J.J., Jackson, R.L., and Green, I., Elasto-plastic hemispherical contact
models for various mechanical properties., IMechE J. of Eng. Trib. -Part J., 2004.
218(4): p. 313-322.
[24] Chaudhri, M.M., Hutchings, I. M., Makin, P. L., Plastic compression of spheres.,
Philosophical Magazine, 1984. 49(4): p. 493-503.
[25] Greenwood, J.A. and Rowe, G.W., Deformation of surface asperities during bulk
plastic  ow., J. Appl. Phys., 1965. 36: p. 667-668.
[26] Marsh, D.M., Plastic  ow in glass., Proc. R. Soc. A, 1964.279(6): p. 420-435.
[27] Johnson, K.L., The correlation of indentation experiments., J. Mech. Phys. Solids,
1970. 18: p. 115-126.
[28] Majumdar, A. and Bhushan, B., Fractal model of elastic-plastic contact between
rough surfaces., ASME J. of Tribol., 1991. 113(1): p. 1-11.
70
[29] Greenwood, J.A. and Williamson, J.B.P, Contact of nominally  at surfaces., Proc.
R. Soc. Lond., 1966. A(295): p. 300-319.
[30] Majumdar, A., and Tien, C. L., Fractal characterization and simulation of rough
surfaces, Wear, 1990. 136: p. 313-327.
[31] Kogut, L., and Jackson, R. L., A comparison of contact modeling utilizing statis-
tical and fractal approaches, ASME J. Tribol., 2006. 128: p. 213-217.
[32] McCool, J. I., Comparison of models for the contact of rough surfaces, Wear, 1986.
107: p. 37-60.
[33] Chang, W.R., Etsion, I., and Bogy, D. B., An elastic-plastic model for the contact
of rough surfaces, ASME J. Tribol., 1987. 109(2): p. 257-263.
[34] Yan, W., and Komvopoulus, K., Contact analysis of elastic-plastic fractal surfaces,
J. Appl. Phys., 1988. 84: p. 3617-3624.
[35] Polycarpou, A., and Etsion, I., Analytical approximations in modeling contact
rough surfaces, ASME J. Tribol., 1999. 121(2): p. 234-239.
[36] Majumdar, A. and Bhushan, B., Role of fractal geometry in roughness characteri-
zation and contact mechanics., ASME J. of Tribol., 1990. 112(2): p. 205-216.
[37] Brizmer, V., Zait, Y., Kligerman, Y., and Etsion, I., The e ect of contact conditions
and material properties on Elastic-Plastic Spherical Contact., J. Mech. Mater.
Struct., 2006. 1(5): p. 865-879.
[38] Kogut, L., and Etsion, I., A  nite element based elastic-plastic model for the con-
tact of rough surfaces, Tribol. Trans., 2003. 46: p. 383-390.
[39] Eid, H. and Adams, G. G., An elastic-plastic  nite element analysis of interacting
asperities in contact with a rigid  at, J. Appl. Phys., 2007. 40: p. 7432-7439.
[40] Archard, J.F., Elastic deformation and the Laws of friction., Proc. R. Soc. Lond.,
1957. A(243): p. 190-205.
71
Appendices
72
Appendix A
Derivation of critical values for sinusoidal contact
The critical or initial interference to cause yielding can be derived independently of
hardness as
!c =
8
:? ?C ?Sy
2?E0
2
R
9
; (A.1)
where C is given by
C = 1:295exp(0:736??) (A.2)
and the Poisson?s ratio of the material that has the lowest value of C ?Sy, and thus yields
 rst should be used. The critical force, Fc is then calculated at the critical interference, !c,
to be
Fc = 43
8
:R
E0
9
;28:C
2 ?? ?Sy
9
;3 (A.3)
Similarly, the critical contact area is
Ac = ?3
8
:CSyR
2E0
9
;2 (A.4)
73
By dividing the critical force (Eq. (A3)) by the critical area (Eq. (A4)) the critical contact
pressure is de ned as
pc =
8
:2CSy
3
9
; (A.5)
These critical values can be modi ed for use in the sinusoidal case considered in this work
by substituting in the relation
R =
8
: 1
4?2 f2
9
; (A.6)
which is the radius of curvature for the tip of a sinusoidal surface. In addition, there are
actually two sphere-like peaks of contact which occur during a  2 area of sinusoidal contact
(see Fig. 1). This results in a modi ed critical force given by
Fc = 16?
8
: 1
 f2E0
9
;28:C ?Sy
2
9
;3 (A.7)
Similarly, the modi ed critical contact area is
Ac = 2?
8
: C ?Sy
8 f2E0
9
;3 (A.8)
Since R does not appear in Eq. (A5), the critical pressure is not in uenced by fi. It
should be noted that Eqs. (A5, A7 and A8) are only valid for the initial stages of con-
tact between sinusoidal surfaces (when the contacting surfaces are geometrically similar to
74
spheres). However, these critical values are useful in determining when a sinusoidal contact
begins deforming in the elasto-plastic range.
75
Appendix B
matlab code for contact area versus frequency iteration
hx = xlsread(0sdata0;1;0c2 : c39010);
hx = hx?10 10;
v = length(hx);
/********Length of sample********/
L = 400e 6;
/********Fourier transform of surface pro le********/
v = v=2;
Y = fft(hx;v)=v;
v = v=2;
YY = abs(Y);
YY = sqrt(Y:?conj(Y));
/********Material properties********/
E = 1=(1 :332)?200?109;
Sy = 1e8;
v1 = 0:33;
/********Johnson area coe?ents********/
w2 = 3=2=pi;
w3 = sqrt(2)?pi?E;
w4 = 3=8=pi;
k = 1;
/********Initial area and Intial load********/
Ari = 3:3?1:5?10(  6);
Fa = 0:0001?E ?Ari;
/********Highest level of applied load********/
76
for(j = 100)
Lo = Fa?(10((j  1)=100 5));
Lo1(j) = Lo=E=Ari;
Atotal(1) = Ari;
Atotal(2) = Atotal(1);
Atotalep(1) = Ari;
Atotalep(2) = Atotalep(1);
/********Frequency level iteration********/
for(i = 2 : v)
I(i 1) = i;
/********Amplitude and Frequency calculation********/
beta = YY(i);
f = (i 1)=L;
l = 1=f;
h = 2?f2;
C3 = Lo=((Atotal(i 1)));
pbar = (w3?beta?f);
prat = C3=pbar;
/********Elastic multiscale model********/
if(C3 > pbar)
Af(i) = Atotal(i 1)=h=Atotal(i 1);
else
AJ1 = (2?pi=f2)?(w4?prat)(2=3);
AJ2 = (1=f2)?(1 (w2)?(1 prat));
/********Area comparison for every frequency level********/
if(prat < 0:8)
Af(i) = AJ1?(1 (prat)(1:51))+AJ2?(prat)(1:04);
else
Af(i) = AJ2;
77
end
end
Atotal(i) = Atotal(i 1)?h?Af(i);
if(Atotal(i 1) > Atotal(i))
a(j) = Atotal(i);
be(i 1) = a(j)=Ari;
else
Atotal(i) = Atotal(i 1);
a(j) = Atotal(i);
be(i 1) = a(j)=Ari;
end
/********Elasto-plastic multiscale model********/
C4 = Lo=Atotalep(i 1);
Fasp = C4=f2;
pbar = (w3?beta?f);
pbarep = pbar?(4:172?Sy=E +0:0173)?sqrt(beta?f);
if(pbarep > pbar)
pbarep = pbar;
end
/********Calculation of constants********/
C = 1:295?exp(0:736?v1);
d = 3:8?(E ?beta=Sy=l)0:11;
/********Critical area and Critical force********/
Fcrit = 1=6=pi?(1=beta=f2=E)2 ?(C=2?Sy)3;
Acrit = (2=pi)?(C ?Sy=8=beta=f2=E)2;
/********Elasto-plastic area comparison********/
78
if(Fasp > Fcrit)
if(C4 > pbarep)
Afep(i) = Atotalep(i 1)=h=Atotalep(i 1);
else
Aep = 2?Acrit(d=(1+d))?(3?C4=f2=4=C=Sy)(d=(1+d));
AJ2ep = (1=f2)?(1 3=2=pi?(1 C4=pbarep));
Afep(i) = Aep?(1 (C4=pbarep)1:51)+AJ2ep?(C4=pbarep)1:04;
end
else
if(C4 > pbar)
Afep(i) = Atotalep(i 1)=h=Atotalep(i 1);
else
A1ep = (2?pi=f2)?(3=8=pi?(C4=pbar))(2=3);
A2ep = (1=f2)?(1 3=2=pi?(1 (C4=pbar)));
if((C4=pbar) < :8)
Afep(i) = A1ep?(1 (C4=pbar)1:51)+A2ep?(C4=pbar)1:04;
else
Afep(i) = A2ep;
end
end
end
Atotalep(i) = Atotalep(i 1)?h?Afep(i);
if(Atotalep(i 1) > Atotalep(i))
bep(i 1) = Atotalep(i);
b1ep(i 1) = bep(i 1)=Ari;
else
Atotalep(i) = Atotalep(i 1);
bep(i 1) = Atotalep(i);
79
b1ep(i 1) = bep(i 1)=Ari;
end
end
end
/********Lowest level of applied load********/
for(j = 10)
Lo = Fa?(10((j  1)=100 5));
Lo1(j) = Lo=E=Ari;
Atotal(1) = Ari;
Atotal(2) = Atotal(1);
Atotalep(1) = Ari;
Atotalep(2) = Atotalep(1);
/********Frequency level iteration********/
for(i = 2 : v)
I(i 1) = i;
/********Amplitude and Frequency calculation********/
beta = YY(i);
f = (i 1)=L;
l = 1=f;
h = 2?f2;
C3 = Lo=((Atotal(i 1)));
pbar = (w3?beta?f);
prat = C3=pbar;
/********Elastic multiscale model********/
if(C3 > pbar)
Af(i) = Atotal(i 1)=h=Atotal(i 1);
else
AJ1 = (2?pi=f2)?(w4?prat)(2=3);
AJ2 = (1=f2)?(1 (w2)?(1 prat));
80
/********Area comparison for every frequency level********/
if(prat < 0:8)
Af(i) = AJ1?(1 (prat)(1:51))+AJ2?(prat)(1:04);
else
Af(i) = AJ2;
end
end
Atotal(i) = Atotal(i 1)?h?Af(i);
if(Atotal(i 1) > Atotal(i))
a(j) = Atotal(i);
be(i 1) = a(j)=Ari;
else
Atotal(i) = Atotal(i 1);
a(j) = Atotal(i);
be(i 1) = a(j)=Ari;
end
/********Elasto-plastic multiscale model********/
C4 = Lo=Atotalep(i 1);
Fasp = C4=f2;
pbar = (w3?beta?f);
pbarep = pbar?(4:172?Sy=E +0:0173)?sqrt(beta?f);
if(pbarep > pbar)
pbarep = pbar;
end
/********Calculation of constants********/
81
C = 1:295?exp(0:736?v1);
d = 3:8?(E ?beta=Sy=l)0:11;
/********Critical area and Critical force********/
Fcrit = 1=6=pi?(1=beta=f2=E)2 ?(C=2?Sy)3;
Acrit = (2=pi)?(C ?Sy=8=beta=f2=E)2;
/********Elasto-plastic area comparison********/
if(Fasp > Fcrit)
if(C4 > pbarep)
Afep(i) = Atotalep(i 1)=h=Atotalep(i 1);
else
Aep = 2?Acrit(d=(1+d))?(3?C4=f2=4=C=Sy)(d=(1+d));
AJ2ep = (1=f2)?(1 3=2=pi?(1 C4=pbarep));
Afep(i) = Aep?(1 (C4=pbarep)1:51)+AJ2ep?(C4=pbarep)1:04;
end
else
if(C4 > pbar)
Afep(i) = Atotalep(i 1)=h=Atotalep(i 1);
else
A1ep = (2?pi=f2)?(3=8=pi?(C4=pbar))(2=3);
A2ep = (1=f2)?(1 3=2=pi?(1 (C4=pbar)));
if((C4=pbar) < :8)
Afep(i) = A1ep?(1 (C4=pbar)1:51)+A2ep?(C4=pbar)1:04;
else
Afep(i) = A2ep;
end
end
end
82
Atotalep(i) = Atotalep(i 1)?h?Afep(i);
if(Atotalep(i 1) > Atotalep(i))
bep(i 1) = Atotalep(i);
b1ep(i 1) = bep(i 1)=Ari;
else
Atotalep(i) = Atotalep(i 1);
bep(i 1) = Atotalep(i);
b1ep(i 1) = bep(i 1)=Ari;
end
end
end
figure(1);
loglog(I;be;0bo0;I;b1ep;0rs0;I;b2e;0bx0;I;b21ep;0rd0)
xlabel(0Frequencyiteration;i0)
ylabel(0Ar=A0n)
set(z;0Interpreter0;0tex0)
83
Appendix C
matlab code for frequency spectrum and load versus
contact area
hx = xlsread(0sdata0;1;0b2 : b39010);
hx = hx?10 10;
/********Length of the sample********/
v = length(hx);
L = 400e 6;
/********Fourier transform of surface pro le********/
v = v=2;
Y = fft(hx;v)=v;
v = v=2;
YY = sqrt(Y:?conj(Y));
/********Material properties********/
E = 1=(1 :332)?200e9;
Sy = 1e9;
v1 = 0:33;
/********Johnson area coe?cients********/
w2 = 3=2=pi;
w3 = sqrt(2)?pi?E;
w4 = 3=8=pi;
k = 1;
/********Nominal area and Initial Load********/
Ari = 3:3?1:5?10(  6);
Fa = 0:0001?E ?Ari;
I = 1 : 1 : v +1;
84
II = log(I);
/********Load iteration procedure from the intial load********/
for(j = 1 : 1 : 1000)
Lo = Fa?(10(5 (j  1)=100));
Lo1(j) = Lo=E=Ari;
Atotal(1) = Ari;
Atotal(2) = Atotal(1);
Atotalep(1) = Ari;
Atotalep(2) = Atotalep(1);
/********Frequency iteration for every level of applied load********/
for(i = 2 : v)
/********Amplitude and Frequency calculation********/
beta = YY(i);
f = (i 1)=L;
l = 1=f;
h = 2?f2;
C3 = Lo=((Atotal(i 1)));
/********Elastic multiscale model********/
/********Critical pressure and pressure ratio for elastic case********/
pbar = (w3?beta?f);
prat = C3=pbar;
/********Area comparison for every frequency level********/
if(C3 > pbar)
Af(i) = Atotal(i 1)=h=Atotal(i 1);
else
AJ1 = (2?pi=f2)?(w4?prat)(2=3);
AJ2 = (1=f2)?(1 (w2)?(1 prat));
if(prat < 0:8)
if(i) = AJ1?(1 (prat)(1:51))+AJ2?(prat)(1:04);
else
85
Af(i) = AJ2;
end
end
Atotal(i) = Atotal(i 1)?h?Af(i);
if(Atotal(i 1) > Atotal(i))
a(j) = Atotal(i);
else
Atotal(i) = Atotal(i 1);v
a(j) = Atotal(i);
end
/********Elasto-plastic multiscale model********/
/********Initial load, asperity force and critical pressure calculation********/
C4 = Lo=Atotalep(i 1);
Fasp = C4=f2;
pbar = (w3?beta?f);
pbarep = pbar?(4:172?Sy=E +0:0173)?sqrt(l=beta);
/********Comparing elastic and elasto-plastic critical pressures********/
if(pbarep > pbar)
pbarep = pbar;
end
/********Calculation of constants for elasto-plastic case********/
C = 1:295?exp(0:736?v1);
d = 3:8?(E ?beta=Sy=l)0:11;
Fcrit = 1=6=pi?(1=beta=f2=E)2 ?(C=2?Sy)3;
Acrit = (2=pi)?(C ?Sy=8=beta=f2=E)2;
/********Comparing asperity and critical force********/
if(Fasp > Fcrit)
if(C4 > pbarep)
Afep(i) = Atotalep(i 1)=h=Atotalep(i 1);
else
86
/********Elasto-plastic area comparison********/
Aep = 2?Acrit(1=(1+d))?(3?C4=f2=4=C=Sy)(d=(1+d));
AJ2ep = (1=f2)?(1 3=2=pi?(1 C4=pbarep));
Afep(i) = Aep?(1 (C4=pbarep)(1:51))+AJ2ep?(C4=pbarep)(1:04);
end
else
if(C4 > pbar)
Afep(i) = Atotalep(i 1)=h=Atotalep(i 1);
else
A1ep = (2?pi=f2)?(3=8=pi?(C4=pbar))(2=3);
A2ep = (1=f2)?(1 3=2=pi?(1 (C4=pbar)));
if((C4=pbar) < :8)
Afep(i) = A1ep?(1 (C4=pbar)(1:51))+A2ep?(C4=pbar)(1:04);
else
Afep(i) = A2ep;
end
end
end
Atotalep(i) = Atotalep(i 1)?h?Afep(i);
if(Atotalep(i 1) > Atotalep(i))
b(j) = Atotalep(i);
else
Atotalep(i) = Atotalep(i 1);
b(j) = Atotalep(i);
end
end
aar(j) = a(j)=Ari;
aarep(j) = b(j)=Ari;
end
xlswrite(0ssdata0;aar0;1;0d10)
xlswrite(0ssdata0;aarep0;1;0e10)
87
/********Plots********/
figure(2);
loglog(Lo1;aar;0k :0;Lo1;aarep;0k x0)
axis([010e 4010e 1])
xlabel(0F=EA0n)
ylabel(0Ar=A0n)
z = legend(0Elastic0;0Elasto plastic0;2);
set(z;0Interpreter0;0tex0)
figure(3);
x = logspace( 1;3);
loglog(I;YY(1 : v +1);0k?0)
xlabel(0i0)
ylabel(0LogAmplitude(m)0)
gridon
88