Friction ISSN 2223-7690
https://doi.org/10.1007/s40544-023-0785-z CN 10-1237/TH
RESEARCH ARTICLE
Approximate contact solutions for non-axisymmetric homogeneous
and power-law graded elastic bodies: A practical tool for design
engineers and tribologists
Valentin L. POPOV*, Qiang LI*, Emanuel WILLERT*
Institute of Mechanics, Technische Universität Berlin, Berlin 10623, Germany
Received: 17 February 2023 / Revised: 13 April 2023 / Accepted: 31 May 2023
© The author(s) 2023.
Abstract: In two recent papers, approximate solutions for compact non-axisymmetric contact problems of
homogeneous and power-law graded elastic bodies have been suggested, which provide explicit analytical
relations for the force–approach relation, the size and the shape of the contact area, as well as for the pressure
distribution therein. These solutions were derived for profiles, which only slightly deviate from the axisymmetric
shape. In the present paper, they undergo an extensive testing and validation by comparison of solutions with
a great variety of profile shapes with numerical solutions obtained by the fast Fourier transform (FFT)-assisted
boundary element method (BEM). Examples are given with quite significant deviations from axial symmetry
and show surprisingly good agreement with numerical solutions.
Keywords: normal contact; non-axisymmetric indenter; extremal principle; generalized method of dimensionality
reduction (MDR); functional elastic grading
1 Introduction
Let us take a look back at the history of contact
mechanics. In 1882, Hertz [1] solved the elastic
normal contact problem of parabolic (not necessarily
axisymmetric) bodies. This solution was expanded
in 1941 by Föppl [2, 3] and in 1942 by Schubert [4],
who found an analytical solution for contacts of
bodies with arbitrary dependence on the polar radius
but restricted to axisymmetric shapes (this theory
became widely known due to Ref. [5] by Sneddon).
Despite the simplicity and elegance of this solution,
the restriction to axisymmetric contacts was too
strong to open the way for practical engineering
applications.
An attempt to overcome the restriction of axial
symmetry was undertaken in 1990 by Barber and
Billings [6]. Their approach is based on Betti’s
reciprocity theorem, as suggested by Shield [7] in
1967 and the extremal principle found by Barber [8].
However, Barber and Billings [6] merely illustrated
their method by examples of contacts with “linear
profiles” (i.e., pyramids with polygonal cross sections),
since an analytical execution of their procedure is
possible only for this case. In Ref. [9] (also in the
correction (Ref. [10])), the extremal principle of
Barber [8] was applied to contacts of profiles, which
are not axially symmetrical, but slightly deviate from
the axial symmetry. In Ref. [9], several examples
were considered, which showed surprisingly good
agreement with more rigorous numerical solutions.
The accuracy of the approximation is comparable,
e.g., to the accuracy of the Cattaneo [11]–Mindlin [12]
approximation for tangential contacts and Fabrikant’s
approximation [13] for the pressure distribution under
a rigid flat punch of arbitrary shape, which are widely
used in contact mechanical applications. Moreover,
in Ref. [14], the analytical approximate solution in
Ref. [9] has been generalized for the application to
power-law graded materials.
* Corresponding authors: Valentin L. POPOV, E-mail: v.[email protected]; Qiang LI, E-mail: qiang.li@tu-berlin.de; Emanuel WILLERT,
E-mail: e.will[email protected]
2 Friction
| https://mc03.manuscriptcentral.com/friction
In the present paper, the solutions found in
Refs. [9, 14] are extensively validated by comparison
with “numerically exact” solutions obtained by the fast
Fourier transform (FFT)-assisted boundary element
method (BEM) for various convex profiles, that are,
in part, deviating strongly from the axisymmetric
shape. Thus, we suggest that the presented analytical
solutions provide a simple and reasonably accurate
method for solving general normal contact problems
with compact contact areas in engineering design and
tribology.
The remainder of the manuscript is organized as
follows: In Section 2, the theoretical foundations of
the approach used in Refs. [9, 14] are summarized.
In Section 3, the solution for non-axisymmetric
contacts of homogeneous materials found in Ref. [9]
is reproduced for the convenience of the reader.
After that, in Section 4, analytical examples for the
homogeneous solution are shown and compared to
rigorous numerical solutions in Section 5. Sections
6–8 repeat the structures of Sections 3–5, for the case
of power-law elastic grading, and a discussion of the
advantages, drawbacks, possible extensions, and scope
of application of the suggested method concludes the
manuscript.
2 Fundamentals
The first basic idea behind the present approach is to
apply Betti’s reciprocity theorem to contact problems.
This idea and its first applications have been described
by Shield [7], who showed that the normal force FN(A)
appearing due to the indentation of the profile f(x,y)
to a depth d (where A is the contact area in this state)
is given by Eq. (1):
*
N( ) (,)( (,))dd
A
FA pxydfxy xy
(1)
where the pressure distribution p*(x, y) is that under
a flat punch with the cross section shape A, which
is indented by a unit distance. Eq. (1) is an exact
statement, which, however, does not allow any
practical application, as neither the correct contact
area nor the pressure distribution under a flat-ended
punch with this cross section is known.
An important step for the practical application of
Eq. (1) was made by Barber [8] in 1974: He proved
that the correct contact area fulfilling the usual contact
conditions (the pressure inside the contact area is
positive, and there is no interpenetration; outside the
contact the pressure is zero, and the distance
between surfaces is positive), corresponds to the
maximum force at a given indentation depth. Thus,
Barber [8] formulated an extremal principle, which
can be used for finding the shape of the contact area,
either by rigorous variation of the contact boundary or
using an approximation in the sense of a Ritz ansatz.
Barber’s proof [8] is based on a harmonic function
representation of the elastic displacements. Let us
briefly demonstrate an alternative approach (which is
easier to generalize for non-homogeneous problems),
based on Mossakovski’s [15] ingenious idea of
understanding a general indentation problem as a
series of incremental flat punch indentations.
We can write Eq. (1) in the form of
N() () ()FA dkA GA (2)
with the flat punch contact stiffness
*
() (,)dd
A
kA p xy xy (3)
and a functional
*
( ) (,)(,)dd
A
GA p xyfxy x y (4)
As the difference between two normal contact
configurations with the contact regions A and A + dA,
following Mossakovski [15], can be thought of as an
incremental indentation dd by a flat punch with the
planform A, the contact stiffness (Eq. (3)) is universal.
Hence,
NNN N
ddd
() ()
ddd
FFF F
A
A
kA kA
d d Ad Ad
(5)
and therefore (the case dA = 0, once again, corresponds
only to the flat punch contact)
N0
F
A
(6)
which is the necessary condition for Barber’s maximum
principle [8]. That “proof”1 applies to any normal
1 Mathematical rigor would require to show that A, indeed, always
corresponds to a maximum of the force; but for reasons of space,
we will skip that step here.
Friction 3
www.Springer.com/journal/40544 | Friction
contact problem, which can be thought of as a
superposition of flat punch indentations, and especially
to normal contacts of (locally isotropic) layered or
functionally graded elastic materials.
Nonetheless, to constructively apply Barber’s
principle [8], the pressure distribution p* under a
flat-ended punch with an arbitrary cross section has
to be known. Barber and Billings [6] proposed to
use Fabrikant’s approximation [13] for this pressure
distribution. This is the last step, which closes the
procedure and reduces it to the solution of a variational
problem for the contact boundary line. Fabrikant’s
hypothesis [13] is that the stress distribution in polar
coordinates (r,
) is given in a good approximation
by Eq. (7):
*
*
022
0
2()
(, ) ( ())
()
Ea
pr r a
Jar
(7)
where ()a
is the equation for the contact boundary
in (r,
) in the contact plane, and J0 is the linear
moment
2
00()dJa
(8)
E* is the composed effective elastic modulus,
22
12
*
12
11
1
EE
E
(9)
where E1 and E2 are the Young’s moduli of the contacting
bodies, and ν1 and ν2 are their Poisson’s ratios.
The motivation for the ansatz (Eq. (7)) is not rigorous,
but robust and convincing. Eq. (7) is known to
be exact for elliptical punches with arbitrary
eccentricity (as was already found by Hertz [1]), and
it provides the correct asymptotic behavior in the
vicinity of the contact boundary if the boundary line
is smooth [15] (for punches with sharp corners, the
corner singularity of the pressure distribution is of
a higher order than the captured singularity by
Fabrikant’s approximation [13], as will be demonstrated
in Section 5).
Note that Fabrikant’s approximation [13]
(Eq. (7)) is still not yet completely defined, as for
its unambiguous definition the origin of the polar
coordinate system has to be chosen. Both Fabrikant
[13] and Barber [16] suggest to choose it at the
centroid of the area A. However, there are no cogent
theoretical reasons for this choice. The basic motivation
for the use of Fabrikant’s approximation [13] does
not depend on the exact position of this point. This
freedom can be used for choosing the position of the
polar center in the most convenient way to simplify
the calculations.
With the approximation (Eq. (7)), Eq. (1) becomes
*2()
N00 22
0
()( (,))dd
2
() ()
aadfrrr
E
FA Jar
(10)
The boundary of the contact area should now be
found by maximizing this functional. We will omit
the derivation of the solution of this variational
problem, which was shown in Ref. [9]. The main
results will, nonetheless, be stated in Section 3.
Before, however, let us consider inhomogeneous
materials with elastic grading of the form
0
() ( 1)
m
Ez Ez m
(11)
i.e., the elastic modulus shall vary with depth z
according to a power-law. The homogenous problem
is always included as the special case (m = 0). The
generalization of Eq. (7) for the indentation of a
power-law graded elastic half-space by a rigid flat
punch of arbitrary cross section has been given in
Ref. [14] and reads (the index “m” corresponds to
the inhomogeneous problem, while the index “0”
shall always denote a solution for the homogenous
problem)
1/2
*2
*
2
2
(, ) 1 ( ())
()
m
m
m
m
Er
pr ra
Ja
(12)
with the effective “modulus” [17]
2
*0
2
11 cos
22
1
33
sin 22 2
11
1
m
E
mm
Em C
mm
C
m
m
(13)
where Γ denotes the Gamma function, and the nonlinear
4 Friction
| https://mc03.manuscriptcentral.com/friction
moment of the contact area is
21
0() d
m
m
Ja
(14)
Note that for the contact of two elastic materials,
a summation of the inverse effective “moduli”—in the
spirit of Eq. (9)—is possible, if both materials have the
same exponent m of elastic grading [18].
Hence, the generalization of Eq. (10) for power-law
graded media is given by Eq. (15):
*2()
N1
00 122
2
2((,))dd
()
() (() )
a
m
m
m
m
Edfr rr
FA Jaar
(15)
Again, to obtain the approximate contact solution, the
variational problem resulting from Barber’s extremal
principle [8] needs to be solved. The complete procedure
is detailed in Ref. [14], and its most important results
will be stated briefly in Section 6.
3 Approximate solution for normal contact
of elastically homogeneous and slightly
non-axisymmetric profiles
Let us consider an “arbitrary” profile
(, )zfr
(16)
underlying the following restrictions:
1) The profile deviates only weakly from an
axisymmetric one.
2) f(r,
) is a monotonously increasing function
of the polar radius r.
Under these conditions, the following simple
analytical procedure can be applied, as shown in
Ref. [9]:
I. The origin of the polar coordinate system is placed
at the lowest point of profile, so that
(0, ) 0f
(17)
II. In the second step, an “equivalent axisymmetric
profile” is determined, which is simply the profile,
averaged over the angles.
2
0
1
() (, )d
2
fr fr
(18)
III. The profile (Eq. (16)) is decomposed into an
axisymmetric part and the deviation
(, ) (, ) ()fr fr fr
(19)
IV. With the equivalent axisymmetric profile
(Eq. (18)), the usual solution procedure of the method
of dimensionality reduction (MDR) is applied [19]. In
particular, the transformed profile is determined as
0022
()
() d
arf r
Ga r
ar
(20)
and
0
0022
d()
() d
d
a
Gfr
ga a r
aar
(21)
V. The relation between the indentation depth and
the effective contact radius a is determined by Eq. (22):
0()dga (22)
VI. The normal force is given by Eq. (23):
*
N0
2( ())FEdaGa (23)
VII. The true non-axisymmetric contact area is given
by Eq. (24):
() ()aaa
(24)
where ()a
is determined as
00
0
(,) (,)
() ()
Ga aga
aag a
(25)
with the prime denoting the first derivative with
respect to the (first or only) functional argument, and
0022
(, )
(, ) d
arfr
Ga r
ar
(26)
and
0
0022
() (, )
(, ) d
a
Gfr
g
aar
aar
(27)
VIII. Finally, the pressure distribution can be
calculated from the integral
Friction 5
www.Springer.com/journal/40544 | Friction
**
**
() *
0
*
** 2 2
() ( )
(, ) d ()
()
()
a
r
aga
Ea
pr a
a
aa r
(28)
where the effective contact radius must be understood
as a function of (),a
as given implicitly by Eqs. (24)
and (25). The upper star on the contact quantities under
the integral in Eq. (28) denote the values during the
indentation procedure, as the pressure distribution
is determined by the integral over all incremental
indentation steps, i.e., from an indentation depth
d* = 0 until the final value d* = d [9].
VIIIa. A much more convenient way to determine
the pressure distribution is scaling the axisymmetric
solution for the contact pressure (axi)
p
*
(axi) 0
22
()d
a
r
g
xx
rE
paxr
(29)
to the actual contact area, i.e.,
(axi)
(, ) ()
r
pr p a
(30)
While this will exactly give the same result as the
evaluation of the more complex Eq. (28) in the case
of power-law indenters, as will be shown below, it
will most probably always provide a satisfactory
approximate solution for the real contact pressure
distribution.
Eq. (30) is a generalization of Fabrikant’s ansatz [13]
for flat-ended punches. In the latter, the pressure
distribution under an arbitrary flat punch is equal to
the one under an axisymmetric punch, but “rescaled”
to the true shape of the contact area. Similarly, Eq. (30)
states that the contact pressure is equal to the one
under an “equivalent axisymmetric indenter”, but
rescaled to the true contact area.
In the present paper, it will be shown that Eq. (30)
is the exact approximation for power-law profiles
(with arbitrary cross section), independently of the
exponent of the power-law. Although it is not proved
in a general case, the independence of the exponent
gives hope that it will be a good approximation for
arbitrary profiles. In that regard, one has to also bear
in mind that the complete procedure is approximate
(albeit yielding very good results, as will be
demonstrated in Sections 5 and 8); therefore, in most
cases, probably not much is gained by evaluating the
seemingly more rigorous Eq. (28), compared to the
scaling idea expressed in Eq. (30).
The whole procedure can be summarized as follows.
First, an axisymmetric profile is produced by averaging
the given profile over the polar angle. After that,
the contact problem is solved for this equivalent
axisymmetric profile. Finally, the true contact area
is found by Eqs. (24) and (25), and the pressure
distribution is calculated from either Eq. (28) or the
more convenient Eq. (30).
4 Analytical examples for homogeneous
materials
4.1 Solution for power-law shapes
Consider a profile having the form
(, ) ()
k
fr r
(31)
which means that all vertical cross sections are
self-similar, differing only by a scaling factor. For
all profiles, which can be written as a product of a
radial and an angular function, the decomposition
(Eq. (19)) is very simple, as it affects only the angular
factor.
()
(, ) ( () )
k
k
fr r
fr r
(32)
with the average value of ()
2
0
1()d
2
(33)
For the MDR-transformed profile, we get
1
00
00
0
() () 1
() ()
1
2
() 1
22
k
k
a
Ga kk
ga ka
k
kk
(34)
Accordingly, for the deviations (Eqs. (26) and (27)),
6 Friction
| https://mc03.manuscriptcentral.com/friction
we have
1
00
00
(, ) () ()
1
(, ) () ()
k
k
a
Ga kk
ga ka
(35)
and the boundary of the contact area is given by
Eq. (36):
22()
() 1 (1)(1)
kk
aa
kk kk
(36)
where the effective contact radius is determined by
Eq. (37):
1/
0()
k
d
ak
(37)
Note that Eq. (36) differs from the corresponding
Eq. (67) in Ref. [9]. The correct one is Eq. [36], while
Eq. [67] in Ref. [9] contains an error, which was
corrected in Ref. [10].
The normal force, according to Eq. (23), equals to
11/
*
N0
2()
1
kk
k
k
FEdk
k
(38)
Hence, for the average contact pressure p
1
*1/
0
2
2()
1
kk
Nk
FkE
pdk
k
a
(39)
For the pressure distribution in the contact area
normalized by the average pressure, we obtain
Eq. (40) according to Eq. (28).
1
1
22
d
(, ) (1)
2
1
()
k
pr k
p
r
a
(40)
It is the same result as that for an axisymmetric contact
(page 78, Ref. [18]), scaled to the non-axisymmetric
contact area ()a
, as expressed in Eq. (30). However,
this phenomenon is only possible due to the self-
similarity property
const.
() ()
aa
aa
(41)
Therefore, it is not generally correct, at least not in the
rigorous (albeit asymptotic) sense, as shown above for
power-law (and hence self-similar) profiles.
4.2 Linear (conical) profiles
In the special case (k = 1), considering that κ0(1) = π/2,
Eqs. (37), (38), (36), and (40) can be simplified as
*2
N
*
2
2
253()
() 22
()
(, ) arcosh (0 ())
2
d
a
Ed
F
d
a
Ea
pr r a
r
(42)
with the area hyperbolic cosine function arcosh(…).
4.3 Parabolic profiles
In the special case k = 2, considering that κ0(2) = 2,
Eqs. (37), (38), (36), and (40) become
3
*
N
2
*
2
2
4
32
52()
() 233
2
(, ) 2 1 ( ())
()
d
a
d
FE
d
a
dr
pr E r a
a
(43)
5 BEM results for elastically homogeneous
bodies with complex shapes
Now, we shall compare the asymptotic analytical
results obtained in the previous section with rigorous
numerical simulations of the corresponding problems,
done with the BEM for an elastic half-space [20],
accelerated by the FFT. All problems are solved under
displacement-controlled conditions, i.e., the indentation
depth shall be prescribed.
Let us introduce two normalized measurements
to quantify the error made by the analytical, but
approximate solution, compared to the numerically
Friction 7
www.Springer.com/journal/40544 | Friction
exact simulation, namely the normalized error of the
total normal force
(th) (sim)
NN
(sim)
N
FF
eF
(44)
and the normalized standard deviation for the contact
boundary
(th) (sim) 2
1
(th)
1()
()
N
ii
i
eN
a
a
(45)
where the superscripts (th) and (sim) represent
theoretical and simulated values. (th)
a is the average
value of (th).a N is the number of discretization points
along the contact boundary, and the theoretical effective
contact radius is determined from the indentation
depth via Eq. (37).
5.1 Indenters with cross sections of regular
polygons with sharp corners
First, we consider regular n-polygons given by Eq. (46):
0
2
( ) cos ( 0,1,..., 1)
22
jj n
n
jj
nn nn
(46)
Hence,
0sin
n
n
(47)
In Fig. 1, the theoretically predicted contact area
(the red line) is compared to the area resulting from
a BEM simulation (the grey filled shape), for the
indentation by conical indenters (k = 1) with n = 3, 4,
5, and 6. In all cases, the theoretical and numerical
results are normalized by the same value, specifically
the theoretical effective contact radius, which can be
determined immediately from the indentation depth,
based on the first equation of Eq. (42). Apparently,
the approximate solution works better for larger
values of n. This is to be expected, as the limit n → ∞
corresponds to the indentation by a perfect cone,
for which the analytical solution is, of course, exact.
Moreover, for all values, except n = 3, the agreement
between the analytical and numerical results is very
satisfactory.
In Fig. 2, the error measurements introduced in
Eqs. (44) and (45) are shown for conical (k = 1) and
parabolic (k = 2) indenters, as well as for different
values of n. It can be seen that the error of the
approximate solution rapidly decreases with the
increasing n.
In Figs. 3 and 4, the pressure distributions in contacts
with pyramidal and parabolic indenters with square
Fig. 1 Comparison of contact areas for pyramid indenters (k = 1)
with regular n-polygon cross sections. The contact size is normalized
by the theoretical value of the effective contact radius, as shown
in the first equation of Eq. (42).
Fig. 2 Normalized errors of calculated normal force and average
deviation of contact boundary between approximate solution and
numerically exact BEM results for conical (k = 1) and parabolic
(k = 2) indenters.
8 Friction
| https://mc03.manuscriptcentral.com/friction
Fig. 3 Pressure distributions in contact of pyramid with square
cross section (k = 1 and n = 4). (a) Pressure at cross sections of
a, b, c, and d corresponding to 0°, 30°, 40°, and 45°, respectively,
as shown in subplots (b) and (c). The symbols are the numerical
results. (b) Contour line diagram of pressure distribution according
to BEM simulation. (c) Contour line diagram of pressure distribution
according to approximate analytical solution.
cross sections are presented—both from the BEM
simulation and the approximate analytical solution.
The numerical results for the pressure distribution
converge rapidly with the mesh size, except for
the sharp edge (i.e., the cross section at 45°), which
leads to the appearance of a weak (logarithmic) stress
singularity. Accordingly, the theoretical asymptotic
prediction is in good agreement with the BEM results
almost everywhere, with the exception of the immediate
vicinity of the sharp edges.
5.2 Indenters with cross sections of regular polygons
with rounded corners
As shown in Section 5.1, the analytical approximate
solution gives very good results, except for the vicinity
of sharp edges of the indenter profile. To better
understand the role of sharp edges, let us consider
the same regular power-law polygonal indenters as
before, but with rounded corners.
Figure 5 shows the results for the same cross sections,
Fig. 4 Pressure distributions of parabolic indenters with square
cross section (k = 2 and n = 4). (a) Pressure at cross sections of
a, b, c, and d corresponding to 0°, 30°, 40,° and 45°, respectively,
as shown in subplots (b) and (c). Symbols are numerical results.
(b) Contour line diagram of pressure distribution according to
BEM simulation. (c) Contour line diagram of pressure distribution
according to approximate analytical solution.
Fig. 5 Comparison of contact areas for pyramid indenters (k = 1),
whose cross sections are polygons with rounded corners. The contact
boundary is normalized by the theoretical value of the effective
contact radius (red: approximate solution, grey: BEM simulation).
Friction 9
www.Springer.com/journal/40544 | Friction
as shown in Fig. 1, but with rounded corners. Figures 6
and 7 show the corresponding results for the pressure
distribution, along the same directions, as shown
in Figs. 3 and 4, respectively, for polygon cross
sections with sharp corners. The pressure distribution
under an indenter with rounded corners has no
singularity, and the simulation results rapidly
converge everywhere with the mesh size. However, a
significant maximum of pressure along the diagonal
direction and correspondingly significant deviations
between numerical and analytical results still remain.
5.3 Power-law indenters with irregular cross
sections
One might argue that the quality of the asymptotic
solution demonstrated in Sections 5.1 and 5.2 is only
due to the fact that regular polygons are still relatively
similar to axisymmetric indenter cross sections (at
least, they exhibit discrete rotational symmetry). To
check this hypothesis, we randomly generated several
irregular horizontal indenter cross sections and
Fig. 6 Pressure distributions at different cross sections of pyramid
with square cross section (k = 1 and n = 4). (a) Pressure at cross
sections of a, b, and c corresponding to 0°, 30°, and 45°, respectively,
as shown in subplot. The symbols are the numerical results.
(b) Contour line diagram of pressure distribution according to BEM
simulation. (c) Contour line diagram of pressure distribution
according to approximate analytical solution.
Fig. 7 Pressure distributions at three cross sections of parabolic
indenters with square cross section (k = 2 and n = 4). (a) Pressure
at cross sections of a, b, and c corresponding to 0°, 30°, and 45°,
respectively, as shown in the subplot. The symbols are the
numerical results. (b) Contour line diagram of pressure distribution
according to BEM simulation. (c) Contour line diagram of pressure
distribution according to approximate analytical solution.
compared the corresponding normal contact solutions
based on the approximate analytical solution and the
BEM simulation. However, self-similarity was retained,
as the radial dependence of the indenter profile was
still chosen in the form of a power-law.
In Fig. 8, the results are shown for the error
measurements introduced in Eqs. (36) and (37), for
conical (k = 1) and parabolic (k = 2) indenters, for
160 irregular cross sections that were generated
randomly. The diagrams give the normalized error
of the approximate solution as a function of the
standard deviation of the cross section from perfect
rotational symmetry. Naturally, the stronger the profiles
deviate from axial symmetry, the larger the error of
the asymptotic solution. However, no qualitative
“misjudgements” of the asymptotic solution can be
detected.
Interestingly, the approximate solution provides
consistently better results for the parabolic indenters
(k = 2) than for the conical ones (k = 1).
In Fig. 9, a detailed comparison of calculated contact
areas is shown for six selected irregular shapes, marked
with the red triangles 1–6 in Fig. 8.
10 Friction
| https://mc03.manuscriptcentral.com/friction
Fig. 8 Normalized errors of calculated normal force and average
deviation of contact boundary between approximate solution and
numerically exact BEM results, depending on standard deviation
of irregular angular profile function ()
. 160 random profile
cross sections were realized. The blue stars correspond to the
pyramid indenters (k = 1), and the orange circles correspond to
the parabolic indenters (k = 2). The examples marked by 1–6 are
shown in Fig. 9.
5.4 Indenters with three-dimensional irregular
shapes
For the application of the approximate solution in
various tribological contexts, it is interesting to know
whether the solution can also be used satisfactorily in
the case of a three-dimensionally irregular indenter
shape, e.g., to analyze local features of rough surfaces.
So, let us turn our attention away from power-law
shapes, and consider a general three-dimensional
indenter profile.
In Fig. 10, a comparison is shown between the
approximate solution and the BEM simulation, for
the relations between the macroscopic contact quantities
(normal force, indentation depth, and contact area) in
the indentation of an elastic half-space by the general
indenter shape, as shown in the subplot of Fig. 10(a). All
relations are in properly normalized variables, and it is
obvious that the agreement between the approximate
and the rigorous numerical solutions is very good.
Figure 11 gives the corresponding comparison for
the pressure distributions at three instances of the
indentation procedure, as marked in Fig. 10(b).
Fig. 9 Comparison of contact areas for six pyramid indenters (k = 1) with irregular cross sections. The orange curves are the
approximate solution. The contact boundary is normalized by the theoretical value of the effective contact radius (red: approximate
solution, grey: BEM simulation).
Friction 11
www.Springer.com/journal/40544 | Friction
Fig. 10 (a) Normal force–indentation depth and (b) contact
area–indentation depth dependencies for indentation by general
three-dimensional indenter. The subplot of (a) is the profile of the
indenter. The red stars are the approximate solution, and the blue
cycles are the BEM simulation with the 1024 × 1024 meshing
grid. Comparisons for the shape of the contact area at three
different indentation depths are shown in (b).
6 Approximate solution for normal contact
of elastically inhomogeneous, slightly
non-axisymmetric profiles
Let us now turn our attention to the mathematically
slightly more general—yet, physically, probably more
specific—problem of the indentation of a power-law
graded elastic half-space by a slightly non-axisymmetric
indenter.
Compared to the algorithm described in Section 3,
the steps I–III, i.e., the definition of the profile and its
separation into axisymmetric and non-axisymmetric
components, remain the same. However, the contact
solution has to be generalized to account for the elastic
grading, as shown in Eq. (11).
Specifically, the generalized auxiliary profiles for
the inhomogeneous problem are given by Eq. (48):
0221
0221
()d
() ()
d
() d
(, )d
(, ) ()
(, )
a
mm
m
m
a
mm
m
m
rf r r
Ga ar
G
ga arfr r
Ga ar
G
ga a
(48)
The relations between the indentation depth, effective
contact radius, and normal force are determined by
the axisymmetric solution (Eq. (49)):
1
*
N
()
2()
1
m
m
m
mm
daga
da
FE Ga
m
(49)
and the deviation of the contact boundary from circular
form can be calculated from Eq. (50).
(1 ) ( , ) ( , )
() () ()
mm
mm
mGa aga
amg a a g a
(50)
Finally, the pressure distribution is found as
*
*
*221
*()
1
**
*
**
1
()
(() )
(, ) () d d()
d
m
a
m
m
r
a
a
ar
E
pr ad
a
aa
(51)
where the upper star on variables under the integral
denotes contact quantities during the indentation
process from d* = 0 to the final value d* = d.
Once again, a more convenient way of determining
the pressure distribution is to scale the respective
axisymmetric result
*
(axi)
221
()d
()
a
mm
rm
Egxx
r
paxr
(52)
to the real asymmetric contact area, as expressed in
Eq. (30). As in the homogeneous case, this will provide
exactly the same solution as the evaluation of Eq. (51)
12 Friction
| https://mc03.manuscriptcentral.com/friction
in the case of power-law indenters with arbitrary
exponents and cross sections—as will be shown
below—and will therefore probably be in close
agreement with that of Eq. (51) in any even more
general case.
7 General approximate inhomogeneous
solution for power-law shapes
As an analytical example, we will consider the
power-law indenter with arbitrary self-similar cross
sections again, as shown in Eq. (31). The auxiliary
MDR profiles are given by Eqs. (53) and (54):
1
() () 1
() ()
1
2
1
() 21
222
km
mm
km
mm
m
a
Ga k
km
ga ka
k
m
kkm
(53)
1
(, ) () ()
1
(, ) () ()
km
mm
km
mm
a
Ga k
km
ga ka
(54)
Hence, the boundary of the contact area equals to
22()
() 1 11
kk
aa
kk m kk m
(55)
where the effective contact radius is determined by
Eq. (56):
1/
()
k
m
d
ak
(56)
Moreover, we obtain for the normal force
*1
N
2
(1)(1)
m
m
k
FEda
km m
(57)
and the pressure distribution in the contact area,
normalized by the average pressure
1
1
221
d
(, ) (1)(1)
2()
k
m
pr km m
p
(58)
with the same definition of ρ, as shown in Eq. (40).
Once again, Eq. (58) is the same result as that for an
axisymmetric contact (page 265, Ref. [18]), scaled to
the non-axisymmetric contour a()
.
Fig. 11 Pressure distributions at three indentation depths, as marked in Fig. 10(b). (a) Approximate solution and (b) BEM simulation.
Friction 13
www.Springer.com/journal/40544 | Friction
8 BEM results for elastically inhomogeneous
bodies with complex shapes
Let us, once more, compare the analytical approximate
results obtained above to rigorous numerical
simulations based on the BEM [21]. We will not
reproduce all the extensive numerical studies shown
in Section 5 for homogeneous materials to account
for the elastic grading, because all main findings are
similar in the inhomogeneous case. To demonstrate
the influence of material grading, some results will
suffice.
In Fig. 12, a comparison is shown for the predictions
of the contact domain, based on the analytical
approximate solution and BEM simulations, for
pyramidal indenters (k = 1) with rounded regular
n-polygon cross sections (n = 3, 4, and 6) and for two
exponents of the power-law elastic grading, specifically
m = 0.5 and −0.3. It is apparent that the approximate
solution works significantly better for the case of
positive elastic grading (m > 0, i.e., a soft surface
with a harder material core). In fact, for m > 0, the
approximate solution provides consistently better
results than that in the homogeneous case (m = 0).
This is to be expected, as it is known—and obvious
from the pressure distribution in Eq. (12)—that
positive elastic grading reduces the order of edge and
corner singularities for the stress distributions, and
therefore the quality of the approximate solution (which
is only compromised at sharp edges or corners, as
shown throughout this manuscript) will be improved
by positive grading (and in turn reduced for negative
elastic grading, i.e., hard surfaces with a softer
material core).
This phenomenon is confirmed again in Fig. 13,
showing the comparison for the predicted pressure
distribution due to the indentation of a graded elastic
half-space with m = 0.5 by a pyramidal indenter (k = 1)
with rounded square cross section (n = 4). As can be
seen, the prediction based on the analytical approximate
solution is in almost perfect agreement with the
rigorous numerical solution.
9 Discussion
The main milestones of contact mechanics, which had
(and still have) a lasting effect on practical tribological
applications in science and engineering, were the
Fig. 12 Comparison of contact areas for pyramid indenters (k= 1) in contact with graded materials. The contact boundary is
normalized by the theoretical value of the effective contact radius (red: approximate solution, grey: BEM simulation).
14 Friction
| https://mc03.manuscriptcentral.com/friction
Fig. 13 Pressure distributions at three cross sections for pyramid
with square cross section (m = 0.5, k = 1, and n = 4). (a) Pressure
at cross sections of a, b, and c corresponding to 0°, 30°, and 45°,
respectively. (b) Contour line diagram of pressure distribution
according to BEM simulation. (c) Contour line diagram of pressure
distribution according to approximate solution.
original work of Hertz [1] from 1882 for elliptical
parabolic profiles, followed by the solution of Föppl [3]
and Schubert from [4] in 1941 and 1942, respectively
(which became widespread thanks to Ref. [5] by
Sneddon in 1965) for axially symmetric but otherwise
arbitrary indenter shapes. The solution presented in
the paper at hand is an extension of this axisymmetric
solution to non-axisymmetric profiles.
As the Boussineq problem (i.e., the frictionless
normal contact problem) for arbitrary contact domains
is too complicated to allow for an exact, general
analytical solution, it has been believed for a very long
time that the scope of analytical contact mechanics is
restricted to very specific geometries, specifically,
problems with either axial or plane symmetry, which
were solved in a general form already in Ref. [4].
However, in the present paper, we have laid out a
general procedure for the solution of the Boussinesq
problem for compact, but otherwise arbitrary, contact
domains, which retains the analytical simplicity of
the axisymmetric solution, but which—despite its
approximate nature and asymptotic character of its
derivation—has proved highly robust in its predictions
even in the case of contacts that are far from axial
symmetry.
Nonetheless, some restrictions must be kept in
mind when applying the obtained procedure to real
engineering contacts; restrictions, which, on the one
hand, originate from the modelling abstractions of
the Boussinesq problem (i.e., linear elasticity and the
half-space approximation), and on the other hand,
from the approximate nature of the solution, which,
at least selectively, should be checked by more rigorous
numerical simulations, e.g., based on the BEM or
finite element method (FEM).
On the other hand, the procedure can be used to
solve other classes of contact problems that reduce to
the Boussinesq problem, e.g., the viscoelastic normal
contact—via the elastic–viscoelastic correspondence
principle [22, 23]—or the tangential contact with friction,
within the framework of the Cattaneo [11]–Mindlin [12]
approximation.
The last point deserves a brief elaboration: It has
been shown that the reduction of the tangential
contact with friction to the frictionless normal
contact problem via the principle of Jäger [24] and
Ciavarella [25], in an approximate sense, is also
possible for general three-dimensional [26] and even
rough contacts [27]. The quality of this approximation
(i.e., of the reduction procedure) is of the same order
as the quality of the approximate normal contact
solution discussed in the present manuscript. In other
words, this approximate contact solution can be
applied straightforwardly to incorporate tangential
(frictional) forces.
The approximate contact solution presented in this
manuscript can be used for different tasks in tribology
and engineering. On the one hand, it can be applied
for the fast analysis of macroscopic contacts with
complex shapes, e.g., within the framework of
indentation testing. On the other hand, it may serve
as a tribological tool for modelling single microcontacts
(“asperities”) with complex (or random) shapes in
the contact of rough surfaces. In that regard, one,
however, has to bear in mind that the solution is
intended (in the first place) for the analysis of single
contacts with simply connected (compact) contact
domains.
Friction 15
www.Springer.com/journal/40544 | Friction
Acknowledgements
This work has been conducted under partial financial
support from Deutsche Forschungsgemeinschaft (DFG)
(Grant Nos. PO 810/66-1 and LI 3064/2-1).
Author contributions
Valentin L. POPOV and Emanuel WILLERT obtained
the analytical results, and Qiang LI executed the
numerical simulations. All authors contributed to the
writing of the manuscript.
Declaration of competing interest
The authors have no competing interests to declare
that are relevant to the content of this article. The
author Valentin L. POPOV is the Editorial Board
Member of this journal.
Open Access This article is licensed under a Creative
Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and
reproduction in any medium or format, as long as
you give appropriate credit to the original author(s)
and the source, provide a link to the Creative Commons
licence, and indicate if changes were made.
The images or other third party material in this
article are included in the article’s Creative Commons
licence, unless indicated otherwise in a credit line to
the material. If material is not included in the article’s
Creative Commons licence and your intended use is
not permitted by statutory regulation or exceeds the
permitted use, you will need to obtain permission
directly from the copyright holder.
To view a copy of this licence, visit
http://creativecommons.org/licenses/by/4.0/.
References
[1] Hertz H. Über die Berührung fester elastischer Körper. J
Reine Angew Math 92: 156–171 (1882) (in German)
[2] Popova E, Popov V L. Ludwig Föppl and Gerhard Schubert:
Unknown classics of contact mechanics. Z Angew Math Me
100(9): e202000203 (2020)
[3] Föppl L. Elastische beanspruchung des erdbodens unter
fundamenten. Forschung auf dem Gebiet des Ingenieurwesens
A 12(1): 31–39 (1941) (in German)
[4] Schubert G. Zur frage der druckverteilung unter elastisch
gelagerten tragwerken. Ing Arch 13(3): 132–147 (1942)
(in German)
[5] Sneddon I N. The relation between load and penetration in
the axisymmetric Boussinesq problem for a punch of arbitrary
profile. Int J Eng Sci 3(1): 47–57 (1965)
[6] Barber J R, Billings D A. An approximate solution for the
contact area and elastic compliance of a smooth punch of
arbitrary shape. Int J Mech Sci 32(12): 991–997 (1990)
[7] Shield R T. Load–displacement relations for elastic bodies.
Z Angew Math Phys 18(5): 682–693 (1967)
[8] Barber J R. Determining the contact area in elastic-indentation
problems. J Strain Anal Eng 9(4): 230–232 (1974)
[9] Popov V L. An approximate solution for the contact problem
of profiles slightly deviating from axial symmetry. Symmetry
14(2): 390 (2022)
[10] Popov V L. Correction: Popov, V.L. An approximate solution
for the contact problem of profiles slightly deviating from
axial symmetry. Symmetry 2022, 14, 390. Symmetry 14(10):
2108 (2022)
[11] Cattaneo C. Sul contatto de due corpi elastici: Distribuzione
locale deglisforzi. Rendiconti dell Accademia Nazionale dei
Lincei 27: 342–348 (434–436, 474–478) (1938) (in Italian)
[12] Mindlin R D. Compliance of elastic bodies in contact. J Appl
Mech 16(3): 259–268 (1949)
[13] Fabrikant V I. Flat punch of arbitrary shape on an elastic
half-space. Int J Eng Sci 24(11): 1731–1740 (1986)
[14] Willert E. On Boussinesq’s problem for a power-law graded
elastic half-space on elliptical and general contact domains.
Materials 16(12): 4364 (2023)
[15] Mossakovskii V I. Compression of elastic bodies under
conditions of adhesion (axisymmetric case). J Appl Math
Mech 27(3): 630–643 (1963)
[16] Barber J R. Contact Mechanics. Cham (Switzerland): Springer
Cham, 2018.
[17] Booker J R, Balaam N P, Davis E H. The behaviour of an
elastic non-homogeneous half-space. Part II—Circular and
strip footings. Int J Numer Anal Met 9(4): 369–381 (1985)
[18] Popov V L, Heß M, Willert E. Handbook of Contact
Mechanics: Exact Solutions of Axisymmetric Contact Problems.
Berlin (Germany): Springer Berlin Heidelberg, 2019.
[19] Popov V L, Heß M. Method of Dimensionality Reduction in
Contact Mechanics and Friction. Berlin (Germany): Springer
Berlin Heidelberg, 2015.
[20] Pohrt R, Li Q. Complete boundary element formulation for
normal and tangential contact problems. Phys Mesomech
17(4): 334–340 (2014)
[21] Li Q, Popov V L. Boundary element method for normal
non-adhesive and adhesive contacts of power-law graded
elastic materials. Comput Mech 61(3): 319–329 (2018)
[22] Lee E H. Stress analysis in visco-elastic bodies. Q Appl
16 Friction
| https://mc03.manuscriptcentral.com/friction
Math 13(2): 183–190 (1955)
[23] Khazanovich L. The elastic-viscoelastic correspondence
principle for non-homogeneous materials with time translation
non-invariant properties. Int J Solids Struct 45(17): 4739–4747
(2008)
[24] Jäger J. A new principle in contact mechanics. J Tribol
120(4): 677–684 (1998)
[25] Ciavarella M. The generalized Cattaneo partial slip plane
contact problem. I—Theory. Int J Solids Struct 35(18):
2349–2362 (1998)
[26] Ciavarella M. Tangential loading of general three-dimensional
contacts. J Appl Mech 65(4): 998–1003 (1998)
[27] Paggi M, Pohrt R, Popov V L. Partial-slip frictional response
of rough surfaces. Sci Rep 4: 5178 (2014)
Valentin L. POPOV. He is a full
professor at Technische Universität
(TU) Berlin, Germany. He studied
physics and obtained his Ph.D.
degree in 1985 from Lomonosov
Moscow State University, Russia.
In 1985–1998, he worked at the
Institute of Strength Physics and Materials Science of
the Russian Academy of Sciences, Russia. After that,
he was a guest professor in the field of theoretical
physics at University of Paderborn, Germany, from
1999 to 2002. Since 2002, he has been the head of
Department of System Dynamics and the Physics
of Friction at TU Berlin, Germany. He has published
over 300 papers in leading international journals
and is the author of the book “Contact Mechanics and
Friction: Physical Principles and Applications”, which
appeared in ten editions in German, English, Chinese,
Russian, Spanish, and Japanese. He is the editor-in-
chief of Frontiers in Mechanical Engineering/Tribology,
member of editorial boards of many international
journals, and organizer of more than 20 international
conferences and workshops over diverse tribological
themes. Prof. Valentin L. POPOV is the Honorary
Professor of Tomsk Polytechnic University, Russia, of
East China University of Science and Technology,
China, and of Changchun University of Science and
Technology, China, and Distinguished Guest Professor
of Tsinghua University, China. His areas of interest
include tribology, nanotribology, tribology at low
temperatures, biotribology, the influence of friction
through ultrasound, as well as the numerical simulation
of contact and friction, earthquakes, or synovial joint
regenerative rehabilitation.
Qiang LI. He is a postdoctoral
researcher at TU Berlin, Germany.
He studied mechanical engineering
at East China University of Science
and Technology, China. In 2014,
he obtained his Ph.D. degree at
TU Berlin, Germany, and now he
works as a scientific researcher at the Department of
System Dynamics and the Physics of Friction at TU
Berlin, Germany, headed by Prof. Valentin L. POPOV.
He has published over 20 papers in international
journals, including Physical Review Letters. His
scientific interests include tribology, elastomer
friction, hydrodynamic lubricated contact, numerical
simulation of frictional behaviors, and the boundary
element method.
Emanuel WILLERT. He studied
engineering science and mechanical
engineering at TU Berlin, Germany
and Tomsk Polytechnic University,
Russia. In 2019, he received his Ph.D.
degree with honors from TU Berlin,
Germany and the prize of the
Dimitris N. CHORAFAS Foundation
for his Ph.D. thesis on the fundamentals and
applications of contact–impact problems. He has
published over 20 papers in international journals,
is a co-author of the “Handbook of Contact Mechanics:
Exact Solutions of Axisymmetric Contact Problems”, and
currently works as a postdoctoral researcher at the
Department of System Dynamics and the Physics
of Friction at TU Berlin, Germany. His main
research interests are dynamic contact problems of
inhomogeneous and inelastic materials as well as the
fundamentals of friction and wear.
