Document [original]

Citation: Burger, H.; Forsbach, F.;

Popov, V.L. Boundary Element

Method for Tangential Contact of a

Coated Elastic Half-Space. Machines

2023,11, 694. https://doi.org/

10.3390/machines11070694

Academic Editor: Matthijn B. De

Rooij

Received: 24 May 2023

Revised: 19 June 2023

Accepted: 24 June 2023

Published: 1 July 2023

Licensee MDPI, Basel, Switzerland.

This article is an open access article

distributed under the terms and

conditions of the Creative Commons

Attribution (CC BY) license (https://

creativecommons.org/licenses/by/

4.0/).

machines

Article

Boundary Element Method for Tangential Contact of a Coated

Elastic Half-Space

Henning Burger * , Fabian Forsbach and Valentin L. Popov *

Department of System Dynamics and Friction Physics, Technische Universität Berlin, 10623 Berlin, Germany

*Correspondence: [email protected] (H.B.); v[email protected] (V.L.P.)

Abstract:

We present a formulation of the boundary element method (BEM) for simulating the tan-

gential contact with an elastic half-space coated with an elastic layer with different elastic properties.

We use the fast Fourier-transform-based formulation of BEM, while the fundamental solution is

determined directly in the Fourier space. Numerical tests are validated by comparison with available

asymptotic analytical solutions for a very thin and a very thick layer, as well as with FEM calculations

for layers with finite thickness.

Keywords:

boundary element method; tangential contact; coatings; contact mechanics; coated surface

1. Introduction

Coatings are widely used to influence and to improve mechanical, electrical, thermal,

adhesive, capillary, and other mechanical properties of surfaces; an overview of related

problems and applications can be found in [

]. One of the most prominent applications

of coatings is wear reduction [

]. The key for understanding and design of coatings is

understanding the contact mechanics of coated materials. Other than for the half-space,

there are no simple analytical solutions for coated materials. However, analytical solutions

have been obtained for limiting cases. Thus, in [

], an analytic asymptotic theory of a normal

contact with a thin elastic layer on a rigid foundation was considered. Further developments

included analytical work related to rough surfaces [

] as well as contact problems with

account of surface tension [

]. Analytical solutions for the tangential contact of coated

surfaces has been derived in [

] for isotropic media, in [

] for a transversely isotropic elastic

layer, and in [

] for the case of a sliding spherical indenter. Numerical solutions of non-

adhesive contact problems for contacts with elastic half-space have been developed since

1990s in the group of Q. Wang (see a review in [

]) and have been extended later to adhesive

contacts [

] and contacts with graded materials [

]. Numerical simulation and calculation

of normal contact of coated surfaces has also been extensively studied. This includes

modelling the contact of specific body shapes [

], or creating contact models, e.g., using

finite element analysis [

], to investigate different contact configurations. The complete

solution for normal contact (both non-adhesive and adhesive) with coated elastic bodies has

been given in [

]. Numerical solutions for tangential contact problems of coated surfaces

are mostly associated with the study of partial slip, as in [

]. In [

], Z. Wang et al.

applied a similar procedure to O’Sullivan and King in [

], using

Papkovich-Neuber

elastic

potentials to derive the corresponding frequency response functions. A semi-analytical

method (SEM) was used to solve the contact problem. Besides SEM and FEM [

], the BEM

is a method that can be used to solve contact problems very efficiently. However, no BEM

solution for the tangential contact of coated surfaces has been presented so far. Therefore,

an extension of the method developed in [

] to tangential contact is described in the

present work; it considers only the tangential part of the contact problem. The FFT-based

formulation of the BEM is used for the solution. Since no points inside the body must be

discretized, but only the points on the surface, this method has a high numerical efficiency

compared to other methods.

Machines 2023,11, 694. https://doi.org/10.3390/machines11070694 https://www.mdpi.com/journal/machines

Machines 2023,11, 694 2 of 15

The remainder of this paper is structured as follows: Section 2, derives in detail the

fundamental solution needed for the BEM formulation. In Section 3, the derived solution is

compared with results from the investigation of limiting cases and FEM simulations. In the

last section a conclusion is drawn.

2. FFT-Based BEM for Tangential Contact of Coated Surfaces

Consider a coated elastic half-space as schematically shown in Figure 1. The layer

having thickness

is assumed to consist of a linearly elastic isotropic material with Young’s

modulus

and Poisson’s ratio

ν1

. The half-space is also an isotropic material with elastic

constants

and

ν2

. The origin of coordinates is placed on the surface of the layer and the

-axis points in the direction of the half-space. The interface between the layer and the

underlying elastic half-space has the coordinate z=h.

Machines 2023, 11, x FOR PEER REVIEW 2 of 17

of the contact problem. The FFT-based formulation of the BEM is used for the solution.

Since no points inside the body must be discretized, but only the points on the surface,

this method has a high numerical eﬃciency compared to other methods.

The remainder of this paper is structured as follows: Section 2, derives in detail the

fundamental solution needed for the BEM formulation. In Section 3, the derived solution

is compared with results from the investigation of limiting cases and FEM simulations. In

the last section a conclusion is drawn.

2. FFT-Based BEM for Tangential Contact of Coated Surfaces

Consider a coated elastic half-space as schematically shown in Figure 1. The layer

having thickness h is assumed to consist of a linearly elastic isotropic material with

Young’s modulus 1

E and Poisson’s ratio

1. The half-space is also an isotropic material

with elastic constants 2

E and

2. The origin of coordinates is placed on the surface of

the layer and the z- axis points in the direction of the half-space. The interface between

the layer and the underlying elastic half-space has the coordinate =zh

Figure 1. A schematic representation of a coated body. The elastic constants of the layer are 1

E and

1 while those of the half-space are 2

E and

2. The thickness of the layer is h and the coordi-

nate origin is on the outer surface of the coating.

For a numerical simulation, a square area with the side length L is considered,

which is discretized with N grid points in each direction. Each square simulation cell

has the same size xyΔ=Δ=Δ

(see Figure 2). For the application of BEM, it is further

assumed that the pressure or stress in each cell is uniform. The usual method for calculat-

ing the tangential displacements u resulting from a tangential stress distribution τ

with the BEM is to perform a direct FFT of the pressure distribution, multiplying it with

the FFT of the fundamental solution and finally performing the inverse FFT as follows [15]

[() ()],=⋅u IFFT FFT FFT

U0τ (1)

where U0 is the fundamental solution, giving the displacement of surface points result-

ing from a single localized tangential force. This procedure is possible for all laterally ho-

mogeneous systems, for which the displacement is represented as a convolution of stress

distribution and fundamental solution. In Fourier-space the convolution transforms to

simple multiplication. A detailed explanation can be found in Ref. [9]. Thus, to calculate

the displacements, both the known fundamental solution and the stress distribution must

be Fourier-transformed. As suggested in [15], it is much easier to derive the fundamental

solution directly in the Fourier-space than in the real space. This also omits one of the

operations in (1). In the following, this fundamental solution is to be derived, that is, in

terms of Equation (1), the factor ()FFT U0.

Figure 1.

A schematic representation of a coated body. The elastic constants of the layer are

and

ν1

while those of the half-space are

and

ν2

. The thickness of the layer is

and the coordinate origin is

on the outer surface of the coating.

For a numerical simulation, a square area with the side length

is considered, which

is discretized with

grid points in each direction. Each square simulation cell has the

same size

∆x=∆y=∆

(see Figure 2). For the application of BEM, it is further assumed

that the pressure or stress in each cell is uniform. The usual method for calculating the

tangential displacements

resulting from a tangential stress distribution

with the BEM

is to perform a direct FFT of the pressure distribution, multiplying it with the FFT of the

fundamental solution and finally performing the inverse FFT as follows [15]

u=IFFT[FFT(U0)·FFT(τ)], (1)

where

is the fundamental solution, giving the displacement of surface points resulting

from a single localized tangential force. This procedure is possible for all laterally homo-

geneous systems, for which the displacement is represented as a convolution of stress

distribution and fundamental solution. In Fourier-space the convolution transforms to

simple multiplication. A detailed explanation can be found in Ref. [

]. Thus, to calculate

the displacements, both the known fundamental solution and the stress distribution must

be Fourier-transformed. As suggested in [

], it is much easier to derive the fundamental

solution directly in the Fourier-space than in the real space. This also omits one of the

operations in (1). In the following, this fundamental solution is to be derived, that is, in

terms of Equation (1), the factor FFT(U0).

Machines 2023,11, 694 3 of 15

Machines 2023, 11, x FOR PEER REVIEW 3 of 17

Figure 2. Representation of the mesh used for the numerical simulation. An area is shown enlarged

to illustrate the discretization and, by way of example, the stress in a cell.

For the derivation of the fundamental solution in Fourier-space, a distribution of tan-

gential stresses

and yz

acting on the surface of the layer is considered in the form

of a plane wave:

xz x x

τττ

kr , (2)

yz y y

τττ

. (3)

x and

are here the amplitudes of the corresponding component,

is the

wave vector and

is the radius vector in the contact plane. In the further text, the symbol

k, which is not printed in bold, denotes the absolute value of the wave vector,

k=k

. For

simplicity, without loss of generality, we can assume that the direction of the wave vector

is given by the

- axis. Thus, Equations (2) and (3) can be written as:

0ikx

ττ

, (4)

0ikx

ττ

. (5)

To obtain equations that contain the displacements and can be evaluated using

boundary conditions, the equilibrium equation of an elastic isotropic medium is used:

grad div (1 2 )

+− ∇ =

uu0

, (6)

where

∇

is the (three-dimensional) gradient operator. The displacements

will in the

- direction also have the form of a plane wave:

000

xyz

() () ()

ikx ikx ikx

x y z x xy yz z

uuuuze uze uze=++= + +ue e e e e e

. (7)

The vectors

y and

are unit vectors pointing in the direction of the coordi-

nate axes.

, y

and

denote the projections of the displacement vector on the cor-

responding directions and the symbols 0

and 0

denote the amplitudes which

depend only on the vertical coordinate

The operators appearing in (6) read:

Figure 2.

Representation of the mesh used for the numerical simulation. An area is shown enlarged

to illustrate the discretization and, by way of example, the stress in a cell.

For the derivation of the fundamental solution in Fourier-space, a distribution of

tangential stresses

τxz

and

τyz

acting on the surface of the layer is considered in the form of

a plane wave:

τxz =τx=τ0

xeikr, (2)

τyz =τy=τ0

yeikr. (3)

τ0

and

τ0

are here the amplitudes of the corresponding component,

is the wave

vector and

is the radius vector in the contact plane. In the further text, the symbol

, which

is not printed in bold, denotes the absolute value of the wave vector,

k=|k|

. For simplicity,

without loss of generality, we can assume that the direction of the wave vector is given by

the x-axis. Thus, Equations (2) and (3) can be written as:

τx=τ0

xeikx, (4)

τy=τ0

yeikx. (5)

To obtain equations that contain the displacements and can be evaluated using bound-

ary conditions, the equilibrium equation of an elastic isotropic medium is used:

grad div u+ (1−2ν1,2)∇2u=0, (6)

where

∇

is the (three-dimensional) gradient operator. The displacements

will in the

x-direction also have the form of a plane wave:

u=uxex+uyey+uzez=u0

x(z)eikxex+u0

y(z)eikxey+u0

z(z)eikxez. (7)

The vectors

and

are unit vectors pointing in the direction of the coordinate

axes.

and

denote the projections of the displacement vector on the corresponding

directions and the symbols

and

denote the amplitudes which depend only on the

vertical coordinate z.

The operators appearing in (6) read:

div u=∂

∂xhu0

x(z)eikxi+∂

∂yhu0

y(z)eikxi+∂

∂zhu0

z(z)eikxi=iku0

x(z)eikx +∂u0

z(z)

∂zeikx, (8)

Machines 2023,11, 694 4 of 15

grad div u=−k2u0

x(z)eikx +ik ∂u0

z(z)

∂zeikxex+ik ∂u0

x(z)

∂zeikx +∂2u0

z(z)

∂z2eikxez, (9)

∇2u=

∂2

∂x2+∂2

∂z2hu0

x(z)eikxiex+∂2

∂x2+∂2

∂z2hu0

y(z)eikxiey+∂2

∂x2+∂2

∂z2hu0

z(z)eikxiez=

h−k2u0

x(z)eikx +∂2u0

x(z)

∂z2eikxiex+−k2u0

y(z)eikx +∂2u0

y(z)

∂z2eikxey+

h−k2u0

z(z)eikx +∂2u0

z(z)

∂z2eikxiez.

(10)

After substitution of these expressions into (6), we obtain:

∂2u0

x(z)

∂z2+ik

1−2ν1,2

∂u0

z(z)

∂z−2(1−ν1,2)k2

1−2ν1,2

x(z) = 0, (11)

∂2u0

y(z)

∂z2−k2u0

y(z) = 0, (12)

∂2u0

z(z)

∂z2−ik

2(ν1,2 −1)

∂u0

x(z)

∂z+(1−2ν1,2)k2

2(ν1,2 −1)u0

z(z) = 0. (13)

We look for solutions of the differential equation system in the form:

x(z)=Aeλz;u0

y(z)=Beλz;u0

z(z)=Ceλz. (14)

Substituting (14) into Equations (11)–(13), we get:

Aλ2+Cik

1−2ν1,2

λ−A2(1−ν1,2)k2

1−2ν1,2

=0, (15)

Bλ2−Bk2=0, (16)

Cλ2−Aik

2(ν1,2 −1)λ+C(1−2ν1,2)k2

2(ν1,2 −1)=0. (17)

Equation (16) is completely decoupled from (15) and (17) and gives

λ1,2 =k

−k

Equations (15) and (17) have non-trivial solutions if their determinant vanishes:



λ2−2(1−ν1,2)k2

1−2ν1,2

1−2ν1,2 λ

−ik

2(ν1,2−1)λ λ2+(1−2ν1,2)k2

2(ν1,2−1)

=0. (18)

The solution of this equation gives the characteristic equation with the four roots

λ3,4,5,6 =k,−k,k,−k. The general solution has the form:

u(1)

x(x,z)=u0

x(z)eikx =A1ekz +A2e−kz +A3zekz +A4ze−kzeikx, (19)

u(1)

y(x,z)=u0

y(z)eikx =B1ekz +B2e−kzeikx, (20)

u(1)

z(x,z)=u0

z(z)eikx =C1ekz +C2e−kz +C3zekz +C4ze−kzeikx. (21)

Machines 2023,11, 694 5 of 15

The superscript

(1)

indicates that the solutions are valid only inside the coating

(0≤z≤h). The general solution inside the half-space has the same form with a different

set of coefficients and the superscript (2):

u(2)

x(x,z)=u0

x(z)eikx =A5ekz +A6e−kz +A7zekz +A8ze−kzeikx, (22)

u(2)

y(x,z)=u0

y(z)eikx =B3ekz +B4e−kzeikx, (23)

u(2)

z(x,z)=u0

z(z)eikx =C5ekz +C6e−kz +C7zekz +C8ze−kzeikx. (24)

Substitution of (19)–(21) and (22)–(24) into Equations (11)–(13) leads to:

C1=−iA1−A3(3−4ν1)

k,

C2=iA2+A4(3−4ν1)

k,

C3=−iA3,C4=iA4.

(25)

C5=−iA5−A7(3−4ν2)

k,

C6=iA6+A8(3−4ν2)

k,

C7=−iA7,C8=iA8.

(26)

We use the following boundary conditions:

1. The displacements of the half-space in infinite depth are zero:

u(2)

x(x,z→∞)=0, u(2)

y(x,z→∞)=0, u(2)

z(x,z→∞)=0. (27)

Continuity of displacements at the interface between the half-space and the coating,

as they are assumed to be bonded together:

u(1)

x(x,h)=u(2)

x(x,h),u(1)

y(x,h)=u(2)

y(x,h),u(1)

z(x,h)=u(2)

z(x,h). (28)

Continuity of stresses at the interface between the half-space and coating, as they are

assumed to be bonded together:

σ(1)

zz (x,h)=σ(2)

zz (x,h),τ(1)

xz (x,h)=τ(2)

xz (x,h),τ(1)

yz (x,h)=τ(2)

yz (x,h). (29)

4. Vanishing of normal stresses at the contact plane:

σ(1)

zz (x,z=0)=0. (30)

5. Given x-component of the tangential stress distribution on the surface:

τ(1)

xz (x,z=0)=τx=τ0

xeikx. (31)

6. Given y-component of the tangential stress distribution on the surface:

τ(1)

yz (x,z=0)=τy=τ0

yeikx. (32)

After their evaluation we obtain:

A5=0, A7=0, B3=0, C5=0, C7=0, (33)

Machines 2023,11, 694 6 of 15

A1e2hk +A2+A3he2hk +A4h−A6−A8h=0, (34)

B1e2hk +B2−B4=0, (35)

A1ke2hk −A2k+A3e2hk(4ν1+hk −3)+A4(4ν1−hk −3)+A6k−A8(4ν2−hk −3)=0, (36)

1+ν1hA1ke2hk +A2k+A3e2hk(2ν1+hk −2)−A4(2ν1−hk −2)i−

1+ν2[A6k−A8(2ν2−hk −2)] =0 ,

(37)

1+ν1hA1ke2hk −A2k+A3e2hk(2ν1+hk −1)+A4(2ν1−hk −1)i+

1+ν2[A6k−A8(2ν2−hk −1)] =0 ,

(38)

E1(1+ν2)B1e2hk −B2+E2(1+ν1)B4=0, (39)

A1k+A2k+2A3(−1+ν1)−2A4(−1+ν1)=0, (40)

1+ν1

[A1k−A2k+A3(−1+2ν1)+A4(−1+2ν1)] =τ0

x, (41)

2(1+ν1)[B1k−B2k]=τ0

y. (42)

For the plain tangential contact, we search for the displacements in

- and

- direction

at the contact surface

u(1)

x(x,z=0)

and

u(1)

y(x,z=0)

. So, we calculate

and

using Equations (34)–(42). The solutions for

are substituted into (19)

and we obtain:

u(1)

x(x,z=0)=2ν2

1−1hAe−4hk +4Bhke−2hk +Di

E1k−Ae−4hk +4Bh2k2e−2hk +2Ce−2hk +Dτ0

xeikx. (43)

Here, the constants A,B,Cand Dare given by the following expressions:

A=E2

2(1+ν1)2(−3+4ν1)+E2

1(1+ν2)2(−3+4ν2)−

2E1E2(1+ν1)(1+ν2)(−3+2ν1+2ν2),(44)

B=E2

2(1+ν1)2+E2

1(1+ν2)2(−3+4ν2)−2E1E2(1+ν1)−1+ν2+2ν2

2, (45)

C=E2

2(1+ν1)25−12ν1+8ν2

1+E2

1(1+ν2)2(−3+4ν2)−

2E1E2−1+ν1+2ν2

1−1+ν2+2ν2

2,(46)

D=−E2

2(1+ν1)2(−3+4ν1)−E2

1(1+ν2)2(−3+4ν2)+

2E1E2(1+ν1)(1+ν2)(5−6ν2+ν1(−6+8ν2)) .(47)

Now we substitute B1and B2into (20) and obtain:

u(1)

y(x,z=0)=2(1+ν1)hFe−2hk −Gi

E1kFe−2hk +Gτ0

yeikx, (48)

where the constants Fand Gare given by the following expressions:

F=E2(1+ν1)−E1(1+ν2), (49)

Machines 2023,11, 694 7 of 15

G=E2(1+ν1)+E1(1+ν2). (50)

We write the solutions (43) and (48) for uxand uyin the following abbreviated form:

ux=u(1)

x(x,z=0)=Φx(k)τ0

xeikx ,

uy=u(1)

y(x,z=0)=Φy(k)τ0

yeikx .(51)

In the above derivation, we have chosen the

-axis along the wave vector. However,

on the back Fourier transformation, integration goes over all possible wave vectors at the

given stress on the surface. To be able to perform this operation, we now write the result

(51) in a coordinate system where both the wave vector and the stress vector have arbitrary

directions. To this end a new coordinate system

(x0,y0)

is considered, which is rotated

relative to (x,y)by angle ϕas shown in Figure 3.

Machines 2023, 11, x FOR PEER REVIEW 7 of 17

()

()( )

()()

222

21 1112 2

12 1 1 2 2

15128 1 34

21212

CE E

EE ,

ννν ν ν

νν νν

=+ −+++−+−

−+ + −+ +

(46)

()( )()( )

()() ( )

()

21 112 2

12 1 2 2 1 2

134134

211 56 68

DE E

EE .

νν νν

νν νν ν

=− + − + − + − + +

++ −+−+

(47)

Now we substitute

and

into (20) and obtain:

()

() ikx

Fe G

ux,z e,

Ek Fe G

−



+−



== 



(48)

where the constants

and G are given by the following expressions:

()()

2112

11FE E

νν

=+−+

, (49)

()()

2112

11GE E

νν

=+++

. (50)

We write the solutions (43) and (48) for

and

in the following abbreviated

form:

()()

() ikx

xx x x

() ikx

yy y y

uux,z ke,

uux,z ke.

===Φ

(51)

In the above derivation, we have chosen the

-axis along the wave vector. However,

on the back Fourier transformation, integration goes over all possible wave vectors at the

given stress on the surface. To be able to perform this operation, we now write the result

(51) in a coordinate system where both the wave vector and the stress vector have arbitrary

directions. To this end a new coordinate system

()

', 'xy

is considered, which is rotated

relative to

()

by angle

as shown in Figure 3.

Figure 3. Representation of the new coordinate system

()

', '

, which is rotated relative to

()

by the angle

The coordinate transformations read:

cos sin

sin cos

x' x y

y' x y

uu u ,

ϕϕ

=−

(52)

Figure 3.

Representation of the new coordinate system

(x0,y0)

, which is rotated relative to

(x,y)

the angle ϕ.

The coordinate transformations read:

ux0=uxcos ϕ−uysin ϕ,

uy0=uxsin ϕ+uycos ϕ,(52)

τ0

x=τ0

x0cos ϕ+τ0

y0sin ϕ,

τ0

y=−τ0

x0sin ϕ+τ0

y0cos ϕ,(53)

x=x0cos ϕ+y0sin ϕ,

y=−x0sin ϕ+y0cos ϕ.(54)

Substitution of (53) and (54) into (51) gives:

ux=Φx(k)τ0

x0cos ϕ+τ0

y0sin ϕeik(x0cos ϕ+y0sin ϕ),

uy=Φy(k)−τ0

x0sin ϕ+τ0

y0cos ϕeik(x0cos ϕ+y0sin ϕ).

(55)

Further substitution of (55) into (52) leads to the result:

ux0=hΦx(k)τ0

x0cos2ϕ+τ0

y0sin ϕcos ϕ−Φy(k)−τ0

x0sin2ϕ+τ0

y0sin ϕcos ϕieik(x0cos ϕ+y0sin ϕ),

uy0=hΦx(k)τ0

x0sin ϕcos ϕ+τ0

y0sin2ϕ+Φy(k)−τ0

x0sin ϕcos ϕ+τ0

y0cos2ϕieik(x0cos ϕ+y0sin ϕ).

(56)

Machines 2023,11, 694 8 of 15

If we assume that the stress vector

is directed along the axis

, then

τ0

y0=

0 and (56)

takes the form:

ux0=τ0

x0Φx(k)cos2ϕ+Φy(k)sin2ϕeik(x0cos ϕ+y0sin ϕ),

uy0=τ0

x0Φx(k)−Φy(k)sin ϕcos ϕeik(x0cos ϕ+y0sin ϕ).(57)

These equations show that the traction along the axis

lead to displacements both

in the direction of traction and perpendicular to it. The reverse is also true: a displace-

ment in the direction

will lead to the appearance of stress component perpendicular to

this direction.

Now it is possible to calculate the displacements in the direction of the loading and

perpendicular to it. To use it in a BEM code, we need to convert it into an inverse Fast

Fourier Transform so that:

ux=IFFThnΦx(k)cos2ϕ+Φy(k)sin2ϕo·FFT(τ)i, (58)

uy=IFFTΦx(k)−Φy(k)sin ϕcos ϕ·FFT(τ). (59)

The operation

h·i

denotes an element-wise multiplication since both terms are 2D

matrices. For the inverse problem, the calculation of the tangential stress distribution from

a given displacement field, the conjugate gradient method can be used. This completes the

formulation of the BEM for the purely tangential contact of coated systems.

3. Comparison with Limiting Cases and FEM Solutions

In order to check the correctness of the derivation and the resulting Equations (58)

and (59), comparisons are made with other solutions. For this purpose, limiting cases are

investigated and results from FEM calculations are used. The limiting cases are, first, the

case of a layer with infinite thickness and, second, the case of a thin layer on a rigid surface.

Comparison solutions can be easily created for these two cases.

To start with the general case, the comparison with FEM solutions is presented first.

More detailed information on the FEM model and how to obtain the FEM results can be

found in Appendix A. For comparison, we consider a boundary value problem in which a

circular contact area of diameter 2

is displaced tangentially by

. The resulting tangential

force

and the resulting tangential contact stiffness

are calculated. We are interested in

the influence of the ratio of elastic moduli

E1/E2

on the contact stiffness at different ratios

a/h. The Poisson’s ratio of the layer and the substrate should be the same (ν1=ν2=0.3).

To work only with dimensionless parameters, the normalized tangential contact stiffness

kT,norm is defined, which can be calculated as follows:

kT,norm =kT

kT,hom

. (60)

The contact stiffness resulting from BEM or FEM simulations is thus divided by the

contact stiffness that would result without the layer. This contact stiffness can be calculated

analytically for our case [19]

kT,hom =2aG∗

2, (61)

with G∗

2=4G2/(2−ν2)and G2as the shear modulus of the substrate.

Figure 4shows the results and the comparison. The normalized tangential contact

stiffness is plotted against the ratio of elastic moduli for different ratios a/h. The different

types of lines represent the BEM solutions, and the markings represent the FEM solutions.

The comparison was made for each case: the contact radius is smaller, larger, or equal to

the layer thickness. In addition, two other ratios a/hwere calculated with the BEM.

Machines 2023,11, 694 9 of 15

Machines 2023, 11, x FOR PEER REVIEW 9 of 17

tangential force

and the resulting tangential contact stiﬀness

are calculated. We

are interested in the influence of the ratio of elastic moduli

on the contact stiﬀness

at diﬀerent ratios

. The Poisson’s ratio of the layer and the substrate should be the

same

()

0.3

νν

. To work only with dimensionless parameters, the normalized tan-

gential contact stiﬀness ,normT

is defined, which can be calculated as follows:

,norm

,hom

(60)

The contact stiﬀness resulting from BEM or FEM simulations is thus divided by the

contact stiﬀness that would result without the layer. This contact stiﬀness can be calcu-

lated analytically for our case [19]

,hom 2

kaG=

(61)

with

()

22 2

42GG

=−

and

as the shear modulus of the substrate.

Figure 4 shows the results and the comparison. The normalized tangential contact

stiﬀness is plotted against the ratio of elastic moduli for diﬀerent ratios

. The diﬀerent

types of lines represent the BEM solutions, and the markings represent the FEM solutions.

The comparison was made for each case: the contact radius is smaller, larger, or equal to

the layer thickness. In addition, two other ratios

were calculated with the BEM.

Figure 4. The normalized tangential contact stiﬀness

,normT

plotted against the ratio of elastic mod-

uli

for diﬀerent ratios

. The BEM solutions are represented by diﬀerent types of lines

for the diﬀerent ratios ah

. The FEM solutions are represented by markers that have the same color

as the corresponding line. The Poisson’s ratio of the layer and the substrate are

0.3

νν

As can be seen, the agreement between the results of the two simulation methods is

very good. Moreover, all curves meet at the point ,norm

k= and

1EE=

, as it is to be

expected, since it is the case of the homogeneous half-space. We can also see that the layer

Figure 4.

The normalized tangential contact stiffness

kT,norm

plotted against the ratio of elastic moduli

E1/E2

for different ratios

a/h

. The BEM solutions are represented by different types of lines for the

different ratios

a/h

. The FEM solutions are represented by markers that have the same color as the

corresponding line. The Poisson’s ratio of the layer and the substrate are ν1=ν2=0.3.

As can be seen, the agreement between the results of the two simulation methods is

very good. Moreover, all curves meet at the point

kT,norm =

1 and

E1/E2=

1, as it is to

be expected, since it is the case of the homogeneous half-space. We can also see that the

layer stiffness has a direct influence on the tangential contact stiffness of the contact system.

When

E1>E2

kT,norm >

1, since the contact stiffness is larger than in the homogeneous

case. The reverse case also applies, so that

kT,norm <

1 if

E1<E2

. The influence of the

layer depends of course on its thickness. For thinner layers

(a/h>1)

kT,norm

tends to

become a constant with a value of 1, since the properties of the substrate dominate the

contact configuration. For thicker layers

(a/h<1)

, the relation between

kT,norm

and the

ratio

E1/E2

becomes linear, as the properties of the layer are more dominant, which can

also be represented analytically. If the layer is thick enough, we can assume it to be a

homogeneous half-space. In this case, (61) can be used for

, but with

G∗

1=4G1/(2−ν1)

Thus the Equation (60) reads (with ν1=ν2=0.3):

kT,norm =2aG∗

2aG∗

=G1

=E1

. (62)

This relation can be seen in Figure 4for a/h=0.05.

3.1. Limiting Cases for Tangential Contact without Slip

The procedure for checking the limiting cases is similar to the comparison with the

FEM solutions. We again consider a circular contact area with diameter 2

, which is

displaced tangentially by

, and calculate the tangential force and the resulting tangential

contact stiffness

. Since in both limiting cases the layer thickness is important, we are

now interested in the influence of the ratio

a/h

on the contact stiffness. To investigate both

cases in an efficient way, the following considerations were made.

For the first case, the infinite layer thickness, very large values for

should be used.

Now we can use the assumption we have already used for the comparison with the FEM

results: If the layer is thick enough, it can be assumed to be an elastic half-space with

the appropriate elastic parameters. Thus, the results of the BEM simulation should be

consistent with those obtained by assuming a homogeneous elastic half-space with the same

parameters as the layer. As written above, the tangential contact stiffness in this case can be

calculated with (61), with the difference that now

G∗

is used. For the second limiting case,

the thin layer on a rigid surface, very small values for

should be used and, in addition,

Machines 2023,11, 694 10 of 15

the value of

must be large. In this case, it is also possible to calculate the tangential

contact stiffness analytically. The derivation of the equation will be briefly illustrated.

Considering a thin elastic layer of thickness

on a rigid surface, the shear strain

γxz

can be approximated as follows:

γxz =∂ux

∂z≈ux

h. (63)

The thin elastic layer acts as a three-dimensional Winkler foundation. Thus, the shear

strain does not depend on

and is constant under the displaced area. The resulting shear

stress is then:

τxz =G1γxz ≈G1

h. (64)

This results in the following equation for the tangential contact stiffness:

kT,ana ≈G1

πa2

h. (65)

The use of the model of a Winkler foundation for a thin elastic layer can also be

found in other literature [

]. To illustrate the comparison, we now use the normalized

tangential contact stiffness again by dividing each calculated contact stiffness by that of

the homogeneous half-space (but now with the parameters of the layer). Since the elastic

modulus of the substrate should not matter in the first case and must be large in the second,

E1/E2=10−9is used to mimic a rigid surface as closely as possible.

In Figure 5the results are shown. The BEM results are represented by solid lines

and the limiting cases by dashed or dotted lines. In addition to the previous theoretical

considerations and for the sake of completeness, the comparison is performed for two

different values of

ν1

, since

kT,hom =

aG∗

depends on this quantity. However, as can

be seen, the influence of the Poisson’s ratio of the layer on the normalized tangential

contact stiffness is not large. The more interesting point is that the limiting cases are very

well reproduced with the derived BEM solution. For very thick layers

(a/h1)

, the

normalized tangential contact stiffness becomes 1, as this is very close to the homogeneous

case. For very thin layers

(a/h1)

there is a linear relation between

kT,norm

and

a/h

, as

can be seen from the equations. To give information about the deviations between BEM

results and the analytical results: For

a/h=

0.1, the deviation between the BEM result and

the analytical result for the homogeneous half-space is 4%. The same deviation results for

the case of the thin layer on a rigid surface for a/h=22.

Machines 2023, 11, x FOR PEER REVIEW 11 of 17

The use of the model of a Winkler foundation for a thin elastic layer can also be found

in other literature [20]. To illustrate the comparison, we now use the normalized tangential

contact stiﬀness again by dividing each calculated contact stiﬀness by that of the homo-

geneous half-space (but now with the parameters of the layer). Since the elastic modulus

of the substrate should not matter in the first case and must be large in the second,

10EE

−

= is used to mimic a rigid surface as closely as possible.

In Figure 5 the results are shown. The BEM results are represented by solid lines and

the limiting cases by dashed or dotted lines. In addition to the previous theoretical con-

siderations and for the sake of completeness, the comparison is performed for two diﬀer-

ent values of

, since

hom 1

kaG= depends on this quantity. However, as can be seen,

the influence of the Poisson’s ratio of the layer on the normalized tangential contact stiﬀ-

ness is not large. The more interesting point is that the limiting cases are very well repro-

duced with the derived BEM solution. For very thick layers

()

1ah<< , the normalized

tangential contact stiﬀness becomes

, as this is very close to the homogeneous case. For

very thin layers

()

1ah>> there is a linear relation between

,normT

k and

, as can be

seen from the equations. To give information about the deviations between BEM results

and the analytical results: For

0.1ah=

, the deviation between the BEM result and the

analytical result for the homogeneous half-space is

. The same deviation results for

the case of the thin layer on a rigid surface for

22ah=

Figure 5. Comparison between the BEM results (solid lines) and analytical results of the limiting

cases (dotted and dashed lines) for diﬀerent ratios ah

and for two diﬀerent

. The elastic mod-

ulus of the substrate is much larger than that of the layer

()

1210EE −

3.2. Limiting Cases for the Tangential Contact of a Parabolic Indenter Considering Slip

The same limiting cases are now tested for the contact of a parabolic indenter, taking

partial slip into account. For this purpose, a parabolic indenter with the radius of curva-

ture

is first pressed in by

and then tangentially displaced by

. For the calcula-

tion, we assume that the normal and the tangential contact can be considered as

Figure 5.

Comparison between the BEM results (solid lines) and analytical results of the limiting

cases (dotted and dashed lines) for different ratios

a/h

and for two different

ν1

. The elastic modulus

of the substrate is much larger than that of the layer E1/E2=10−9.

Machines 2023,11, 694 11 of 15

3.2. Limiting Cases for the Tangential Contact of a Parabolic Indenter Considering Slip

The same limiting cases are now tested for the contact of a parabolic indenter, taking

partial slip into account. For this purpose, a parabolic indenter with the radius of curvature

is first pressed in by

and then tangentially displaced by

. For the calculation, we

assume that the normal and the tangential contact can be considered as decoupled. We

also assume Coulomb friction with the coefficient of friction

. With these assumptions, an

iterative procedure can be used to investigate the partial slip and calculate the

stick-slip

regions for each incremental displacement

∆ux

. To do this, the no-slip condition

|τ|<µp

is checked for each contact point after each incremental displacement. The pressure distri-

bution

is calculated using the BEM formulation presented in [

]. As can be seen, only

the absolute value of the tangential stress

is considered in the no-slip condition. Thus, the

Cattaneo-Mindlin assumption is used, and the exact directions of the displacements and

the friction force are not considered [21,22].

For the comparison between the calculations and the limiting cases, we are interested

in the ratio of the radii of the stick region

and the total contact area

for each

∆ux

To derive asymptotic analytical solutions, the Ciavarella-Jäger assumption is used, with

which the tangential contact can be determined by the corresponding frictionless normal

contact [23,24]. Thus, among others, the following equation can be used:

τ(r)=µ[p(r;a)−p(r;c)] . (66)

For the case of infinite layer thickness, the solution of the homogeneous half-space can

be used as before, which has the same elastic parameters as the layer. The derivation can

be found for example in [20] and the resulting equation is:

c

ahom =1−Fx

µFN1/3

. (67)

For the case of the thin layer on a rigid surface, the model of a three-dimensional

Winkler foundation is again used. Thus, Equation (64) can be used for the shear stress in

the stick region. The derivation of the pressure distribution

p(r;e

is similar to that of the

shear stress. The strain can be approximated as:

εzz =duz

dz≈ −1

hδ−r2

2R, (68)

where

r2/2R

describes the shape of the indenter. The known relation between the indenta-

tion depth

and the corresponding contact radius

δ=e

a2/2R

, can be used and thus the

following pressure distribution is obtained:

p(r;e

a)=−σzz =G1

2(1−ν1)

(1−2ν1)e

a2−r2

2R. (69)

The representation with

serves to facilitate the transformation of the equation after

(64) and (69) have been substituted into (66). In this way, the following result is obtained:

c

aana =1−(1−2ν1)

2(1−ν1)

µd1/2

. (70)

The comparison is shown in Figure 6. The tangential displacement is normalized to

ux,max

, i.e., the maximum tangential displacement at which the stick region just vanishes in

the homogeneous case. It can be calculated as follows:

ux,max =µ2−ν1

2(1−ν1)δ. (71)

Machines 2023,11, 694 12 of 15

Machines 2023, 11, x FOR PEER REVIEW 13 of 17

The asymptotic analytical results are represented by the solid and dashed lines and

the BEM results for diﬀerent ratios

theo

by markers. The contact radius

theo

2aR

is used to have a fixed value for comparison since the real contact radius changes for dif-

ferent layer thicknesses. In addition, the contact configuration with the stick-slip regions

05ca .=

for

theo 1ah=

is shown as an inset image.

Figure 6. The ratio of radii ca

plotted against the ratio

,maxxx

. The analytical results of the

limiting cases (solid and dashed lines) are compared with the BEM results (markers) for diﬀerent

ratios

theo

. The elastic parameter data are given and the contact configuration at

05ca .=

for

theo

1ah=

is shown as an inset image.

As can be seen, the agreement between the BEM results and the limiting cases is very

good. For very thick layers

()

theo

002ah.= the course of the BEM result is equal to the

course described by Equation (67). For both results, the ratio

max

xx,

uu = be-

comes zero due to normalization. For thinner layers, larger tangential displacements are

required to reach the full slip state, as can be seen for

theo 1ah=

. The ratio of the radii

becomes zero for

max

xx,

uu >. The agreement between the BEM result and Equation (70)

for the limiting case of the thin layer can be seen for

theo

100ah=

. Although the curves

do not overlap as well as in the infinite layer thickness case, the deviation between both

results is only about 1%.

For the same values of

theo

, the tangential stress distributions normalized to their

maximum are plotted against the radial coordinate normalized to the contact radius

and represented by markers (see Figure 7). In addition, the homogeneous case is plotted

and represented by a solid line. As before, the homogeneous case can be represented very

well by the result of

theo

002ah.=

. Even for

theo

1ah=

there is no big diﬀerence to the

homogeneous case with the used normalizations, especially in the slip region. The course

for a very thin layer

()

theo

100ah= looks quite diﬀerent. In the stick region the tangential

stress is constant, which is in very good agreement with Equation (64), in which no

Figure 6.

The ratio of radii

c/a

plotted against the ratio

ux/ux,max

. The analytical results of the

limiting cases (solid and dashed lines) are compared with the BEM results (markers) for different

ratios

atheo/h

. The elastic parameter data are given and the contact configuration at

c/a=

0.5 for

atheo/h=1 is shown as an inset image.

The asymptotic analytical results are represented by the solid and dashed lines and

the BEM results for different ratios

atheo/h

by markers. The contact radius

atheo =√2Rδ

is used to have a fixed value for comparison since the real contact radius changes for

different layer thicknesses. In addition, the contact configuration with the stick-slip regions

at c/a=0.5 for atheo/h=1 is shown as an inset image.

As can be seen, the agreement between the BEM results and the limiting cases is very

good. For very thick layers

(atheo/h=0.02)

the course of the BEM result is equal to the

course described by Equation (67). For both results, the ratio

c/a

ux/ux,max =

1 becomes

zero due to normalization. For thinner layers, larger tangential displacements are required

to reach the full slip state, as can be seen for

atheo/h=

1. The ratio of the radii becomes

zero for

ux/ux,max >

1. The agreement between the BEM result and Equation (70) for the

limiting case of the thin layer can be seen for

atheo/h=

100. Although the curves do not

overlap as well as in the infinite layer thickness case, the deviation between both results is

only about 1%.

For the same values of

atheo/h

, the tangential stress distributions normalized to their

maximum are plotted against the radial coordinate normalized to the contact radius

and

represented by markers (see Figure 7). In addition, the homogeneous case is plotted and

represented by a solid line. As before, the homogeneous case can be represented very

well by the result of

atheo/h=

0.02. Even for

atheo/h=

1 there is no big difference to

the homogeneous case with the used normalizations, especially in the slip region. The

course for a very thin layer

(atheo/h=100)

looks quite different. In the stick region the

tangential stress is constant, which is in very good agreement with Equation (64), in

which no dependency on the ratio

r/a

is given. Thus, the tangential stress distributions

of the BEM results and the limiting cases can also be successfully compared using the

above assumptions.

Machines 2023,11, 694 13 of 15

Machines 2023, 11, x FOR PEER REVIEW 14 of 17

dependency on the ratio

is given. Thus, the tangential stress distributions of the BEM

results and the limiting cases can also be successfully compared using the above assump-

tions.

Figure 7. Tangential stress distributions at

05ca .=

for diﬀerent ratios

theo

, plotted against

the ratio

, normalized to its maximum and represented by markers. The homogeneous case is

represented by a solid line. Information on the elastic parameters is given.

4. Conclusions

The derivation of a fundamental solution for the calculation of the tangential contact

of a coated elastic half-space and its integration into the boundary element method are

presented. To make the integration as simple as possible, the solution was derived in Fou-

rier space. It is assumed that the layer is bonded to the substrate and there is no normal

load to consider. The main conclusions are as follows:

•

The derived solution works very well, as comparison with asymptotic analytical so-

lutions for a very thin and a very thick layer, as well as with FEM results for finite

layer thicknesses, shows very good agreement.

•

With the new formulation, tangential contact problems between an arbitrarily shaped

indenter and an elastic half-space coated with a layer with diﬀerent elastic properties

can be simulated and computed. In addition, arbitrary tangential loads can be con-

sidered.

•

It might be interesting to take a closer look at the findings that have been pointed out

and briefly discussed.

•

Some strong assumptions, such as decoupling, Cattaneo-Mindlin and Ciavarella-Jä-

ger, are used to study the stick-slip regions of a parabolic indenter. A study without

these assumptions, that is with accounting of coupling terms and considering the

direction of the displacements and the frictional force will be presented in a separate

publication.

Author Contributions: conceptualization, H.B. and V.L.P.; analytical solution, H.B. and V.L.P.; soft-

ware, H.B.; validation, H.B. and F.F.; formal analysis, H.B.; writing—original draft preparation, H.B.;

Figure 7.

Tangential stress distributions at

c/a=

0.5 for different ratios

atheo/h

, plotted against

the ratio

r/a

, normalized to its maximum and represented by markers. The homogeneous case is

represented by a solid line. Information on the elastic parameters is given.

4. Conclusions

The derivation of a fundamental solution for the calculation of the tangential contact

of a coated elastic half-space and its integration into the boundary element method are

presented. To make the integration as simple as possible, the solution was derived in

Fourier space. It is assumed that the layer is bonded to the substrate and there is no normal

load to consider. The main conclusions are as follows:

•

The derived solution works very well, as comparison with asymptotic analytical

solutions for a very thin and a very thick layer, as well as with FEM results for finite

layer thicknesses, shows very good agreement.

•

With the new formulation, tangential contact problems between an arbitrarily shaped

indenter and an elastic half-space coated with a layer with different elastic prop-

erties can be simulated and computed. In addition, arbitrary tangential loads can

be considered.

•

It might be interesting to take a closer look at the findings that have been pointed out

and briefly discussed.

•

Some strong assumptions, such as decoupling, Cattaneo-Mindlin and Ciavarella-Jäger,

are used to study the stick-slip regions of a parabolic indenter. A study without these

assumptions, that is with accounting of coupling terms and considering the direction

of the displacements and the frictional force will be presented in a separate publication.

Author Contributions:

Conceptualization, H.B. and V.L.P.; analytical solution, H.B. and V.L.P.; software,

H.B.; validation, H.B. and F.F.; formal analysis, H.B.; writing—original draft preparation, H.B.; editing

and final text, V.L.P.; visualization, H.B.; supervision, V.L.P.; project administration, V.L.P.; funding

acquisition, V.L.P. All authors have read and agreed to the published version of the manuscript.

Funding:

We acknowledge support by the German Research Foundation and the Open Access

Publication Fund of TU Berlin. This research was funded by the Deutsche Forschungsgemeinschaft

(DFG, German Research Foundation) under project number PO 810/69-1.

Institutional Review Board Statement: Not applicable.

Data Availability Statement: Not applicable.

Acknowledgments: The authors acknowledge valuable discussions with Q. Li and I. Lyashenko.

Machines 2023,11, 694 14 of 15

Conflicts of Interest:

The authors declare no conflict of interest. The funder had no role in the design

of the study, in the collection, analyses, or interpretation of data, in the writing of the manuscript, or

in the decision to publish the results.

Appendix A FE-Model of the Tangentially Loaded Coated Half-Space

The finite element simulations were carried out using the commercial software Abaqus

Standard. As shown in Figure A1, the tangential loading was applied by prescribing a

tangential displacement to a circular area on the coating with radius

. We modelled the

coating of thickness

as a cylinder and the half-space below as a hemisphere of the same

radius. Due to symmetry in the

x−z

plane, it is sufficient to use a half model with the

appropriate symmetry conditions. We meshed the half of the cylinder with approximately

25 k–100 k C3D8R elements (depending on the ratio

a/h

) and the half of the hemisphere

using approximately 300k C3D10 elements. Especially for the case of a softer half-space

G2<G1

, any regular boundary conditions (fixed or free) on the outside influence the results

falsely. Thus, we used infinite elements CIN3D8 that are connected to the outer boundary

of the cylinder and hemisphere that assume a linear far field solution. As described in the

Abaqus Documentation [

], the radial dimension of these elements must be the distance

of their inner point to the pole (radius

in Figure A1). For each geometry

a/h

, both

the mesh size and the dimension of the region with regular mesh were determined by

studying the convergence of the tangential contact stiffness. The ratio

R/a=

20 showed

good convergence for all studied cases. For the homogeneous half-space (

G2=G1

), the

agreement with the available analytical solution in Equation (61) is very good with less

than 1% error for all geometries (see also Figure 4).

Machines 2023, 11, x FOR PEER REVIEW 16 of 17

Figure A1. Finite element model and mesh of the tangentially loaded coated half-space.

References

1. Bhushan, B.; Peng, W. Contact Mechanics of Multilayered Rough Surfaces. Appl. Mech. Rev. 2002, 55, 435–480.

https://doi.org/10.1115/1.1488931.

2. Khadem, M.; Penkov, O.V.; Yang, H.; Kim, D. Tribology of multilayer coatings for wear reduction: A review. Friction 2017, 5,

248–262. https://doi.org/10.1007/S4054401701817.

3. Jaﬀar, M.J. Asymptotic behaviour of thin elastic layers bonded and unbonded to a rigid foundation. Int. J. Mech. Sci. 1989, 31,

229–235. https://doi.org/10.1016/0020740389901136.

4. Persson, B.N.J. Contact Mechanics for layered materials with randomly rough surfaces. J. Phys. Condens. Matter 2012, 24, 10.

https://doi.org/10.1088/0953-8984/24/9/095008.

5. Li, S.; Yuan, W.; Ding, Y.; Wang, G. Indentation load-depth relation for an elastic layer with surface tension. Math. Mech. Solids

2019, 24, 1147–1160. https://doi.org/10.1177/1081286518774090.

6. Westmann, R. Layered Systems Subjected to Asymmetric Surface Shears. Proc. R. Soc. Edinburgh Sect. A Math. Phys. Sci. 1963, 66,

140–149. https://doi.org/10.1017/S0080454100007792.

7. Fabrikant, V. Tangential contact problem for a transversely isotropic elastic layer bonded to an elastic foundation. J. Eng. Math.

2011, 70, 363–388. https://doi.org/10.1007/S1066501094184.

8. O’Sullivan, T.C.; King, R.B. Sliding Contact Stress Field Due to a Spherical Indenter on a Layered Elastic Half-Space. J. Tribol.

1988, 110, 235–240. https://doi.org/10.1115/1.3261591.

9. Wang, Q.; Sun, L.; Zhang, X.; Liu, S.; Zhu, D. FFT-Based Methods for Computational Contact Mechanics. Front. Mech. Eng. 2020,

6, 61. https://doi.org/10.3389/fmech.2020.00061.

10. Pohrt, R.; Popov, V.L. Adhesive Contact Simulation of Elastic Solids Using Local Mesh-Dependent Detachment Criterion in

Boundary Elements Method. Facta Univ. Ser. Mech. Eng. 2015, 13, 3–10.

11. Li, Q.; Popov, V.L. Boundary element method for normal non-adhesive and adhesive contacts of power-law graded elastic

materials. Comput. Mech. 2018, 61, 319–329. https://doi.org/10.1007/s00466-017-1461-9.

12. Chen, Z.; Etsion, I. Recent Development in Modeling of Coated Spherical Contact. Materials 2020, 13, 460.

https://doi.org/10.3390/ma13020460.

13. Shisode, M.P.; Hazrati, J.; Mishra, T.; de Rooij, M.B.; van den Boogaard, A.H. Semi-analytical contact model to determine the

flattening behavior of coated sheets under normal load. Tribol. Int. 2020, 146, 106182. https://doi.org/10.1016/j.tri-

boint.2020.106182.

14. Chen, Z.; Goltsberg, R.; Etsion, I. A universal model for a frictionless elastic-plastic coated spherical normal contact with mod-

erate to large coating thicknesses. Tribol. Int. 2017, 114, 485–493. https://doi.org/10.1016/j.triboint.2017.05.020.

15. Li, Q.; Pohrt, R.; Lyashenko, I.A.; Popov, V.L. Boundary element method for nonadhesive and adhesive contacts of a coated

elastic half-space. Eng. Tribol. 2019, 234, 73–83. https://doi.org/10.1177/1350650119854250.

16. Wang, Z.J.; Wang, W.Z.; Wang, H.; Zhu, D.; Hu, Y.Z. Partial Slip Contact Analysis on Three-Dimensional Elastic Layered Half

Space. ASME J. Tribol. 2010, 132, 021403. https://doi.org/10.1115/1.4001011.

Figure A1. Finite element model and mesh of the tangentially loaded coated half-space.

References

1. Bhushan, B.; Peng, W. Contact Mechanics of Multilayered Rough Surfaces. Appl. Mech. Rev. 2002,55, 435–480. [CrossRef]

Khadem, M.; Penkov, O.V.; Yang, H.; Kim, D. Tribology of multilayer coatings for wear reduction: A review. Friction

2017

5, 248–262. [CrossRef]

Jaffar, M.J. Asymptotic behaviour of thin elastic layers bonded and unbonded to a rigid foundation. Int. J. Mech. Sci.

1989

31, 229–235. [CrossRef]

Persson, B.N.J. Contact Mechanics for layered materials with randomly rough surfaces. J. Phys. Condens. Matter

2012

,24, 10.

[CrossRef]

Machines 2023,11, 694 15 of 15

Li, S.; Yuan, W.; Ding, Y.; Wang, G. Indentation load-depth relation for an elastic layer with surface tension. Math. Mech. Solids

2019,24, 1147–1160. [CrossRef]

Westmann, R. Layered Systems Subjected to Asymmetric Surface Shears. Proc. R. Soc. Edinburgh Sect. A Math. Phys. Sci.

1963

66, 140–149. [CrossRef]

Fabrikant, V. Tangential contact problem for a transversely isotropic elastic layer bonded to an elastic foundation. J. Eng. Math.

2011,70, 363–388. [CrossRef]

O’Sullivan, T.C.; King, R.B. Sliding Contact Stress Field Due to a Spherical Indenter on a Layered Elastic Half-Space. J. Tribol.

1988,110, 235–240. [CrossRef]

Wang, Q.; Sun, L.; Zhang, X.; Liu, S.; Zhu, D. FFT-Based Methods for Computational Contact Mechanics. Front. Mech. Eng.

2020

6, 61. [CrossRef]

10.

Pohrt, R.; Popov, V.L. Adhesive Contact Simulation of Elastic Solids Using Local Mesh-Dependent Detachment Criterion in

Boundary Elements Method. Facta Univ. Ser. Mech. Eng. 2015,13, 3–10.

11.

Li, Q.; Popov, V.L. Boundary element method for normal non-adhesive and adhesive contacts of power-law graded elastic

materials. Comput. Mech. 2018,61, 319–329. [CrossRef]

12. Chen, Z.; Etsion, I. Recent Development in Modeling of Coated Spherical Contact. Materials 2020,13, 460. [CrossRef]

13.

Shisode, M.P.; Hazrati, J.; Mishra, T.; de Rooij, M.B.; van den Boogaard, A.H. Semi-analytical contact model to determine the

flattening behavior of coated sheets under normal load. Tribol. Int. 2020,146, 106182. [CrossRef]

14.

Chen, Z.; Goltsberg, R.; Etsion, I. A universal model for a frictionless elastic-plastic coated spherical normal contact with moderate

to large coating thicknesses. Tribol. Int. 2017,114, 485–493. [CrossRef]

15.

Li, Q.; Pohrt, R.; Lyashenko, I.A.; Popov, V.L. Boundary element method for nonadhesive and adhesive contacts of a coated elastic

half-space. Eng. Tribol. 2019,234, 73–83. [CrossRef]

16.

Wang, Z.J.; Wang, W.Z.; Wang, H.; Zhu, D.; Hu, Y.Z. Partial Slip Contact Analysis on Three-Dimensional Elastic Layered Half

Space. ASME J. Tribol. 2010,132, 021403. [CrossRef]

17.

Jin, F.; Wan, Q.; Guo, X. Plane Contact and Partial Slip Behaviors of Elastic Layers with Randomly Rough Surfaces. ASME J. Appl.

Mech. 2015,82, 091006. [CrossRef]

18.

Komvopoulos, K. Finite Element Analysis of a Layered Elastic Solid in Normal Contact with a Rigid Surface. ASME J. Tribol.

1988

110, 477–485. [CrossRef]

19.

Popov, V.L.; Heß, M.; Willert, E. Handbook of Contact Mechanics—Exact Solutions of Axisymmetric Contact Problems; Springer-Verlag

GmbH: Heidelberg, Germany, 2018; pp. 127–141.

20.

Zhao, Y.; Liu, H.C.; Morales-Espejel, G.E.; Venner, C.H. Response of elastic/viscoelastic layers on an elastic half-space in rolling

contacts: Towards a new modeling approach for elastohydrodynamic lubrication. Tribol. Int. 2023,186, 108545. [CrossRef]

21. Cattaneo, C. Sul Contatto di due Corpore Elastici: Distribuzione degli sforzi. Rend. Dell’ Acad. Naz. Dei Lincei 1938,27, 342–348.

22. Mindlin, R.D. Compliance of elastic bodies in contact. J. Appl. Mech. 1949,16, 259–268. [CrossRef]

23. Jäger, J. A New Principle in Contact Mechanics. J. Tribol. 1998,120, 677–684. [CrossRef]

24. Ciavarella, M. Tangential Loading of General Three-Dimensional Contacts. J. Appl. Mech. 1998,65, 998–1003. [CrossRef]

25. Dassault Systems Simulia Corp. SIMULIA User Assistance; Dassault Systems Simulia Corp.: Providence, RI, USA, 2021.

Disclaimer/Publisher’s Note:

The statements, opinions and data contained in all publications are solely those of the individual

author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to

people or property resulting from any ideas, methods, instructions or products referred to in the content.