An Analysis of Elastic Rough Contact Models
by
Yang Xu
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
Aug 04, 2012
Keywords: Rough Surfaces Contact, Elastic, Statistical Model, Boundary Element Model,
Finite Element Model
Copyright 2012 by Yang Xu
Approved by
Dan Marghitu, Co-Chair, Professor of Mechanical Engineering
Robert Jackson, Co-Chair, Associate Professor of Mechanical Engineering
Hareesh Tippur, McWane Professor of Mechanical Engineering
George Flowers, Dean, Graduate School
Abstract
The problem of elastic contact between one nominally  at rough surface and one rigid
 at is studied. A variety of contact models have been applied to this problem since the
pioneering work of Archard in 1956. Those contact models can be divided into 4 categories
namely, the statistical, multi-scale, semi-numerical and deterministic numerical models. In
Chapter 1, the methods of describing and generating rough surfaces are introduced. Addi-
tionally, the Hertzian contact theory is brie y illustrated. In Chapter 2, the developments
of the statistical, multi-scale and semi-numerical models are reviewed brie y by introducing
several important models, e.g., the Greenwood-Williamson (GW) model. The contact area
to load relations predicted by those models are compared to that by deterministic numerical
model (the  nite element model) which is treated as the "exact" solution. In Chapter 3
and 4, the developments and applications of two popular deterministic numerical methods
namely, the boundary element method (BEM) and the  nite element method (FEM), for the
elastic rough surface contact are reviewed in detail. The contact area to load relations pre-
dicted by the BEM and the FEM are compared. The contact areas and the contact pressure
distributions at di erent stages of deformation are also explored. Conclusions of this study
can be drawn in the following aspects: (1) the asperity-based models, e.g., the statistical and
semi-numerical models, are only valid for low load condition, (2) the simpli ed multi-scale
model is good for both early contact and complete contact, (3) the contact area to load
relations predicted by the FEM and the BEM have very good agreement, at least for the
case of isotropic surfaces.
ii
Table of Contents
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
List of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Rough surface geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Statistical description of rough surface . . . . . . . . . . . . . . . . . 1
1.1.2 Method of generating rough surface . . . . . . . . . . . . . . . . . . . 5
1.1.3 Rough surfaces data used in this study . . . . . . . . . . . . . . . . . 7
1.2 Single asperity contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Elastic asperity contact . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2 Statistical, multi-scale and semi-numerical rough surfaces contact model . . . . 13
2.1 Statistical rough surfaces contact model . . . . . . . . . . . . . . . . . . . . 13
2.1.1 Greenwood-Williamson (GW) model . . . . . . . . . . . . . . . . . . 14
2.1.2 Advanced statistical models . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Multi-scale model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1 Majumdar and Bhushan (MB) fractal contact model . . . . . . . . . 18
2.2.2 Jackson and Streator (JS) multi-scale model . . . . . . . . . . . . . . 19
2.3 Semi-numerical deterministic model . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 Ciavarella et al?s model . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3 Finite element method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
iii
3.1 Surface modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2 Boundary condition and e ect of thickness . . . . . . . . . . . . . . . . . . . 35
3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4 Boundary element method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1 Boundary integral equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2 Formulation of rough surfaces contact problem . . . . . . . . . . . . . . . . . 51
4.3 Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.4 Fast algorithm for deformation integral . . . . . . . . . . . . . . . . . . . . . 56
4.5 Convergence criterion and surface separation . . . . . . . . . . . . . . . . . . 58
4.6 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
A In uence Coe cient Kijkl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
iv
List of Figures
1.1 Real rough surface measured by pro lometer. . . . . . . . . . . . . . . . . . . . 2
1.2 Real rough surface pro le measured by pro lometer. . . . . . . . . . . . . . . . 3
1.3 Details of a 2D asperity pro le in rough surface pro le in Fig. 1.2 . . . . . . . . 4
1.4 Details of a 2D rough surface pro le with 15 sampling points from Fig. 1.2 . . . 5
1.5 3 dimensional pro le (a) and colored contour (b) of surface 1, i.e., fractal surface. 7
1.6 3 dimensional pro le (a) and colored contour (b) of surface 2, i.e., real rough
surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.7 Asperity contact model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1 Contact area-load relation of surface 1 predicted by statistical, multi-scale and
FEM models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 Contact area-load relation of surface 2 predicted by statistical, multi-scale and
FEM models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.3 Contact area-load relation of surface 1 predicted by simpli ed multi-scale models
[34] with di erent de nition of Bmax. . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4 Contact area-load relation of surface 2 predicted by simpli ed multi-scale models
[34] with di erent de nition of Bmax. . . . . . . . . . . . . . . . . . . . . . . . . 29
v
3.1 Schematic representation of uniformly meshed elastic solid with smooth surface
(only the mesh on the nominal contact surface are visible). . . . . . . . . . . . . 32
3.2 Schematic representation of a badly distorted rough surface element. . . . . . . 34
3.3 Details of surface mesh of surface 1. . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.4 Details of surface mesh of surface 2. . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.5 Schematic representation of cyclic symmetry and displacement coupling . . . . . 37
3.6 Schematic representation of boundary conditions: (a) cyclic symmetry, (b) con-
strained side surfaces and (c) free side surfaces . . . . . . . . . . . . . . . . . . . 38
3.7 Dimensionless contact area A to dimensionless load P relation of surface 1. The
FE models with di erent boundary conditions illustrated in Fig. 3.6. . . . . . . 39
3.8 Dimensionless contact area A to dimensionless load P relation of surface 2. . . 40
3.9 Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.10 Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.11 Dimensionless contact area A to dimensionless load P relation of surface 1.
Value of thickness: H = 0:5Lx; Lx; 5Lx. . . . . . . . . . . . . . . . . . . . . . . 41
3.12 Dimensionless contact area A to dimensionless load P relation of surface 2.
Value of thickness: H = 0:5Lx; Lx; 5Lx. . . . . . . . . . . . . . . . . . . . . . . 42
3.13 Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 1. Value of thickness: H = 0:5Lx; Lx; 5Lx. . . . . . . . . . . . . . . . 42
vi
3.14 Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 2. Value of thickness: H = 0:5Lx; Lx; 5Lx. . . . . . . . . . . . . . . . 43
3.15 Contact area contours and contact pressure distribution of surface 1 when contact
load is light(a-b), medium(c-d) and high(e-f). Light load: P = 4:49 10 4,
A = 0:0489. Medium load: P = 0:00218, A = 0:283. High load: P = 0:00438,
A = 0:948. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.16 Contact area contours and contact pressure distributions of surface 2 when con-
tact load is light(a-b), medium(c-d) and high(e-f). Light load: P = 1:91 10 5,
A = 0:0173. Medium load: P = 9:62  10 4, A = 0:419. High load:
P = 0:00438, A = 0:894. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.17 Correlation of rough surface geometry (a) and pressure distribution (b) of surface
1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.18 Correlation of rough surface geometry (a) and pressure distribution (b) of surface
2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1 Semi-in nite contact body and its rough boundary. . . . . . . . . . . . . . . . . 49
4.2 Schematic representation of displacement due to uniform pressure distribution. . 51
4.3 Contact area to load relation of surface 1 predicted by the boundary element
method and the  nite element method. . . . . . . . . . . . . . . . . . . . . . . . 61
4.4 Contact area to load relation of surface 2 predicted by the boundary element
method and the  nite element method. . . . . . . . . . . . . . . . . . . . . . . . 62
4.5 Schematic representation of rough surfaces deformation corresponding to dif-
ferent FFT-based fast integral method: NR+DC-FFT before (a) and after (c)
deformation; VM+FFTSK before (b) and after (d) deformation. . . . . . . . . . 62
vii
List of Tables
1.1 Statistical and fractal parameters of two surfaces . . . . . . . . . . . . . . . . . 8
viii
List of Abbreviations
(d0) Dimensionless surface separation, (d0) = d0= 
 A Contact area on the single asperity
 P Contact load on the single asperity
 De ection (or penetration, indentation) of the single asperity
 ,  Amplitude and wavelength
 c Critical asperity de ection
 x,  y Resolutions in the x and y directions
 s Asperity density
 1,  2 Poisson?s ratio of contact body 1 and 2
 Probability density function
 s Root mean square of asperity height
A Total contact area
a Radius of contact area on the single asperity
A Dimensionless contact area, A = A=An
An Nominal contact area
Ar Real contact area
B Ratio of amplitude to wavelength, B =  = 
ix
Bmax Maximum value of B
d Surface separation between the rigid  at and the mean level of undeformed rough
surface
d0 Surface separation between the rigid  at and the mean level of deformed rough
surface
d Dimensionless surface separation, d = d= 
E E ective material modulus, E = [(1  21)=E1 + (1  22)=E2] 1
E1, E2 Young?s modulus of contact body 1 and 2
g Gap between the rigid  at and the deformed rough surface
G;D Fractal parameters
H Thickness of the substrate of the rough surface
h Rough surface height
H0 Hurst exponent
h0x, h00x Slope and curvature in the x direction
h0y, h00y Slope and curvature in the y direction
Kijkl(Dms) In uence function
Lx, Ly Sampling length in the x and y directions
m0;m2;m4 Spectral moments
P Total contact load
p Pressure distribution
x
P Dimensionless contact load, P = P=(E An)
R Average radius of asperities or e ective radius in the Hertzian contact theory, R =
(R 11 +R 12 ) 1
R( ) Autocorrelation function
Ra Center-line average
Rq( ) Root mean square (RMS) roughness
Rx, Ry Radius of asperity in the x and y directions
S(w) Power spectrum density
Sy Yield stress
w Rough surface deformation
xi
Chapter 1
Introduction
1.1 Rough surface geometry
1.1.1 Statistical description of rough surface
In this section, methods of describing rough surface geometries are introduced. Two
popular ways of generating the rough surface by computer code will be discussed. At the
end of this section, details of two rough surfaces data, one generated by computer code and
another obtained from pro lometer, will be illustrated. These two rough surface data will
be used as the inputs of all the elastic rough surfaces contact models discussed in the rest of
the thesis.
Figs. 1.1 and 1.2 illustrate the typical rough surface and rough surface pro le along the
x direction at certain y coordinate measured by pro lometer. Fig. 1.3 shows a 2D asperity
pro le from the rough surface pro le in Fig. 1.2. The peaks on the rough surface or rough
surface pro le are called asperities. Sampling points, i 2;i 1;i;i+ 1 and i+ 2, and their
connecting lines form the pro le of a 2D asperity in Fig. 1.3. Sampling point on the asperity
is called summit, i in Fig. 1.3, when its height is larger than all of its neighboring sampling
points, i 1 and i+ 1.
If we neglect the surface roughness when two elastic bodies are in contact, the corre-
sponding contact area is called nominal contact area, An. When the contact interfaces are
rough, the real contact area, Ar, is a fraction of nominal contact area, An. The ratio of Ar
to An is called the contact ratio.
According to the size, rough surface features can be divided into 3 categories [1] namely,
macrodeviation, waviness and roughness. Surfaces of road and landscape are belonging to
1
Figure 1.1: Real rough surface measured by pro lometer.
macrodeviation because the asperity height of those surfaces is within macroscale and usually
the ratio of amplitude to wavelength is considerably large. For waviness surfaces, the ratio
of distance between neighboring asperities to the asperity height is usually more than 40
[1]. However, the waviness pro le is still visible. Roughness are de ned as peaks and valleys
which have very small value of height deviations comparing to dimensions of sampling length
(length of nominal contact area) and are generally invisible to the unassisted eye.
It is di cult to describe rough surface deterministically due to its multi-scale structure,
which will be discussed later, and random nature. Therefore, simplifying mathematics such
as statistical parameters are used. Here we only illustrate the statistical parameters of a
rough surface pro le (2D). Extension to rough surface (3D) is very similar to the 2D case.
hi;i = 1;:::;N is the height of rough surface pro le obtained from a pro lometer with the
resolution  x and one small portion of a rough surface pro le with 15 sampling points is
2
0 100 200 300 400 500 600 700 800 900 1000
?0.04
?0.03
?0.02
?0.01
0
0.01
0.02
0.03
0.04
 x [10
?6
 m]
 z [mm]
Figure 1.2: Real rough surface pro le measured by pro lometer.
plotted in Fig. 1.4. N is number of sampling points. Ra, center-line average, describes the
roughness of rough surface pro le
Ra = 1N
NX
i=1
jhij
and level at the average height of surface pro le is called mean level. Another similar param-
eter, the root mean square (RMS) roughness,  or Rq, is also widely used as a measurement
of average roughness
 =
vu
ut 1
N
NX
i=1
h2i
and  is more sensitive to large deviations from the mean level [2]. However, the center-line
average, Ra, and root mean square,  , cannot uniquely describe the rough surface pro le
3
x
z
i-1
i+1
i
i-2
i+2
Figure 1.3: Details of a 2D asperity pro le in rough surface pro le in Fig. 1.2
statistically, thus the spectral moments method are introduced [11]
m0 = 1N
NX
i=1
hi2
m2 = 1N 2
N 1X
i=2
h0i2
m4 = 1N 2
N 1X
i=2
h00i 2:
h0i and h00i are the slope and curvature of sampling point i, respectively. h0i and h00i can be
numerically determined by using central di erence
h0i = hi+1 hi 12 
x
h00i = hi+1 2hi +hi 1 2
x
; i = 2;:::;N 1: (1.1)
4
x
z
i
K
L
mean level
mean aseprity level
d0
?[
rigid flat
d
Figure 1.4: Details of a 2D rough surface pro le with 15 sampling points from Fig. 1.2
m0 is the same as the square of  , i.e., m0 =  2. m2 and m4 are called skewness (SK) and
Kurtosis (K) and commonly normalized by  2 and  4.
In many elastic rough surfaces models, only the heights of the summits are used as
the input due to the fact that contacts dominantly occur on asperities. Greenwood and
Williamson [3] found that asperity height in many real surfaces can be described by a Gaus-
sian distribution
 (h) = 1p2  
s
exp
  h2
2 2s
 
where surface height is centered about the asperity mean level, see Fig. 1.4.  s is the root
mean square of the asperity height ( is the root mean square of the rough surface pro le).
Asperity mean level is determined by averaging the asperity heights.
1.1.2 Method of generating rough surface
Commonly, rough surface data is obtained using a pro lometer. For model evaluation,
an alternative is to generate fractal surfaces, which will be discussed later, by the com-
puter code. Two popular algorithms of generating fractal surfaces namely, the Weierstrass-
Mandelbrot function and the Random Midpoint Algorithm, can be found in the past litera-
tures.
Fractal geometry
5
Fractal geometry [22, 23] can be used to characterize the multi-scale nature of the rough
surfaces. A fractal geometry should be continuous, non-di erentiable and statistically self-
a ne, in a mathematical sense. Real rough surfaces obtained from pro lometer show exactly
the continuous and non-di erentiable properties. If a similar surface pro le is obtained as
the sampling resolution  x decreases, it can be described as statistically self-a ne. Due to
the above reasons, fractal geometries can be used to approximate real rough surfaces with a
multi-scale structure. It is found that the Weierstrass-Mandelbrot (WM) function meets all
the properties of a fractal geometry. The following WM function [23] can be used
h(x) = G(D 1)
nmaxX
n=n1
cos(2  nx)
 (2 D)n ; 1 D 2 (1.2)
to approximate the real 2D rough surface pro le. G is the length scale of the rough surface
pro le. D is the fractal dimension.  n describes the frequency spectrum of the surface
roughness. n1 and nmax relate to the cut-o frequencies which are dependent on sampling
length Lx and sampling resolution  x. Since a rough surface pro le has a  nite length,
frequency cannot be less than the lowest frequency  n1  1=Lx. Also, the rough surface
pro le does not have any surface height information beyond the highest frequency level
 nmax  1= x. For a 3D rough surface [24, 25], rough surface height can be written as
h(x;y) = L
 G
L
 (D 2) ln 
M
 1=2 MX
m=1
n0maxX
n=1
 (D 2)nfcos( m;n) 
cos
 2  n(x2 + y2)1=2
L cos
 
tan 1
 y
x
 
  mM
 
+ m;n
 
g; 2 D 3 (1.3)
where n0max = int
 
log(L= )
log 
 
is the maximum frequency level. L and  are the sampling
length and sampling resolution (same in the both directions) of the 3D fractal surface, re-
spectively. M denotes the number of superposed ridges used to construct the surface.  m;n
is a random variable which can be achieved by any random generator. Majumdar and Tien
[27] showed  = 1:5 is a suitable value for high spectral density and phase randomization.
6
Fractal parameters G and D can be extracted from rough surface data. A fourier transfor-
mation of the autocorrelation function R( ) of fractal surfaces generated by Eq. 1.2 and 1.3,
i.e., the power spectrum density S(w), may result in the following power law form
S(w) = G
2(D 1)
2ln 
1
w(5 2D) (1.4)
where w is the frequency of roughness. If S(w) is plotted in a log-log plot, D is related
to the slope of S(w) and G can be extracted from S(w = 0). It is obvious that G and D
are scale-independent parameters comparing to the scale-dependent surface parameters, e.g.,
m0, m2 and m4 [27, 26].
Random Midpoint Algorithm
Another way of generating self-a ne fractal surface, which is called random midpoint
algorithm, is developed by Voss [28]. In stead of scale-invariant fractal parameters used
in WM function, Eqs. 1.2 and 1.3, only one parameter H0, Hurst exponent, is needed.
The relationship between Hurst exponent and Fractal dimension is: H0 = 3 H(3D) and
H0 = 2 H(2D).
1.1.3 Rough surfaces data used in this study
0
0.5
1
1.5
2
x 10
?5
0
0.5
1
1.5
2
x 10
?5
0
1
2
3
4
5
x 10
?8
 
x[m]
y[m]
 
z[m]
x[m]
y[m]
 
 
2 4 6 8 10 12 14 16
x 10
?6
2
4
6
8
10
12
14
16
x 10
?6
0.5
1
1.5
2
2.5
3
3.5
4
x 10
?8
[m]
(a)
(b)
Figure 1.5: 3 dimensional pro le (a) and colored contour (b) of surface 1, i.e., fractal surface.
7
Table 1.1: Statistical and fractal parameters of two surfaces
Parameter Surface 1 Surface 2 Parameter Surface 1 Surface 2
m0 (m2) 6:7081 10 16 5:9061 10 14 m2x 3:6506 10 4 2:6607 10 5
m4x (m 2) 8:2388 109 8:7907 105 m2y 3:9087 10 4 9:7553 10 5
m4y (m 2) 9:1915 109 2:3517 107  s(m) 3:4070 10 8 2:4331 10 7
 s (m 2) 8:8281 1011 1:3600 109 R(m) 1:2098 10 5 0:0013
G (m) 9:4600 10 14 D 2:4400
Lx (m) 1:6000 10 5 9:9200 10 4 Ly (m) 1:6000 10 5 9:9200 10 4
 x (m) 1:25 10 7 8 10 6  y (m) 1:25 10 7 8 10 6
0
0.2
0.4
0.6
0.8
1
x 10
?3
0
0.5
1
x 10
?3
?2
?1.5
?1
?0.5
0
0.5
1
x 10
?6
x[m]
y[m]
z[m]
x[m]
y[m]
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
1
2
3
4
5
6
7
8
9
10
x 10
?4
?15
?10
?5
0
5
x 10
?7
[m]
(a)
(b)
Figure 1.6: 3 dimensional pro le (a) and colored contour (b) of surface 2, i.e., real rough
surface.
Two rough surfaces data namely, surface 1 and surface 2, are used throughout the rest of
the thesis as the inputs of all the elastic rough surfaces contact model. Surface 1 is generated
by the WM function, Eq. 1.3, and surface 2 is measured by pro lometer from a standard
rough surface sample. The main statistical and fractal parameters are illustrated in Table.
1.1.  s and R are the asperity density and average radius of asperity, respectively.  s is the
ratio of number of asperities to the nominal area An. The radius of curvature, Ri, of each
8
asperity i can be de ned as
Ri = (Rx +Ry)=2 (1.5)
Rx = (1 (h0x)2)1:5=h00x
Ry = (1 (h0y)2)1:5=h00y;
h0x, h0y, h00x and h00y are slopes and curvatures of asperity i along x and y directions (Eq. 1.1).
Lx and Ly are sampling lengths in x and y directions, An = Lx Ly.  x and  y are the
resolutions along x and y directions. Surface 1 consists of 129 129 sampling points and
surface 2 is 125 125.
m2 and m4 along both x and y are m2x, m2y, m4x and m4y and are listed in Table 1.1.
Di erences between spectral moments of surface 1 in the x and y directions are negligible
comparing to that of surface 2. However, they di er greatly for surface 2. Consequently,
surface 1 can be treated as an isotropic surface and surface 2 as an anisotropic surface.
1.2 Single asperity contact
Between two approaching rough surfaces, contacts initially occur at the summits of
asperities. Thus, at low load conditions, the majority of contact areas grow independently
(without interactions between them). Thus the resultant rough surface contact load P and
contact area A are the superpositions of contact load  P, and the area of single contacting
asperities  A.
P =
NcX
i=1
 Pi
A =
NcX
i=1
 Ai (1.6)
Note that Nc is the number of contacting asperities at a given surface separation d. Com-
monly, rough surface contact models are called asperity-based contact models if Eq. 1.6 is
9
used to obtain the total contact force P and area A. As a result, it is necessary to derive an
accurate single asperity contact model for the elastic contact accounting for di erent shapes
of contact interfaces.
It is unlikely to accurately approximate the complex shapes of asperities using simple
geometrical entities. At the low load condition, nearly all the contacts occur near the tips
of the asperities and many e orts have been put on modeling the geometry of the tip of
asperity and its corresponding asperity contact model. The most popular geometrical entity
used to approximate the shape of tip is a hemisphere, see Fig. 1.7.
1.2.1 Elastic asperity contact
Deformed
Undeformed
R
Rigid Flat
a
Figure 1.7: Asperity contact model
As we mentioned before that spherical cap is the most popular geometrical entity to
approximate the tip of asperities where contacts are more likely to occur, thus the well-
known Hertzian spherical contact theory can be used here to model the contact between two
hemispherical caps of radius R1 and R2, respectively. Contact between two hemispherical
caps can be further simpli ed to the contact between an equivalent hemispherical cap and
one rigid  at, see Fig. 1.7. The e ective radius, R, can be written as
R = (R 11 +R 12 ) 1
and the Young?s modulus of the equivalent hemispherical cap (e ective material modulus,
E ) is
E = [(1  21)=E1 + (1  22)=E2] 1
10
where Ei and  i, i = 1;2, are Young?s modulus and Poisson?s ratio of the two hemispherical
caps, respectively. In the rest of the article, the contact between two asperities will always
be simpli ed as the contact between one asperity of e ective radius R and one rigid  at.
Thus the contact area,  A, and contact load,  P, due to contact between two hemispherical
caps, are described by Hertz [32]
 A =  a2 =  R (1.7)
 P = 4
3E
 R12 32: (1.8)
where a is the radius of contact area.  is the de ection (or penetration, indentation) as
illustrated in Fig. 1.7. The pressure distribution de ned in polar coordinate, p(r), within
the contact area is
p(r) = p0
p
1 r2=a2; r2[0;a]
where p0 is maximum contact pressure.
Commonly, critical deformation,  c, is used to predict the onset of yielding in the sphere.
Yielding  rst occurs beneath or on the contact interface between two asperities with hemi-
spherical tips. When deformation is larger than  c, i.e.,  >  c, deformation beneath or on
the contact interface may become elastoplastic. For the ductile material, the onset of  rst
yield is governed by the von-Mises criterion. Chang et al [4] and Jackson and Green [5]
derived similar expressions of  c based on the von-Mises criterion. The one from Jackson
and Green is
 c =
  p0
2E 
 2
R
where
p0=Sy = 1:295 exp(0:736 )
where Sy is the yield stress of the material in unidirectional tension (or compression). Sy
and  are the values corresponding to the minimum p0. For the brittle materials, onset of
11
 rst yield is governed by Tresca maximum tensile stress criterion and the corresponding  c
is derived numerically by Brizmer et al[6].
The above derivation of critical deformation  c is only valid for the perfect slip condition.
 c for full stick conditions can be found in [6].  c for the partial slip condition is currently
not available.
12
Chapter 2
Statistical, multi-scale and semi-numerical rough surfaces contact model
2.1 Statistical rough surfaces contact model
The statistical type of contact model is still the most popular model used in rough
surfaces contact. In stead of using the complete roughness data, only probability density
function  (h) is used.  (h) means the probability of the asperity with the height between
h and h+dh. Assume we know the exact expression of  (h), then the number of asperities
with height h is
N(h) =  (h)  N
where  N is the total number of asperities in the rough surface. Thus, the total contact load
(P) and area (A) as a function of d can be written as
P(d) =  N
+1Z
d
 P (h)dh (2.1)
A(d) =  N
+1Z
d
 A (h)dh (2.2)
d is the surface separation between rigid  at and mean asperity level, see Fig. 1.4. Defor-
mation  of each contacting asperity is
 = h d
Commonly, two rough surfaces contact can be simpli ed as one equivalent rough surface
contacting a rigid  at. The height of equivalent surface, h(x), can be written as h(x) =
13
h1(x)+h2(x). h1 and h2 are two rough surfaces. Additionally, the spectral moments, m0;m2
and m4, of the equivalent rough surface pro le are
m0 = (m0)1 + (m0)2 (2.3)
m2 = (m2)1 + (m2)2 (2.4)
m4 = (m4)1 + (m4)2 (2.5)
(m0)i;(m2)i and (m4)i;i = 1;2, are spectral moments of two rough surfaces. Eq. 2.1 was
 rst developed by Greenwood and Williamson [3].
2.1.1 Greenwood-Williamson (GW) model
In 1966, Greenwood and Williamson combined random process with classical Hertzian
contact to deal with rough surfaces contact. In their pioneering theory, they adopted the
following assumptions:
1. The asperity height distribution is Gaussian, i.e., the probability density function  (h) is
 (h) = 1p2  
s
exp
 
 h
2
2 2s
 
(2.6)
where  s is the standard deviation of asperity height. The asperity height h is centered
about mean asperity level.
2. Asperity contact is modeled by the Hertzian spherical contact theory (Eq.1.7);
3. The radius of all asperities are assumed constant (R);
4. Ignore adhesion, shoulder-shoulder contact, neighboring asperity interaction and asperity
coalescing.
By using dimensionless asperity height h = h= s and dimensionless surface separation
d = d= s. Asperity probability density function  (h ), total contact load P(d ) and area
14
A(d ) in the GW model can be written as
 (h ) = 1p2 exp
 
 (h
 )2
2
 
(2.7)
P(d ) = An s
+1Z
d 
 P (h )dh (2.8)
A(d ) = An s
+1Z
d 
 A (h )dh (2.9)
An and  s are nominal contact area and asperity density, respectively. Multiplication of
An and  is the total number of asperities in the rough surface sample (  N).  P and  A are
calculated by Eq. 1.7.
2.1.2 Advanced statistical models
Whitehouse and Archard (WA) model
Whitehouse and Archard [9] (WA) abandoned the assumption of a constant asperity ra-
dius. Rough surfaces obtained from a pro lometer show that higher asperities have sharper
tips. Whitehouse and Archard incorporated this correlation into the surface statistics and
assumed an exponential-like autocorrelation function, R( ). They derived the joint prob-
ability density function,  (h ;C), of both summit height and summit curvature. C is a
non-dimensional summit curvature. Their joint probability density function was used by
Onions and Archard [10] in deriving elastic contact load (P) and area (A) following the
framework of the GW model, Eq. 2.7.
Nayak Model
Nayak [11] (1971) derived more generalized joint probability functions for both surface
pro le (2D) and rough surfaces (3D) based on 2 and 3 dimensional joint probability function:
 (h;@h@x;@2h@x2 ) and  (h;@h@x;@h@y;@2h@x2; @2h@x@y;@2h@y2 ).
Bush et al (BGT) Model
15
For simplicity, the curvatures of the contacting summits used in Nayak?s statistical
model is mean curvature  m = ( @2h@x2 + @2h@y2 )=2 and Hertzian spherical contact theory is
used. Bush et al [12] (BGT) (1975) extended the Nayak statistical model by considering
two principle curvatures of asperities along the x and y directions. Asperity contact model
in the BGT model becomes more generalized Hertzian elliptic contact model. Due to its
complexity, only the asymptotic solution of the BGT model [14] is written here
A=An =
r  
m2
P
E An (2.10)
which has a linear relation between contact area A and contact load P.
McCool Model
The above mentioned advanced statistical models all introduced more than one random
variables in their joint probability density function  which cause their models to be more
complex (multiple integrals) than the GW model and inconvenient to apply. The GW model
has a clean-cut expression but the inaccurate inputs of asperity density  s, asperity r.m.s
 s and constant radius R may cause the inaccurate contact load (P) and area (A). McCool
[13] (1987) derived the following closed-form expression of inputs for GW model
 s = 16 p3
 m
4
m2
 
(2.11)
R = 0:375
  
m4
 0:5
(2.12)
 s =
 
1 0:8968 
 0:5
m0:50 (2.13)
based on the Nayak statistical model for isotropic rough surfaces.  = m0m4m2
2
is called the
bandwidth parameter introduced by Nayak [11].
Zhao and Chang (ZC) Model and Ciavarella et al (CGP) Model
In the above statistical models, the substrate beneath the rough surface is assumed to be
rigid, thus the interaction between the neighboring asperities due to the elastic deformation
16
of substrate is neglected and the local asperity deformation is calculated as:  = h d. Due
to the asperity interaction, the contact load may cause the redistribution of asperity height.
Consequently, the probability density function  of undeformed asperity distribution may
no longer be the same. Zhao and Chang (ZC)[15] derived a modi ed asperity deformation
equation
 = h d+ 1:12
pw
lpm
E 
by introducing the displacement of the mean asperity level caused by contact load. wl is the
contact load on the single asperity and pm is the average pressure over the nominal contact
area An, i.e., pm = P=An. When modi ed asperity deformation  is applied in the GW
model, Eq. 2.7, the modi ed GW model becomes nonlinear and the iterative method has
to be used in order to obtained a convergent solution. Ciavarella et al (CGP) [16] derived a
similar modi ed asperity deformation equation.
 = h d pmA1=2n =E 
A comprehensive comparison between GW+McCool (GW model with McCool?s inputs),
full BGT and Nayak models have recently been performed by Carbone and Bottiglione [14].
Two important conclusions have been drawn:
1. The contact area-load curve of all the models quickly deviate from the linear relation
obtained at the light load condition.
2. The slope of contact area-load curve predicted by BGT and Nayak models converges to
unity as surface separation d reaches a very high value with vanishing contact load and area.
However, GW+McCool doesn?t have such unity slope at high values of d.
Even though the advanced statistical models, e.g., Nayak and BGT models, are derived
from a complete joint probability density function with multiple random variables, they are
proved to be inaccurate by numerical simulations [34] at medium and heavy load conditions.
The major reason that causes the inaccuracy is that the statistical model is based on single
17
asperity contact model. At the low load condition (small contact area and large surface
separation), the majority of contact areas of single asperities grow independently, i.e. no
asperity coalescing occurs. Thus the contact load and area are the summation of contact
loads and areas on single asperities. The contact load-area curve shows a linear relation
only at vanishing load condition. When rough surface contact reaches medium and even
heavy load conditions, asperity coalescing becomes more and more severe. Usually statistical
models may therefore overestimate both the contact area and load [34].
2.2 Multi-scale model
Experimental results [2, 23, 26] show that the statistical surface parameters, e.g., spec-
tral moments m0;m2 and m4, are dependent on sampling resolution. Consequently, results of
statistical model may deviate from each other if the same rough surface is measured with dif-
ferent resolutions. Due to the multi-scale nature of rough surfaces, Archard [21] developed
the  rst multi-scale elastic rough surfaces contact model (perhaps the  rst rough surface
contact model). The structure of rough surfaces used in Archard?s model is described as
"protuberances on protuberances" type. Each protuberance has the hemispherical shape
and more smaller protuberances are uniformly distributing on its surface. As more scales
are being involved, the contact area-load relation approaches linearity. However, Archard?s
multi-scale model cannot be practically used on a real rough surface because of its unrealistic
rough surface structure.
2.2.1 Majumdar and Bhushan (MB) fractal contact model
Majumdar and Bhushan (MB)[23] developed the  rst elastic multi-scale contact model
based on fractal geometry. The MB model can only deal with a 2-D rough surface pro le
with the fractal dimension D lying between 1 and 2. They assumed that the contact spot
distribution n(a) is
n(a) = D2 a
D=2
l
aD=2+1 (2.14)
18
where al is the maximum contact spot area and can be calculated by al = A2 DD . A is the
contact area and can be obtained by the relation
A=An = 1p2  
+1Z
d
exp( h2=2)dh (2.15)
which assumes that the surface pro le follows a Gaussian distribution.  is the r.m.s of
rough surface pro le. Actually, the fractal surface can be treated as the stack of sinusoidal
waves with random phases based on the mathematical expression (Eq. 1.2), and contact load
P calculated in the MB model is the superposition of contact load  P of di erent frequency
domain.  P is calculated based on Hertzian spherical contact theory and asperity deformation
is assumed as the amplitude of each sinusoidal wave. Then we have
P =
alZ
0
4p E GD 1a(3 D)=2
3 n(a)da: (2.16)
Eq. 2.15 and Eq. 2.16 show a nonlinear relation between A and P. The MB fractal model
has been criticized due to the following reasons:
(1) The contact spot distribution, Eq. 2.14, is only geometrically-dependent and it has
nothing to do with the elastic properties and contact condition;
(2) The real contact area function, Eq. 2.15, is only valid for the surface follows Gaussian
distribution, which is unrealistic for the fractal surface.
(3) Using amplitude of each sinusoidal wave as the asperity deformation is too simpli ed.
2.2.2 Jackson and Streator (JS) multi-scale model
Jackson and Streator (JS)[29] developed an isotropic 3-D multi-scale contact model
based on 3D sinusoidal contact model ([30, 31]). The main ideal of the JS model is that a
2-D rough surface pro le can be decomposed into stacks of sinusoidal waves with di erent
frequencies fi (or wavelength  i) and amplitude  i. Asperity density  i and asperity contact
19
radius Ri of each sinusoidal waves can be easily determined:
 i = 2fi
Ri = 14 2 
if2i
on each frequency level. Johnson et al [31] provided two asymptotic solutions about contact
between a sinusoidal surface and a rigid base. When  p p ,
(  AJGH)1 =  f2
i
3
8 
 p
p 
and  p approaches p 
(  AJGH)2 = 12f2
i
(1 32 [1  p=p ])
where  p is the average contact pressure over nominal contact area ( 2i) and p is the average
pressure to cause complete contact
p =p2 E  i= i:
Jackson and Streator curve- tted the two asymptotic solutions to bridge the gap between
them:
 Ai = (  AJGH)1
 
1 
  p
p 
 1:51!
+ (  AJGH)2
  p
p 
 1:04
(2.17)
 Ai = (  AJGH)2 (2.18)
The above empirical functions have been validated by both  nite element method [33] and
FFT-based deterministic model [34]. The assumptions made in the JS model are:
1. Higher frequency domains are superimposed upon lower frequency domain, which is
similar to Archard?s model [21]. The di erence is that the shape of protuberance used in
20
Archard?s model is spherical and while it is sinusoidal in the JS model.
2. Each scale carries the same total load
3. A given scale cannot increase the contact area beyond what is experienced by the frequency
level below it. Assumption 3 ensures real contact area, Ar, is bounded (Ar  An) and
guarantees the convergence of Ar as frequency level increases.
4. The load at each frequency level is shared equally among all the asperities at that level.
Thus the following contact model can be gained
A =
 i
maxY
i=1
 Ai i
!
An (2.19)
P =  Pi iAi 1 (2.20)
By using fast Fourier transformation (FFT), we can decompose the rough surface into imax
frequency levels. For a given contact load P, contact load  Pi =  pi( 2i) on each asperity
within frequency domain i can be determined by Eq. 2.20. Then corresponding contact area
 Ai can be calculated. These steps are iteratively performed until frequency level (imax) is
reached. Based on Eq. 2.19, the real contact area (A) can be determined.
Jackson [35] further simpli ed the JS model by using
A = Pp2 E0B
max
; Bmax = max(B) (2.21)
and this simpli ed multi-scale model, Eq. 2.21, can be applied for both 2D and 3D rough
surface data. For 3D rough surface data, B =  i= i is derived by taking the FFT of each 2D
row of 3D surface and then averaging the results [34].
2.3 Semi-numerical deterministic model
In the above sections, we brie y introduced the statistical models, e.g., the GW model,
and multi-scale models, e.g., the MB and JS models. In the statistical models, the concept
21
of asperity contact is essential in determining the contact load and area. Additionally, only
limited information of roughness, e.g., spectral moments (m0, m2, and m4) in the GW
model, is utilized. This is also true for the multi-scale. Only frequencies and amplitudes of
each sinusoidal components are used as the inputs. In the next chapter, the application of
deterministic numerical methods for rough surface contact will be introduced in detail. Since
a majority of contacts occur at the vicinity of asperities, coordinates of asperities (strictly
speaking the summits) are essential in predicting the contact load-area relation. This ideal
can be achieved by using the coordinates of summits and corresponding radius of curvatures
as the inputs. Because the real surface roughness data is still partially utilized, we call the
rough surfaces contact models, which adopt the above ideal, the semi-numerical models in
this thesis.
Considerable e orts are been placed on the development of semi-numerical deterministic
models in the following aspects:
(1) The Radius of curvature distribution are considered and more complex methods of eval-
uating the radius of curvature at each contact asperity are developed [72, 73, 74];
(2) More accurate asperity interaction models are included [75, 76, 65, 77];
(3) Asperity coalescing is included [78, 79, 80, 81, 77].
The concept of asperity contact is retained which makes these models di erent from
numerical deterministic models. However, in order to model more accurate asperity interac-
tion, asperity coalescing and asperity radius of curvature, coordinates of summits are utilized
instead of statistical parameters.
In this section, only one typical semi-numerical deterministic model, developed by
Ciavarella et al [75],are introduced.
2.3.1 Ciavarella et al?s model
Contact model studied by Ciaverella et al. [75] consists of one rigid rough surface
and one elastic half-space. Since only asperities are essential for the rough surface contact
22
according to the GW model, the rough surface is characterized by the vectors of asperity
coordinates x, y and z and radius of curvature R. Asperities are determined based on the
criterion that it should be the maxima comparing to the neighboring 4 sampling points. The
radius of curvature R can be determined by
Rx = (1 h
0
x
2)1:5
h00x
Ry = (1 h
0
y
2)1:5
h00y
R = (Rx +Ry)=2;
where h0i and h00i , i = x and y, are calculated by Eq. 1.1.
As the rigid rough surface approaches the elastic half-space, the contact areas grow
independently on the top surface of the elastic half-space. The shape of each contact area
can be determined by Hertzian spherical contact theory:
p(r) = p0
p
1 r2=a2 r2[0;a]
 = a2=R
a =
 3
4
 P R
E 
 1=3
: (2.22)
p0 is maximum contact pressure. a is contact area radius.  is the deformation of asperity.
 P is external load acting on asperity. For a typical Hertzian pressure distribution within the
contact area aj, the corresponding vertical displacement of asperity i, wi, may be written as
[31]
wi = 1R
j
(a2j r2=2); r aj (2.23)
wi = 1 R
j
h
(2a2j r2)asin
 aj
r
 
+ aj
q
r2 a2j
i
; r>aj (2.24)
23
where r is the distance between asperity i and asperity j: r = p(xi xj)2 + (yi yj)2.
Total displacement of asperity i due to all the Hertzian contact is
wi =
NcX
j=1;i6=j
 1
 Rj
h
(2a2j r2)asin
 aj
r
 
+ aj
q
r2 a2j
i 
+ 1R
i
a2i (2.25)
where Nc is the number of contacting asperities. Eq. 2.25 assumes that only one asperity,
itself, lies within the contact area of radius ai. This assumption is reasonable when contact
areas grow independently. However, the prediction of contact load and area at high load
conditions may not be accurate due to the fact of asperity coalescing. For those contacting
asperities, the following governing equation should be obeyed:
hi = d+wi; i = 1;:::;Nc
where hi is asperity height measured from the mean level. d is surface separation.
Because the contacting asperities are unknown in advance, an iterative scheme needs to
be adopted. Ciavarella et al. [75] introduced the following correction scheme
 ai = Ri2a
i
(hi d)
ai+1 = ai +  ai (2.26)
by di erentiating Eq. 2.22. The contact area ai are updated iteratively until  ai is su -
ciently small.
2.4 Results and discussion
In the above sections, we have brie y introduced some important models in the de-
velopment of statistical (GW, Nayak, BGT, ZC and CGP) and multi-scale (Archard, MB
and JS) models. In this section, statistical (GW, Nayak, asymptotic BGT and CGP) and
24
the simpli ed multi-scale, Eq. 2.21, models are compared to the  nite element model, the
solution of which is treated as "exact" and will be discussed later.
Surface 1 and surface 2, discussed in the Section 1.1.2, are used as the inputs of the
elastic rough surfaces contact models. Surfaces are compressed by the rigid  at. The essential
statistical parameters of both rough surfaces needed as inputs for all of the statistical models,
are listed in Table 1.1. m0, m2 and m4 are the average values of that along the x and
y directions. Contact load (P) contact area (A) are in the dimensionless forms: P =
P=(AnE ), A = A=An.
10
?3
10
?1
10
0
Dimensionless load P/(E
*
A
n
)
Dimensionless area A/A
n
 
 
GW
Nayak
BGT Asymptotic
CGP
Ciavarella et al model
The simplified multi?scale model
FEM
Surface 1
Figure 2.1: Contact area-load relation of surface 1 predicted by statistical, multi-scale and
FEM models.
Figs. 2.1 and 2.2 illustrate the relations betweenA andP predicted by the GW, Nayak,
asymptotic BGT, CGP, Ciavarella et al?s semi-numerical model, the simpli ed multi-scale
model and the FE model. Value of Bmax used in the simpli ed multi-scale model is the
average value of Bmax of each row of rough surface data.
25
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
Dimensionless load P/(E
*
A
n
)
Dimensionless area A/A
n
 
 
GW
Nayak
BGT Asymptotic
CGP
Ciavarella et al model
The simplified multi?scale model
FEM
Surface 2
Figure 2.2: Contact area-load relation of surface 2 predicted by statistical, multi-scale and
FEM models.
Fig. 2.1 shows the relations between A and P of (nearly) isotropic surface, e.g., surface
1. All the statistical models discussed above are based on the assumption that the rough
surface is isotropic. The asymptotic BGT model has a better accuracy than the GW and
Nayak models since it includes both the joint probability function and more generalized
elliptic contact model. The multi-scale model is very similar to the asymptotic BGT model.
As a results, the Nayak model has a better accuracy than the GW model. At the low load
condition, the CGP and GW models have the same prediction because the e ect of asperity
interaction is not severe. As the load P increases, the contact area A predicted by the
CGP models deviates from that predicted by the GW model. In CGP model, however, only
displacement of mean asperity level is used to show the redistribution of asperity height
caused by asperity interaction. Relative changes of asperity heights between neighboring
asperities are not included. Under the frame work of the statistical models, it is impossible to
26
include a more accurate asperity interaction model because asperity coordinates are unknown
and that is one of the main reasons to develop the semi-numerical deterministic models. Even
though Ciavarella et al?s semi-numerical deterministic model includes more accurate asperity
interaction, its accuracy is only better than that of the Nayak models and worse than that
of the BGT model, the simpli ed multi-scale model and the FE model. Two reasons can
explain the defects of Ciavarella et al?s semi-numerical deterministic model: (1) Surface 1 is
not perfectly isotropic. The spherical asperity contact model used in Ciavarella et al?s semi-
numerical model may not accurately describe the asperity contact than the more generalized
elliptic asperity contact model used in the asymptotic BGT model. (2) Asperity coalescing is
not included. It is surprising to see that both the simpli ed multi-scale and the asymptotic
BGT models has the the same and most accurate prediction of the relation between A and
P . This con rms that Archard?s type of structure of rough surface combined with sinusoidal
protuberances is reasonable to be used in predicting the A to P relation, although the
selection of Bmax may be important.
Fig. 2.2 illustrates the relations between A and P of the anisotropic surface, i.e.,
surface 2. Spectral moments, m4x and m4y, are not in the same order of magnitude, see
Table. 1.1, which causes surface 2 to show a strong anisotropic nature. As the results of
all the statistical models based upon the isotropic assumption, they may not have good
predictions. Ciavarella et al?s semi-numerical model uses the spherical asperity contact
model and cannot accurately capture the asperity contacts of the anisotropic rough surface
(a majority of them are elliptic contacts). The gap between the predictions of the simpli ed
multi-scale and FE models becomes smaller as the load P increases within the small load
range. The reason of the high accuracy of the simpli ed multi-scale model is because the
parameter Bmax is the average value of the Bmax of each row of the rough surface data.
This average process may more or less capture the strong e ect of anisotropic nature on the
relation between A and P .
27
10
?4
10
?3
10
?2
10
?2
10
?1
10
0
Dimensionless load P/(E
*
A
n
)
Dimensionless area A/A
n
 
 
JS ave(B
max
)
JS max(B
max
)
JS min( B
max
)
FEM
Surface 1
Figure 2.3: Contact area-load relation of surface 1 predicted by simpli ed multi-scale models
[34] with di erent de nition of Bmax.
Figs. 2.3 and 2.4 illustrate the relations between A and P predicted by the simpli ed
multi-scale and FEM models. Di erent de nitions of Bmax:
\ave(Bmax)", Bmax value used in Eq. 2.21 is the average of Bmax of each row of rough surface;
\max(Bmax)", Bmax value used in Eq. 2.21 is the maximum of Bmax of each row of rough
surface;
\min(Bmax)", Bmax value used in Eq. 2.21 is the minimum of Bmax of each row of rough
surface.
are used in the simpli ed multi-scale model in order to test the e ect of Bmax on the relations
between A and P .
Solutions of the FE model of both isotropic and anisotropic surfaces lie between two limit
solutions, i.e., min(Bmax) and max(Bmax). As the load P (or contact area A ) increases, the
FE solution transits from min(Bmax) curve to max(Bmax) curve. This may be explained by
28
10
?5
10
?4
10
?3
10
?2
10
?1
10
0
Dimensionless load P/(E
*
A
n
)
Dimensionless area A/A
n
 
 
JS ave(B
max
)
JS max(B
max
)
JS min( B )
FEM
Surface 2
max
Figure 2.4: Contact area-load relation of surface 2 predicted by simpli ed multi-scale models
[34] with di erent de nition of Bmax.
the following reasons: The de nition of B is the ratio of the amplitude  to the wavelength
 . After the Fourier transformation of the 3D isotropic rough surface heights, we may have
stacks of isotropic sinusoidal surfaces with di erent values of B. Assume we have two cross-
sections of 3D isotropic sinusoidal surfaces with di erent values of B, i.e., B1 and B2. The
wavelength,  , of both sinusoidal waves are the same. The total load, needed to  atten those
surfaces, is dependent on the B. Thus, the sinusoidal surfaces with the larger value of B
have larger resistance to the deformation. Consequently, at the beginning of rough surface
contact, i.e., low load conditions, a majority of deformation is due to the components of
sinusoidal surfaces with smaller value of B. That is why the relation between A and P with
min(Bmax) matches the solution of the FE model at the low load condition for both isotropic
and anisotropic cases. As the load P increases, more and more sinusoidal components have
been  attened. At the point of complete contact, the sinusoidal component with the largest
29
B, i.e., max(Bmax), is  nally  attened. That is why when the A = 1, prediction of P by
the simpli ed multi-scale model with max(Bmax) has a good agreement with that of the FE
model. The two limiting solutions, i.e., min(Bmax) and max(Bmax), can exactly predict the
asymptotic solutions at the early and complete contact when rough surfaces are isotropic,
i.e., surface 1. For the anisotropic surfaces, i.e., surface 2, small mismatches can be found.
This may due to the fact that the simpli ed multi-scale model is based on the isotropic
assumption.
2.5 Conclusion
Some important statistical, multi-scale and semi-numerical models are discussed in this
chapter. The statistical (the GW, Nayak, asymptotic BGT and CGP), the simpli ed multi-
scale, Eq. 2.21, models and Ciavarella et al?s model are compared to the  nite element
model. Conclusions drawn from the above comparison can be generalized as:
(1) Statistical models are only valid for extremely small load conditions which has already
been accepted. Numerical results shown in Figs. 2.1 and 2.2 cannot prove it because the
FEM solution is not accurate at low load conditions. The limited capability of the FEM
model will be discussed in the following chapter.
(2) Except for adjacent asperity interaction, asperity coalescing may have even larger e ect
on the relation between A and P . The di erence between Ciavarella et al.?s model with and
without asperity coalescing can be found in A errante et al?s paper [77]. Clearly, the accuracy
of the one including the e ect of asperity coalescing increases dramatically. Unfortunately,
A errante et al?s solution is not included in the above comparison because it is too di cult
to recreate within the scope of this work
(3) The e ect of the anisotropic nature of rough surfaces considerably in uences the accuracy
of all the elastic rough surface contact models, especially the statistical models.
30
Chapter 3
Finite element method
As the computational capabilities have become cheaper and more powerful, more and
more e orts have been placed on the application of deterministic numerical methods in
the simulations of rough surface contact. Due to the fact that it is not easy (sometimes
impossible) to obtain the accurate contact area data from the experimental results, numerical
solutions have always been treated as \exact"even though it de nitely deviates from real
solutions because of discretized and round o error. The  nite element method (FEM) and
the boundary element method (BEM) are two popular methods used in rough surfaces contact
simulation.
Comparing to the massive usage of the BEM in the research of rough surface contact, the
application of the FEM (especially 3 dimensional rough surfaces analysis) starts to appear
just at the beginning of this century. Komvopoulos and Choi [59] may be the  rst to applied
the FEM to rough surface contact. However, their  nite element (FE) model was restricted
to plane-strain problem. Rough surface is modeled by uniformly distributed cylindrical
asperities and it is not real rough surface. As the storage capacity and computation speed
increase, solving a million nodes model is no longer as much of a challenge. Thanks to the
parallel computation technique, the computational time is dramatically saved by running the
simulation in super computers with multiple cores. Hyun et al. [60, 61] and Pei et al. [62]
applied the FEM in elastic and elastic-plastic rough surface contact with nearly one million
nodes, respectively. Similar analysis have also been done by Sahoo and Ghosh [63], Walter
and Mitterer [64], Yastrebov et al [65], Bryant et al [66], Olshevskiy et al [67] and Megalingam
and Mayuram [68]. Additionally, Thompson [69] studied thermal contact resistance between
contacting rough surfaces using the commercial software ANSYSTM. Timzer [70] studied the
31
thermomechanical contact between randomly rough surfaces. A complete literature review
of the development and application of the FEM in rough surfaces contact is provided by
Thompson and Thompson [36].
In the current work, a FE model of rough surface contact is built and excuted using
the commercial software ANSYSTM. The FE modeling in ANSYSTMis achieved by using the
APDL language [71]. The advantage of using APDL language instead of the GUI operation
is that the input parameter can be easily adjusted and data and solutions of the FE model
can be easily accessed and exported.
3.1 Surface modeling
x
z
y
Lx
Ly
H
Node
Element
A
?x
?y
Figure 3.1: Schematic representation of uniformly meshed elastic solid with smooth surface
(only the mesh on the nominal contact surface are visible).
Generally, the rough surfaces used in FEM analysis can be divided into computer-
simulated surfaces [60, 62, 63, 70, 68] and real surfaces[64, 65, 66, 67, 69]. No matter what
kind of rough surfaces data are used, meshing the rough surfaces and its sub-surface is always
32
a challenging task. Currently, two methods have been widely used in the above references,
namely, moving nodes and direct surface generation.
Moving Nodes [60]
Assume we have a rectangular volume with a  at top surface of Lx Ly and a thickness,
H (see Fig. 3.1). This volume has been discretized using a uniform mesh. Node A is
locating on the boundary or inside the volume, with the coordinates (x;y;z). Coordinate z
is measured from bottom surface. In order to obtain the rough surface h(x;y) on the top
surface and to prevent the moving procedure from distorting the shape of the elements, we
move the node A by
 z = h(x;y)
 z
H
 a
only along the z direction. a is a constant de ned by users. H is the dimension measured
between the top and bottom surface. Nodes on the top surface will move by the amount of
h(x;y) and the coordinate change along the z direction drops exponentially to zero until the
changes reach the nodes on the bottom surface.
Direct Surface Generation
The idea of this method is to directly generate the rough surface through the sampling
point cloud, then the bottom-up technique [69] is used to create a volume for meshing. Four
neighboring sampling points can form a Coon?s patch [71] (those sampling points are not co-
planar). By using the Boolean operation, we can create the entire rough surface. However,
great care should be placed when using Boolean operation, especially in ANSYSTM, to avoid
failure.
Due to the fact that direct surface generation is not robust in ANSYSTM, we choose
the method of moving nodes as the strategy for generating the rough surface mesh in the
current work. Usually, meshing of a solid with a rough surface results in more than 200;000
nodes, and therefore moving all the nodes may cost considerable operation time during
pre-processing. In the current work, only nodes on the top surface are moving.
33
x
z
y
A
B
C
D
H
G
F
E
A1
B1
C1
D1
A2
C2
D2
B2
(xi, yj, H)
hij
?x
?y
Figure 3.2: Schematic representation of a badly distorted rough surface element.
The elastic solid, see Fig. 3.1, is meshed uniformly at intervals  x and  y along the
x and y directions.  x and  y are the resolutions of the pro lometer along the x and
y directions. Since only the nodes on top surface are moved, the value of  z(x;y) equals
h(x;y). If the amount of moving (hij) at node (xi;yj) is larger than the dimension of interval
 x and  y, the shape of the element may be badly distorted (see element A2B2C2D2EFGH
in Fig. 3.2). In order to avoid distorted elements, a rough surface is selected for which the
max(jhj), the maximum absolute value of the rough surface height, is less than the mesh
interval (min( x; y)). Fig. 3.3 and 3.4 illustrates the details of the mesh of surface 1 and
surface 2 when the procedure of the moving nodes is achieved in ANSYSTM. When max(jhj)
of the rough surface is larger than the min( x; y), all the nodes need to be moved, as was
done in [60], in order to avoid element distortion.
Size of the rigid  at is chose to be much larger than that of the nominal contact area of
rough surface in this study.
34
Figure 3.3: Details of surface mesh of surface 1.
3.2 Boundary condition and e ect of thickness
It is easy to achieve the assumptions namely, cyclic symmetry and in nite thickness,
declared in the beginning of this chapter in ANSYSTM.
For the former one, i.e., cyclic symmetry, we should couple all the degrees of freedom
(D.O.F) of nodes on the opposite side surfaces. Fig. 3.5 shows the top view of the in nite
rough surface consists of  nite rough surface cells. All the rough surface cells are the same
and only one of them has been modeled in ANSYSTM. Points A and B are two nodes on
the opposite side surfaces 1 and 2, respectively, and share the same coordinates on cyclic
symmetric boundary. During each iteration, ANSYSTMcalculates the nodal displacement at
A in the x, y and z directions, i.e., Ux(A);Uy(A), Uz(A) and assigns the same Ux, Uy and
Uz values to node B. Similarly, the values of Ux, Uy and Uz are enforced to be the same at
35
Figure 3.4: Details of surface mesh of surface 2.
nodes C and D. In the current work, only weak-coupling is applied in the FE model, i.e.,
Ux(A) = Ux(B) and Uy(C) = Uy(D).
Besides the cyclic symmetric boundary condition (Fig. 3.6(a)), other two boundary
conditions namely, the constrained side surfaces (Fig. 3.6(b)) and free side surfaces (Fig.
3.6(c)) are also studied in the current work. All the nodal displacements on side surfaces are
constrained in Fig. 3.6(b) and are free in Fig. 3.6(c).
The contact area to load (A  P ) and the contact area to surface separation (A  (d0) )
relations from the FE models with the above three boundary conditions are obtained in
ANSYSTMand compared in Figs. 3.7-3.10. It is surprising to  nd that the A to P and
A to (d0) relations of the three di erent boundary conditions are almost identical. This
phenomenon may result from the negligible displacements along the in-plane directions, x
and y, compared to that in the vertical direction during the rough surfaces contact.
36
cyclic symmetric boundary
rough surface cell
x
y
Ly
Lx
rough surface cell
rough surface cell
D
C
AB
D
C
B
A
BA
C
D
(1) Neighboring rough surface cells are sharing 
the cyclic symmetric boundary
(2) Gap between two neighboring rough 
surface cell is exaggerated. 
1
2
3
4
1
2
3
4
1
2
3
4
Figure 3.5: Schematic representation of cyclic symmetry and displacement coupling
The second assumption is in nite thickness and in other words: the contact results, at
least contact load, P , area, A , and surface separation, (d0) , relations, should be indepen-
dent of thickness, H. We tested the FE models with di erent values of thickness namely,
H = 0:5Lx, H = Lx and H = 5Lx. If we represent the thickness in the function of RMS
roughness  , we have: surface 1, H = 308:88 , H = 617:76 and H = 3088:80 , surface
2, H = 2041:6 , H = 4082:6 and H = 2041:6 . The corresponding contact area to load
relations are plotted in Figs. 3.11-3.14. We can conclude from Figs. 3.11-3.14 that the
contact load, P , area, A , and surface separation, (d0) , relations are independent of the
thickness H, at least, when thickness is larger than 0:5Lx.
3.3 Results and discussion
In the above section, we have shown that the A to P and A to (d0) relations predicted
by the FE models with three di erent boundary conditions have good agreement (see Figs.
37
Infinite elastic body
Rough surface cell
(a)
Load
(b)
(c)
Load
Load
Displacement 
Constraint
Figure 3.6: Schematic representation of boundary conditions: (a) cyclic symmetry, (b) con-
strained side surfaces and (c) free side surfaces
3.7-3.10). Additionally, the A to P and A to (d0) relations are independent of thickness,
H, when H > 0:5Lx (see Figs. 3.11-3.14). The A to P solution used in the comparison
in Chapter 2 is obtained by the FE models with free side surfaces boundary condition and
H = Lx.
Fig. 3.15 and 3.16 illustrate the contact area and pressure distribution of surface 1 and
surface 2 at di erent contact stages and are the solutions of the FE model with free side
surfaces boundary condition and H = Lx. At the low contact force (see Fig. 3.15(a-b) and
3.16(a-b)), contact spots of asperities grow independently. At the medium load condition (see
Fig. 3.15(c-d) and 3.16(c-d)), asperity coalescing becomes severe and independent contact
spots form into contact patches with irregular shapes. At the high load condition (see Fig.
3.15(e-f) and 3.16(e-f)), only the not in contact spots grow independently. For the nearly
isotropic surface, i.e., surface 1, shapes of the contact spots at the low load condition (a)
and the not in contact spots at high load condition (e) are close to circular and ellipsoid.
For surface 2, which shows a strong anisotropic nature, the shapes of contact spots at the
low load condition (a) and the not in contact spots at high load condition (e) are more like
long strips and ellipsoids.
38
The correlation between rough surface height distribution and pressure distribution can
be found in Figs. 3.17 and 3.18 when complete contact is nearly achieved. Because the
contact has not been reached complete, the correlation of surface 2 is not obvious.
0 0.005 0.01 0.015 0.02 0.025 0.03
0
0.2
0.4
0.6
0.8
1
Surface 1
Dimensionless contact load P/(E* A
n
)
Dimensionless contact area A/A
n
 
 
Free side surf
D.O.F coupling on side surf
Constrain normal D.O.F on side surf
Figure 3.7: Dimensionless contact area A to dimensionless load P relation of surface 1.
The FE models with di erent boundary conditions illustrated in Fig. 3.6.
3.4 Conclusion
Application of the  nite element method is discussed in this chapter. The following
conclusions may be drawn from the above discussion:
1. Pre-processing in the FEM is tedious and time consuming, especially during the meshing
procedure.
2. According to the contact area to load results with di erent boundary condition, the FE
model in this study is valid for periodic domain problem.
3. Global contact results (contact load, area and surface separation) relations are indepen-
dent of the thickness of the substrate of the rough surface, at least when H > 0:5Lx.
39
0 1 2 3 4 5 6
x 10
?3
0
0.2
0.4
0.6
0.8
1
Surface 2
Dimensionless contact load P/(E* A
n
)
Dimensionless contact area A/A
n
 
 
Free side surf
D.O.F coupling on side surf
Constrain normal D.O.F on side surf
Figure 3.8: Dimensionless contact area A to dimensionless load P relation of surface 2.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
0
0.2
0.4
0.6
0.8
1
Surface 1
Dimensionless surface separation d?/?
Dimensionless contact area A/A
n
 
 
Free side surf
D.O.F coupling on side surf
Constrain normal D.O.F on side surf
Figure 3.9: Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 1.
40
0 0.5 1 1.5 2 2.5 3
0
0.2
0.4
0.6
0.8
1
Surface 2
Dimensionless surface separation d?
Dimensionless contact area A/A
n
 
 
Free side surf
D.O.F coupling on side surf
Constrain normal D.O.F on side surf
/?
Figure 3.10: Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 2.
0 0.005 0.01 0.015 0.02 0.025 0.03
0
0.2
0.4
0.6
0.8
1
Surface 1
Dimensionless contact load P/(E* A
n
)
Dimensionless contact area A/A
n
 
 
H = L = 617.76
H = 5L = 3088.80
H = 0.5L = 308.88
?
?
?
Figure 3.11: Dimensionless contact area A to dimensionless load P relation of surface 1.
Value of thickness: H = 0:5Lx; Lx; 5Lx.
41
0 1 2 3 4 5 6
x 10
?3
0
0.2
0.4
0.6
0.8
1
Surface 2
Dimensionless contact area A/A
n
Dimensionless contact load P/(E* A
n
)
 
 
H = L = 617.76
H = 5L = 3088.80
H = 0.5L = 308.88
?
?
?
Figure 3.12: Dimensionless contact area A to dimensionless load P relation of surface 2.
Value of thickness: H = 0:5Lx; Lx; 5Lx.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
0
0.2
0.4
0.6
0.8
1
Surface 1
Dimensionless surface separation d?/?
Dimensionless contact area A/A
n
 
H = L = 617.76
H = 5L = 3088.80
H = 0.5L = 308.88
?
?
?
Figure 3.13: Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 1. Value of thickness: H = 0:5Lx; Lx; 5Lx.
42
0 0.5 1 1.5 2 2.5 3
0
0.2
0.4
0.6
0.8
1
Surface 2
Dimensionless surface separation d?/ ?
Dimensionless contact area A/A
n
 
H = L = 617.76
H = 5L = 3088.80
H = 0.5L = 308.88
?
?
?
Figure 3.14: Dimensionless contact area A to dimensionless surface separation (d0) relation
of surface 2. Value of thickness: H = 0:5Lx; Lx; 5Lx.
43
 
 
2 4 6 8 10 12 14 16
x 10
?6
2
4
6
8
10
12
14
16
x 10
?6
 
 
2 4 6 8 10 12 14 16
x 10
?6
2
4
6
8
10
12
14
16
x 10
?6
x[m]
y[m]
 
 
2 4 6 8 10 12 14 16
x 10
?6
2
4
6
8
10
12
14
16
x 10
?6
Contact
No Contact
(a) (b)
(e)
 
 
2 4 6 8 10 12 14 16
x 10
?6
0
1
2
3
4
5
6
7
 
 
2 4 6 8 10 12 14 16
x 10
?6
0
1
2
3
4
5
6
7
8
9
10
x[m]
 
 
2 4 6 8 10 12 14 16
x 10
?6
0
2
4
6
8
10
12
14
16
(c) (d)
(f)
[Gpa]
[Gpa]
[Gpa]
Figure 3.15: Contact area contours and contact pressure distribution of surface 1 when
contact load is light(a-b), medium(c-d) and high(e-f). Light load: P = 4:49 10 4, A =
0:0489. Medium load: P = 0:00218, A = 0:283. High load: P = 0:00438, A = 0:948.
44
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
1
2
3
4
5
6
7
8
9
x 10
?4
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
1
2
3
4
5
6
7
8
9
x 10
?4
x[m]
y[m]
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
1
2
3
4
5
6
7
8
9
x 10
?4
Contact No Contact
(e)
(a)
(c)
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
0
0.5
1
1.5
2
2.5
3
x[m]
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
[Gpa]
[Gpa]
[Gpa]
(b)
(d)
(f)
Figure 3.16: Contact area contours and contact pressure distributions of surface 2 when
contact load is light(a-b), medium(c-d) and high(e-f). Light load: P = 1:91 10 5, A =
0:0173. Medium load: P = 9:62 10 4, A = 0:419. High load: P = 0:00438, A = 0:894.
45
x[m]
y[m]
 
 
2 4 6 8 10 12 14 16
x 10
?6
2
4
6
8
10
12
14
16
x 10
?6
0.5
1
1.5
2
2.5
3
3.5
4
x 10
?8
[m]
(b)
x[m]
y[m]
 
 
2 4 6 8 10 12 14 16
x 10
?6
2
4
6
8
10
12
14
16
x 10
?6
0
2
4
6
8
10
12
14
16
(a)
Figure 3.17: Correlation of rough surface geometry (a) and pressure distribution (b) of
surface 1.
x[m]
y[m]
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
1
2
3
4
5
6
7
8
9
10
x 10
?4
?15
?10
?5
0
5
x 10
?7
[m]
x[m]
y[m]
 
 
1 2 3 4 5 6 7 8 9
x 10
?4
1
2
3
4
5
6
7
8
9
x 10
?4
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
(b)(a)
Figure 3.18: Correlation of rough surface geometry (a) and pressure distribution (b) of
surface 2.
46
Chapter 4
Boundary element method
Rough surface modeling is always the biggest problem when the FEM is chosen to
simulate the contact. The accuracy of the FEM can be directly a ected by the mesh quality.
Each potential contact geometrical entity, i.e., asperities, and its subsurface needs to be
discretized using a highly re ned mesh. Usually, for a small or medium size rough surface
contact problem, the ultimate number of nodes required to obtain an \accurate "result is on
the order of millions or even more [36].
The BEM, however, has a natural advantage in solving rough surface contact problem.
Based on the help of the in uence function, only the boundary of the problem domain needs
to be discretized, i.e., more nodes are saved to discretize the rough surface model with a
 ner mesh. Additionally, the solution of stress and displacement within the domain are
more accurate than that on the boundary [37]. While the BEM has long been criticized
for its dense, full rank, non-symmetrical coe cient matrix (symmetrical and sparse in the
FEM) created during the formulation, thus the storage of the coe cient matrix and solving of
the whole system becomes a challenging task for the personal computer (PC) with a limited
storage capacity or even the workstation. Looking back into the history of the development of
the BEM, there is a long period, during which only thousands of D.O.Fs can be solved by the
BEM due to the limited storage capacity of typical computational device [37, 38]. However,
as the source point is closer to the boundary of the domain, the solution may become less
accurate due to the singular integration. Thus the subsurface stress distribution, which is
close to the rough surface, may not be as accurate as expected. The  rst disadvantage, dense
coe cient matrix, can be overcome by using a fast multipole boundary element method [38].
47
Even though the BEM can handle more rough surface sampling points than the FEM,
it is always unfeasible to include an entire fully detailed sampled surface in any numerical
method. Generally, only a part of the surface data is collected by a pro lometer with
relatively  ne resolution down to the nano scale. The sampled area may result in millions
of sampling points and not to mention the entire rough surface. In order to overcome both
computation and storage capacity limits, we assume:
1. Rough surface contact pair, including rough surface and rigid  at, is cyclicly symmetric
along x and y(Fig. 3.5). Each rough surface cell in Fig. 3.5 represents the same sampled
rough surface which cyclicly distributes along the x and y directions.
2. Nominally elastic  at rough surface has an in nite depth.
Assumption 1 enables the analysis of rough surface contact to restrict within a small
sampled area with  ne resolution and surface data to reveal the multi-scale structure of
rough surfaces as much as possible. Additionally, assumption 1 dramatically decreases the
computation time and storage memory. Assumption 2 actually is a byproduct of assumption
1, since the in-plane dimensions of the sampled surface is negligible compared to the in-plane
dimensions of entire rough surface. The combination of assumptions 1 and 2 are equivalent
to the assumption of an elastic half-space.
In the rest of this chapter, applications of the BEM in modeling the rough surface contact
are discussed thoroughly. A special kind of the BEM is derived by adopting some essential
assumptions. Methods of solving the nonlinear boundary integral equation are discussed.
Additionally, fast algorithms for deformation integral are also introduced. Solutions, the
relations of contact area and load, of the BEM are compared to that of the FE model.
4.1 Boundary integral equation
Boundary integral equation of three dimensional elastostatic problem is [37]
ui(p) +
ZZ
 0
Tij(p;Q)uj(Q)dS(Q) =
ZZ
 0
Uij(p;Q)tj(Q)dS(Q); p;Q2 0 (4.1)
48
where u and t indicate the surface displacement and surface traction on boundary  0. Tij is
the traction kernel denotes the traction at point p which lies on a certain plane due to the
unit point load at point Q. Uij is the displacement kernel denotes the displacement of point
p due to the unit point load. Tij and Uij are derived from Kelvin fundamental solution of
an in nite solid. i, j = x, y and z.
z
y
Lx
Ly
x
Roughness
Elastic half-space
Rigid Flat
z
x
h(x,y)
Figure 4.1: Semi-in nite contact body and its rough boundary.
Let us consider a special case, in which contact occurs between one nominally  at rough
surface and one rigid  at, see Fig. 4.1. The body with the rough surface is assumed to be a
semi-in nite body with the boundary z = h(x;y) which is measured from mean level. The
rigid  at has a  nite dimension in the x (Lx) and y(Ly) direction, see Fig. 4.1. Consequently,
boundary  0 in Eq. 4.1 now indicates the rough surface which is the only boundary of the
semi-in nite contact body. Commonly, the height of the rough surface, h, is negligible when
comparing to the in-plane dimension of the contact body, e.g., Lx and Ly. Thus we can
further assume that z = 0 plane can approximate the boundary of semi-in nite contact
body.
De ne domain  as part of domain  0 which coincides with the in-plane region of rigid
 at: f(x;y)2 j0 x Lx;0 y Lyg. The traction components within the region  0  
is zero, i.e., tx(Q) = ty(Q) = tz(Q) = 0;Q2  0  . Because the contact is frictionless,
thus tx(Q) = ty(Q) = 0;Q2 . ux(Q) and uy(Q) are negligible comparing to the dominant
49
deformation uz(Q) within the contact region  . Based on above assumptions, Eq. 4.1 can
be rewritten as
uz(p) +
ZZ
 0
T33(p;Q)uz(Q)dS(Q) =
ZZ
 
U33(p;Q)t3(Q)dS(Q); p;Q2 (4.2)
T33 is the traction kernel indicating the value of traction in z direction at point p, i.e.,
t3(p), due to the unit displacement in z direction at point Q, i.e., u3(Q). Point p and Q are
two points inside the in nite body according to Kelvin fundamental solution. In our special
case, Eq. 4.2, p and Q are all on the boundary  . From a physical perspective, point Q
cannot \feel "any load due to the unit displacement at point p, because of lack of continuum
medium to transfer the traction. As a results of above explanation, T33 = 0 within the
boundary  and Eq. 4.2 can be further simpli ed as
uz(p) =
ZZ
 
U33(p;Q)t3(Q)dS(Q); p;Q2 (4.3)
Eq. 4.3 is the original equation which is wildly used in modeling the contact problem in
tribology, especially contact mechanics with smooth interfaces [32]. The main assumptions
used in deriving the Eq. 4.3, the reduced form of three dimensional boundary integral
equation of elastostatic problem, are:
(1) One contact solid with a rough boundary is assumed to be a semi-in nite body with the
simpli ed boundary z = 0;
(2) Rough surface displacements along x and y are negligible comparing to that in z, i.e.,
ux uzand uy  uz;
(3) The contact is frictionless.
50
4.2 Formulation of rough surfaces contact problem
x-y plane of coordinate system coincides with the mean level of rough surface. d is the
surface separation measured between rigid  at and mean level. The origin of the coordinate
system is placed at the bottom left corner of the  domain.  is discretized a into uniform
2D grid such that
f(x;y)2 jxi = (i 1) x;yj = (j 1) y;i = 1;:::;nx;j = 1;:::;nyg
with element size  x  y.  x and  y are also the sampling resolutions of rough surface
data in x and y directions. nx and ny are number of the sampling points along x and y
directions.
2a
2b
(xk, yl)
p=pkl
(xi, yj)
z
y
x
wij
Elastic half-space
?
Figure 4.2: Schematic representation of displacement due to uniform pressure distribution.
Assume a piece-wise constant contact pressure (same as traction t3) distribution p within
the domain  , see Fig. 4.2, i.e., the pressure p within a rectangular cell (shaded area) , 2a 2b
(same as  x  y) centered about point (xk;yl), has a constant value pkl.
Rewriting Eq. 4.3 in a discretized form with commonly used symbols in contact me-
chanics, we have
wij =
X
(xi;yj)2 
X
(xk;yl)2 
Kijkl pkl (4.4)
51
where
Kijkl =
ZZ
 c
U33(xi;yj;xk;yl)dxdy; (xk;yl)2 c
wij and pkl are the displacement and contact pressure along z direction at at boundary point
(xi;yj), respectively.  c is the domain of shaded cell centered about point (xk;yl). Kijkl is
called in uence coe cient which derived by Love [39] and expression of Kijkl can be found
in Appendix B.
De ne two sub-domains: Ic and Inc. Points within domain Ic are in contact with rigid
 at and are not in contact within Inc.
The discretized governing equation of modeling the contact between a nominally  at
rough surface and one rigid  at is
gij =
X
(xi;yj)2 
X
(xk;yl)2 
Kijkl pkl +d hij; (xi;yj)2 ; (xk;yl)2Ic (4.5)
where gij is the gap in z direction between rough surface and rigid  at at point (xi;yj). The
boundary conditions for Eq. 4.5 are
pij6= 0; gij = 0; 8(xi;yj)2Ic (4.6)
pij = 0; gij 0; 8(xi;yj)2Inc
4.3 Solver
Because the unknown in Eq. 4.5 is pressure p and the value of p within the non-contact
region Inc is zero according to the free-traction boundary condition 4.6. Eq. 4.5 only needs
to be solved within region Ic. Venner and Lubrecht [40] suggested to use the matrix D
Dms = Kijkl; m =ji kj+ 1;s =jj lj+ 1
52
to replace the fourth order tensor Kijkl due to the fact that in uence coe cients Kijkl is only
dependent on the relative distance between (xi;yj) and (xk;yl). Then
nx 1X
m=0
ny 1X
s=0
Dms pkl = hij d; (xi;yj);(xk;yl)2Ic (4.7)
Because domain Ic is unknown before solving Eq. 4.7, an iterative method needs to be
carried out and an initial guess of domain Ic, Ic(0), is needed. An easy way of predicting the
Ic(0) is by  nding the points (xi;yj), which is above the rigid  at.
Several numerical methods have been developed in the past 40 years in order to solve
Eq. 4.7, namely, matrix inversion [32, 41] Gauss-Seidel relaxation [40], Newton-Raphson
relaxation [40] and variational method [49, 50].
Matrix Inversion
The summation term in Eq. 4.7 can be treated as multiplication between the matrix A
and unknown pressure vector p:
p = A 1(h d)
Matrix A is the in uence coe cients assembly from the Kijkl. Pressure vector, p, can be
solved by inverting A. The disadvantage of the matrix inversion method is (1) unrealistic
memory spaces are needed to save matrix A, and (2) solving unknown pressure vector p
by the direct solver, such as Gaussian elimination, requires O(N3) operations because A is
dense and non-symmetric (here, N is the number of contacting points). Consequently, matrix
inversion can only be used in the small-scale contact problem with less than thousands of
points, even though it is robust. Matrix Inversion was a popular solver of rough surfaces
contact model during the 1980s [41, 42, 43, 44].
Gauss-Seidel Relaxation
53
During one iteration, say ith iteration, contact domain Ic(i) is  xed and Eq. 4.7 is a linear
equation which can be solved by Gauss-Seidel relaxation
 pij = ~pij +w1 ij (4.8)
where
 ij = hij d 
nx 1X
m=0
ny 1X
s=0
Dms ~pkl
where  pij and ~pij are the pressure value after and before the Gauss-Seidel relaxation at point
(xi, yj). D11 is the in uence coe cient when (i, j) coincides with (k, l). w1 is the under-
relaxation coe cient used to stabilize and in some cases accelerate the iteration process.
However, Eq. 4.8, is not stable because changing the value of pressure at (xi, yj) may
change the deformation of all the points within the domain Ic more or less [40]. Venner and
Lubrecht [40] suggested to use the following correction scheme:
 pij = ~pij + ij 14( i+1 j + i 1 j + i j+1 + i j 1) (4.9)
where
 ij = (hij d 
nx 1X
m=0
ny 1X
s=0
Dms ~pkl)=(D11 D12=2 D21=2)
for 2 dimensional relaxation in order to stabilize the iteration scheme. The Gauss-Seidel
relaxation scheme with Eq. 4.9 has been wildly used in both elastohydrodynamic lubrication
[45, 40] and dry rough surface contact[46, 47, 48] problems.
Newton-Raphson Relaxation
Rewriting Eq. 4.5 here
gij =
X
(xi;yj)2Ic
X
(xk;yl)2Ic
Dms pkl +d hij
54
Expanding the above equation in Taylor series and only keeping the  rst direvative term, we
have
gij + @gij@p
kl
 pkl = 0; (xi;yj);(xk;yl)2Ic (4.10)
where @gij@p
kl
= Dms. As the distance between (xi;yj) and (xk;yl) increases, value of Dms
decreases dramatically. This implies that those pressure changes far away from point (i,j)
may have little e ect on the equilibrium of Eq. 4.5. Thus only 5 points: (xi;yj), (xi+1;yj),
(xi 1;yj), (xi;yj+1), (xi;yj 1), are picked
D11 pij +D21 pi+1 j +D21 pi 1 j +D12 pi j+1 D12 pi j 1 = gij
to join the correction scheme. Line relaxation can be used to perform more robust iteration:
D11  pij +D21  pi+1 j +D21  pi 1 j = gij D12 ~pi j+1 D12  pi j 1 (4.11)
  pij and  ~pi j are the residual pressure at (xi;yj) after and before one complete Newton-
Raphson line iteration. Eq. 4.11 is a non-linear equation and an extra iteration is needed
to ensure the accuracy of   pij satis es the required convergence criterion. An even simpler
method can be achieved by only using 3 points along the same dominant direction, e.g.,
D11  pij +D21  pi+1 j +D21  pi 1 j = gij (4.12)
to avoid the nonlinearity. Pressure solution p is then updated by Eq. 4.8.
Variational Method
Kalker and van Randen [49] established the variational principle for the frictionless con-
tact problem and applied it to frictionless non-Hertzian contact problems. The unique pres-
sure solution causes the total complementary energy to take the minimum value. However,
variational method cannot guarantee the solution satisfying exactly the boundary condition
55
in Eq. 4.6 [50]. Application of the variational method can be found in rough surface contact
[50] and rough contact with coated layer [51].
4.4 Fast algorithm for deformation integral
Let us rewrite the deformation integral, Eq. 4.4, in the governing equation Eq. 4.5
wij =
nx 1X
m=0
ny 1X
s=0
Dms pkl; (xi;yj);(xk;yl)2 :
If iterative methods, e.g., Gauss-Seidel, Netwon-Raphson and variational methods, are used
in solving Eq. 4.5, then the above integral equation needs to be evaluated in each iteration
(see Eq. 4.8 and 4.12). For the current 3D rough surface contact, O(N2) operations are
required to compute the whole surface deformation. As more and more sampling points
are used in modeling the rough surface, computational time for each iteration increases
dramatically to an unrealistic value. Currently, the maximum number of sampling points
used in rough surface contact modeling with the application of BEM is 16385 16385 [48].
In order to speed up the evaluation of surface displacement, two fast algorithms, namely
Multi level multi integration (MLMI) and fast Fourier transformation (FFT), are developed
and wildly used in tribology.
Multi Level Multi Integration (MLMI)
Multi level multi integration (MLMI) method was developed by Brandt and Lubrecht
[52] based on the idea of the multi-grid method which was originally used in solving elasto-
hydrodynamic lubrication problem. Evaluation of wij is done on the coarsest level with the
minimum number of points. Pressure distribution pij on the  nest grid is transformed to the
coarsest level by a restriction technique. After the evaluation of wij on the coarsest level,
wij is transformed back to  nest grid by a prolongation technique. By the help of MLMI,
only O(NlogN) operations are required to solve the deformation integral. The application
of the MLMI method can be found in the area of elastohydrodynamic lubrication [45, 40]
56
and contact mechanics [46, 53, 57, 47, 48].
Fast Fourier Integration
The fast Fourier transformation method was  rst introduced by Stanley and Kato [50]
for solving deformation of rough surface contact. Johnson et al. [31], however, may be the
 rst to apply the Fourier transformation in waviness contact. In the Westergarrd solution
of sinusoidal surfaces contact [30], we can  nd the correlation between the pressure p and
surface displacement u when complete contact occurs:
p(x;y) = cos( x)cos( y)
w(x;y) = 2E0( 2 + 2)1=2 cos( x)cos( y)
 = 2 nL
x
;  = 2 mL
y
:
A similar conclusion can also be found in sinusoidal pro le (2D) contact. By decomposing
the known pressure distribution p into the frequency domain by FFT [58], we can obtain
the corresponding amplitudes of surface displacement. Then surface displacement wmn cor-
responding to the certain frequency pair (m;n) can be found by an inverse FFT. Then the
resultant surface displacement w can be evaluated by the superposition of wmn. Again,
O(N2) operations are reduced to O(NlogN). The FFT-based method developed by Stanley
and Kato, however, is strictly applicable to the periodic discrete domain only. Consequently,
this makes the FFT-based method ideal for solving contact problems for nominally  at rough
surfaces which can be modeled periodically, see Fig. 3.5. Rough surfaces in concentrated
contact problems, i.e.,  nite nominal contact areas, are non-periodic and application of this
FFT-based method in such cases may introduce error [53, 55].
To address this issue, Nogi and Kato [54] developed a di erent FFT-based method.
Their ideal is to transfer the deformation integral: w = D p to frequency domain. By
multiplying the Fourier transforms of D and p, we can obtain the value of wij through
57
inverse Fourier transformation
w = i t2( t2(D)  t2(p))
where i t2() and  t2() are 2 dimensional inverse FFT and FFT, respectively. Di erent
from the FFT-based method developed by Stanley and Kato, this method is only valid for
non-periodic domain, i.e., valid for concentrated contact problems only.
However, the FFT-based method developed by Nogi and Kato may bring errors at
boundary of nominal area [55]. Liu et al. [55] developed a FFT-based method (DC-FFT
method) and can help to avoid the errors at the borders of the domain.
4.5 Convergence criterion and surface separation
Di erent numerical methods, e.g., matrix inversion, Gauss-Seidel, Newton-Raphson and
variational methods, are discussed in details in Section. 4.3. As the ith iteration is  nished,
convergence of the predicted pressure  p(i), is checked by
jj p(i)  p(i 1)jj2=jj p(i)jj2 < p (4.13)
where jj( )jj2 is the 2-norm. In order to make sure that contact area, A, converges and has
similar accuracy, the contact area solution after each iteration is also checked via
jj A(i)  A(i)jj2=jj A(i)jj2 < a (4.14)
where  a = 1=(nx ny). The iteration process will be terminated until the above two criterions,
Eq. 4.13 and Eq. 4.14, is satis ed simultaneously.
As the total contact load, P, increases, the deformation of the surfaces becomes more
severe and this may bring numerical di culties to the iterative method. In order to gain a
convergence solution, the under-relaxation coe cient, w1, in Eq. 4.8 needs to be adaptively
58
adjusted by w1 = w1=a, 1:2 is selected as the value of a in this study, if solution is found to
be divergent. 0:1 is chosen as the initial value of w1. Also, in order to increase the convergent
rate, the load step is divided intoN substeps and the converged pressure solutionp(i) obtained
at substep i is used as the initial pressure guess at substep i+ 1.
According to the assumption of a semi-in nite body, if the surface separation d is mea-
sured from the mean level of the undeformed rough surface, the surface separation may
 nally reach negative in nity as contact load increases. It is possible, however, to de ned a
new surface separation, d0, measured from the mean level of deformed rough surface
d0 = d (
nxX
i=1
nyX
j=1
wij)=(nxny) (4.15)
which excludes the rigid body displacement of the mean level. wij is the displacement on
z = 0 plane along the z direction at (xi;yj). Evidently, d0 varies from an initial value of d
at  rst contact to 0 when complete contact is reached. A similar approach is considered by
Wilson et al. [56]
4.6 Results and discussion
In this study, the governing Eq. 4.7 is solved by two di erent methods. One uses the
Newton-Raphson iterative method and the FFT-based fast integral method by Liu et al. [54]
(NR+DC-FFT). Another one uses the variational method and the FFT-based fast integral
method by Stanley and Kato [50] (VM+FFTSK). The same dimensionless group are used
for the results in order to compare with the FE model: P = P=(AnE ), A = A=An.
Fig. 4.3 shows a good agreement of the contact load to area relation of surface 1
predicted by the VM+FFTSK and the FE model. The FE Model is proved earlier that
it is suitable for modeling the nominally  at rough surface contact with periodic domain.
NR+DC-FFT overestimates the contact load for a given contact area and the di erence
increases as contact area increases. A possible reason causing the solution of NR+DC-FFT
59
to be di erent from that of VM+FFTSK and the FE model can be found in Fig. 4.5.
When the DC-FFT is selected as the fast integral method, the contact problem is de ned
as shown in Fig. 4.5(a). As shown, a  nite rigid  at is contacting a semi-in nite body with
a rough surface. When the load P is applied on the rigid  at, deformation of rough surface
is dominant comparing to that of the domain out side the rough surface contact region, see
Fig. 4.5(c). Tension on the top surface of the semi-in nite body can be found and this may
cause more material to slip from the borders of the contact region in stead of occupying the
gap between the rough surface and the rigid  at. This tension force may dramatically cause
a drop of contact area. When the FFTSK is solved using the fast integral method, rough
surface contact pairs distributes periodically (cyclically), as shown in Fig. 4.5(b). As the
load P is applied on the rigid  at, the top surface of the semi-in nite body deforms uniformly,
see Fig. 4.5(d). According to mass conservation, surface points within the contact region
have no where to  ow except to occupy the gap between the rough surface and the rigid  at.
In Fig. 4.5(c), the tangential displacement of surface points is considerable due to surface
tension. Those tangential displacements may delay the formation of a larger contact area.
That is why the contact area A predicted by NR+DC-FFT is underestimated for a given
contact load P . However, the gap can be found between the solutions of surface 2 predicted
by VM+FFTSK and the FE model (Fig. 4.4). Di erence may be caused by anisotropic
nature of surface data.
4.7 Conclusion
In this chapter, application of the boundary element method is discussed in details. The
contact load to area relations obtained from the BEM and the FEM are in good agreement
if the same boundary conditions are chosen. For elastic rough surface contact, the boundary
element method has a great advantage in saving computation time. The rough surface
meshing process in the commercial FEM software, e.g., ANSYSTM, is tedious and a large
number of sampling points (> 400 400) is not allowed due to unrealistic computation time
60
0 0.005 0.01 0.015 0.02 0.025
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Dimensionless load P/(E
*
A
n
)
Dimensionless area A/A
n
 
Newton?Raphson+DC-FFT 
Variational method+FFT
FEM
Figure 4.3: Contact area to load relation of surface 1 predicted by the boundary element
method and the  nite element method.
and lack of storage capacity. The FEM, however, may have advantage when deformation
between rough surfaces is nonlinear, e.g., elastic-plastic, viscous-elastic and creep.
61
0 1 2 3 4 5 6 7
x 10
?3
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Dimensionless contact area A/A
n
Dimensionless contact load P/(E* A
n
)
 
 
Newton?Raphson+DC-FFT 
Variational method+FFT
FEM
Figure 4.4: Contact area to load relation of surface 2 predicted by the boundary element
method and the  nite element method.
P
P
Rigid Flat
Surface Tension
Rough Surface
 Rigid Flat
Rough Surface is periodically distributing
(a)
(b)
(c) (d)
P
P
Figure 4.5: Schematic representation of rough surfaces deformation corresponding to di er-
ent FFT-based fast integral method: NR+DC-FFT before (a) and after (c) deformation;
VM+FFTSK before (b) and after (d) deformation.
62
Chapter 5
Conclusion and future work
In the current thesis, developments and applications of statistical, multi-scale, semi-
numerical and deterministic numerical models are discussed. According to the conclusions
of each chapter, we can conclude that the BEM has the advantages in modeling the rough
surface contact among all the numerical models.
In the future work, current rough surface contact model built by the BEM will include
the e ect of friction for both full stick and partial slip conditions. Additionally, elastoplastic
and adhesion can also be easily embedded in the framework of BEM.
The two asymptotic solutions with min(Bmax) and max(Bmax) in the conclusion of
Chapter 2 can be used to curve- t the empirical solution of elastic rough contact for isotropic
surface. For anisotropic surface, the accurate model of anisotropic sinusoidal contact should
be derived.
63
Bibliography
[1] Goryacheva, I.G., 1998 "Contact Mechanics in Tribology," Kluwer, pp.11-12
[2] Thomas, T.R., 1982 "Rough Surfaces," Longman, pp.92
[3] Greenwood, J.A., and Williamson, J.B.P., 1966, "Contact of Nominally Flat Surfaces,"
Proc. R. Soc. London, Ser. A, 295, pp. 300-319.
[4] Chang, W.R., Etsion, I., and Bogy, D.B., 1987, "Elastic Plastic Model for the Contact
of Rough Surfaces," ASME J. Tribol., 109, pp. 257-262.
[5] Jackson, R.L., and Green, I., 2005, "A Finite Eleemnt Study of Elasto-Plastic Hemi-
spherical Contact Against a Rigid Flat," ASME J. Tribol., 127, pp. 343-354.
[6] Brizmer, V., Kligerman, Y., and Etsion, I., 2006, "The E ect of Contact Conditions
and Material Properties on the Elasticity Terminus of a Spherical Contact," J. Mech.
Mater. Struc., 1, pp. 865-879.
[7] Zhao, Y., Maletta, D.M., and Chang, L., 2000, "An Asperity Microcontact Model
Incorporating the Transition From Elastic Deformation to Fully Plastic Flow," ASME
J. Tribol., 122, pp. 86-93.
[8] Kogut, L., and Etsion, I., 2002, "Elastic-Plastic Contact Analysis of a Sphere and a
Rigid Flat," ASME J. Appl. Mech., 69(5), pp. 657-662
[9] Whitehouse, D.J., and Archard, J.F., 1970, "The Properties of Random Surfaces of
Signi cance in Their Contact," Proc. Roy. Soc. Lond A., 316, pp. 97-121.
[10] Onions, R.A., and Archard, J.F., 1973, "The Contact of Surfaces Having a Random
Structure," J. Phys. D: Appl. Phys., 6, pp. 289-304.
[11] Nayak, P.R., 1971, "Random Process Model of Rough Surfaces," AMSE J. Lub. Tech.,
93, pp. 398-407.
[12] Bush, A.W., Gibson, R.D., and Thomas, T.R., 1975, "The Elastic Contact of a Rough
Surface," Wear, 19, pp. 163-168.
[13] McCool, J.I., 1987, "Relating Pro le Instrument Measurements to the Functional Per-
formace of Rough Surfaces," ASME, J. Tribol., 109, pp. 264-270.
[14] Carbone, G., and Bottiglione, F., 2008, "Asperity Contact Theories: Do They Predict
Linearity Between Contact Area and Load?," J. Mech. Phys.Solids 56, pp. 2555C2572.
64
[15] Zhao, Y., and Chang, L., 2001, "A Model of Asperity Interactions in Elastic-Plastic
contact of Rough Surfaces," ASEM, J. Tribol., 123, pp. 857-864.
[16] Ciavarella, M., Greenwood, J.A., and Paggi, M., 2008, "Inclusion of "Interaction" in
the Greenwood and Williamson Contact Theory," Wear, 265, pp. 729-734.
[17] Nayak, P.R., 1973, "Random Process Model of Rough Surfaces in Plastic Contact,"
Wear, 26, pp. 305-333.
[18] Greenwood, J.A., 2007, "A Note on Nayak?s Third Paper," Wear, 262, pp. 225-227.
[19] Ma, X., Mahtijn, R., and Dik, S., 2010, "A Load Dependent Friction Model for Fully
Plastic Contact Conditions," Wear, 269, pp. 790-796.
[20] Sepehri, A., and Farhang, K., 2008, "On Elastic Interaction of Nominally Flat Rough
Surfaces," ASME J. Tribol., 130, pp. 011014-1
[21] Archard, J.F., 1957, "Elastic Deformation and the Laws of Friction," Proc. Roy. Soc.
A., 243, pp. 190-205.
[22] Mandelbrot, B.B., 1982, The Fractal Geometry of Nature, W.H.Freeman, New York.
[23] Majumdar, A., and Bhushan, B., 1991, "Fractal Model of Elastic-Plastic Contact Be-
tween Rough Surface," ASME J. Tribol., 113, pp. 1-11.
[24] Ausloos, M., and Berman, D.H., 1985, A multivariate Weierstrass-Mandelbrot function,
Proc. Roy. Soc. A., 400, pp. 331-350.
[25] Yan, W., and Komvopoulos, K., 1998, Contact Analysis of Elastic-Plastic Fractal Sur-
faces., J. Appl. Phys., 84(7), pp. 3617-3624.
[26] Kogut, L., and Jackson, R.L., 2006, "A Comparison of Contact Modeling Utilizing
Statistical and Fractal Approaches," ASME J. Tribol., 128, pp. 213-217.
[27] Majumdar, A., and Tien, C.L., 1990, "Fractal Characterization and Simulation of Rough
Surfaces," Wear, 136, pp. 313-327.
[28] Voss, R.F., 1985, "Fractals in nature: from characterization to simulation," Springer-
Verlag, Berlin, pp. 805.
[29] Jackson, R.L., and Streator, J.L., 2006, "A Multi-Scale Model for Contact Between
Rough Surfaces," Wear, 261, pp. 1337-1347.
[30] Westergaard, H.M., 1933, Bearing Pressures and Cracks, ASME J. Appl. Mech., 6, pp.
49-53.
[31] Johnson, K.L., Greenwood, J.A., and Higginson, J.G., The Contact of Elastic Regular
Wavy Surfaces, Int. J. Mech. Sci., 27(6), pp. 383-396.
[32] Johnson, K.L., 1985, Contact Mechanics, Cambridge Univeristy Press.
65
[33] Krithivasan, V., and Jackson, R. L., 2007 "An Analysis of Three-Dimensional Elasto-
Plastic Sinusoidal Contact," Tribol Lett, 27(1), pp. 31-43.
[34] Jackson, R.L., and Green, I., 2011, "On the Modeling of Elastic Contact between Rough
Surfaces," Tribol Trans, 54, pp. 300-314.
[35] Jackson, R.L., 2010, "An Analytical Solution to an Archard-Type Fractal Rough Surface
Contact Model," Tribol Trans, 53(4), pp. 543-553.
[36] Thomson, M.K., 2010, "Considerations for the Incorporation of Measured Surface in
Finite Element Models," SCANNING, 32, pp. 183-198.
[37] Becker, A.A., 1992, The Boundary Element Method in Engineering, McGrall-Hill, Lon-
don.
[38] Liu, Y.J., 2009, Fast Multipole Boundary Element Method - Theory and Applications
in Engineering, Cambridge University Press, Cambridge.
[39] Love, A.E.H., 1952, A Treatise on the Mathematical Theory of Elasticity, 4th Edn,
Cambridge University Press, Cambridge.
[40] Venner, C.H., 2000, Multilevel Methods in Lubrication, Elsevier, Amsterdam.
[41] Webster, M.N., and Sayles, R.S., 1986, A Numerical Model for the Elastic Frictionless
Contact of Real Rough Surfaces, ASME J. Tribol, 108, pp. 314-320
[42] Lei, W.T. and Cheng, H.S., 1985, Computer Simulation of Elatic Rough Surface Con-
tacts, Trans. ASLE, 1985, 28, pp. 172-180.
[43] Seabra, J., and Berthe, D., 1987, In uence of Surface Rougness and Waviness on the
Normal Pressure Distribution in the Hertzian Contact, ASME J. Tribol, 109, pp. 462-
470.
[44] West, M.A., Webster, M.N., and Sayles, R.S., 1987, A New Method for Rough Surface
Contact Analysis, Proc. I. Mech. E., C161/87, pp.945-955.
[45] Venner, C.H., 1991, Multigrid solution to the EHL line multigrid techniques, Ph.d
Thesis, University of Twente, Enschede.
[46] Lubrecht, A.A., and Ioannides, E., 1991, "Problem and the Associated Subsurface Stress
Field, Using Multilevel Techniques," J. Tribol, 113, pp. 128-133.
[47] Colin, F., 2001, "Comparison of FFT-MLMI for Elastic Deformation Calculations," J.
Tribol, 123, pp. 884-887.
[48] Sainsot, P., and Lubrecht, A.A., 2010, "E cient Solution of the Dry Contact of Rough
Surfaces: A Comparison of Fast Fourier Transform and Multigrid Methods," Proc.
IMechE Part J: J. Engineering Tribology., 225, pp. 441- 448.
66
[49] Kalker, J.J., and van Randen, Y., 1972, "A Minimum Principle for Frictionless Elastic
Contact with Application to Non-Herzian Half-Space Contact Problems," J. Eng. Math.,
6, pp. 193.
[50] Stanley, H.M., and Kato, T., 1997, "An FFT-Based Method for Rough Surface Con-
tact," ASME J. Tribol, 119, pp. 481-485.
[51] Bhushan, B., and Peng, W., 2001, "Contact Mechanics of Multiplayered Rough Sur-
faces," ASME Appl Mech Rev, 55(5), pp. 435-480.
[52] Brandt, A., and Lubrecht, A.A., 1990, "Multilevel Matrix Multiplication and Fast so-
lution of Integral Equations," J. Comp. Phys., 90, pp. 348-370.
[53] Polonsky, I.A., and Keer, L.M., 1999, "A Numerical Method for Solving Rough Contact
Problems Based on the Muti-Level Multi-Summation and Conjugate Gradient Tech-
niques," Wear, 231, pp. 206-219.
[54] Nogi, T., and Kato, T., 1997, In uence of A Hard Surface Layer on the Limit of Elastic
Contact:Part I. Anlysis Using A Real Surface Model, ASME J. Tribol., 119, pp. 493-500.
[55] Liu, S.B., Wang, Q., and Liu, G., 2000, "A Versatile Method of Discrete Convolution
and FFT(DC-FFT) for Contact Analyses," Wear, 243, pp. 101-111.
[56] Wilson, W.E., Angadi, S.V., and Jackson, R.L., 2010, "Surface Separation and Contact
Resistance Considering Sinusoidal Elastic-Plastic Multi-Scale Rough Surface Contact,"
Wear, 268, pp. 190-201.
[57] Polonsky, I.A., and Keer, L.M., 2000, "Fast Methods for Solving Rough Contact Prob-
lems: A Comparative Study," ASME J. Tribo., 122, pp. 36-41.
[58] Cooly, C.W., and tukey, J.W., 1965, "An Algorithm for the Machine Calculation of
Complex Fourier Series," Math. Comput., 19, pp. 297-301.
[59] Komvopoulos, K., and Choi, D.H., 1992, "Elastic Finite Element Analysis of Multi-
Asperity Contacts," ASME J. Tribo, 114, pp. 823-831.
[60] Hyun, S., Pei, L., Molinari, J.F., and Robbins, M., 2004, "Finite Element Analysis of
Contact Between Elastic Self-A ne Surfaces," Phy. Rev. E, 70, 026117.
[61] Hyun, S., and Robbins, M., 2007, "Elastic Contact Between Rough Surface: E ect of
Roughness at Large and Small Wavelengths," Tribology International, 40, pp. 1413-
1422.
[62] Pei, L., Hyun, S., Molinari, J.F., and Robbins, M., 2005, "Finite Element Modeling
of Elasto-Plastic Contact Between Rough Surfaces," J. Mech. Phys. Solids, 53, pp.
2385-2409.
[63] Sahoo, P., and Ghosh, N., 2007, "Finite Element Contact Analysis of Fractal Surfaces,"
J. Phys. D: Appl. Phys. 40, pp. 4245-4252.
67
[64] Walter, C., and Mitterer, C., 2009, "3D versus 2D Finite Element Simulation of the Ef-
fect of Surface Roughness on Nanoindentation of Hard Coatings," Surface and Coatings
Technology, 203, pp. 3286-3290.
[65] Yastrebov, V.A., Durand, J., Proudhon, H., and Cailletaud, G., 2011, "Rough Sur-
face Contact Analysis by Means of the Finite Element Method and of a New Reduced
Model," C.R.Mecanique, 339, pp. 473-490.
[66] Bryant, M.J., Evans, H.P., and Snidle, R.W., 2011, "Plastic Deformation in Rough
Surface Line Contacts-A Finite Element Study," Tribology International, 46(1), pp.
269-278.
[67] Olshevskiy, A., Yang, H.I., and Kim, C.W., 2011, "Finite Element Simulation of Inelas-
tic Contact for Arbitrarily Shaped Rough Bodies," Proc. IMechE. Part C:J. Mechanical
Engineering Science, 0, pp. 1-12.
[68] Megalingam, A., and Mayuram, M.M., 2012, "Comparative Contact Analysis Study of
Finite Element Method Based Deterministic, Simpli ed Multi-Asperity and Modi ed
Statistical Contact Models," ASME J. Tribo, 134, 014503.
[69] Thompson, M.K., 2007, A Multi-Scale Iterative Approach for Finite Element Modeling
of Thermal Contact Resistance, Ph.D. Thesis, Massachusetts Institute of Technology.
[70] Temizer, I., 2010, "Thermomechanical Contact Homogenization with Random Rough
Surfaces and Microscopic Contact Resistance," Tribology International, 44, pp. 114-124.
[71] ANSYS Help V13.0, ANSYS Inc
[72] Jamari, J., de Rooij, M.B., and Schipper, D.J., 2007, "Plastic Deterministic Contact of
Rough Surfaces," ASME J. Tribo, 129, pp. 957- 962.
[73] Pasaribu, H.R., and Schipper, D.J., 2005, "Application of a Deterministic Contact
Model to Analyze the Contact of a Rough Surface Against a Flat Layered Surface,"
ASME J. Tribo, 127, pp. 451- 455.
[74] Hariri, A., Zu, J.W., and Ben Mrad, R., 2006, "n-Point Asperity Model for Contact
Between Nominally Flat Surfaces," ASME J. Tribo, 128, pp. 505-514.
[75] Ciaverella, M., Del ne, V., and Demelio, G., 2006, "A "re-vitalized" Greenwood and
Williamson Model of Elastic Contact Between Fractal Surfaces," J. Mech. Phys. Solids,
54, pp. 2569-2591.
[76] Yeo, C.D., Katta, R.R., Lee, J.K., and Polycarpou, A.A., 2010, "E ect of Asperity Inter-
actions on the Rough Surface Elastic Contact Behavior: Hard Film on Soft Substrate,"
Tribology International, 43, pp. 1438-1448.
[77] A errante, L, Carbone, G., and Demelio, G., 2012, "Interacting and Coalescing Hertzian
Asperities: A New Multiasperity Contact Model," Wear, 278-279, pp. 28-33.
68
[78] Nayak, P.R., 1973, "Random Process Model of Rough Surfaces in Plastic Contact,"
Wear, 26, pp. 305-333.
[79] Eid, H., and Adams, G.G., 2007, "An Elastic-Plastic Finite Element Analysis of In-
teracting Asperities in Contact with a Rigid Flat," J. Phys. D: Appl. Phys., 38, pp.
2841-2847.
[80] Greenwood, J.A., 2007, "A Note on Nayak?s third paper," Wear, 262, pp. 225-227.
[81] Ma, X., de Rooij, M.B., and Schipper, D.J., 2010, "A Load Dependent Friction Model
for Fully Plastic Contact Conditions," Wear, 269, pp. 790-796.
69
Appendix A
In uence Coe cient Kijkl
Finding the stresses and displacements in an elastic half-space due to the point load
belonging to the Boussinesq problem [32]. Love [39] derived the vertical displacement of an
arbitrary surface point (xi;yj) due to the uniform pressure distribution p = pkl within the
rectangular area 2a 2b with the center (xk;yl), see Fig. 2.1.1.
 E wij
p = Kijkl =(x+a) ln
 (y + b) +f(y + b)2 + (x + a)2g1=2
(y b) +f(y b)2 + (x + a)2g1=2
 
(A.1)
(y +b) ln
 (x + a) +f(y + b)2 + (x + a)2g1=2
(x a) +f(y + b)2 + (x a)2g1=2
 
(x a) ln
 (y b) +f(y b)2 + (x a)2g1=2
(y + b) +f(y + b)2 + (x a)2g1=2
 
(y b) ln
 (x a) +f(y b)2 + (x a)2g1=2
(x + a) +f(y b)2 + (x + a)2g1=2
 
where
x =jxi xkj; y =jyj ylj
70