# Contact analyses for anisotropic half space: effect of the anisotropy on the pressure distribution and contact area

Caroline Bagault, Daniel Nelias, Marie-Christine Baietto

► To cite this version:

Caroline Bagault, Daniel Nelias, Marie-Christine Baietto. Contact analyses for anisotropic half space: effect of the anisotropy on the pressure distribution and contact area. Journal of Tribology, 2012, 134, pp.031401-1-8. 10.1115/1.4006747. hal-00938302

HAL Id: hal-00938302

https://hal.science/hal-00938302

Submitted on 31 Aug 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
# Contact Analyses for Anisotropic Half Space: Effect of the Anisotropy on the Pressure Distribution and Contact Area

Caroline Bagault, Daniel Nélias, Marie-Christine Baietto

Université de Lyon, CNRS, INSA-Lyon, LaMCoS UMR5259, Villeurbanne, F69621 France

A contact model using semi-analytical methods, relying on elementary analytical solutions, has been developed. It is based on numerical techniques adapted to contact mechanics, with strong potential for inelastic, inhomogeneous or anisotropic materials. Recent developments aim to quantify displacements and stresses of an anisotropic mate-rial contacting both an isotropic or anisotropic material. The influence of symmetry axes on the contact solution will be more specifically analyzed.

## 1 Introduction

Engineering problems are becoming more complicated when trying to reduce the gap between the model and the actual application. It means that less restrictive assumptions should be made, or in other words more physics should be implemented in the model. Among the challenges to succeed in it, the material properties should be considered accurately. Supposing the material is isotropic is not enough. For most composite and mono-crystal materials their compositions or the elaboration and manufacturing processes imply that it exists one or two main directions or even a general anisotropy. The influence of the anisotropy orientations has to be taken into account to better predict the stress state in order to optimize the service life of industrial components.

The Green's functions of a point force applied in an isotropic infinite space was first solved by Kelvin [1]. Boussinesq [2] derived surface Green's functions for a force normal to the free surface, in isotropic solids. Then, Mindlin [3], by adding a complementary part of the solution, reduced the Kelvin's infinite space functions to the half space Green's functions. Ting and Lee [4] developed an explicit expression of the Green's functions in terms of the Stroh eigenvalues [5,6], for an anisotropic infinite space. Also using the Stroh formalism, Ting [7] obtained the anisotropic Green's functions of a point force in a half space in the Fourier transformed domain. With inverse fast Fourier transforms, solutions can be finally obtained in the physical domain. Pan and Yuan studied the three-dimensional Green's functions of point forces in anisotropic bimaterials [8], in anisotropic trimaterials [9] and also at the interface of an anisotropic bimaterial [10]. Ciavar-ella et al. [11] presented a method for solving the 3D contact problem for anisotropic materials by using the standard Hertzian solution. Li and Wang [12] worked also with the Hertzian solution in order to analyze the contact problem for two anisotropic piezoelectric bodies pressed together. Borodich [13] solved some contact problems of anisotropic elastodynamics by applying his method of integral characteristics of solutions to boundary-initial value problems. Swanson [14] used the procedure outlined by Willis [15], a numerical contour integration to determine the contact area and the pressure distribution, combined with the Pagano solution [16] to obtain detailed stress fields. Aizikovich et al. [17] worked on the influence of the elastic properties variation along the depth on the contact area and stresses. Gao and Pharr [18] were interested on the effective moduli of elastically anisotropic

solids under normal and tangential contacts. In their books, Rand and Rovenski [19] focused on the procedure of obtaining analytical solutions in anisotropic elasticity and Galin [20] dealt with contact problems and got onto contact problems for an anisotropic half-plane. Lin and Ovaert [21] studied the rough surface contact for anisotropic materials in 2D. Then He and Ovaert [22] developed a 3D contact model between a rough rigid sphere and a semi-infinite anisotropic elastic body, by applying the line integral of Barnett-Lothe tensors [23] on oblique planes. More recently, the anisotropic elastic layer was the subject of many articles. Kahya et al. [24] worked on the plane receding contact problem, Batra and Jiang [25] studied an analytical solution, and Argatov [26] focused on the depth-sensing indentation. Brock and Georgiadis [27] presented a class of multiple-zone sliding contact problems, including frictional and thermal effects on an anisotropic half space. Clements and Ang [28] solved some contact problems for inhomogeneous anisotropic elastic materials. Another method, the boundary element method, is used by Blázquez et al. [29] for generalized plane problems and by Rodríguez-Tembleque et al. [30] in 3D to study contact problems in anisotropic solids.

In this paper, starting from three dimensional Green's functions for anisotropic bimaterials [8], these functions are derived for anisotropic half space, as a sum of infinite space Green's functions and a complementary part, similar to the Mindlin's superposition method. The explicit expression of the Green's functions for anisotropic infinite space by Ting and Lee [4] are recalled and used as starting point to express the half space anisotropic contribution. The complementary part that should be added to account for the free surface of the half-space is finally given.

Several methods can be used for contact simulation of anisotropic materials. Semi-analytical methods have proven their efficiency in contact mechanics and are developed here to account for anisotropy of materials. The main advantage is here the small computing time compared to the finite element (FE) method that is however widely used for many contact problems. The contact model with its specificities when at least one of the bodies in contact behaves anisotropically will be first presented. Then, the influence of some parameters on the contact pressure distribution will be studied. The effect of the anisotropy orientation on the contact solution will be also investigated.

## 2 Contact Resolution

The semi-analytical method (SAM) consists of the summation of elementary solutions known analytically. One of the difficulties is the derivation or the identification of these elementary
analytical solutions, such as the well-known Boussinesq and Cerruti solutions in isotropic elasticity. The framework of the three-dimensional problem is simplified here by assuming the contact between one anisotropic elastic half space and a rigid body. The contact area is small in comparison to the dimensions of bodies justifying the assumption of half spaces. Each point of the surface is assigned a value of the pressure corresponding to the total load divided by the surface area. Analytical solutions giving the contributions of normal and tangential loading assumed uniform over a single rectangular element will be used. By summation the elastic deflection at each point within and near the contact area will be derived.

The elastic displacements are expressed by a double discrete convolution product between influence coefficients and the pressure or shear at the contact surface. The normal problem and the tangential problem in partial slip are therefore solved. The solution is performed by minimizing the complementary energy, so the contact pressure is constrained to be positive everywhere and there is no interpenetration. An algorithm is developed based on the conjugated gradient method (CGM). To accelerate the calculation, the fast Fourier transforms (FFT) are used to perform the double convolution product between the pressure and the influence coefficient matrix, at each iteration of the CGM. Once the contact problem solved, the strains in the half space are calculated.

The semi-analytical contact solver is based on the pioneering work of Jacq et al. [31] for elastic-plastic contacts. The solver has since been developed and improved in several ways. Boucly et al. [32] and Chen et al. [33] introduced thermal aspects. Fulleringer and Nelias [34] focused on the influence of a cuboid of uniform plastic strain in a half space, on the normal and tangential displacements of surface points, and derived the corresponding Green's functions in an analytical form. Then Leroux et al. [35] and Zhou et al. [36] studied the effects of the presence of inhomogeneous inclusions within a half-space, both on the contact pressure distribution and the subsurface stress field. Discrete convolution (DC) and 3D FFT have been since then widely used in the contact solvers. Leroux and Nelias [37] worked also on the stick-slip problem for a sphere in contact with a flat half-space containing unidirectional cylindrical fibers. Gallego et al. [38] proposed an algorithm based on the conjugate gradient method to account for the coupling between normal and tangential effects, which is required for frictional contact problems between elastically dissimilar materials. The same group of researchers also improved the contact algorithm in several ways so that it becomes affordable to simulate wear for 3D contact problems, cycle after cycle [39–41]. Note that SAM could be also linked with the FE method and used as a zoom on the contact to account for the effect of roughness or simulate wear [42]. More recently Chaise and Nelias [43] improved the numerical model to account for kinematic and isotropic hardening and analyzed the problem of a rolling load versus indentation. They also proposed a method to predict the coefficient of restitution when an elastic or rigid sphere is impacting an elastic-plastic half-space [44].

### 3 Anisotropic Influence Coefficients

Anisotropic materials are defined by the elastic stiffness tensor $C_{ijkl}$, which satisfies the full symmetry $C_{ijkl} = C_{jkl} = C_{klij}$. These materials can be divided into several parts, cubic, orthotropic, and anisotropic. Depending on the family classes of materials, three, nine or twenty one parameters are necessary for defining completely the elastic tensor. Note that in this work the main directions or axes of symmetry are not necessarily the same than those of the contact ones. In such a situation a base change has to be carried out. The elastic stiffness tensor links stresses to strains with this following relation:

$$ \sigma_{ij} = C_{ijkl} \varepsilon_{kl}. \qquad (1) $$

The influence coefficients, which link the loading to the displacements or to the strains, for an anisotropic material, are obtained with the Green's functions. The calculation of the displacements $u$ and the stresses $\sigma$ anywhere in the material is made in two stages. An infinite part [4] determines these elements for the anisotropic infinite space and a complementary part [8] is used to limit them to the half space,

$$ \begin{cases} u = u_{\infty} + u_{\text{comp}} \\ \sigma = \sigma_{\infty} + \sigma_{\text{comp}} \end{cases} \qquad (2) $$

These complementary functions can be described explicitly in the Fourier domain, whereas it is more complicated to obtain their formulation in the physical domain due to the general anisotropy of the material. They are therefore obtained in the physical domain by applying the inverse Fourier transform to the explicit formulas of the Fourier domain. In this paper, intermediate results are not all presented for conciseness.

A concentrated force $\vec{f} = (f_1, f_2, f_3)$ is applied at the origin of the reference frame, which is located on the surface of the anisotropic elastic half space ($x_3 \ge 0$). In the absence of body forces, the equations of equilibrium in terms of displacements $u_k$ are written as $C_{ijkl}u_{k,lj} = 0$. Three matrices $3 \times 3$ are defined with the tensor $C_{ijkl}$ and the vectors $n$ and $m$, which form a right handed triad with the position vector $x$. $Q, R$, and $T$ are the double projections of the elastic stiffness tensor,

$$ \vec{x} = \begin{pmatrix} r \cos \theta \sin \phi \\ r \sin \theta \sin \phi \\ r \cos \phi \end{pmatrix} \qquad (3) $$

$$ Q_{ik} = C_{ijkl} n_j n_s, \quad R_{ik} = C_{ijkl} n_j m_s, \quad T_{ik} = C_{ijkl} m_j m_s \qquad (4) $$

Six distinct eigenvalues $p$ are obtained by calculating the roots of

$$ \det(\Gamma(p)) = \det(Q + p(R + R^T) + p^2 T) = 0 \qquad (5) $$

The roots are three pairs of conjugate complex numbers

$$ p_j = \alpha_j + i\beta_j, \quad \beta_j > 0 \quad (j = 1, 2, 3) \qquad (6) $$

with $\alpha_j$ and $\beta_j$ real coefficients.

The complex eigenvectors $a$, Eq. (7), are not trivial solutions. The eigenvectors $b$ are derived by Eq. (8)

$$ [Q + p(R + R^T) + p^2 T]a = 0 \qquad (7) $$

$$ b = (R^T + pT)a = -\frac{1}{p}(Q + pR)a \qquad (8) $$

If $p_j$, $a_j$, and $b_j$ ($j = 1, 2, ..., 6$) are the eigenvalues and the associated eigenvectors, it is defined that

$$ \begin{align} & \operatorname{Im}(p_j) > 0 \quad (j = 1, 2, 3) \notag \\ & p_{j+3} = \bar{p}_j, \quad a_{j+3} = \bar{a}_j, \quad b_{j+3} = \bar{b}_j \tag{9} \\ & A = [a_1, a_2, a_3], \quad B = [b_1, b_2, b_3] \notag \end{align} $$

By supposing that the eigenvalues are distinct and the eigenvectors satisfy the normalization relation

$$ b_i^T a_j + a_i^T b_j = \delta_{ij} \qquad (10) $$

with $\delta_{ij}$ the Kronecker $\delta$, then the Green's functions are obtained by superposing the six eigensolutions. The calculated displacements and stresses are real. The determination of the roots' sextic equation from which the complex eigenvectors are derived is made numerically.
$\hat{\Gamma}$ is the adjoint of $\Gamma$ such that

$$ \Gamma(p)\hat{\Gamma}(p) = \det(\Gamma(p))I \quad (11) $$

The adjoint matrix $\hat{\Gamma}(p)$ is a polynomial of degree four in p

$$ \hat{\Gamma}(p) = \sum_{n=0}^{4} p^n \hat{\Gamma}^{(n)} \quad (12) $$

with $\hat{\Gamma}^{(n)}$ the real matrices which are independent of p. The infinite part of the displacement can be written as [4]

$$ \vec{u}_{\infty} = \frac{1}{4\pi r} \frac{1}{\det(T)} \sum_{n=0}^{4} q_n \vec{f} \hat{\Gamma}^{(n)} \quad (13) $$

The real coefficients $q_n$ are defined by

$$ q_n = \begin{cases} \frac{-1}{2\beta_1\beta_2\beta_3} \left[ \text{Re} \left\{ \frac{p_1^n}{(p_1 - \bar{p}_2)(p_1 - \bar{p}_3)} + \dots \right\} - \delta_{n2} \right] \\ \frac{-1}{2\beta_1\beta_2\beta_3} \left[ \text{Re} \left\{ \frac{p_1^{n-2}\bar{p}_2\bar{p}_3}{(p_1 - \bar{p}_2)(p_1 - \bar{p}_3)} + \dots \right\} - \delta_{n2} \right] \\ \text{for } n = \begin{cases} 0, 1, 2 \\ 2, 3, 4 \end{cases} \end{cases} \quad (14) $$

where $\delta_{n2} = 1$ if $n=2$ and $\delta_{n2} = 0$ if $n \neq 2$. $q_2$ has two alternative expressions. The infinite part of the stresses is deduced from

$$ \begin{aligned} \sigma_{ij} &= C_{ijkls} \varepsilon_{ks} \\ &= \frac{1}{2} C_{ijkls} (u_{k,s} + u_{s,k}) \end{aligned} \quad (15) $$

To calculate the complementary part of the displacements and the stresses, a vector c is needed

$$ c = Da \quad (16) $$

The matrix D is defined by

$$ \begin{aligned} D(1,j) &= C_{11jx}n_x + pC_{11j3} \\ D(2,j) &= C_{12jx}n_x + pC_{12j3} \\ D(3,j) &= C_{22jx}n_x + pC_{22j3} \end{aligned} \quad (17) $$
with $j = 1, 2, 3$ and $\vec{n} = (\cos \theta, \sin \theta, 0)^T$

and

$$ c_{j+3} = \bar{c}_j (j = 1,2,3), \quad C = [c_1, c_2, c_3] \quad (18) $$

A matrix G is determined by

$$ G = \bar{B}^{-1}B \quad (19) $$

The traction vector $\vec{t}$ on the plane $x_3$ = constant and the in-plane stress vector $\vec{s}$ are deduced from

$$ \vec{t} \equiv (\sigma_{13}, \sigma_{23}, \sigma_{33}) = (C_{13ks}u_{k,s}, C_{23ks}u_{k,s}, C_{33ks}u_{k,s}) \quad (20) $$

$$ \vec{s} \equiv (\sigma_{11}, \sigma_{12}, \sigma_{22}) = (C_{11ks}u_{k,s}, C_{12ks}u_{k,s}, C_{22ks}u_{k,s}) \quad (21) $$

Assuming $x_3 \neq 0$, the displacements are reduced to [8]

$$ \vec{u}(x_1, x_2, x_3) = \vec{u}_{\infty}(x_1, x_2, x_3) + \frac{1}{4\pi^2} \left[ \int_0^{2\pi} \bar{A}G_u A^T d\theta \right] \vec{f} \quad (22) $$

where the matrix $G_u$ is defined by

$$ (G_u)_{ij} = \frac{(G)_{ij}}{-\bar{p}_i x_3 - (x_1 \cos \theta + x_2 \sin \theta)} \quad (23) $$

Similarly, the Green's stresses can be written as

$$ \begin{aligned} \vec{t}(x_1, x_2, x_3) &= \vec{t}_{\infty}(x_1, x_2, x_3) + \frac{1}{4\pi^2} \left[ \int_0^{2\pi} \bar{B}G_t A^T d\theta \right] \vec{f} \\ \vec{s}(x_1, x_2, x_3) &= \vec{s}_{\infty}(x_1, x_2, x_3) + \frac{1}{4\pi^2} \left[ \int_0^{2\pi} \bar{C}G_t A^T d\theta \right] \vec{f} \end{aligned} \quad (24) $$

with the matrix $G_t$ as

$$ (G_t)_{ij} = \frac{(G)_{ij}}{[-\bar{p}_i x_3 - (x_1 \cos \theta + x_2 \sin \theta)]^2} \quad (25) $$

If the studied point is at the surface of the anisotropic half-space ($x_3 = 0$), the formulation is slightly different because of a singularity in the complementary part. The Green's displacements then become

$$ \vec{u}(x_1, x_2, 0) = \vec{u}_{\infty}(x_1, x_2, 0) - \frac{1}{4\pi r} \times \left\{ \frac{1}{\pi} \int_0^{2\pi} \frac{\bar{A}GA^T}{\cos(\theta - \theta_0)} d\theta + i[\bar{A}GA^T]_{\theta=\theta_0 \pm \pi/2} \right\} \vec{f} \quad (26) $$

by using the Cauchy principal value for calculating the integral. Green's stresses, where the integral is defined by the Hadamard finite part, can be derived as

$$ \begin{aligned} \vec{t}(x_1, x_2, 0) &= \vec{t}_{\infty}(x_1, x_2, 0) + \frac{1}{4\pi r^2} \times \left\{ \frac{1}{\pi} \int_0^{2\pi} \frac{\bar{B}GA^T}{\cos^2(\theta - \theta_0)} d\theta \mp i \frac{d[\bar{B}GA^T]}{d\theta} \Big|_{\theta=\theta_0 \pm \pi/2} \right\} \vec{f} \end{aligned} \quad (27) $$

$$ \begin{aligned} \vec{s}(x_1, x_2, 0) &= \vec{s}_{\infty}(x_1, x_2, 0) + \frac{1}{4\pi r^2} \times \left\{ \frac{1}{\pi} \int_0^{2\pi} \frac{\bar{C}GA^T}{\cos^2(\theta - \theta_0)} d\theta \mp i \frac{d[\bar{C}GA^T]}{d\theta} \Big|_{\theta=\theta_0 \pm \pi/2} \right\} \vec{f} \end{aligned} \quad (28) $$

# 4 Application

## 4.1 Validation for an Isotropic Case.
The first validation is a comparison between the results obtained using the anisotropic model when dealing with an isotropic half space submitted to an Hertzian loading and the analytical solution given in Eqs. (29) and (30). Assuming a rigid indenter, with a spherical tip of radius $R_{\text{indenter}}$, in contact with an isotropic half space, the Hertz pressure $P_{\text{Hertz}}$ and the Hertz contact radius $a_{\text{Hertz}}$ can be written as

$$ P_{\text{Hertz}} = \frac{3P}{2\pi a_{\text{Hertz}}^2} \quad (29) $$

$$ a_{\text{Hertz}} = \left( \frac{3PR_{\text{indenter}}(1 - \nu_{\text{solid}}^2)}{4E_{\text{solid}}} \right)^{(1/3)} \quad (30) $$

$E$ and $\nu$ are the Young's modulus and the Poisson's ratio of the isotropic half space, respectively, $P$ represents the load. The material is assumed cubic, almost isotropic, because isotropic materials lead to double roots which create singularities. The Coulomb's modulus is equal to the isotropic Coulomb's modulus plus 0.7%. Two meshes are studied, the first one 47 × 47 × 25 elements (the
<table><caption>**Table 1** The influence of the mesh and the model on the CPU time</caption><tr><td>Element number in the contact area</td><td>23 × 23</td><td>71 × 71</td></tr><tr><td>Mesh</td><td>47 × 47 × 25</td><td>95 × 95 × 50</td></tr><tr><td>CPU isotropic model</td><td>15 s</td><td>1 min 49 s</td></tr><tr><td>CPU anisotropic model</td><td>2 min 37 s</td><td>20 min 25 s</td></tr></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Pressure profile for an isotropic case</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Finite element model with a detailed view of the contact area</figcaption></figure>

two first numbers refer to the directions parallel to the surface and the third one to the depth) with 23 × 23 elements within the contact area, and the second one, a refined mesh 95 × 95 × 50 elements, with 71 × 71 elements within the contact area. Note that for contact problems it is usually admitted that 10 points are usually sufficient to accurately describe a specific wavelength.

The gap between the Hertz solution and the numerical solution is 0.6%, and between the standard mesh and refined mesh is 1.6%, for the pressure. The (anisotropic) numerical solution tends towards the (isotropic) analytical solution for isotropic properties.

The difference between computing time with isotropic and anisotropic models stems from the calculation of the influence

<table><caption>**Table 2** Number and type of elements in the FE model</caption><thead><tr><td>Area</td><td>Type</td><td>Quantity</td></tr></thead><tbody><tr><td>Sphere</td><td>C3D10</td><td>22276</td></tr><tr><td>Contact area of the half space</td><td>C3D20</td><td>54000</td></tr><tr><td>The area closest to the contact area</td><td>C3D4</td><td>49656</td></tr><tr><td>The largest area</td><td>C3D8R</td><td>22942</td></tr><tr><td>Edges of the half space</td><td>CIN3D8</td><td>3448</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Pressure profile for an orthotropic material</figcaption></figure>

<table><caption>**Table 3** Elastic parameters of studied materials</caption><thead><tr><th>Figs. 3(a), 6, 8<br>Figs. 9, 10 (for θ = 0 deg)</th><th>Figs. 3(b), 4, 5, 7</th><th>Figs. 11, 12</th></tr></thead><tbody><tr><td>E<sub>1</sub>/E<sub>2</sub> = 1</td><td>E<sub>1</sub>/E<sub>2</sub> ≠ 1</td><td>E</td></tr><tr><td>E<sub>3</sub>/E<sub>2</sub> ≠ 1</td><td>E<sub>3</sub>/E<sub>2</sub> = 1</td><td></td></tr><tr><td>ν<sub>12</sub> = ν<sub>13</sub> = ν<sub>23</sub>, ν<sub>ij</sub> = ν<sub>ji</sub>.E<sub>i</sub>/E<sub>j</sub></td><td></td><td>ν</td></tr><tr><td>G<sub>12</sub> = G<sub>13</sub> = G<sub>23</sub> = E<sub>2</sub>/(2(1 + ν<sub>12</sub>))</td><td></td><td>G</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig. 4 Influence of $E_1$ on the contact pressure</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5 Influence of $E_1$ on the contact area</figcaption></figure>

coefficients, whose calculation is more complex and lengthy for the latter model. Note that the influence coefficients should be computed only once for a given anisotropic material, and stored to be used for further simulations. By using a finite element model, presented afterwards, results for an anisotropic half space are obtained with a CPU time of 144 h on the same class of the personal computer (Intel® Core™ Duo CPU T9600 @ 2.80 GHz 1.59 GHz, 3.48 Go RAM). According to Table 1 and Fig. 1, the standard mesh is a good compromise, it gives accurate results in a reasonable computing time. It will be used in the following parts.

## 4.2 Validation for an Anisotropic Case.
A finite element model has been developed with the commercial FE packageABAqus (version 6.9), in order to validate the anisotropic semi-analytical method.

An isotropic sphere, with Young's modulus $10^6$ times larger than the half space Young's modulus and identical Poisson's ratio, is in contact with a semi-infinite half space. The mesh is shown in Fig. 2 and described in Table 2. It includes semi-infinite elements

<figure><img src="image_7.png" /><figcaption>Fig. 6 Influence of $E_3$ on the contact pressure</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Influence of $E_1$ on the maximum pressure</figcaption></figure>

CIN3D8 on the edges of the body, and quadratic elements with 20 nodes and 27 integration points (C3D20) in the contact area. In order to respect the Hertz conditions, the sphere radius is 30 times larger than the contact radius. To optimize the computation time, properties of symmetry along the direction 1 are used.

The first step consists of validating the FE mesh through a comparison of the contact solution obtained for a Hertzian loading. The error for the pressure distribution (i.e., on the maximum pressure and the contact radius) is smaller than 0.4%.

Then, two anisotropic cases are studied. The half space is considered orthotropic, with the same Poisson's ratio and the same Coulomb's modulus. For the first case, the Young's modulus along direction 3 (the depth) is twice as large as the Young's modulus in directions 1 and 2. The FE and SAM pressure distributions are plotted in Fig. 3(a). The error for the maximal value is lower than 0.2%. For the second case, the Young's modulus along direction 1 (on the surface) is twice as large as those in directions 2 and 3. A difference of 2.2% for the maximum pressure is observed in Fig. 3(b).

## 4.3 Parametric Studies.
The influence of the properties of an elastic anisotropic half space which is in contact with a rigid indenter, with a spherical tip, is studied here. The contact parameters, i.e., the contact area and the pressure distribution, will be more specifically investigated. The depth corresponds to direction
<figure><img src="image_9.png" /><figcaption>Fig. 8 Influence of $E_3$ on the maximum pressure</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9 Influence of the material's orientation on the contact pressure (for $\theta_m = 0°$, $E_3/E_2 = 3$)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10 Influence of the material's orientation on the maximum pressure (Caption is valid for $\theta_m = 0°$)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11 Influence of the Poisson's ratio on the contact area</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12 Influence of the Coulomb's modulus on the contact area</figcaption></figure>

3, which means that the surface is defined by directions 1 and 2. In these examples, materials are orthotropic, with the same Poisson's ratio and the same Coulomb's modulus. The material data used in the numerical simulations are synthetized in Table 3.

The influence of $E_1$ and $E_3$ on the contact pressure distribution is shown in Figs. 4–6. It is observed that a change of the Young's modulus along a direction parallel to the surface ($E_1$ here) has a moderate effect (Fig. 4). The maximum contact pressure is increased by 3% when $E_1$ is doubled. Moreover it can be observed that the contact area is no more circular (Fig. 5). The effect of the Young's modulus along the depth ($E_3$) is more pronounced, an increase of 32% of the maximum contact pressure and decrease of the contact radius by 11% when the Young's modulus is increased by a factor 2 (Fig. 6). Figures 7 and 8 summarize the influence of the Young's modulus on the maximum pressure. A slight modification of $E_1$ around $E_{isotropic}$ does not have a significant effect,
<figure><img src="image_14.png" /><figcaption>Fig. 13 Influence of the Young's modulus along axis 1 (parallel to the surface) on the indentation curve</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14 Influence of the Young's modulus along axis 3 (normal to the surface) on the indentation curve</figcaption></figure>

conversely to $E_3$ that does have a pronounced effect. Nevertheless, a monotonic increase of the maximum contact pressure is observed when $E_1$ increases (see Fig. 7). The maximum pressure drops quickly when $E_3$ decreases (see Fig. 8).

The effect of the material's orientation relative to the contact is shown in Fig. 9. Here the Young's modulus of the material in direction 3 is equal to $E_3 = 3 E_{\text{isotropic}}$, but the material main direction is different from that of the contact. When the orientation angle $\theta_m$ around the 1-axis increases up to 90°, the numerical solution converges progressively to the solution where $E_3 = E_{\text{isotropic}}$ and $E_2 = 3 E_{\text{isotropic}}$. This is confirmed in Fig. 10, where the maximum pressure is plotted versus the angle, for different ratios between $E_3$ and $E_{\text{isotropic}}$. Hence, the method is valid whatever the orientation angle compared to the surface.

In Figs. 11 and 12, the influence of the Poisson's ratio $\nu$ and the Coulomb's modulus $G$ are observed more precisely. The studied material is cubic. Dividing (multiplying) by 2 the Poisson's ratio increases by 6% (decreases by 14%) the contact radius. In terms of contact area it corresponds to an increase of 12% or a decrease of 30% for the contact area, respectively. A modification of the Coulomb's modulus by a factor 2 changes the contact radius by 6% if $G > G_{\text{isotropic}}$ or 11% if $G < G_{\text{isotropic}}$. Therefore $\nu$ and $G$ have a significant influence.

Finally it is also interesting to have a look on the effect of the anisotropy on the normal load versus displacement response, as monitored during an indentation test. It is confirmed in Fig. 13 that a change of the Young's modulus along a direction parallel to the surface hardly affects the load-displacement curve. Conversely one may observe a very significant modification of the indentation curve when the Young's modulus along the direction corresponding to depth is varying, see Fig. 14. Note that the use of the equivalent Young's modulus as proposed by Swanson [14] does not allow us to reproduce these indentation curves.

## 5 Conclusions

A semi-analytical method has been developed for the contact problem of anisotropic elastic materials. The model has been validated by comparison with the Hertz solution for isotropic materials and by comparison with a FE model for anisotropic materials. It is found that the stiffness along the normal to the contact has a strong influence on the contact solution in terms of pressure distribution and contact size; an increase in $E_3$ leads to a higher maximum contact pressure and a smaller contact radius. Conversely a change of the Young's modulus along a direction parallel to the surface (plane (1, 2)) does not significantly affect the contact pressure distribution but the contact area is no more circular. The performance of the method is highlighted by analyzing the effect of the orientation of the material main directions compared to the surface normal.

## Acknowledgment

This work was supported by the French FUI Project Innolub.

## Nomenclature

$a_{\text{Hertz}}$ = Hertzian contact radius

$a_i, b_i, c_i$ = eigenvectors associated with $p_i$

$A, B, C$ = matrices of the eigenvectors

$C_{ijkl}$ = elastic stiffness tensor

$E$ = Young's modulus

$E_i$ = Young's modulus in direction $i$

$G$ = Coulomb's modulus

$\bar{M}$ = conjugate of $M$

$M^{-1}$ = inverse of $M$

$M^T$ = transpose of $M$

$\nu$ = Poisson's ratio

$P$ = contact pressure distribution

$P_{\text{Hertz}}$ = maximum Hertzian pressure

$p_i$ = eigenvalues

$Q, R, T$ = double projections of the elastic stiffness tensor

$r, \theta, \phi$ = spherical coordinates of the studied point

$\sigma_{ij}$ = stress tensor

$\vec{s}$ = in-plane stress vector

$\vec{l}$ = traction vector

$\theta_m$ = orientation angle of the material

$\vec{u}$ = displacement vector

$u_{\text{comp}}$ = relative to the complementary part

$u_{\infty}$ = relative to the infinite part

$u_{k,s}$ = partial derivative of $u_k$ by $x_s$

$x_i$ = Cartesian coordinates of the studied point

## References

[1] Kelvin, L., 1848, "Note on the Integration of the Equations of Equilibrium of an Elastic Solid," Cambridge and Dublin Math. J., **3**, pp. 87–89.

[2] Boussinesq, J., 1885, *Application des Potentiels a l'étude de l'équilibre et du Mouvement des Solides Elastiques*, Gauthier-Villars.

[3] Mindlin, R., 1936, "Force at a Point in the Interior of a Semi-Infinite Solid," Physics (N.Y.), **7**(5), pp. 195–202.

[4] Ting, T., and Lee, V., 1997, "The Three-Dimensional Elastostatic Green's Function for General Anisotropic Linear Elastic Solids," Q. J. Mech. Appl. Math., **50**(3), pp. 407–426.

[5] Stroh, A., 1958, "Dislocations and Cracks in Anisotropic Elasticity," Philos. Mag., **3**(30), pp. 625–646.
[6] Stroh, A., 1962, "Steady State Problems in Anisotropic Elasticity," J. Math. Phys., **41**(2), pp. 77-103.

[7] Ting, T., 1996, *Anisotropic Elasticity: Theory and Applications*, Oxford University Press, NY.

[8] Pan, E., and Yuan, F., 2000, "Three-Dimensional Green's Functions in Anisotropic Bimaterials," Int. J. Solids Struct., **37**(38), pp. 5329-5351.

[9] Yang, B., and Pan, E., 2002, "Three-Dimensional Green's Functions in Anisotropic Trimaterials," Int. J. Solids Struct., **39**(8), pp. 2235-2255.

[10] Pan, E., and Yang, B., 2003, "Three-Dimensional Interfacial Green's Functions in Anisotropic Bimaterials," Appl. Math. Model., **27**(4), pp. 307-326.

[11] Ciavarella, M., Demelio, G., Schino, M., and Vlassak, J., 2001, "The General 3D Hertzian Contact Problem for Anisotropic Materials," Key Eng. Mater., **221**, pp. 281-292.

[12] Li, X., and Wang, M., 2006, "Hertzian Contact Of Anisotropic Piezoelectric Bodies," J. Elast., **84**(2), pp. 153-166.

[13] Borodich, F., 2000, "Some Contact Problems of Anisotropic Elastodynamics: Integral Characteristics and Exact Solutions," Int. J. Solids Struct., **37**(24), pp. 3345-3373.

[14] Swanson, S., 2004, "Hertzian Contact of Orthotropic Materials," Int. J. Solids Struct., **41**, pp. 1945-1959.

[15] Willis, J., 1966, "Hertzian Contact of Anisotropic Bodies," J. Mech. Phys. Solids, **14**(3), pp. 163-176.

[16] Pagano, N., 1970, "Exact Solutions for Rectangular Bidirectional Composites and Sandwich Plates," J. Compos. Mater., **4**(1), pp. 20-34.

[17] Aizikovich, S., Alexandrov, V., Kalker, J., Krenev, L., and Trubchik, I., 2002, "Analytical Solution of the Spherical Indentation Problem for a Half-Space With Gradients With the Depth Elastic Properties," Int. J. Solids Struct., **39**(10), pp. 2745-2772.

[18] Gao, Y., and Pharr, G., 2007, "Multidimensional Contact Moduli of Elastically Anisotropic Solids," Scr. Mater., **57**(1), pp. 13-16.

[19] Rand, O., and Rovenskii, V., 2005, *Analytical Methods in Anisotropic Elasticity: With Symbolic Computational Tools*, Birkhauser.

[20] Galin, L., 2008, *Contact Problems: The Legacy of LA Galin*, Springer Verlag.

[21] Lin, Y., and Ovaert, T., 2004, "A Rough Surface Contact Model For General Anisotropic Materials," J. Tribol., **126**(1), pp. 41-49.

[22] He, L., and Ovaert, T., 2008, "Three-Dimensional Rough Surface Contact Model for Anisotropic Materials," J. Tribol., **130**(2), pp. 021402.

[23] Barnett, D., and Lothe, J., 1975, "Line Force Loadings on Anisotropic Half-Spaces and Wedges," Phys. Norv., **8**(1), pp. 13-22.

[24] Kahya, V., Ozsahin, T., Birinci, A., and Erdol, R., 2007, "A Receding Contact Problem for an Anisotropic Elastic Medium Consisting of a Layer and a Half Plane," Int. J. Solids Struct., **44**(17), pp. 5695-5710.

[25] Batra, R., and Jiang, W., 2008, "Analytical Solution of the Contact Problem of a Rigid Indenter and an Anisotropic Linear Elastic Layer," Int. J. Solids Struct., **45**(22), pp. 5814-5830.

[26] Argatov, I., 2011, "Depth-Sensing Indentation of a Transversely Isotropic Elastic Layer: Second-Order Asymptotic Models for Canonical Indenters," Int. J. Solids Struct., **48**(25-26), pp. 3444-3452.

[27] Brock, L., and Georgiadis, H., 2007, "Multiple-Zone Sliding Contact With Friction on an Anisotropic Thermoelastic Half-Space," Int. J. Solids Struct., **44**(9), pp. 2820-2836.

[28] Clements, D., and Ang, W., 2009, "On Some Contact Problems For Inhomogeneous Anisotropic Elastic Materials," Int. J. Eng. Sci., **47**(11-12), pp. 1149-1162.

[29] Blázquez, A., Mantić, V., and Pars, F., 2006, "Application of BEM to Generalized Plane Problems for Anisotropic Elastic Materials in Presence of Contact," Eng. Anal. Boundary Elem., **30**(6), pp. 489-502.

[30] Rodriguez-Tembleque, L., Buroni, F., Abascal, R., and Sáez, A., 2011, "3D Frictional Contact of Anisotropic Solids Using BEM," Eur. J. Mech. A/Solids, **30**(2), pp. 95-104.

[31] Jacq, C., Nélias, D., Lormand, G., and Girodin, D., 2002, "Development of a Three-Dimensional Semi-Analytical Elastic-Plastic Contact Code," J. Tribol., **124**(4), pp. 653-667.

[32] Boucly, V., Nelias, D., Liu, S., Wang, Q., and Keer, L., 2005, "Contact Analyses for Bodies with Frictional Heating and Plastic Behavior," J. Tribol., **127**(2), pp. 355-364.

[33] Chen, W., and Wang, Q., 2008, "Thermomechanical Analysis of Elastoplastic Bodies in a Sliding Spherical Contact and the Effects of Sliding Speed, Heat Partition, and Thermal Softening," J. Tribol., **130**(4), pp. 041402.

[34] Fulleringer, B., and Nélias, D., 2010, "On the Tangential Displacement of a Surface Point Due to a Cuboid of Uniform Plastic Strain in a Half-Space," J. Appl. Mech., **77**(2), pp. 021014.

[35] Leroux, J., Fulleringer, B., and Nelias, D., 2010, "Contact Analysis in Presence of Spherical Inhomogeneities Within a Half-Space," Int. J. Solids Struct., **47**(22-23), pp. 3034-3049.

[36] Zhou, K., Keer, L., and Wang, Q., 2011, "Semi-Analytic Solution for Multiple Interacting Three-Dimensional Inhomogeneous Inclusions of Arbitrary Shape in an Infinite Space," Int. J. Numer. Methods Eng., **87**(7), pp. 617-638.

[37] Leroux, J., and Nélias, D., 2011, "Stick-Slip Analysis of a Circular Point Contact Between a Rigid Sphere and a Flat Unidirectional Composite With Cylindrical Fibers," Int. J. Solids Struct., **48**, pp. 3510-3520.

[38] Gallego, L., Nelias, D., and Deyber, S., 2010, "A Fast and Efficient Contact Algorithm for Fretting Problems Applied to Fretting Modes i, ii, and iii," Wear, **268**(1-2), pp. 208-222.

[39] Nélias, D., Boucly, V., and Brunet, M., 2006, "Elastic-Plastic Contact Between Rough Surfaces: Proposal for a Wear or Running-in Model," J. Tribol., **128**(2), pp. 236-244.

[40] Gallego, L., Nélias, D., and Jacq, C., 2006, "A Comprehensive Method to Predict Wear and to Define the Optimum Geometry of Fretting Surfaces," J. Tribol., **128**(3), pp. 476-485.

[41] Gallego, L., and Nélias, D., 2007, "Modeling of Fretting Wear Under Gross Slip and Partial Slip Conditions," J. Tribol., **129**(3), pp. 528-535.

[42] Gallego, L., Fulleringer, B., Deyber, S., and Nelias, D., 2010, "Multiscale Computation of Fretting Wear at the Blade/Disk Interface," Tribol. Int., **43**(4), pp. 708-718.

[43] Chaise, T., and Nélias, D., 2011, "Contact Pressure and Residual Strain in 3D Elasto-Plastic Rolling Contact for a Circular or Elliptical Point Contact," J. Tribol., **133**(4), pp. 041402.

[44] Chaise, T., Nélias, D., and Sadeghi, F., 2011, "On the Effect of Isotropic Hardening on the Coefficient of Restitution for Single or Repeated Impacts Using a Semi-Analytical Method," Tribol. Trans., **54**(5), pp. 714-722.