Document [original]

Original Article

Boundary element method for

nonadhesive and adhesive contacts

of a coated elastic half-space

Qiang Li

, Roman Pohrt

, Iakov A Lyashenko

1,2

and Valentin L Popov

1,3,4

Abstract

We present a new formulation of the boundary element method for simulating the nonadhesive and adhesive contact

between an indenter of arbitrary shape and an elastic half-space coated with an elastic layer of different material. We use

the Fast Fourier Transform-based formulation of boundary element method, while the fundamental solution is deter-

mined directly in the Fourier space. Numerical tests are validated by comparison with available asymptotic analytical

solutions for axisymmetric flat and spherical indenter shapes.

Keywords

Boundary element method, coated elastic body, adhesive contact, contact mechanics, coatings

Date received: 30 July 2018; accepted: 10 May 2019

Introduction

Layered systems of materials having diﬀerent mech-

anical properties have attracted a lot of scientiﬁc

interest over the last decades.

1,2

A well-chosen coating

can improve the structural, mechanical, optical, or

thermal properties at the surface of a bulk material.

Layered structure can be created by ion implantation,

vacuum deposition, nanostructure burnishing, laser

implantation, and other manufacturing technologies.

Coatings are widely used e.g. for reducing wear,

increasing corrosion resistance, controlling friction,

inﬂuencing adhesion properties, or manipulating

thermal. Due to the signiﬁcant inﬂuence surface

layers have on mechanical properties, a multitude of

experimental techniques has been developed for the

characterization of coatings, in particular measuring

their elastic properties.

For the case of nonadhesive elastic contact with

coated systems, theoretical solutions of various inden-

ters have been obtained. These solutions range from

the early asymptotic solution for line and circular con-

tacts on a single layer coating on an otherwise rigid

foundation,

4,5

to the semi-analytical solution for axi-

symmetric contacts on a general multilayer substrate.

Usually the integral transform method or images

method are used for achieving an analytical solu-

tion.

7,8

The adhesive contact between a rigid sphere

and an elastic multilayer coated half-space was inves-

tigated by use of an integral transform formulation

and Maugis-type adhesion model.

Solutions for axi-

symmetric contacts on a single layer were found using

a JKR-type (Johnson–Kendall–Roberts) adhesion

model.

For a randomly rough surface in contact

with a coated half-space only approximate analytical

theory is available.

Numerical methods have also been intensively

developed to study the contact and tribological

behavior of layered materials. The ﬁnite element

method (FEM) is most commonly used. It is very

versatile and can be applied for various structures

without the restriction of linear material behavior.

In contrast to FEM, the Fast Fourier Transforms

(FFT)-based boundary element method (BEM) is

suitable for all problems where the elastic and geomet-

rical behavior is linear. In part because of its much

higher numerical eﬃciency for contact problems, the

BEM evolved to be the standard method in research

and development. The boundary element formulation

was presented for a contact of an arbitrary shaped

indenter with a homogeneous half-space.

Later,

Berlin University of Technology, Berlin, Germany

Sumy State University, Sumy, Ukraine

National Research Tomsk State University, Tomsk, Russia

National Research Tomsk Polytechnic University, Tomsk, Russia

Corresponding author:

Valentin L Popov, Technische Universita

¨t Berlin, Str. des 17. Juni, Berlin

10623, Germany.

Email: v[email protected]

Proc IMechE Part J:

J Engineering Tribology

2020, Vol. 234(1) 73–83

!IMechE 2019

Article reuse guidelines:

sagepub.com/journals-permissions

DOI: 10.1177/1350650119854250

journals.sagepub.com/home/pij

it was extended to include JKR-type adhesion.

The method was validated by available analytical solu-

tions including parabolic contacts (classical JKR solu-

tion), toroidal indenters,

and ﬂat elliptical

indenters.

Recently, it was also applied to contacts

between ﬂat-ended indenters of complicated shape and

a ﬂat soft body.

The BEM was extended also to

include contacts with power-law graded materials.

In the present paper, we propose a further gener-

alization of BEM for the case of a coated half-space.

Several numerical tests will be carried out and the

results will be compared with the known analytical

solutions.

FFT-based BEM for layered half-space

We consider a half-space with a single elastic coating

of thickness h, elastic modulus E1and Poisson ratio

v1. The corresponding elastic constants of the half-

space are E2and v2(Figure 1). The origin of coord-

inates is placed at the surface of the layer so that the

interface between two media is located at z¼h.

In previous versions of the BEM for contact of

homogeneous and power-law graded materials,

12,17

we proceeded from the fundamental solution in

coordinate space and the corresponding integral for-

mulation of the stress–displacement relation. This

integral relation has the form of a convolution of

the surface pressure distribution with the fundamental

solution U

. This fundamental solution represents the

deformation resulting from a single localized normal

force. For the numerical solution of the contact prob-

lem, we consider a square region on the surface of the

body with the size LL, discretized with Ncells in

each direction. The size of each of the N2square cells

is x¼y¼. Pressure is assumed to be uniform in

each cell (see Figure 2).

In the discretized form, the pressure–displacement

relation can be written as

u¼Kp ð1Þ

where pis stress distribution acting on the surface

(vector of the length N2with values of pressure in

the corresponding discrete cells), uis the normal dis-

placement of surface elements due to applied pressure,

and Kis the compliance matrix having the size N4.

The contact problem is solved in BEM iteratively. In

each step, the displacements for a given (assumed)

pressure distribution are determined by evaluation

of equation (1). Because of the convolution-type

integral equation and the resulting structure of

matrix K, the operation is performed using direct

and inverse FFT

u¼IFFT FFT U0ÞFFTðpÞð½ð2Þ

The number of operation for performing this oper-

ation is on the order of OðN2logNÞ(as compared with

OðN4Þfor direct evaluation of equation (1).

The compliance matrix Kis basically a long version

of the discretized fundamental solution of the prob-

lem. To be used in equation (2), a known fundamental

solution U0must ﬁrst be Fourier-transformed.

However, if the fundamental solution can be found

directly in the Fourier space, this step could be

omitted. In the present paper, we determine the fun-

damental solution for a layered system directly in the

Fourier space in analytic form thus saving both one

Fourier transform and memory space.

For the derivation of the fundamental solution in

Fourier space, we consider a pressure distribution

acting on the surface of the layer in the form of a

plane wave with wave vector kand amplitude p0

p¼p0eikr ð3Þ

Here ris the (two-dimensional) radius vector in the

contact plane. Here and further in the text non-bold

symbol, kdenotes the absolute values of vector,

k¼k

. For simplicity, without loss of generality, let

us use the direction of wave vector kas the x-axis,

thus eikr ¼eikx and equation (3) can be simpliﬁed as

p¼p0eikr ¼p0eikx.

Figure 1. Scheme of the system under consideration. An

elastic layer with thickness h, elastic modulus E

, and Poisson

ratio 

is located on top of an elastic half-space with elastic

parameters E

and 

.Figure 2. Discretization of the simulation area.

74 Proc IMechE Part J: J Engineering Tribology 234(1)

In the homogeneous case, with E1¼E2¼Eand

v1¼v2¼v, the vertical displacement of the surface

(z¼0) is given by

uzx,yðÞ¼uzkðÞeikx with uzkðÞ¼2p0

Ekð4Þ

where Eis the reduced elastic modulus

E¼E=1v2



. In the following, we provide the cor-

responding solution for a layered system.

The equilibrium equation of an elastic isotropic

medium reads

grad div uþ12v1,2



r2u¼0ð5Þ

where v1is the Poisson number of the layer and v2

that of the half-space. The displacement vector uwill

also have the form of a plane wave

u¼uxexþuyeyþuzez¼u0

xðzÞeikxexþu0

zðzÞeikxez

ð6Þ

Here exand ezare unit vectors in directions of the

wave vector and perpendicular to the contact plane

correspondingly. Symbols ux,uyand uzdenote projec-

tions of the displacement vector on the corresponding

directions. The projection in ydirection is zero,

uy¼0. The amplitudes u0

xand u0

zare only functions

of the vertical coordinate z.

Operators appearing in equation (5) read

div u¼@

@xu0

xzðÞeikx



þ@

@zu0

zðzÞeikx



¼iku0

xzðÞeikx þ@u0

zðzÞ

@zeikx ð7Þ

grad div u¼k2u0

xzðÞeikx þik @u0

zðzÞ

@zeikx



þik @u0

xðzÞ

@zeikx þ@2u0

zðzÞ

@z2eikx



ð8Þ

r2u¼@2

@x2þ@2

@z2



xzðÞeikx



þ@2

@x2þ@2

@z2



zðzÞeikx



¼k2u0

xzðÞeikx þ@2u0

xðzÞ

@z2eikx



þ@2u0

zðzÞ

@z2eikx k2u0

zzðÞeikx



ð9Þ

After substitution of these expressions into equa-

tion (5), we have the following relations

@2u0

xzðÞ

@z2þik

12v1,2

@u0

zzðÞ

@z2ð1v1,2Þk2

12v1,2

xzðÞ¼0ð10Þ

@2u0

zzðÞ

@z2ik

2ðv1,2 1Þ

@u0

xzðÞ

@zþð12v1,2Þk2

2ðv1,2 1Þu0

zzðÞ¼0

ð11Þ

We search solution of this system in the form

xzðÞ¼Aelz;u0

zzðÞ¼Belzð12Þ

Substitution of equation (12) into equations (10)

and (11) gives

Al2þBik

12v1,2

l2ð1v1,2Þk2

12v1,2

A¼0ð13Þ

Bl2ik

2ðv1,2 1ÞAlþð12v1,2Þk2

2ðv1,2 1ÞB¼0ð14Þ

The systems (13) and (14) have only trivial solution

if its determinant vanishes

l22ð1v1,2Þk2

12v1,2

ik

2ðv1,2 1Þll

2þð12v1,2Þk2

2ðv1,2 1Þ



¼0ð15Þ

This characteristic equation has four roots l1,2,3,4 ¼

k,k,k,k. Thus the general solution has the form

u1ðÞ

zx,zðÞ¼u0

zzðÞeikx

¼A1ekz þA2ekz þA3zekz þA4zekz



eikx ð16Þ

u1ðÞ

xx,zðÞ¼u0

xzðÞeikx

¼B1ekz þB2ekz þB3zekz þB4zekz



eikx ð17Þ

inside the coating (04z4h) and the same general

form with another set of coeﬃcients inside the half-

space (z4h)

u2ðÞ

zx,zðÞ¼u0

zzðÞeikx

¼A5ekz þA6ekz þA7zekz þA8zekz



eikx ð18Þ

u2ðÞ

xx,zðÞ¼u0

xzðÞeikx

¼B5ekz þB6ekz þB7zekz þB8zekz



eikx ð19Þ

The superscripts (1) and (2) indicate the coating

and half-space respectively. There are 16 coeﬃcients

to be determined. Substitution of u0

zzðÞand u0

xzðÞin

equations (16)–(19) into the diﬀerential equations (10)

and (11) generates

B1¼iA

1þA334v1

ðÞ



B2¼iA

2A434v1

ðÞ



,B3¼iA3,B4¼iA4

ð20Þ

Li et al. 75

B5¼iA

5þA734v2

ðÞ



B6¼iA

6A834v2

ðÞ



,B7¼iA7,B8¼iA8

ð21Þ

We use the following boundary conditions:

(a) Displacements of half-space at inﬁnite depth are

zero: uð2Þ

xx,z!1ðÞ¼0 and uð2Þ

zx,z!1ðÞ¼0;

(b) Continuity of displacements at the interface

between the half-space and coating: uð1Þ

zx,hðÞ¼

uð2Þ

zx,hðÞ,uð1Þ

xx,hðÞ¼uð2Þ

xx,hðÞ;

plane (frictionless problem): ð1Þ

zx x,0ðÞ¼0;

(d) Given normal stress distribution at the surface,

equation (3): ð1Þ

zz x,z¼0ðÞ¼p0eikx;

(e) Continuity of stresses and strains at the interface

between the half-space and coating: ð1Þ

zz x,hðÞ¼

ð2Þ

zz x,hðÞ,ð1Þ

zx x,hðÞ¼ð2Þ

zx x,hðÞ.

The ﬁrst boundary condition (a) leads to

A5¼0, A7¼0, B5¼0, B7¼0ð22Þ

and the others lead to the system of linear algebraic

equations

E1A1k12v1

ðÞþA2k1þ2v1

ðÞ½

þA314v1þ4v2



þA414v1þ4v2







1þv1

ðÞ12v1

ðÞ ¼p0

ð23Þ

A1e2kh þA2þA3he2kh þA4hA6A8h¼0

ð24Þ

kA1þA334v1þkhðÞ½e2kh kA2

þA434v1khðÞþkA6A834v2khðÞ¼0

ð25Þ

1þv1

ðÞ12v1

ðÞ

e2kh A1k12v1

ðÞð



þA31þkh 4v12v1kh þ4v2



þA2k1þ2v1

ðÞ

þA41kh 4v1þ2v1kh þ4v2



E2

1þv2

ðÞ12v2

ðÞ

A6k1þ2v2

ðÞ½

þA81kh 4v2þ2v2kh þ4v2



¼0

ð26Þ

1þv1

ðÞ

A1kþA322v1þkhðÞðÞe2kh



þA2kþA42þ2v1þkhðÞ

E2

1þv2

ðÞ

A6kþA82þ2v2þkhðÞ½¼0

ð27Þ

A1kþA2kþ2A31v1

ðÞþ2A41þv1

ðÞ¼0ð28Þ

For the plain normal contact problem, we only

need normal displacements at the contact surface.

The solution of the systems (23)–(28) can be substi-

tuted into equation (16) and we obtain

uzx,z¼0ðÞ

¼2p01v2



Ae4kh þBkhe2kh þD



kE1Ae4kh Bk2h2e2kh þ2Ce2kh þDðÞ

eikx

ð29Þ

where the constants A, B, C, D are given by the fol-

lowing expressions

A¼E234v1

ðÞ1þv1

ðÞE134v2

ðÞ1þv2

ðÞ½

E11þv2

ðÞE21þv1

ðÞ½

B¼4E21þv1

ðÞþE134v2

ðÞ1þv2

ðÞ½

E11þv2

ðÞE21þv1

ðÞ½

C¼E2

14v23ðÞ1þv2

ðÞ

2E1E2v1þ1ðÞ2v11ðÞv2þ1ðÞ2v21ðÞ

þE2

28v2

112v1þ5



1þv1

ðÞ

D¼E21þv1

ðÞþE134v2

ðÞ1þv2

ðÞ½

E234v1

ðÞ1þv1

ðÞþE11þv2

ðÞ½

ð30Þ

For any given pressure distribution p, the vertical

displacement of the surface at z¼0 can now be calcu-

lated explicitly using equation (29)

u¼IFFT 21v2



Ae4kh þBkhe2kh þD

kAe4kh Bk2h2e2kh

þ2Ce2kh þDÞ

()

FFTðpÞ

ð31Þ

Note that the middle item in the square bracket is

the function of k. For some certain k, it is constant

and one can consider it as the coeﬃcient of the term

FFT(p) for the corresponding value of k. Numerically

both terms are 2D matrices and the operation

<>will then be element-wise multiplication.

Similar to equation (4), this procedure only gives

results for k6¼ 0, in other words, the pressure distri-

bution must have no DC component. The usual BEM

procedure reduces to performing the FFT of pressure

distribution, multiplying the result with the analytical

fundamental solution (equation (29)) and performing

inverse FFT to ﬁnd the displacement ﬁeld. The inverse

problem of ﬁnding pressure for producing given

deformations can be solved by the conjugate graded

method.

Conjugate graded method is a widely used

numerical algorithm for solution of systems of linear

equations, and it has been used in contact problems

76 Proc IMechE Part J: J Engineering Tribology 234(1)

frequently. A detailed discussion of this method and

its extension in contact mechanics can be found in the

paper.

Proper preconditioning

allows keeping the

number of iteration steps of the procedure bounded

by approximately 10 independently of the complexity

of contact conﬁguration and mesh size. These two

steps complete the formulation of BEM for nonadhe-

sive contacts with layered systems. For adhesive con-

tacts, an additional detachment criterion is needed

which is discussed in the remainder of this section.

In each step of an adhesive BEM simulation, the

pressure distribution is calculated in all discretized

grid cells and it has to be decided whether each point

remains in contact or detaches. In Pohrt and Popov

and Popov et al.,

it was suggested to make the deci-

sion based on the energy balance criterion of Griﬃth.

For a nonperiodic system of a homogenous elastic

half-space and a rigid indenter, this leads to a local

mesh-dependent detachment criterion: A surface elem-

ent at the boundary of the contact loses its contact as

soon as tensile stress in this element exceeds the critical

value given by

c¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

E

1

0:473201 

rð32Þ

Here is the speciﬁc work of adhesion between

the indenter and substrate, and E

1¼E1=1v2



Note that this criterion contains only elastic proper-

ties of the coating. This criterion applies also to

layered systems, as long as the size of the discrete

cell is much smaller than the thickness of the layer.

Under this assumption, the elastic energy released due

to the detachment of an element is completely ‘‘con-

ﬁned’’ in the coating, thus the detachment criterion

has the same form as in the case of the homogeneous

material

9,12

with elastic properties of the coating.

The calculation procedure for numerical simula-

tion of an adhesive contact is basically the same as

for nonadhesive contacts. The main diﬀerence is in the

condition for the loss of contact. Instead of requiring

all normal stresses to satisfy p40, the condition

p4cis imposed. Because the adhesive solution

is potentially not unique, we can only approach it

from the state of full contact.

If the entire detachment process is considered, then

starting from full contact, the indenter is moved

upwards by a distance d(displacement-controlled

pull oﬀ) in each step. First it is assumed that the con-

tact area does not change, so that all displacements of

contact points are augmented by d. In the second

stage, the new stress distribution p0is calculated,

which satisﬁes the new displacement ﬁeld (inverse

problem). In the third stage, stresses are checked in

all elements at the boundary of the contact area. If the

tensile stress in an element is larger than the critical

value (32), this element detaches (the stress is set zero),

and a reduced contact area A0is obtained. The stress

distribution is calculated again with the new contact

area and the contact criterion is checked again. This

iteration procedure is continued until the tensile stress

in each element is smaller than cin (32) which means

that the correct contact zone and stresses have been

found. Then the simulation continues with the next

pull-oﬀ step.

Numerical results and comparison with

analytic solutions

In recent studies,

10,22,23

it was shown that the solution

for axis-symmetric adhesive contact problems can be

deduced from the solution of the nonadhesive contact

problem: The critical separation distance dcin adhe-

sive contact is determined by the equation

dksaðÞ

2¼2að33Þ

where ksaðÞis the dependency of the incremental stiﬀ-

ness on the contact radius afor the nonadhesive con-

tact problem. Let us apply this relation to the limiting

case of an elastic layer bonded to a rigid substrate.

The asymptotically exact result in this case with the

additional condition ahis given by

ks¼a2

E1ð34Þ

E1¼E1

1v1

1þv1

ðÞ12v1

ðÞ ð35Þ

Substitution into equation (33) provides the follow-

ing critical values for the indentation depth and the

pull-oﬀ force.

dc¼ ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2h

s,Fc¼ksdc¼a2ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

E1

sð36Þ

In the frame of the proposed BEM formulation, we

can simulate this limiting case by assuming a very large

ratio of E2=E1. Simulation results of the pull oﬀ of a

ﬂat cylindrical indenter are shown in Figure 3 for two

ratios E2=E1¼105and E2=E1¼102with E1¼

2109Pa, v1¼v2¼0:3anda¼50h. The displacement

and the force are normalized to the critical values (36):

F¼F=Fc

vs ~

d¼d=dc

. As expected, the force–

displacement relation is linear up to the moment of

sudden complete detachment. In the case of E2=E1¼

102, this instability point does not match the critical

values (36) ~

F¼1and ~

d¼1. This is due to the fact

that the two ratios E2=E1and a=hare of the same order

of magnitude. Therefore, the asymptotic solution is not

suitable for predicting the detachment. For E2=E1¼

105, the condition for the asymptotical solution is satis-

ﬁed, thus it can be used to validate the correctness of

our numerical criterion. Indeed, at the point of sudden

detachment, the normalized numerical values approach

(1, 1) very closely.

Li et al. 77

Indentation of a parabolic indenter

For a thin elastic layer, an asymptotically exact ana-

lytic solution exists for arbitrary indenter shapes pro-

vided the condition ahis satisﬁed. In the limiting

case, E2!1, displayed in Figure 4, the solution for

a parabolic indenter reads

10,25

F¼~

E1a2

4Rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2h

! ð37Þ

d¼a2

2Rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2h

sð38Þ

The critical values of force, separation, and contact

radius are given by

dcrit ¼ ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2h

s,acrit ¼8R2h



1=4

,Fcrit ¼2R

ð39Þ

With dimensionless variables, ~

a¼a

acrit,~

d¼d

dcrit

F¼F

Fcrit

, the dependencies of the normal force on

indentation depth and contact radius (37) and (38)

can be written in the form

F¼~

a42~

a2ð40Þ

d¼~

a21ð41Þ

These relations are plotted with a black dashed line

in Figure 5(a) and (b). The corresponding numerical

results are shown by the curve with ‘‘plus’’ symbols

for the case of ¼15, which we will deﬁne later

in equation (47), now indicating the case of ah.

The other parameters can be found in the following

discussion. Numerical and asymptotical solutions are

a very close match. The small discrepancy may again

be due to ﬁnite values of E2=E1and a=hused in the

numerical simulation.

In the opposite limiting case of the contact radius

being small compared to the thickness of the layer,

analytical solutions exist in form of asymptotic

series

in dimensionless parameter "¼a=h1

F¼4E

1a3

3R1"38a1

3



13Rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2E

1a

p2E

1a2

ð42Þ

d¼a2

R1"4a0

3"316a1

5þ"432a0a1

92



Rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2E

1a

pE

1a21"2a0

"316a1

3þ"416a0a1

32



ð43Þ

The coeﬃcients in equations (42) and (43) are

given by

am¼1ðÞ

22mm!ðÞ

2Z1

uðÞu2mdu ð44Þ

uðÞ¼2KLe4uLþKþ4uK þ4u2K



e2u

1LþKþ4u2KðÞe2uþKLe4uð45Þ

with coeﬃcients

K¼1n

1þn34v1

ðÞ

,L¼34v2

ðÞn34v1

ðÞ

34v2

ðÞþn,

n¼E21þv1

ðÞ

E11þv2

ðÞ ð46Þ

The dependencies (42) and (43) for Fand ddepend

on the following adhesion parameter

24,26

¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2R2

E

1h3

sð47Þ

Figure 3. Dependencies of the dimensionless elastic force vs.

dimensionless distance for the adhesive contact of a cylindrical

indenter on an elastic layer.

Figure 4. Adhesive contact between a rigid indenter and an

elastic layer bound to a rigid foundation.

78 Proc IMechE Part J: J Engineering Tribology 234(1)

Numerically we carried out the pull-oﬀ simulation

with ﬁve diﬀerent adhesion parameters ranging

from 0.1 to 15, where the constant parameters are

set as E1¼109Pa, E2¼e100 Pa (which means 1, cor-

responding to the limiting case of rigid foundation),

1¼0:3, h¼2 mm, ¼100 J=m2and is varied by

changing the radius of curvature Rof the indenter.

The results are shown in Figure 5 in the same dimen-

sionless coordinates as given by equations (40) and

(41). The curves for ¼0:1 and 0.2 corresponding

to small contact radii are compared with the asymp-

totic relations (42) and (43) while that of large param-

eter are compared with the asymptotic relations (40)

and (41). In both limiting cases, we see very good

agreement between numerical and analytical results.

For intermediate values of an interesting behavior

can be observed. For example in the case of ¼0:5,

for the small indentation depth when the contact

radius is much smaller than the layer thickness

ah, a good coincidence can be observed between

numerical (triangles) and analytical results (dash-dot

line). With an increasing indentation depth, the ana-

lytical solution is not valid any more. At large inden-

tation depths, the numerical results approach the

dashed line (the other limiting case ah) due to a

large contact radius. We thus conclude that numerical

results coincide with all available analytical results in

region of their validity.

Instead of varying the adhesive parameter, let us

now look at the inﬂuence of the elastic parameters of

the foundation on the detachment process. The

approximate equations (42) and (43) provide results

with a high accuracy (see Figure 5) for low values of

. We choose ﬁxed values ¼0:1 and R21 mm

with the other parameters identical as in the above

case but at varying values of E2=E1. Results for

Figure 5. Dependencies of dimensionless contact force on dimensionless indentation depth and contact radius, for the adhesive

indentation of a spherical indenter into a layered counter body. Curves for different values of adhesion parameter aare shown.

Dashed lines depict the dependencies (40) and (41) for the case of a>>h. Dash-dot and solid lines are given by expressions (42) and

(43) for the case of a<<h. Symbols are numerical results obtained by the BEM presented in this paper.

Figure 6. (a) Force–displacement and (b) force–contact radius relation for different ratios of E

Li et al. 79

Fover ~

dand ~

aare shown in Figure 6 and coincide

with analytical approximation very well for diﬀerent

values of E2=E1.

From Figure 6 it can be seen that with increasing

ratio of elastic moduli E2=E1the dependencies

become universal, as the limit of a very rigid half-

space is approached (see Figure 5). Note that in

the case of E2=E1¼0:05, analytical approximation

(solid line) gives wrong result in the range of a large

contact radius ( ~

a41:25) due to the limitation ah,

(see upper limit of Figure 6(b)).

Case studies of flat-ended indenters

Here we present two applications of the above-

validated numerical method: The indentation test of

a square punch on stiﬀ coatings and the adhesive pull

oﬀ of a star-shaped indenter. Both are non-axisym-

metric contact problems and analytical solutions are

not available.

Indentation of a square punch on hard coatings

We consider a nonadhesive contact between a square

indenter with length L0and an elastic half space

coated by a stiﬀ layer (E2=E151). We put special

focus on the contact area A. For the homogeneous

contact, it is known that the contact area is simply

the area of the full square, A0. However, when we

introduce a stiﬀ coating, an interesting behavior can

be observed in the indentation test. If the stiﬀ layer is

thin and the foundation is relatively soft, then the

indentation leads to the loss of contact in the middle

of the square (Figure 7(a)). Because of linearity,

the contact area is independent of the indentation

depth. Figure 7(a) shows the dependency of A=A0

on the ratio of E2=E1for the case of h=L0¼0:5.

It can be seen that the contact area decreases with

decreasing ratio E2=E1(stiﬀer layer) and ﬁnally

approaches a constant plateau. At the plateau, the

contact zone is limited by the edges of the initial

square. There is a critical value of E2=E1dividing

the cases of partial and full contact, e.g. E2=E1

0:15, for the case of the layer thickness h=L0¼0:5.

Figure 7(b) shows that this critical value is dependent

of the relative layer thickness: for a thinner layer,

full contact can be achieved with a softer layer mater-

ial. For a thick layer h=L0¼0:9, the contact is always

complete because the deformation occurs mainly in

the layer. We also ﬁnd that in the case of a very

thin layer h=L0¼0:1 another behavior can be

observed: contact exists both in the center and at

edges (see inset picture in Figure 7(b)).

The described multiple contact behavior has been

observed in the analytical

27,28

and in numerical stu-

dies

for indentation problem of a layered medium

by a rigid one-dimensional (inﬁnite long) ﬂat punch. It

was found that the multiple contact appears only in

the case of stiﬀ layer and soft half space, and the

solutions to the full contact, border contact, and

border center contact depends on the material and

geometrical properties, not on the indentation

depth. For the complicated form of indenters (e.g.,

the studied case) an analytical solution is not possible,

but the BEM shows that behavior is similar to that of

simpler geometries.

Adhesive contact of a star-shaped indenter

The adhesive detachment of a ﬂat indenter with odd

shape from a homogeneous elastic half space has been

studied in Popov et al.

It was found that ﬂat

Figure 7. (a) Dependency of contact area on the ratio of elastic moduli E

. Inset pictures show the contact areas (black color) at

the given ratio. The color picture in the upper left shows a sample contact configuration of the square indenter and the deformed

surface of the layer, with a noncontact bump in the center. (b) Same as (a) but for different values of the relative layer thickness h.

80 Proc IMechE Part J: J Engineering Tribology 234(1)

indenters tend to detach at their outer edges ﬁrst, in

particular at sharp corners. It was also found that for

most shapes Kendall’s solution for the detachment of

a cylindrical punch

applies to both force Fc,homo ¼

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

8Ea3

pand displacement dc,homo ¼ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

2a=E

at the point of ﬁnal detachment. In Kendall’s solution

ais radius of the cylinder, while for odd shapes the

incircle of the shape should be chosen. In the follow-

ing, we will normalize with these two quantities and

also set a¼aincircle.

We would like to exemplify the application of the

above method by simulation of contacts of star-

shaped indenters. The force–displacement relation

for the homogeneous contact can be found in

Figure 8. Similarly to Figures 5 and 6, the force–

displacement dependencies for the star-shaped inden-

ter were obtained by pull-oﬀ BEM simulations. For a

layered system, we set E¼E1=12



. Figure 8

shows the relations for diﬀerent layer thicknesses

and elastic moduli of the layers. It can be seen that

both the thickness of the layer and the elastic proper-

ties have strong inﬂuence on the detachment behavior.

While the dependencies of adhesive force and dis-

placement vary signiﬁcantly, the evolution of the con-

tact area is almost universal and resembles the

homogenous case: The detachment begins from

points which are far from the center and at sharp

corners. When the contact area shrinks to the incircle

of the shape, the two surfaces separate suddenly and

completely. Figure 8(b) right shows the evolution of

the contact area corresponding to the six positions in

all ﬁve curves.

We ﬁnd that deviations from this universal evolu-

tion occur only for extreme values of the stiﬀness

ratio. If the layer is very thin (e.g., h=aincircle ¼0:01)

and soft, then we observe a diﬀerent contact behavior

shown in Figure 9. Here the contact area at the ﬁnal

detachment is signiﬁcantly bigger than the incircle.

The softer the layer, the larger the contact area at

detachment will be.

Conclusions

We generalized the BEM proposed in papers

12,20

for

normal nonadhesive and adhesive contacts of an elas-

tic half-space coated with a layer having a diﬀerent

elastic modulus. It is assumed that the layer is bonded

Figure 8. Example of adhesive detachment of a star-shaped flat indenter from a coated elastic half space. (a) Normalized tensile force

vs. normalized lifting displacement. (b) Evolution of the contact zone during pull off (obtained for homogenous material).

Figure 9. Example of adhesive detachment of a star-shaped

flat indenter as in Figure 8(a) but for very thin and soft layers.

Here the evolution of the contact zone deviates from the

universal behavior of moderate parameters and sudden

detachment occurs at larger contact area.

Li et al. 81

to the elastic half-space and that the contact between

the layer and the indenter is frictionless. The method

is based on the fundamental solution of load–

displacement of the one-layer system in the Fourier

domain and is valid for linear elastic contact

problems. For nonadhesive contact, it is applicable

for arbitrary loading. We presented the simulation

procedure of the adhesive pull oﬀ. The opposite case

of approaching adhesive bodies was not considered in

the present paper. With the suggested BEM formula-

tion, we carried out simulations of adhesive contacts

with cylindrical ﬂat-ended and parabolic indenters

and compared the results with available asymptotic

analytic solutions. We found that numerical results

coincide with all available analytical results in the

regions of their validity. We displayed two sample

applications of non-axisymmetric contacts—the

indentation of a square punch and the pull oﬀ of a

star-shaped ﬂat-ended indenter. In both cases, nontri-

vial behavior has been observed.

Declaration of Conflicting Interests

The author(s) declared no potential conﬂicts of interest with

respect to the research, authorship, and/or publication of

this article.

Funding

The author(s) disclosed receipt of the following ﬁnancial

support for the research, authorship, and/or publication

of this article: Authors acknowledge ﬁnancial support of

the Deutsche Forschungsgemeinschaft (DFG PO 810-55-1)

and the German ministry for research and education

BMBF, grant No. 13NKE011A. This research was also par-

tially supported by the ‘‘Tomsk State University competi-

tiveness improvement program.’’

ORCID iD

Qiang Li https://orcid.org/0000-0001-7458-9450

References

1. Khadem M, Penkov OV, Yang HK, et al. Tribology of

multilayer coatings for wear reduction: a review. Friction

2017; 5: 248–262.

2. Pogrebnjak AD, Bratushka SN, Beresnev VM, et al.

Shape memory effect and superelasticity of titanium

nickelide alloys implanted with high ion doses. Russian

Chem Rev 2013; 82: 1135.

3. Bec S, Tonck A and Loubet JL. A simple guide to deter-

mine elastic properties of films on substrate from nanoin-

dentation experiments. Philosophic Mag 2006; 86: 33–35.

Taylor & Francis: 5347–5358.

4. Meijers P. The contact problem of a rigid cylinder on an

elastic layer. Appl Scient Res 1968; 18: 353–383.

5. Jaffar MJ. Asymptotic behaviour of thin elastic

layers bonded and unbonded to a rigid foundation.

Int J Mech Sci 1989; 31: 229–235.

6. Constantinescu A, Korsunsky AM, Pison O, et al.

Symbolic and numerical solution of the axisymmetric

indentation problem for a multilayered elastic coating.

Int J Solids Struct 2013; 50: 2798–2807.

7. Gupta PK and Walowit JA. Contact stresses between

an elastic cylinder and a layered elastic solid. J Lubricat

Technol 1974; 96: 250–257.

8. Fabrikant VI. Elementary solution of contact prob-

lems for a transversely isotropic elastic layer bonded

to a rigid foundation. J Appl Math Phys 2006; 57:

464–490.

9. Stan G and Adams GG. Adhesive contact between a

rigid spherical indenter and an elastic multi-layer coated

substrate. Int J Solids Struct 2016; 87: 1–10.

10. Papangelo A. Adhesion between a power-law indenter

and a thin layer coated on a rigid substrate. Facta

Universitatis Series Mech Eng 2018; 16: 19–28.

11. Persson BNJ. Contact mechanics for layered materials

with randomly rough surfaces. J Phys Condens Matter

2012; 24: 095008.

12. Pohrt R and Li Q. Complete boundary element formu-

lation for normal and tangential contact problems. Phys

Mesomech 2014; 17: 334–340.

13. Pohrt R and Popov VL. Adhesive contact simulation of

elastic solids using local mesh-dependent detachment

criterion in boundary elements method. Facta

Universitatis Series Mech Eng 2015; 13: 3–10.

14. Argatov I, Li Q, Pohrt R, et al. Johnson-Kendall-

Roberts adhesive contact for a toroidal indenter. Proc

R Soc Lond A Math Phys Eng Sci 2016; 472(2191).

DOI: 10.1098/rspa.2016.0218.

15. Li Q and Argatov I. Onset of detachment in adhesive

contact of an elastic half-space and flat-ended punches

with non-circular shape : analytic estimates and com-

parison with numeric analysis. J Phys D Appl Phys

2018; 145601. DOI: 10.1088/1361-6463/aab28b.

16. Popov VL, Pohrt R and Li Q. Strength of adhesive

contacts: influence of contact geometry and material

gradients. Friction 2017; 5: 308–325.

17. Li Q and Popov VL. Boundary element method for

normal non-adhesive and adhesive contacts of power-

law graded elastic materials. Computat Mech 2017; 61:

319–329.

18. Sneddon IN. Fourier-transform solution of a

Boussinesq problem for a hexagonally aeolotropic elas-

tic half-space. J Mech Appl Math 1992; 45: 607–616.

19. Landau LD and Lifshitz EM. Theory of elasticity.

Course of theoretical physics. 2nd ed. New York:

Pergamon Press, 1970.

20. Vollebregt EAH. A new solver for the elastic normal

contact problem using conjugate gradients, deflation,

and an FFT-based preconditioner. J Computat Phys

2014; 257: 333–351.

21. Griffith AA. VI. The phenomena of rupture and flow in

solids. Philosophic Trans R Soc Lond A Math Phys Eng

Sci 1921; 221: 163–198.

22. Popov VL. Solution of adhesive contact problem on the

basis of the known solution for non-adhesive one. Facta

Universitatis Series Mech Eng 2018; 16: 93–98.

23. Ciavarella M. An approximate JKR solution for a gen-

eral contact, including rough contacts. J Mech Phys

Solids 2018; 114: 209–218.

24. Yang F. Asymptotic solution to axisymmetric indenta-

tion of a compressible elastic thin film. Thin Solid Films

2006; 515: 2274–2283.

25. Li Q and Popov VL. Adhesive contact between a rigid

body of arbitrary shape and a thin elastic coating.

Acta Mechanica 2019; 1: 1–7.

82 Proc IMechE Part J: J Engineering Tribology 234(1)

26. Argatov II, Borodich FM and Popov VL. JKR adhe-

sive contact for a transversely isotropic layer of finite

thickness. J Phys D Appl Phys 2015; 49: 45307.

27. Shield TW and Bogy DB. Multiple region contact solu-

tions for a flat indenter on a layered elastic half space:

plane-strain case. J Appl Mech 1989; 56: 251–262.

28. Urquhart EE and Pindera MJ. Incipient separation

between a frictionless flat punch and an anisotropic

multilayered half plane. Int J Solids Struct 1994; 31:

2445–2461.

29. Chen LF, Urquhart EE and Pindera MJ.

Microstructural effects in multilayers with large

moduli contrast loaded by flat punch. AIAA J 2005;

43: 962–973.

30. Kendall K. The adhesion and surface energy of elastic

solids. J Phys D Appl Phys 1971; 4: 1186.

Li et al. 83