HAL Id: hal-00761201
https://hal.archives-ouvertes.fr/hal-00761201
Submitted on 16 Jan 2019
HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entic research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.
L’archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diusion de documents
scientiques de niveau recherche, publiés ou non,
émanant des établissements d’enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
An Ecient 3D Model of Heterogeneous Materials for
Elastic Contact Applications Using Multigrid Methods
H. Boy, Marie-Christine Baietto, P. Sainsot, Ton Lubrecht
To cite this version:
H. Boy, Marie-Christine Baietto, P. Sainsot, Ton Lubrecht. An Ecient 3D Model of Heterogeneous
Materials for Elastic Contact Applications Using Multigrid Methods. Journal of Tribology, American
Society of Mechanical Engineers, 2012, 134 (2), pp.021401. �10.1115/1.4006296�. �hal-00761201�
Hugo Boffy
1
Ph.D. Student
LaMCoS, CNRS UMR 5259,
Universite´ de Lyon, INSA-Lyon,
Villeurbanne, F69621, France
e-mail: hugo.boffy@insa-lyon.fr
Marie-Christine Baietto
CNRS Research Director
LaMCoS, CNRS UMR 5259,
Universite´ de Lyon, INSA-Lyon,
Villeurbanne, F69621, France
e-mail: marie-christine.baietto@insa-lyon.fr
Philippe Sainsot
Assistant Professor
LaMCoS, CNRS UMR 5259,
Universite´ de Lyon, INSA-Lyon,
Villeurbanne, F69621, France
e-mail: philippe.sainsot@insa-lyon.fr
Antonius A. Lubrecht
Professor
LaMCoS, CNRS UMR 5259,
Universite´ de Lyon, INSA-Lyon,
Villeurbanne, F69621, France
e-mail: ton.lubrecht@insa-lyon.fr
An Efficient 3D Model
of Heterogeneous Materials
for Elastic Contact Applications
Using Multigrid Methods
A 3D graded coating/substrate model based on multigrid techniques within a finite differ-
ence frame work is presented. Localized refinement is implemented to optimize memory
requirement and computing time. Validation of the solver is performed through a com-
parison with analytical results for (i) a homogeneous material and (ii) a graded material.
The algorithm performance is analyzed through a parametric study describing the influ-
ence of layer thickness (0.01 < t/a < 10) and mechanical properties (0.005 < E
c
/E
s
< 10)
of the coating on the contact parameters (P
h
, a). Three-dimensional examples are then
presented to illustrate the efficiency and the large range of possibilities of the model. The
influence of different gradations of Young’s modulus, constant, linear and sinusoidal,
through the coating thickness on the maximum tensile stress is analyzed, showing that the
sinusoidal gradation best accommodates the property mismatch of two successive layers.
A final case is designed to show that full 3D spatial property variations can be accounted
for. Two spherical inclusions of different size made from elastic solids with Young’s mod-
ulus and Poisson’s ratio are embedded within an elastically mismatched finite domain and
the stress field is computed.
1 Introduction
The actual contact between solid surfaces is generally rough
and the stresses induced can only be correctly described using a
detailed 3D model. Current manufacturing processes allow one to
obtain functional graded materials, i.e., materials with 3D varying
properties. A traditional application is a coated material. Never-
theless, the rough coated contact problem is strongly multi-scale:
the characteristic dimensions of the contact, the coating and the
roughness range from the millimeter to the nanometer. A straight-
forward discretization of this multi-scale problem would exceed
the memory and CPU capacity of current (and next generation)
computers. This paper proposes an efficient numerical model that
can handle such multi-scale problem requiring O(10
9
) points and
locally refined grids.
The development and coating selection becomes a complex and
costly task for designers. Numerical simulations can thus guide
the choice of the designer and reduce the cost of an experimental
approach. Substrate and coating are generally elastically dissimi-
lar; thus the contact stresses and the internal stresses will differ
from those generated in the uncoated case and quantification is
necessary. The actual stress distribution is mostly calculated using
numerical or semi-numerical techniques for layered materials.
These numerical simulations can become very difficult when
varying properties and thin coatings are involved. These computa-
tions can aid in fine tuning the material choice to the specific oper-
ating conditions, such as layer thickness, property variations, etc.
Several studies have been performed on this topic and will be out-
lined later on.
Most of the models consider discrete isotropic layers with con-
stant properties and a sharp transition at the interface. Moreover,
perfect bounds at interfaces are assumed and this leads to inherent
stress discontinuities at the interface. The methods used are gener-
ally based on integral transform combined with a Fast Fourier
Transform (FFT) algorithm or on Finite Element (FE) methods.
Ju and Liu [
1] and Leroy et al. [2] used the Fourier transform to
study the thermo-mechanical deformation of a 2D layered elastic
half space subjected to a moving heat source. Plumet and Dubourg
[
3] used these techniques to study the 3D mechanical case. FFT
methods introduce a numerical periodicity error which can vanish
if one extends the surface grid sufficiently far beyond the contact
area. In their work, Polonsky and Keer [
4,5] remedy the periodic-
ity error by using a special correction term. Nogi and Kato [
6]
combined the conjugate gradient method with the FFT technique
for solving the rough contact problem for both homogeneous and
layered solids. FFT methods can deal with coated materials but
are not appropriated for graded ones when property variations do
not follow a single space direction.
Kesler et al. [
7,8] worked on experimental methods to deter-
mine the residual stresses and the intrinsic stresses resulting from
the deposition process. They also studied the Young’s modulus of
a layered or graded material as a function of depth in the deposit.
Giannakopoulos and Suresh [
9] developed semi-analytical models
to adress graded materials with properties varying either as expo-
nential or power laws. More recently, Prasad et al. [
10] studied
the steady-state frictional sliding contact on surfaces of plastically
graded materials using FE techniques.
In 3D mechanical studies, simulations are mostly performed
using FE techniques [
11,12]. Holmberg et al. [11] used these
methods to calculate the real stress components and fracture
toughness of a coated surface and to investigate the effect of resid-
ual stress in the coating. This approach presents a major disad-
vantage for thin layers, which is the large system of equations
arising from the fine mesh. Conventional numerical methods take
unacceptably long times even on modern computers. The compu-
tation speed can be increased by using advanced numerical techni-
ques such as multigrid methods [
13].
Semi-analytical methods have also been used in the case of
graded and coated materials to determine displacement and stress
fields. Fretigny et al. [
14] developed an analytical approach to the
mechanics of adhesive contacts between coated substrates and
1
rigid axisymmetric probes as an extension of a model developed
by Perriot [
15] for the non-adhesive case. Giannakopoulos and
Suresh [
9] studied the stresses in the case of a graded material
with an exponential variation of its Young’s modulus in depth.
The methods used cannot deal with all types of variations.
Watremetz et al. [
16] developed a 2D thermomechanical nu-
merical model based on multigrid methods to solve mechanical
problems including both coated and graded materials. His model
is able to solve problems with any type of property variation. He
particularly showed how a graded variation avoids the discontinu-
ity in the stress components.
Starting from Ref. [
16], extensions were realized to propose a
3D multigrid model for homogeneous bodies submitted to a mov-
ing heat source [
17]. Then a 3D multigrid elastic contact model
for coated and graded materials with 3D spatial property varia-
tions has been developed. This work is based on the finite differ-
ence method which allows one to account for property variations
in all space dimensions. A local refinement strategy is used to
reduce computation time.
The model is presented in the first section. The validation is
outlined in Sec.
2. Comparisons are performed with classical
Hertzian solutions given by Johnson [
18] and analytical solutions
proposed by Giannakopoulos and Suresh [
9]. In Sec. 3, a paramet-
ric study of a coated substrate is conducted. A large range of coat-
ing thickness and E
c
/E
s
ratio is considered to illustrate the
efficiency and robustness of the model. Three-dimensional numer-
ical examples are then presented in Sec.
4 to illustrate the effi-
ciency and the large range of possibilities of the model. Coating/
substrate systems with different property gradations through the
coating thickness, constant, linear and sinusoidal, are submitted to
contact loading. The resulting stress fields within the coating and
the substrate and at the interface are analyzed. A final case is then
designed to show that full 3D spatial property variations can be
accounted for. Two spherical inclusions of different size made
from different materials (a compliant and a stiff one) are embed-
ded within an elastically mismatched finite domain. The stress
field is computed.
2 3D Model and Numerical Method
3D Numerical tools dealing with coating problems are gener-
ally based on Finite Element Methods (FEM) or on semi-
analytical techniques combined with Fast Fourier Transform
(FFT). FEM demand a huge computational effort (CPU and mem-
ory) if one requires locally very fine meshes. FFT methods cannot
deal with three-dimensional varying mechanical properties. A fi-
nite difference discretization can handle such three-dimensional
varying properties, but requires rather important computational
efforts. Using multigrid methods (MG) the convergence speed can
be accelerated and large scale problems can be solved with a lim-
ited computational effort. Moreover, these methods can easily be
coupled with mesh refinement strategies allowing one to deal with
the different contact scales [
13,19].
The first to apply these techniques to a 2D elastic model was
Watremetz et al. [
16]. He studied a material with graded trans-
verse elastic properties.
2.1 Equations. The work is based on a second order finite
difference formulation of 3D Lame´ equations, written for a finite
rectangular domain. The Lame´ coefficients can vary as a function
of space, i.e., k ¼f(x,y,z), l ¼g(x,y,z). Hence the Young’s modu-
lus and Poisson’s ratio of this elastic solid are obviously 3D spa-
tial varying functions (E(x,y,z),(x,y,z)). Thus any property
gradation can be accounted for through the depth, but also along x
and y directions. Applications like embedded inclusions within a
mismatched elastic solid can be addressed as it will be shown in
Sec. 4. These equations can be generalized in terms of displace-
ments (u(x,y,z)) to:
ðku
j; j
Þ
; i
þðlu
i; j
Þ
; j
þðlu
j; i
Þ
; j
¼ 0 i; j ¼ 1; 2; 3
(1)
The partial differential equations PDE are discretized using a sec-
ond order finite difference scheme in the bulk except on the boun-
daries. One of the originalities of this work is that the cross
derivative terms are taken into account as
@
@x
k
@u
@x
due to the prop-
erty variations along all spatial dimensions.
Displacements are determined using boundary conditions (BC).
Dirichlet BC (displacements) are imposed on the five boundaries
and Neuman BC prescribe the load on the top free surface accord-
ing to the following equation:
r
zz
¼ k
@w
@z
þ
@v
@y
þ
@u
@x
r
zy
¼ l
@w
@y
þ
@v
@z
r
zx
¼ l
@w
@x
þ
@u
@z
8
>
>
>
>
>
>
>
<
>
>
>
>
>
>
>
:
(2)
2.2 Multigrid Methods. Conventional iterative methods
converge slowly as the long wavelength error present in the solu-
tion slowly diminishes with successive relaxation sweeps. Multi-
grid techniques overcome this difficulty by using a sequence of
grids. The grids communicate with one another through restriction
and prolongation operators [
20,21]. Efficient multigrid algorithms
for the solution of PDE require good ellipticity, which implies
that non-smooth error components can be solved by local process-
ing. The elasticity equations are elliptic and one can expect a
good multigrid efficiency. The treatment of the Dirichlet BC is
straight forward, however, Neumann BC can seriously degrade
the solver efficiency.
2.3 Boundary Conditions. The boundary conditions should
not modify the efficiency of the multigrid techniques. The Dirich-
let conditions are easily implemented. Concerning the Neumann
ones, in order to maintain the multigrid solver efficiency, the BC
have to fulfil two criteria:
•
consistency on the various grids
•
no introduction of local errors.
Venner and Lubrecht [
13] relax the interior problem on the sur-
face using the Neumann condition to fulfill both criteria. Some
pre-relaxations are necessary to decrease the residual near the
boundaries before solving the entire problem, as described by
Brandt [
21].
2.4 Refinement Strategy. Coated and graded systems are
multi-scale problems as a large difference between coating thick-
ness, loaded zone and the substrate dimensions exists. A uniform
level of resolution is not required everywhere and a proper spatial
resolution, fine for instance around the coating and the interface
and coarse elsewhere, allows a considerable reduction of the com-
putational work. This is especially true if one considers a coated
semi-infinite bulk (bulk size coating thickness and/or contact
width). MG techniques are easily combined with local mesh
refinements techniques. In such cases, the finest grids can be re-
stricted to smaller and smaller subdomains, whereas the coarse
grids cover the entire domain.
The refinement strategy consists in defining successive subdo-
mains corresponding to successive levels where grid refinement is
performed only locally (Fig.
1). The fine grid residuals and the
coarse grid correction are only carried out in the local part where
the fine grid exists [
19]. The refinement strategy is centered on M
(0,0,0). A refined level corresponds to a reduction of the coarsen
volume with a ratio of 3/4 in the x and y directions and 1/2 in the
z direction. The 1/2 z-ratio allows one to maintain a constant num-
ber of grid-points along the z direction. In order to show the effi-
ciency of the refinement strategy, a simple case is studied. A
Hertzian loading is imposed at the surface of a rectangular do-
main. Numerical simulations are conducted with a different
2
number of levels of local refinement, from a one-level to six-level
local refinements. The case considered as the reference includes
six grid levels and no local refinement (number of grid-points on
the finest grid around 67M and u
ref
the reference solution). This
number of grid-points reduces drastically, down to 14K for the
five-level local refinement case. The difference between the refer-
ence solution u
ref
and the solution u
app
obtained using p-level
local grid refinement and the corresponding reduction in computa-
tional time, are summarized in Table
1. A maximum and a mean
displacement error norm are defined respectively as:
u
ref
u
app
1
¼ max
u
ref
u
app
u
ref
u
ref
u
app
1
¼
1
N
X
N
0
u
ref
u
app
u
ref
The maximum error increases with increasing number of local
refined levels, but the accuracy remains pretty good. The time
reduction increases very significantly with the increasing number
of local refined levels and the calculation can be performed up to
300 times faster with almost the same accuracy in this case. These
results show the ability of the model to deal with local refinement
techniques combined with MG. This was required to address large
3D problems requiring high resolution sub-domains.
3 Validation
The next step consists in validating the model and illustrating
its accuracy. Three problems are considered and the numerical
results are compared with analytical solutions to validate (i) the
multigrid model, (ii) the contact solution within the multigrid
framework, and (iii) the description of coatings with elastic prop-
erty gradients through the depth. Calculations based on multigrid
techniques were performed using different number of levels
without local refinement, ranging from 1 (monogrid) to 7 grid
levels. A rectangular domain is considered of length L
x
/a ¼L
y
/
a ¼2
*
L
z
/a ¼16. The coarse grid involves 9*9*5 points and the
finest one 513*513*257 points with a mesh size of 1/32.
3.1 Multigrid Model Validation. A homogeneous bulk sub-
mitted to an Hertzian loading is considered. The solution can be
found [
18] in terms of surface displacement (i.e., z ¼0).
In the contact zone (
r 1):
wð
rÞ¼
1
4
ð2
r
2
Þ
(3)
uð
rÞ¼
ð1 2Þ
ð1 Þ
1
3p
r
ð1 ð1
r
2
Þ
3=2
Þ
x
r
(4)
vð
rÞ¼
ð1 2Þ
ð1 Þ
1
3p
r
ð1 ð1
r
2
Þ
3=2
Þ
y
r
(5)
In the center of the contact (
r ¼ 0):
r
xx
¼ r
yy
¼ð1 þ 2Þ=2
(6)
On the edge of the contact (
r ¼ 1):
r
xx
¼ð1 2Þ=3
(7)
r
yy
¼
r
xx
(8)
With :
x ¼ x=a
0
;
y ¼ y=a
0
;
r ¼
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
x
2
þ
y
2
p
,
u ¼ u=d
0
;
v ¼ v=d
0
;
w ¼ w=d
0
,
r
xx
¼ r
xx
=P
h0
;
r
yy
¼ r
yy
=P
h0
.
The mean error is defined as g ¼ w
ana
w
app
1
, where w
ana
and w
app
are, respectively, the analytical and the approximate nu-
merical solutions of the vertical displacement calculated at each
surface point. This error is plotted versus the mesh size in Fig.
2.
This curve shows a slope of 3.5 which is close to the expected
second order convergence speed value of 4. The difference
Fig. 1 Local refinement technique: (a) coarsen grid, (b) one-level local refinement,
(c) two-level local refinement
Table 1 Computational time reduction and accuracy of the
solver versus the number of p-level local refinements
p u
ref
u
app
1
u
ref
u
app
1
Time reduction
1 6.8 * 10
5
3.6 * 10
5
3.0
2 2.8 * 10
4
1.6 * 10
5
10.0
3 6.2 * 10
4
3.2 * 10
4
34.5
4 1.5 * 10
3
5.0 * 10
4
142.8
5 2.1 * 10
3
6.5 * 10
4
384.6
3
is probably due to the stress singularity at the edge of the
contact. Our results are in a very good agreement with the exact
solution. An accuracy better than 10
–3
is obtained using a mesh
size of 1/32.
3.2 Contact Solution Using MG. The numerical relation
between the displacement w at point M due to a point loading at a
different point N gives the influence coefficients which can be
applied in a contact algorithm [
22]. This algorithm has success-
fully addressed rough homogeneous contact problems. Combining
this algorithm with MG numerical coefficients allows us to
numerically solve contact problems involving materials with vary-
ing properties throughout the depth. A spherical indenter pressed
normally against a half-space is considered. The contact pressure
is numerically determined and compared to the analytical Hertz-
ian solution. The relative error w between the maximum Hertzian
pressure P
h0
and the numerical one is plotted as a function of the
number of grid points in Fig.
3. An accuracy of 10
3
is reached
for a grid spacing smaller than 1/32. This result validates the use
of the numerical influence coefficients within the MG framework.
3.3 Contact Solution for Coatings With Graded Elastic
Properties Through the Depth. Giannakopoulos and Suresh [
9]
gave approximate analytical solutions for both the dimensionless
pressure distribution and contact half-width for a spherical rigid
(parabolic) indenter pressed against a material with a specific
Young’s modulus variation with depth. E is varying according
to an exponential law: EðzÞ¼pe
az
, with
a ¼ a
0
a and
z ¼ z=a
0
;¼ 0.
pðr Þ
ffiffiffiffiffiffiffiffiffiffiffiffiffi
1
r
2
p
for 0 r a
(9)
a
3
1 þ
C
1
a
3
3
1
(10)
With
pðr Þ¼pðrÞ=p
h0
and
a ¼ a=a
0
. Note that a ¼0 corresponds
to the Hertzian case. The C
1
value is determined and stability
issues concerning the exponential law case are described in
Ref. [
23].
For these conditions, MG calculations were performed for
0:25
a 0:25. The maximum contact pressure and contact
half-width are compared to the reference values given by equa-
tions
(9) and (10) and plotted in Figs. 4 and 5.
The error is of the order of one percent independently of the
value of
a. It is thus shown that the agreement between analytical
and numerical solutions is excellent, validating, therefore, the
ability of the MG model to address materials with property grada-
tions versus depth. The next section is devoted to the analysis of
the influence of the in-depth gradation on the internal stress field
with a particular interest to the coating/substrate.
4 Results
In Fig.
6, a half cube of finite dimensions L
x
¼L
y
¼2
*
L
z
is
shown. The Poisson’s ratio is
s
¼
c
in the substrate and in the
coating. The Young’s modulus in the substrate is also assumed
constant and its value is E
s
. The layer is either considered with
constant Young’s modulus throughout the thickness (called a
coating) leading to a discontinuity in material constant across the
coating/substrate interface, or with a specific gradient in E
c
(lin-
ear, sinusoidal…) resulting in a continuous variation in material
constant across the interface. A spherical indenter is pressed on
the domain with a certain load Q. The contact problem is solved
numerically using the MG numerical influence coefficients and
the contact algorithm. In the first subsection, coated materials are
Fig. 2 Surface displacement error
Fig. 3 Relative error in the maximum Hertzian pressure
Fig. 4 Comparison between analytical (GS) and numerical
(MG) dimensionless central pressure
p (0,0) using the exponen-
tial law
Fig. 5 Comparison between analytical (GS) and numerical
(MG) dimensionless contact radius
a using the exponential law
4
