# A 3-D model for a multilayered body loaded normally and tangentially against a rigid body: Application to specific coatings

Sylvie Plumet, Marie-Christine Baietto

► To cite this version:

Sylvie Plumet, Marie-Christine Baietto. A 3-D model for a multilayered body loaded normally and tangentially against a rigid body: Application to specific coatings. Journal of Tribology, 1998, 120 (4), pp.668-676. 10.1115/1.2833764. hal-01951927

HAL Id: hal-01951927

https://hal.science/hal-01951927

Submitted on 7 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific 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 diffusion de documents
scientifiques 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A 3-D Model for a Multilayered Body Loaded Normally and Tangentially Against a Rigid Body: Application to Specific Coatings

S. Plumet

M.-C. Dubourg

Laboratoire de Mécanique des Contacts,
UMR CNRS 5514, INSA de Lyon, Bat 113,
20 Avenue Albert Einstein,
69621 Villeurbanne Cedex, France

Coatings are increasingly used to improve the mechanical and tribological behavior of surfaces. It is necessary to develop models to guide the initial choice of coating/substrate combinations that can withstand the applied loads. A three-dimensional model of an elastic multilayered body, loaded both normally and tangentially against an elliptical rigid body (partial sliding, rolling/sliding conditions), is presented here. This model is based on linear elasticity theory, integral transforms, Fast Fourier Transform, and unilateral contact analysis with friction. Normal and tangential contact conditions between the two bodies are first determined and then used to calculate the multilayered body stress field. One application is given here: The influence of the mechanical properties of coating and substrate, as well as coating thickness, is studied on contact conditions, internal stresses, and potential failure mechanisms.

## 1 Introduction

In many engineering applications, protective coatings are increasingly used to extend the fatigue life of mechanical components in contact (Chang et al., 1990) and to provide low friction coefficients and wear resistance of tribological surface properties (Komvopoulos, 1987). A great deal of research is underway to evaluate the performance of coatings and their failure mechanisms. The analysis of the contact and internal stresses in a multilayered body is therefore of great practical and analytical importance. The purpose of this paper is to present a model capable of solving the normal and tangential contact problem and to determine the internal stresses for a 3-D multilayered body configuration.

At a relatively early stage in the history of continuum mechanics, Hertz (1880) was the first to solve the contact problem of frictionless elastic homogeneous bodies under normal loading. Using the same model, Cattaneo (1938) and Mindlin (1949) treated the case in which the bodies are subjected to both tangential and normal forces. Kalker (1982) proposed a method based on a Simplex algorithm (Conry and Seireg, 1971) to solve the 3-D contact problem under rolling conditions. However, none of methods take into account the finite thickness of layered bodies, nor can then be used to guide the choice of coatings.

Only the use of numerical solutions is possible for layered configurations. The case of an elastic layer overlaying a rigid substrate was treated by Hannah (1951) and Bentall and Johnson (1968), using, respectively, Coker and Filon general stress functions, and a Fourier transform of Sneddon (1951) stress functions. Chen (1971) was one of the first to solve the problem of an elastic layer perfectly bonded to an elastic substrate. Gupta and Walowitt (1973) proposed a general solution for the plane-strain contact problem of layered solids. The latter authors showed that contact pressure distribution may deviate significantly from that obtained previously by Hertz, depending on

the finite thickness and the elastic properties of the layer and the substrate. In 1988, Komvopoulos using the Finite Element Method (FEM), presented a 2-D contact solution under complete sliding conditions for a stiff layer over an elastic substrate. He studied the influence of friction coefficient on contact and internal stresses for different coating thickness and also defined a minimal value of this thickness to minimize Von Mises's maximal value within the substrate. Leroy et al. (1989) used an original half analytical and half numerical approach to analyze the thermomechanical behavior of a multilayered half plane. They used a Fast Fourier Transform algorithm (Brigham, 1974) that avoids singularity problems encountered during numerical integration of the inverse transform, and reduces calculation time. Their application focused on an elastic layer bonded to an elastic substrate, subjected to both load and temperature at the contact surface.

Naturally, the problem has been extended to 3-D layered configurations. Goodman and Keer (1975) proposed a solution for an elastic layer overlaying a rigid substrate. Chiu and Hartnett (1983) solved the normal problem of a sphere over an elastic substrate using Hankel transforms. O'Sullivan and King (1988) presented a new solution for a layered half space under sliding conditions. More recently, Kuo and Keer (1992) investigated a method based on Hankel transforms to solve the contact problem between a spherical indenter and a multilayered transversely isotropic structure bonded to an elastic half space. Very recently, Nogi and Kato (1997), using Fast Fourier Transforms, have extended O'Sullivan and King's method to formulate a contact solution in the Fourier domain for rough surfaces of a layered half space. Ju and Farris (1997) use the Fourier transform to study thermomechanical problems in a spatial Fourier transform domain.

Although these theories take into account a layer over an elastic substrate, there is a lack of models for solving the elliptical contact problem under normal and tangential loading conditions, and furthermore, to analyze multilayered body behavior. In this paper, a rigid ellipsoid over a smooth elastic multilayered body is considered. The method is based on integral transforms and a Fourier algorithm to avoid singularity problems and to reduce calculation time. This model has the advantage, com-
pared to previous 3-D methods, of accounting for multicoated bodies subjected to any loading conditions and not to only complete sliding conditions. As a first step, results under normal loading conditions are presented for an ellipsoid over a layered half-space. They are compared to previous 3-D methods encountered in the literature.

## 2 Theoretical Model

The aim is to solve the normal and tangential contact problems between a rigid ellipsoid and the first layer of an elastic multilayered body (Fig. 1). Pressure and traction distributions at the two body contact surfaces are then used as data to calculate the subsurface stress field. The ellipsoid is defined by two radii of curvatures ($Rx_1$ and $Rx_2$), and the multilayered body by the number of layers ($n$) (Fig. 1). Each layer ($i$) is isotropic, linearly elastic, smooth and homogeneous, and is characterized by its mechanical properties ($E_i$, $v_i$) and finite thickness ($e_i$). The interfacial boundary conditions between layers $i$ and $i+1$ are either a perfect bonding (Fig. 2(a)), an imperfect bonding or a complete sliding condition (Fig. 2(b)). Boundary conditions in terms of nil stresses or nil displacements (Fig. 2(c)) can be applied at the lower border of the last layer (layer $n$).

### 2.1 Contact Problem and Stress Field Formulation.

The contact solution consists in determining normal and tangential pressure distributions that satisfy the boundary conditions inside and outside the a priori unknown contact area. First, a potential contact area is defined. It is discretized into cells where unit and uniform pressure and tractions may act. A new approach for the influence coefficient formulation is proposed here. It is based on the original method founded on Fourier techniques, already used by Leroy et al. (1989) for 2-D solution.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Geometry of 3-D multilayered body</figcaption></figure>

The contact problem is then solved classically. Finally the internal stress field is calculated.

#### 2.1.1 Influence Coefficient Formulation.
The first step in the contact solution is to establish the surface displacements at a point M due to uniform pressure and/or traction centered at point N (see Fig. 3). The basic concept of this 3-D influence coefficient formulation is based on the Fourier transform with respect to space variables $x_1$ and $x_2$. In practice, this is done by the Fast Fourier Transform (FFT) technique (Brigham, 1974). The FFT algorithm presents the advantage of computing the discrete Fourier transform much more rapidly than other avail-

### Nomenclature

$a_{pm}(i, j, k, l)$ = influence coefficient term representing deflection along direction $x_p$ at node $(k, 1, 0)$ due to unit pressure or traction along direction $x_m$ centered at node $(i, j, 0)$ (mm·MPa⁻¹)

$a, b$ = semi-axes of the contact area along directions $x_1$ and $x_2$, respectively, (mm)

$A_j, B_j, C_j, \bar{D}_j, \bar{E}_j, F_j$ = integral coefficients of layer $j$

$d(k, 1, 0)$ = distance between ellipsoid and layer 1 along axis $x_3$ (mm)

$\text{DEP}_i^+, \text{DEP}_i^-$ = displacement field at node $(k, l, m)$ at the lower $(-, x_3 = \sum e_k)$ or upper $(+, x_3 = \sum e_k)$ face of layer $i$ (mm)

$[DS_j], [SD_j], [DD_j], [SS_j]$ = matrices $[3 \times 3]$ relating displacement and stress fields at the lower and upper faces of layer $j$ ($[DS]$, MPa⁻¹, $[SD]$, MPa)

$e_j$ = finite thickness of layer $j$ (mm)

$f_r, g_s$ = Fourier transform frequencies (along $x_1$ and $x_2$ axes) (mm⁻¹)

$[\text{Infl}]$ = influence coefficient matrix $[3 \times 3]$ (mm·MPa⁻¹)

$K_1$ = calculation constant ($K_1 = (\lambda + 3\mu)/(\epsilon_1(\lambda + \mu))$) (mm)

$P$ = normal force (N)

$p(i, j, 0)$ = pressure distribution at node $(i, j, 0)$ (MPa)

$Rx_1, Rx_2$ = ellipsoid radii (mm)

$s\dot{x}_i(i, j)$ = velocity at node $(i, j, 0)$ along direction $x_i$ (mm/s)

$\text{SIG}_i^+, \text{SIG}_i^-$ = stress field at node $(k, l, m)$ at the lower $(-)$ or upper $(+)$ border of layer $i$ (MPa)

$tx_1(i, j, 0)$,

$tx_2(i, j, 0)$ = tangential pressure distribution at node $(i, j, 0)$ along $x_1$ or $x_2$ axis (MPa)

$Tx_1, Tx_2$ = tangential forces along $x_1$ or $x_2$ axis (N)

$u_i(k, l, m) =$

$u_i(m)$ = displacement in layer $j$ along $x_i$ axis at node $(k, l, m)$ (mm)

$u_i^{\text{ell}}(k, l, m)$ = displacement of the ellipsoid along $x_i$ axis at node $(k, l, m)$ (mm)

$V$ = rolling velocity

$\bar{x}$ = $x$ transformed value

$\lambda, \mu$ = Lamé's coefficient

$\beta$ = calculation constant ($\beta = \sqrt{(\epsilon_1^2 + \epsilon_2^2)}$) (mm⁻¹)

$\sigma_{ij}^m(x_3)$ = stress field within layer $m$ at node $(k, l, x_3)$ (MPa)

$\epsilon_1, \epsilon_2$ = Fourier transform parameters ($\epsilon_1 = 2i*\pi f_r, \epsilon_2 = 2i*\pi g_s$) (mm⁻¹)

$\delta_m = \delta_m^{\text{ell}} - \delta_m^1$ = rigid displacement related to ellipsoid and layer 1 displacements along direction $x_m$ (mm)

$\mu_f$ = coefficient of friction at the interface of the two contacting bodies (ellipsoid and layer 1).

$\xi x_i$ = creep ratio along direction $x_i$

$\Delta = \Delta = (\partial^2/\partial x_1^2) + (\partial^2/\partial x_2^2) + (\partial^2/\partial x_3^2)$
<figure><img src="image_3.png" /><figcaption>Fig. 2 Interfacial boundary conditions</figcaption></figure>

able algorithms, without singularity problems in the inverse transform.

Equilibrium equations (Eq. (1)), stress-displacement relations (Eqs. (2-3)) and biharmonic equations (Eq. (4)) are written for the layer k (Timoshenko and Goodier, 1970):

$$ (\lambda + 2\mu) \frac{\partial^2 u_i^k}{\partial x_i^2} + \mu \left( \sum_{j \neq i} \frac{\partial^2 u_j^k}{\partial x_j^2} \right) + (\lambda + \mu) \left( \sum_{j \neq i} \frac{\partial^2 u_i^k}{\partial x_j \partial x_i} \right) = 0 \quad \text{for } i=1,3; j=1,3 \quad (1) $$

$$ \sigma_{ii}^k = (\lambda + 2\mu) \frac{\partial u_i^k}{\partial x_j} + \lambda \sum_{j \neq i} \frac{\partial u_j^k}{\partial x_j} \quad \text{for } i=1,3; j=1,3 \quad (2) $$

$$ \sigma_{ij}^k = \mu \left( \frac{\partial u_i^k}{\partial x_j} + \frac{\partial u_j^k}{\partial x_i} \right) \quad \text{for } i \neq j \quad (3) $$

$$ \Delta\Delta u_i^k = 0 \quad \text{for } i=1,3 \quad (4) $$

The application of an integral transform to (1-4) has the effect of removing partial derivatives with respect to the space variables considered. The space variables considered here are $x_1$ and $x_2$. The Fourier integral transform associated with the $x_1$ and $x_2$ variables for a function $h(x_1, x_2, x_3)$ is defined by:

$$ \overline{h(f_r, g_s, x_3)} = \iint_{-\infty \to \infty} h(x_1, x_2, x_3) \times \exp(-i2\pi(f_r x_1 + g_s x_2)) dx_1 dx_2 \quad (5) $$

Equations (1) and (4) are Fourier transformed to obtain Eqs. (6)-(9) (i is a pure imaginary, $\epsilon_i$, the frequency in the Fourier domain along direction $x_i$).

$$ -(\lambda_k \epsilon_1^2 + \mu_k(2\epsilon_1^2 + \epsilon_2^2)) \overline{u_1^k} + \mu_k \frac{d^2 \overline{u_1^k}}{dx_3^2} \\ -(\lambda_k + \mu_k) \epsilon_1 \epsilon_2 \overline{u_2^k} + i\epsilon_1(\lambda_k + \mu_k) \frac{du_3^k}{dx_3} = 0 \quad (6) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3 Contact problem formulation</figcaption></figure>

$$ -(\lambda_k \epsilon_2^2 + \mu_k(2\epsilon_2^2 + \epsilon_1^2)) \overline{u_2^k} + \mu_k \frac{d^2 \overline{u_2^k}}{dx_3^2} \\ -(\lambda_k + \mu_k) \epsilon_1 \epsilon_2 \overline{u_1^k} + i\epsilon_2(\lambda_k + \mu_k) \frac{du_3^k}{dx_3} = 0 \quad (7) $$

$$ i(\lambda_k + \mu_k) \left( \epsilon_1 \frac{du_1^k}{dx_3} + \epsilon_2 \frac{du_2^k}{dx_3} \right) \\ - \mu_k \beta^2 \overline{u_3^k} + (\lambda_k + 2\mu_k) \frac{d^2 \overline{u_3^k}}{dx_3^2} = 0 \quad (8) $$

$$ \frac{d^4 \overline{u_1^k}}{dx_3^4} - 2\beta^2 \frac{d^2 \overline{u_1^k}}{dx_3^2} + \beta^4 \overline{u_1^k} = 0 \quad \text{for } i=1,3 \quad (9) $$

Equations (6)-(9) are then solved to determine general displacement expressions (Eqs. (10)-(12)) in the frequency domain (i pure imaginary):

$$ \overline{u_1^k}(\epsilon_1, \epsilon_2, x_3) = i((A_k + B_k x_3) \exp(\beta x_3) + (C_k + D_k x_3) \exp(-\beta x_3)) \quad (10) $$

$$ \overline{u_2^k}(\epsilon_1, \epsilon_2, x_3) = i \left( \left( E_k + \frac{\epsilon_1}{\epsilon_2} B_k x_3 \right) \exp(\beta x_3) + \left( F_k + \frac{\epsilon_1}{\epsilon_2} D_k x_3 \right) \exp(-\beta x_3) \right) \quad (11) $$

$$ \overline{u_3^k}(\epsilon_1, \epsilon_2, x_3) = \left( \frac{A_k \epsilon_1}{\beta} + \frac{E_k \epsilon_2}{\beta} + K_1 B_k + \frac{B_k \beta x_3}{\epsilon_1} \right) \exp(\beta x_3) \\ + \left( -\frac{C_k \epsilon_1}{\beta} - \frac{F_k \epsilon_2}{\beta} + K_1 D_k - \frac{D_k \beta x_3}{\epsilon_1} \right) \exp(-\beta x_3) \quad (12) $$

where $\beta = \sqrt{(\epsilon_1^2 + \epsilon_2^2)}$

These expressions depend on six integral constants ($A_k$, $B_k$, $C_k$, $D_k$, $E_k$, and $F_k$) determined by the boundary conditions at the lower and upper faces of layer k (Fig. 4). Displacement and stress expressions are then obtained in matrix form (relations 13 and 14), using the stress-displacement relations and displacement expressions in the frequency domain.

$$ \overline{\text{DEP}}_k^+ = [DS_k] \overline{\text{SIG}}_k^+ + [DD_k] \overline{\text{DEP}}_k^- \quad (13) $$

$$ \overline{\text{SIG}}_k^- = [SS_k] \overline{\text{SIG}}_k^+ + [SD_k] \overline{\text{DEP}}_k^- \quad (14) $$

where:

$$ \overline{\text{DEP}}_k^+ = \left[ \frac{\overline{u_1^k}(\epsilon_1, \epsilon_2, 0)}{i}, \frac{\overline{u_2^k}(\epsilon_1, \epsilon_2, 0)}{i}, \frac{\overline{u_3^k}(\epsilon_1, \epsilon_2, 0)}{i} \right]^T, $$

$$ \overline{\text{DEP}}_k^- = \left[ \frac{\overline{u_1^k}(\epsilon_1, \epsilon_2, e_k)}{i}, \frac{\overline{u_2^k}(\epsilon_1, \epsilon_2, e_k)}{i}, \frac{\overline{u_3^k}(\epsilon_1, \epsilon_2, e_k)}{i} \right]^T, $$

<figure><img src="image_5.png" /><figcaption>Fig. 4 Boundary conditions of layer k</figcaption></figure>
$$
\overline{\text{SIG}}_{k}^{+} = \left[ \frac{\overline{\sigma}_{13}^{k}(\epsilon_{1}, \epsilon_{2}, 0)}{i}, \frac{\overline{\sigma}_{23}^{k}(\epsilon_{1}, \epsilon_{2}, 0)}{i}, \frac{\overline{\sigma}_{33}^{k}(\epsilon_{1}, \epsilon_{2}, 0)}{i} \right]^{T},
$$

$$
\overline{\text{SIG}}_{k}^{-} = \left[ \frac{\overline{\sigma}_{13}^{k}(\epsilon_{1}, \epsilon_{2}, e_{k})}{i}, \frac{\overline{\sigma}_{23}^{k}(\epsilon_{1}, \epsilon_{2}, e_{k})}{i}, \frac{\overline{\sigma}_{33}^{k}(\epsilon_{1}, \epsilon_{2}, e_{k})}{i} \right]^{T}
$$

Exponents + and – refer, respectively, to the upper and lower
faces of layer *k*. Matrices [DS*<sub>*k*</sub>], [DD*<sub>*k*</sub>], [SS*<sub>*k*</sub>], and [SD*<sub>*k*</sub>]
depend on the finite thickness of layer *k* (*e*<sub>*k*</sub>) and on the fre-
quency couple (*ε*<sub>1</sub>, *ε*<sub>2</sub>).

Relations 13 and 14 are obtained for each layer. A matrix
assembly between layers *k* and *k* + 1 is then performed to
remove the vectors (DEP<sub>*k*</sub><sup>-</sup>, SIG<sub>*k*</sub><sup>-</sup>, DEP<sub>*k*+1</sub><sup>+</sup>, and SIG<sub>*k*+1</sub><sup>+</sup>). This
assembly is obtained by taking into account the boundary condi-
tions between the layers (see appendix 2). The case of two
adhesive layers is detailed in Appendix 3. The assembly of
two layers is then considered as only one layer with boundary
conditions at the upper and lower faces respectively (DEP<sub>*k*</sub><sup>+</sup>,
SIG<sub>*k*</sub><sup>+</sup>, DEP<sub>*k*+1</sub><sup>-</sup>, SIG<sub>*k*+1</sub><sup>-</sup>). The method of layer assembly
is applied to the (*n* − 1) interfaces and yields a system of
equations (Eqs. (15) and (16)) which links the transforms of
the displacements and stresses of the surface (DEP<sub>*k*</sub><sup>+</sup>, SIG<sub>*k*</sub><sup>+</sup>) to
those of the lower face of layer *n* (DEP<sub>*n*</sub><sup>-</sup>, SIG<sub>*n*</sub><sup>-</sup>).

$$
\overline{DEP}_1^+ = [DS_a]\overline{SIG}_1^+ + [DD_a]\overline{DEP}_n^- \quad (15)
$$

$$
\overline{\text{SIG}}_{-}^{n} = [\overline{\text{SS}}_{a}] \overline{\text{SIG}}_{1}^{+} + [\overline{\text{SD}}_{a}] \overline{\text{DEP}}_{-}^{n} \quad (16)
$$

Relation (17) is finally obtained straightforwardly from rela-
tions 15 and 16 and the boundary conditions at the lower face
of layer *n*. It is written here in matrix form. It links surface
displacements and unit contact pressure and/or traction, taking
into account the boundary conditions at both the layer interface
and the lower face of layer *n*.

$$
\mathrm{DEP}_{1}^{+} = [\mathrm{Infl}] \mathrm{SIG}_{1}^{+} \quad (17)
$$

Displacements due to any distribution of normal and tangen-
tial tractions are obtained from (17) by superposition. This
relation will be used for the contact solution. The influence
coefficient matrix depends on frequencies, number of layers,
their thicknesses and mechanical properties, and interfacial
boundary conditions. It needs to be computed only once during
the contact simulation for each frequency couple. Symmetries
leading to reduction in matrix size are taken into account. An
inverse Fast Fourier Transform algorithm is used to obtain in-
fluence coefficients in the real domain. This algorithm avoids
singularity problems encountered during numerical integration
of the inverse transform and is inexpensive in computer time
when compared to the direct integration method in the inverse
transform. This is essentially due to the fact that the periodicity
of trigonometric functions is taken into account. The next step
is now concerned with the contact solution.

2.1.2 *Contact Solution.* The contact problem solution is split into two parts: a normal and a tangential problem solution. The direct method is used. The two bodies are pressed together by a force *P* over an area of initially unknown semi-axes *a* and *b*. Sliding and rolling conditions can be considered:

* two tangential loads along x₁ and x₂ axes are applied to the ellipsoid. Coulomb's law of friction is used.

* rolling contact conditions (spin, creep coefficient along x₁ and x₂ axes, etc...) hold (Kalker, 1990).

The classical contact solution is used: normal and tangential problems are uncoupled and solved in turn. The potential contact area is discretized into regular rectangular cells (*n*₁*n*₂) on which the pressure and tractions are assumed to be constant (Fig. 3). The solution of a discrete contact problem is therefore a set of stresses and displacements which satisfies the elasticity

equations of contacting bodies and boundary equations (18, 21, 22, 24, 26, 27, 28, 29) and inequalities (19, 20, 23, 25) in and outside the contact area for normal (Fig. 5) and tangential contact conditions. Thus, at the central point of an area (i, j):

Normal conditions:

$$
\begin{align}
u_3^{\text{ell}}(i, j, 0) + u_3^1(i, j, 0) &= \delta_3 - d(i, j, 0) \tag{18} \\
&\qquad \text{within the contact area}
\end{align}
$$

$$
\begin{align}
u_3^{\text{ell}}(i, j, 0) + u_3^1(i, j, 0) &> \delta_3 - d(i, j, 0) \tag{19} \\
&\quad \text{outside the contact area}
\end{align}
$$

$$
p(i, j, 0) > 0 \quad \text{within the contact area} \tag{20}
$$

$$
p(i, j, 0) = 0 \quad \text{outside the contact area} \tag{21}
$$

$$
\int_{\Gamma_c} \sigma_{33}(i, j, 0) dS = P \quad \text{traction equilibrium} \quad (22)
$$

Tangential conditions (generalities)

$$
|t(i, j, 0)| < \mu_f p(i, j, 0) \quad \text{stick zone} \quad (23)
$$

$$
|t(i, j, 0)| = \mu_f p(i, j, 0) \quad \text{slip zone} \quad (24)
$$

$$
\begin{equation}
\begin{split}
& (u_1^1(i, j, 0) - u_1^\text{ell}(i, j, 0))tx_1(i, j, 0) + \\
& \phantom{(u_1^1(i, j, 0) - u_1^\text{ell}(i, j, 0))tx_1(i, j, 0) +} \text{slip zone}
\end{split}
\tag{25}
\end{equation}
$$

$$
(u_2^1(i, j, 0) - u_2^\text{ell}(i, j, 0))tx_2(i, j, 0) < 0
$$

$$
\begin{equation}
\begin{split}
& (u_1^1(i, j, 0) - u_1^\text{ell}(i, j, 0))tx_2(i, j, 0) + \\
& \quad \text{slip zone} \\
& (u_2^1(i, j, 0) - u_2^\text{ell}(i, j, 0))tx_1(i, j, 0) = 0
\end{split}
\tag{26}
\end{equation}
$$

Partial slip (for *i* = 1, 2)

$$
(u_i^{\text{ell}}(i, j, 0) - u_i^1(i, j, 0)) = \delta_i \quad \text{stick zone} \quad (27)
$$

$$
\int_{\Gamma_c} \sigma_{i3}(i, j, 0) dS = Txi \quad \text{traction equilibrium} \quad (28)
$$

Rolling conditions (for *i* = 1, 2)

$$
\begin{equation}
\begin{split}
\mathrm{s}\dot{x}_i/V = \xi x_i - \varphi \frac{x_2}{c} &+ \left( \frac{\partial \mathrm{u}_j^{\mathrm{ell}}(i, j, 0)}{\partial x_1} - \frac{\partial \mathrm{u}_j^1(i, j, 0)}{\partial x_1} \right)
\end{split}
\tag{29}
\end{equation}
$$

Classical formulation (Eqs. (30)–(31)) of elastic deflection is
performed. The influence coefficient, *a*<sub>*pm*</sub>(*i*, *j*, *k*, *l*), is obtained
in [Infl] after this latter is transformed inversely.

$$
u_3^1(i, j, 0) = \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} a_{33}(i, j, k, l)p(i, j, 0) \quad (30)
$$

$$
\begin{equation}
\begin{split}
u_1^1(i, j, 0) = \sum_{p=1}^{2} & \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} a_{lp}(i, j, k, l)tx_p(i, j, 0) \\
& 1 = 1, 2 \quad (31)
\end{split}
\end{equation}
$$

Normal and tangential contact problems are solved using
classical unilateral analysis with friction, following Kalker
(1982), Mindlin (1949), and Gattina (1987). Normal and tan-
gential traction distributions, contact area, stick and slip zone

<figure><img src="image_6.png" /><figcaption>Fig. 5 Definition of parameters used in a normal contact problem solution</figcaption></figure>
repartition satisfying equations and inequalities are determined. The tangential contact problem formulation does not present convergence difficulties whatever the loading conditions, are due to the use of the Newton-Raphson method (Gattina, 1987). Note further that the maximum number of iterations is not excessive (≈10).

### 2.1.3 Subsurface Stress Field Determination.
Pressure and traction distributions are used as data for the calculation of the total stress field at node (i, j, l). Subsurface stresses and displacements produced by unit pressure and tractions at the contact surface are expressed in the Fourier domain. Therefore partial assembly of the layers is performed to obtain $SIG_k(i, j, l)$, depending on $SIG_k^+(k, m, 0)$ (relation 32), by using relations 13 and 14, the boundary conditions at interfaces and at the lower border of layer n (see Appendix 2).

$$ \overline{SIG}_j(i, j, l) = [\text{Mat}] \overline{SIG}_j^+(k, m, 0) \quad (32) $$

An inverse transform is then performed to obtain the stress field in the spatial domain due to unit pressure or traction. The total stress field at (i, j, l) is then built up by superposition of stresses induced by all pressures and tractions at the contact surface.

## 2.2 Tests.
Different configurations were investigated to test the validity and the accuracy of the model.

1. The matrix assembly was tested. A layer of thickness, $e$, was modeled, on the one hand, as a single layer and, on the other hand, as a multilayered body composed of 10 layers with identical mechanical properties. Perfectly bonded conditions were considered at each internal interface. No difference was obtained in the contact solution and the subsurface stress field (no discontinuity at internal interfaces) between both models.

2. Additional comparisons were made between those results and those:

* of Kalker (1990) and Johnson (1985) for a half space under rolling contact and partial slip conditions for testing the contact problem solution and the stress field determination.

* obtained by using the Finite Element Method (Nastran-Ideas) for a 3-D single layer stress field under one normal loading condition.

* of O'Sullivan and King (1988) concerning layered half space under sliding conditions.

Differences do not exceed 5 percent for all comparisons.

## 2.3 Computer Time.
This 3-D influence coefficient formulation based on Fourier integral transforms and FFT permits considering complex sliding and rolling contact conditions, various layer interfacial boundary conditions and multilayered body configurations. The contact solution and stress field determination are, moreover, performed accurately thanks to grid refinement and do not require considerable memory size and computer time.

For instance, a single layer subjected to one normal load and having large dimensions compared to those of the final contact region, was modeled according to FEM using ABAQUS software by Sassi et al. (1996). These authors solved the normal contact problem and determined stress field within the layer described by 8153 nodes: 144 contact nodes are localised within the contact area. Subsurface stresses are determined along 10 planes parallel to the contact surface. Their calculation takes 110 minutes on an HP 9000-720. This layered configuration is modeled here according to a grid specified by 128 × 128 × 17 nodes along $x_1, x_2$ and $x_3$ axes. This grid was also defined in order to obtain 144 contact nodes. The contact solution and internal stresses require 28 CPU minutes on an HP 9000-715. Our computer system is twice as fast; although we used a more refined mesh, this leads to a calculation twice as fast.

## 2.4 General Remarks About the Model.
A 3-D influence coefficient formulation has been presented in this paper. Some points can be pointed out:

1. A half numerical and half analytical method has been presented to determine influence coefficients in a multilayered half-space or body configuration. No limits on layer thickness are formulated with regard to accuracy: layers of thickness ranging from a micrometer to several centimeters can be considered, as well as successive thin-thick layer configuration.

2. The FFT technique is applicable strictly to periodic contacts. For this reason, the grid must extend sufficiently far beyond the contact area in order to eliminate periodic effects in the solution. For instance, the grid size along direction $x_i$ ($i = 1, 2$), must be five times greater than the contact half width along direction $x_i$ so that the contact solution is not disturbed. But contact periodicity has a more significant effect on tensile stresses than on contact stresses. The computation of the internal stresses consequently needs a grid size along direction $x_i$ about 30 times the contact area half width along direction $x_i$. As demonstrated in part 2.3, a large grid does not increase computer time compared to FEM.

3. The direct method requires $O(N_c^2)$ operations to solve the contact problem, where $N_c$ represents the number of contact points. The method developed by Nogi and Kato (1997) solves the problem of $O(N \log_2 N)$ complexity, where $N$ is the total number of points considered at the contact surface. The difference in computational time used by both methods depends on $N_c$ and $N$. An increase in $N_c$ may lead to longer computer time using the direct method than Nogi and Kato's formulation. Brigham (1974) observed the time required to compute 1-D summation by both the direct and FFT approaches, as a function of $N$. He showed that both methods are competitive if $N$ does not exceed 64. But as the aim of our model is to solve a smooth contact problem, a great number of contact points is not necessary to ensure good accuracy of 1 percent (≈200 points). No significant difference in computer time between Nogi and Kato's method and the one presented here is obtained for this range of contact node number.

4. The method presented here is adapted to further developments, such as for instance, taking into account hysteresis effects between two bodies.

# 3 Application to Specific Coatings

Coating thickness and mechanical properties have a great influence on surface contact conditions (in comparison to uncoated cases) and on internal stresses. Designing against the risk of debonding requires knowledge of these variables.

As a first study, the influence of one coating over an elastic substrate on the normal contact solution and stress field is analyzed. Results are then compared with those obtained using 2-D models (Gupta and Walowitt, 1974; Leroy, 1989; Mao et al., 1995; and Komvopoulos, 1988) and 3-D models (O'Sullivan and King, 1988 and Chiu and Hartnett, 1983).

The main difference between this study and previous ones is that it considers a rigid ellipsoid ($Rx_1 = 120$ mm and $Rx_2 = 200$ mm) pressed on an elastic coated medium. The coating is defined by its thickness ($e_1$) and material properties ($E_1$ and $\nu_1$). It is perfectly bonded to the substrate whose Young modulus ($E_2$), Poisson's coefficient ($\nu_2$), and thickness ($e_2$) are constant and equal, respectively, to 200 GPa, 0.3, and 50 mm. The contact surface is described by $128 \times 128$ cells; the lengths along directions $x_1$ and $x_2$ in which the Fourier transform is applied, depend on the coating/substrate configuration. One normal load, $P$ of 200 N is considered.
Results obtained for coating/substrate configurations are compared to a reference case, corresponding to the uncoated substrate. Results obtained in the reference case are indicated by the subscript 0.

**3.1 Influence of Coating Parameters on Normal Contact Solution.** Maximal pressure Pmax and contact area size variations (A) (see Fig. 5) are studied for different coating/substrate configurations. Contact area size (A) is determined following Eq. (33). This numerical expression was preferred to a mathematical expression because of an imperfect elliptical form of contact area for the layered configuration,

$$A = \sqrt{N_c \Delta x_1 \Delta x_2} \qquad (33)$$

where $N_c$: number of contact points, $\Delta x_i$: length of the cell along direction $x_i$.

A configuration is defined by the coating finite thickness, $e_1$ and the Young modulus ratio between the coating and the substrate: $E1/E2$. These variations are compared to those ($Pmax_0, a_0, b_0, A_0$) for the reference case. Variations of $Pmax/Pmax_0$ and $A/A_0$ are presented versus $e_1/(A_0)^{1/2}$ in Figs. 6 and 7, respectively.

As reported in these figures, variations of $Pmax/Pmax_0$ and $A/A_0$ depend on ratio $E1/E2$. When a stiff coating overlays the substrate, Pmax increases and A decreases with the increase of $e_1/(A_0)^{1/2}$, while reciprocally, an increase in $e_1/(A_0)^{1/2}$ leads to a decrease in Pmax and a decrease in A for a soft layer bonded to a substrate. O'Sullivan and King (1988) and Chiu and Hartnett (1983) described the similar variations of Pmax and the contact area for a spherical layered configuration.

Furthermore, for a 2-D coating/substrate configuration, Gupta and Walowitt (1974), Leroy (1989), Mao et al. (1995), and Komvopoulos (1988) described three kinds of behavior depending on the layer thickness and contact area size. Their results can be reformulated for a 3-D coating/substrate configuration under normal elliptical contact conditions, as follows:

1. $e_1/(A_0)^{1/2} < 0.05$, when the coating thickness is very thin, no significant influence of the coating on the contact solution is noted. Pmax and A are approximately equal to $Pmax_0$ and $A_0$.

2. $0.05 < e_1/(A_0)^{1/2} < 1.5$, when the coating becomes thicker, Pmax increases or decreases depending on ratio $E_1/E_2$. $E1/E2 > 1$ leads to a stiffer coating/substrate configuration than the reference case. As in a half-plane approximation, the configuration with the stiffest mechanical properties has the greatest Pmax values and the smallest contact area. This leads to $Pmax/Pmax_0 > 1$ and $A/A_0 < 1$. Reciprocally, $E1/E2 < 1$ defines a softer coating/substrate configuration than the reference case. $Pmax/Pmax_0$ and $A/A_0$ are respectively smaller and greater than 1.

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Maximal normal pressure variations depending on the finite thickness of the coating</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Contact area size variations depending on the finite thickness of the coating</figcaption></figure>

3. $e_1/((A_0)^{1/2}) > 1.5$, the contact conditions are governed by the coating behavior independently of the substrate. $Pmax/Pmax_0$ and $A/A_0$ tend to values obtained without the substrate.

**3.2 Stress Field Analysis.** Attention is focused on stresses at the coating/substrate interface and on Von Mises stresses in coating and substrate. Attempts at relating stresses to phenomena such as risk of debonding and delamination are presented in first subpart.

**3.2.1 Study of the Interface Between Coating and Substrate.** Stress field analysis at the interface gives information about potential phenomena such as the risk of debonding and delamination. O'Sullivan and King (1988) considered a spherical indenter over an elastic substrate under complete sliding conditions. The interface behavior is analyzed in terms of shear stresses. These authors showed that maintaining low interfacial shear stress is important, and for a layered medium the use of a compliant layer can be beneficial compared to a stiff layer. Komvopoulos (1988) and Leroy (1989) showed that the debonding risk in 2-D layer/substrate configuration, can moreover be linked to the progression of tensile stresses at the interface.

Here, a normally loaded elliptical indenter is considered. Attention is focused on the progression of tensile stresses at the interface, as low values of shear stresses are obtained at the interface under normal loading. The coating is considered perfectly bonded to the substrate. Normal and tangential stresses ($\sigma_{x1x3}$, $\sigma_{x2x3}$, and $\sigma_{x3x3}$) are continuous at the interface. As no equality exists concerning other stresses ($\sigma_{x1x1}$, $\sigma_{x2x2}$, and $\sigma_{x1x2}$), discontinuities appear. Internal tensile stresses, $\sigma_{x1x1}$ and $\sigma_{x2x2}$, behave differently versus directions $x_1$ and $x_2$ respectively due to elliptical contact conditions. As $\sigma_{x2x2}$ differences at the interface layer between coating and substrate are small compared to $\sigma_{x1x1}$ ones, they will not be studied here. $\sigma_{x1x1}$ variations at the coating ($\sigma_{x1x1}^{-1}$) and substrate ($\sigma_{x1x1}^{2+}$) interface for different coating thicknesses ($e_1 = 5, 250$ and $800$ µm) and material properties ($E1 = 100, 400$, and $600$ GPa), are presented in a meridian plane of contact where the severest state of stress occurs (Fig. 8). Figures 9 and 10 represent, respectively, the influence of interface depth and material property on interface behavior.

The coating is subjected to different effects:

* a compression due to normal contact; this compression decreases within the coating thickness (Fig. 9)

* a lateral effect due to the difference of coating and substrate mechanical properties, generating compressive and tractive stresses (Figs. 9 and 10). Further tractive stresses are never obtained for $E1/E2 < 1$.

A coating interface situated close to the contact surface is mainly subjected to the first effect. As the coating thickness
<figure><img src="image_9.png" /><figcaption>Fig. 8 Definition of the meridian plane where the stress field is analyzed</figcaption></figure>

increases, (Figs. 9(b) and 9(c)), the latter effect becomes predominant; tractive $\sigma_{x1x1}$ stresses appear. Furthermore, interface depth is not the only governing parameter of interface behavior. In Fig. 10, the coating thickness is constant and is equal to 250 µm. The interface behavior also depends on the material property of the coating. When $E1/E2 < 1$ (Fig. 10(a)), no significant discontinuity in $\sigma_{x1x1}$ stresses appears at the interface. But for an $E1/E2$ greater than 1 (Figs. 10(b)–(c)), $\sigma_{x1x1}^{1-}$ variations versus $x_1$ are modified, tend to 0 (Fig. 10(b)) or become tractive (Fig. 10(c)). This phenomenon induced by $E1/E2$ was mentioned in 2-D analyses performed by Komvopoulos (1988) and Leroy (1989). According to the 3-D model, the interfacial condition of displacement continuity versus $x_1$ leads to discontinuities in $\sigma_{x1x1}$. For $e_1 = 250$ µm, substrate mainly governs the multilayered body behavior. Therefore the strain field at the interface ($\epsilon_{x1x1}$) is approximately the same as those obtained at a depth of 250 µm in a half space having the same substrate material properties.

Further, $\sigma_{x1x1}$ stresses at the interface may be expressed as a first approximation by $\sigma_{x1x1}^j = E_j\epsilon_{x1x1}$. It gives the following relation between $\sigma_{x1x1}$ interface stresses in the two layers and Young's modulus of each layer:

$$ \frac{\sigma_{x1x1}^{1-}}{\sigma_{x1x1}^{2+}} = \frac{E_1}{E_2} \quad (34) $$

When $E1/E2$ is greater than 1, the ratio between interfacial $\sigma_{x1x1}$ stresses becomes greater than 1 and leads to $\sigma_{x1x1}^{1-} > \sigma_{x1x1}^{2+}$. Reciproically for $E1/E2 < 1$, $\sigma_{x1x1}^{1-}$ is smaller than $\sigma_{x1x1}^{2+}$ as observed in Fig. 10. Furthermore, the amplitude of discontinuities increases with the increase of $E1/E2$ (comparison of Figs. 10(b) and 10(c)).

In Fig. 10(c), tractive $\sigma_{x1x1}$ values are obtained below the contact surface at the interface. The interface will also be subjected to tractive $\sigma_{x1x1}$ stresses in the coating and reciprocally, to compressive stresses in the substrate. No tractive

<figure><img src="image_10.png" /><figcaption>Fig. 10 Interface study, influence of material properties E1 and E2</figcaption></figure>

$\sigma_{x2x2}$ stresses in the coating are obtained at the interface due to elliptical contact. This phenomenon induces flexion mechanisms at interface only along direction $x_1$ (Fig. 11). If the difference between the energy induced by stresses at the interface is greater than the interface energy limit, a debonding risk may appear and give rise to cracking. Designing against this debonding risk implies a change of interface depth through ratio $e_1/\sqrt{A}$. A decrease of this ratio will lead to a reduction in tractive $\sigma_{x1x1}$ and, furthermore, to negative $\sigma_{x1x1}$ stresses in the coating due to an interface depth closer to the contact surface. In addition, a reduction in the ratio $e_1/\sqrt{A}$ may be obtained by a coating thickness decrease or a change in the normal loading condition. For example, an increase in normal load up to 500 N will increase the contact area and also reduce ratio $e_1/\sqrt{A}$.

As already shown by Komvopoulos (1988) and Leroy (1989) for 2-D configurations, tensile stresses at the layer interface leading to delamination depend on both the layer configuration (stiff-soft) and the interface depth location. For 3-D configurations, this remark still holds and, furthermore, can be encountered in two directions, $x_1$ and $x_2$ and in only one direction for, respectively, a spherical and an elliptical contact problem.

### 3.2.2 Von Mises Stress Field Analysis.
Coatings are generally used to unload the substrate. Previous studies (Chiu and Hartnett, 1983; O'Sullivan and King, 1988; Komvopoulos, 1988) showed that a coating/substrate configuration does not necessarily lead to better behavior than an uncoated substrate configuration. This application is presented to confirm that similar conclusions can be drawn for elliptical contact geometry.

$E1/E2$ ranges from 0.5 to 3. Three coating thicknesses are considered. The results are reported in Table 1.

In the reference case, $(\sigma_{\text{umis}})_{\text{MAX}}$ is equal to 368 MPa. As reported in Table 1, $(\sigma_{\text{umis}})_{\text{MAX}}$ in the substrate is not always smaller than $((\sigma_{\text{umis}})_{\text{MAX}})_0$. A compliant coating of intermediate thickness unloads the substrate, but a stiff coating ($E1/E2 > 1$) will overload the substrate for $e_1 < 800$ µm when $E1/E2 = 2$. The increase in $(\sigma_{\text{umis}})_{\text{MAX}}$ in coating with $e_1$, may be ex-

<figure><img src="image_11.png" /><figcaption>Fig. 9 Interface study, influence of interface depth ($E1/E2 = 2$)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 Debonding phenomenon</figcaption></figure>
<table><caption>Table 1 Von Mises stress field analysis</caption><thead><tr><th rowspan="2">Variable</th><th colspan="3">e<sub>1</sub> = 250 µm, variation of E1/E2</th><th colspan="2">E1/E2 = 2, variation of e<sub>1</sub></th></tr><tr><th>E1/E2 = 0.5</th><th>E1/E2 = 2</th><th>E1/E2 = 3</th><th>e<sub>1</sub> = 5 µm</th><th>e<sub>1</sub> = 800 µm</th></tr></thead><tbody><tr><td>(σ<sub>vmis</sub>)<sub>MAX</sub> in coating (MPa)</td><td>250</td><td>380</td><td>454</td><td>371</td><td>518</td></tr><tr><td>depth (mm)</td><td>0.15</td><td>0.15</td><td>0.15</td><td>0.005</td><td>0.13</td></tr><tr><td>(σ<sub>vmis</sub>)<sub>MAX</sub> in substrate (MPa)</td><td>237</td><td>371</td><td>430</td><td>236</td><td>126</td></tr><tr><td>depth (mm)</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.2</td><td>0.8</td></tr></tbody></table>

plained by considering the existence of tractive $\sigma_{x1x1}$ stresses at the interface. Indeed, Von Mises stress is related to other stresses by the relation:

$$ \sigma_{\text{vmis}} = \sqrt{\frac{1}{6}\left[\sum_{i \neq j=1}^{3} (\sigma_{xixi} - \sigma_{xjxj})^2 + \sum_{i \neq j=1}^{3} (\sigma_{ixj})^2\right]} \quad (35) $$

$\sigma_{x1x1}$, $\sigma_{x2x2}$, and $\sigma_{x3x3}$ stresses are mostly compressive within the multilayered body, only tractive $\sigma_{x1x1}$ stresses appear near the interface, thus leading to greater $(\sigma_{x1x1} - \sigma_{x3x3})$ and $(\sigma_{x1x1} - \sigma_{x2x2})$ terms and to an increase in the $(\sigma_{\text{vmis}})$ value at the interface. A stiff coating may also lead to a greater value in $(\sigma_{\text{vmis}})_{\text{MAX}}$ and to debonding risk at the interface depending on its own thickness. When its thickness increases up to 800 µm, the coating mainly governs the behavior of the multilayered configuration, and satisfies the objective of the coating: to unload the substrate. As the highly stressed zones are now situated in the coating, the next aim is thus to estimate the severity of these stresses through a comparison with the coating stress limits.

## 4 Conclusion

A numerical 3-D contact model between the elastic multilayered body and rigid bodies was presented to guide choice among coating/substrate combinations that can withstand the applied loads. This model is able to solve the contact problem under partial slip, rolling/sliding contact conditions and to determine the subsurface stress field. It is based on linear elasticity theory, unilateral contact with friction and Fourier transforms. The application of Fourier techniques avoids singularity problems during inverse integration and reduces computer time. This model further has the advantage of solving the problem thanks to grid refinement and does not require considerable memory size.

As a first step, one application to a specific coating/substrate configuration was analyzed in this paper. Analysis of normal elliptical contact pressure distribution and subsurface stress field was performed. It led to similar conclusions to those obtained previously for 2-D and 3-D coating/substrate configurations: an increase or decrease in maximal pressure and contact area depending on the mechanical properties of the coating, the existence of discontinuity in stresses at the interface between the coating and the substrate, tractive stresses which may lead to debonding risk and an increase of Von Mises stresses in the substrate for a stiff coating of medium thickness. The 3-D stress field approach is interesting in debonding risk analysis, as this risk may exist along one or two directions (for spherical loading conditions) ($x_1$ and $x_2$). A 2-D coating/substrate analysis may hide one of these directions, such as a 3-D analysis restricted to spherical contact geometry.

Further analyses are now required to understand the behavior of a multilayered body in order to eliminate incompatible layer combinations and to minimize fatigue phenomena during the life of the contact under given loading conditions (tangential loading conditions).

Studies are underway to use this model in biomechanics (Plumet and Dubourg, 1997). Indeed, a gamma irradiation effect leads to different mechanical properties versus depth in sterilized tibial inserts and may be responsible for delamination in

prostheses. It requires accounting for several layers within finite thickness of tibial insert.

## Acknowledgments

The authors would like to acknowledge the support of the BRITE-EURAM project BE-7928 which sponsored this work.

## References

Bentall, R. H., and Johnson, K. L., 1968, "An Elastic Strip in Plane Rolling Contact," *International Journal of Mechanical Sciences*, Vol. 10, Mar., pp. 637-663.

Brigham, E. O., 1974, *The Fast Fourier Transform*, Prentice-Hall NJ, NJ, 252 pp.

Cattaneo, C., 1938, "Sul contatto di due corpi elastici," *Acc; dei Lincei, Rend. Ser.*, 6, 27, pp. 342-348.

Chang, T. P., Cheng, H. S., Chou, W. A., and Sproul, W. D., 1990, "A Comparison of Fatigue Failure Morphology Between TiN Coated and Uncoated Lubricated Rollers," *ASME JOURNAL OF TRIBOLOGY*, Vol. 112.

Chiu, Y. P., and Hartnett, M. J., 1983, "A Numerical Solution for Layered Solid Contact Problems With Application to Bearings," *ASME JOURNAL OF LUBRICATION TECHNOLOGY*, Vol. 105, pp. 585-590.

Gattina, J. C., 1987, "Contacts de corps élastiques, effets tangentiels et normaux, formulation et résolution des problèmes inverses et directs," thèse de Doctorat d'état: Insa de Lyon, UCB Lyon I, 180 pp.

Goodman, L. E., and Keer, L. M., 1975, "Influence of an Elastic Layer on the Tangential Compliance of Bodies in Contact," *The Mechanics of the Contact Between Deformable Bodies*, De Pater, A. D., and Kalker, J. J., eds., Delft University Press, pp. 127-151.

Gupta, P. K., Walowitt, J. A., Finkin, E. F., October 1973, "Stress Distributions in Plane Strain Layered Elastic Solids Subjected to Arbitrary Boundary Loading," *ASME JOURNAL OF LUBRICATION TECHNOLOGY*, Vol. 95, pp. 427-433.

Gupta, P. K., Walowitt, J. A., 1974, "Contact Stress Between an Elastic Cylinder and a Layered Elastic Solid," *ASME JOURNAL OF LUBRICATION TECHNOLOGY*, Vol. 96, pp. 250-257.

Hannah, M., July 1950, "Contact Stress and Deformation in a Thin Elastic Layer," *Journal of Mech. and Applied Mathematics*, Vol. 4, pp. 94-105.

Hertz, H., 1882, "Über die beruhrung fester elastischer korper," *J. Reine und angewandte mathematik*, Vol. 92, pp. 156-171.

Johnson, K. L., 1985, *Contact Mechanics*, Edited by Cambridge University Press, 451 pp.

Ju, Y., and Farris, T. N., 1997, "FFT Thermoelastic Solutions for Moving Heat Sources," *ASME JOURNAL OF TRIBOLOGY*, Vol. 119, pp. 156-162.

Kalker, J. J., 1972, "On Elastic Line," *ASME Journal of Applied Mechanics*, pp. 1125-1132.

Kalker, J. J., 1982, "The Contact Between Wheel and Rail," *Report of the Department of Mathematics and Informatics, No. 82-87, Delft: Delft University of Technology*, 36 pp.

Kalker, J. J., 1990, "Three Dimensional Elastic Bodies in Rolling Contact," Kluwer Academic Publishers, eds., 314 pp.

Komvopoulos, K., 1988, "Finite Element Analysis of a Layered Elastic Solid in Normal Contact With a Rigid Surface," *ASME JOURNAL OF TRIBOLOGY*, Vol. 110, pp. 477-485.

Kuo, C. H., Keer, L. M., 1992, "Contact Stress Analysis of a Layered Transversely Isotropic Halspace," *ASME JOURNAL OF TRIBOLOGY*, Vol. 114, pp. 253-262.

Leroy, J. M., Floquet, A., and Villechaise, B., 1989, "Thermomechanical Behavior of Multilayered Media: Theory," *ASME JOURNAL OF TRIBOLOGY*, Vol. 111, pp. 538-544.

Leroy, J. M., 1989, "Modélisation thermoelastique des revetements de surface utilises dans les contacts non lubrifies," thèse de Doctorat, INSA de Lyon, 209 pp.

Mao, K., Sun, Y., and Bell, T., 1995, "A Numerical Model for the Dry Sliding Contact of Layered Elastic Bodies with Rough Surfaces," STLE Preprint No. 95-TC-6A-1, STLE/ASME Tribology Conference in Kissimmee, FL, Oct. 11, 1995, 9 pp.

Mindlin, R. D., 1949, "Compliance of Elastic Bodies in Contact," *ASME Journal of Applied Mechanics*, Vol. 16, pp. 259-268.

Nogi, T., and Kato, T., 1997, "Influence of a Hard Surface Layer on the Limit of Elastic Contact—Part I: Analysis Using a Real Surface Model," *ASME Journal of Tribology*, Vol. 119, pp. 493-500.
O'Sullivan, T. C., and King, R. B., 1988, "Sliding Contact Stress Field Due to a Spherical Indenter on a Layered Elastic Half-Space," ASME JOURNAL OF TRIBOLOGY, Vol. 110, pp. 235-240.

Plumet, S., and Dubourg, M.-C., 1997, "Analysis of Gamma Irradiation Effects on Tibial Insert: Stress Field Analysis Within Sterilised Tibial Insert Model Depending on the Multicoated Body Configuration," Technical report 2.1 number 4 of BRITE/EURAM Contract, BRE2-CT-94-7928.

Rabinovitch, V. L., Sipcic, S. R., and Sarim, V. K., 1994, "Three-Dimensional Unilateral Frictionless Contact Problem for Finite Bodies," ASME Journal of Applied Mechanics, Vol. 61, pp. 54-59.

Timoshenko, S. P., and Goodhier, J. N., 1970, Theory of Elasticity, Third Edition, McGraw-Hill, New York, pp. 409-420.

Sassi, M., Desvignes, and M., Foulon, M., 1996, "Etude tridimensionnelle des contact élastiques sans frottement par la méthode des éléments finis," Mécanique Industrielle et Matériaux, Vol. 49, No. 1, pp. 23-27.

Sneddon, I. N., 1951, "Computation of Stress and Displacements in a Layered Elastic Medium," International Journal of Engineering Science, Vol. 9, pp. 775-800.

# APPENDIX 1

## Fourier Transform

The Fourier transform used in this model is the discrete Fourier transform. It is based on an initial definition of the finite length discrete transform which approximates the continuous Fourier transform (Relation 5). The inverse transform is:

$$h(x_1, x_2, x_3) = \iint_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \overline{h(f_r, g_s, x_3)} \exp(i2\pi(f_r x_1 + g_s x_2)) df dg_s$$

By considering this definition, a derivative property is obtained:

$$\frac{d^n \overline{h(f_r, g_s, x_3)}}{dx_1^n} = (2i\pi f_r)^n \overline{h(f_r, g_s, x_3)}$$

where

$f_r$ represents Fourier frequency along direction $x_1$
$g_s$, Fourier frequency along direction $x_2$

# APPENDIX 2

## Boundary Conditions

Different interfacial boundary conditions can be applied:

1. perfectly bonded:

$$\overline{\text{DEP}_k} = \overline{\text{DEP}_{k+1}} \\ \overline{\text{SIG}_k} = \overline{\text{SIG}_{k+1}}$$

2. imperfectly bonded

$$\overline{\sigma_{33}^k(e_k)} = \overline{\sigma_{33}^{k+1}(e_k)} \\ \overline{u_3^k(e_k)} = \overline{u_3^{k+1}(e_k)} \\ f \cdot \overline{\sigma_{33}^k(e_k)} = \overline{\sigma_{13}^{k+1}(e_k)}$$

Boundary conditions at the lower face of layer n

1. No stresses

$$\overline{\text{SIG}_n} = 0$$

2. No displacements

$$\overline{\text{DEP}_n} = 0$$

# APPENDIX 3

## Layer Assembly in the Case of Two Adhesive Layers

The transfer-matrices and vectors for each layer are:

layer k

$$\overline{\text{DEP}_k^+}, \overline{\text{SIG}_k^+}, \overline{\text{DEP}_k^-} = [\text{DS}_k]\overline{\text{SIG}_k^+} + [\text{DD}_k]\overline{\text{DEP}_k^-} \\ \overline{\text{DEP}_k^-}, \overline{\text{SIG}_k^-}, \overline{\text{SIG}_k^+} = [\text{SS}_k]\overline{\text{SIG}_k^+} + [\text{SD}_k]\overline{\text{DEP}_k^-}$$

layer k + 1

$$\overline{\text{DEP}_{k+1}^-}, \overline{\text{SIG}_{k+1}^+}, \overline{\text{DEP}_{k+1}^+} = [\text{DS}_{k+1}]\overline{\text{SIG}_{k+1}^+} + [\text{DD}_{k+1}]\overline{\text{DEP}_{k+1}^-} \\ \overline{\text{DEP}_{k+1}^-}, \overline{\text{SIG}_{k+1}^-}, \overline{\text{SIG}_{k+1}^+} = [\text{SS}_{k+1}]\overline{\text{SIG}_{k+1}^+} + [\text{SD}_{k+1}]\overline{\text{DEP}_{k+1}^-}$$

The displacements and stresses are continuous at the interface, so:

$$\overline{\text{DEP}_k} = \overline{\text{DEP}_{k+1}^+} \\ \overline{\text{SIG}_k} = \overline{\text{SIG}_{k+1}^-}$$

The transfer-matrices and vectors of the assembly are noted $[DD^*]$, $[SS^*]$, $[SD^*]$, and $[DS^*]$ with:

$$\overline{\text{DEP}_k^+} = [\text{DS}^*]\overline{\text{SIG}_k^+} + [\text{DD}^*]\overline{\text{DEP}_{k+1}^-} \\ \overline{\text{SIG}_{k+1}^-} = [\text{SS}^*]\overline{\text{SIG}_k^+} + [\text{SD}^*]\overline{\text{DEP}_{k+1}^-}$$

They are calculated as follows:

$$[A]^{-1} = [I] - [\text{SD}_k][\text{DS}_{k+1}] \quad \text{with} \quad [I] = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

$$[DS^*] = [DS_k] + [DD_k][DS_{k+1}][A][SS_k] \\ [DD^*] = [DD_k][A][DD_{k+1}] \\ [SS^*] = [SS_{k+1}][A][SS_k] \\ [SD^*] = [SS_{k+1}][A][SD_k][DD_{k+1}] + [SD_{k+1}]$$