Contact mechanics with a multiscale approach

Guillaume Anciaux, Jean-François Molinari

► To cite this version:

Guillaume Anciaux, Jean-François Molinari. Contact mechanics with a multiscale approach. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422155

HAL Id: hal-01422155

https://hal.science/hal-01422155

Submitted on 23 Dec 2016

**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 mechanics with a multiscale approach

Anciaux Guillaume¹*, Molinari Jean-François¹

¹ Ecole Polytechnique Fédérale de Lausanne (EPFL)
Faculté ENAC-IS, Laboratoire de Simulation en Mécanique des Solides (http://lsms.epfl.ch/)
CH-1015 Lausanne, Switzerland

* corresponding author : guillaume.anciaux@epfl.ch

**Abstract** — We present a multiscale coupling method to address contact problems. The components of the model are a molecular dynamics engine, a finite element program and a coupling scheme. We validate the approach on an Hertzian contact. Various measures are provided to highlight limitations and new opportunities in conducting large scale simulations of contact brought by the proposed multiscale approach.

**Keywords** — Hertz contact, MultiScale, Concurrent coupling, Molecular dynamics.

## 1 Introduction

Be they at the macroscopic or at the nano scale, tribological problems are particularly difficult to comprehend as different physical mechanisms interact at disparate scales. While contact mechanics is traditionally described as a macroscopic phenomena with empirical or semi-empirical phenomenological laws (see for instance Coulomb friction law, or Archard wear law), recent efforts are increasingly attempting to capture the statistics of contact forces at contacting asperities all the while incorporating the influence of large-scale elastic fields. This separation of scale presents a unique challenge that has limited, so far, the accuracy of numerical simulations of friction and wear phenomena. Nonetheless, as experiments go down in size with a recent focus on nanotribology, as computational power expands, and with the emergence of multiscale coupling schemes, numerical simulations become increasingly relevant to experimental work. This paper presents a novel approach for handling 3D contact at the nanoscale. It is composed of a direct multiscale platform in which a molecular dynamics package is used for interactions of nano-asperities at the free surface, and a continuum model (Finite Elements) for long range elastic interactions.

Finite Element Methods (FEM) are standard techniques for solving continuum equations. But continuum mechanics is unable to capture the details of force profiles at this scale. In order to represent efficiently the atomic organization and forces at contact clusters, one has often recourse to Molecular Dynamics (MD).

Molecular Dynamics (MD) is a well-established single scale approach to simulate materials at the atomic scale [1], atomic origin of adhesive forces [2] and of friction [3, 4, 5, 6]. The main problem concerning MD approaches lies in computer manipulation of the number of atoms needed to create simulations of consistent sizes.

The limitations of purely atomistic or purely continuum simulations, have motivated research in multi-scale simulations that bridge atomistic simulation and continuum modeling [7]. In this direct multiscale approach, atoms can be used at contacting asperities to capture in great details contact forces, whereas finite elements are used away from the interface to accurately represent
<table><thead><tr><th></th><th>Microscopic units</th><th>Macroscopic units</th></tr></thead><tbody><tr><td>Elastic Modulus [100]</td><td>0.6628 eV</td><td>106.19 GPa</td></tr><tr><td>Poisson ratio</td><td>0.3468</td><td>0.3468</td></tr><tr><td>Mass density</td><td>8.93357e-27 Kg/³</td><td>63.546 g/mol</td></tr></tbody></table>

FIG. 1 – Table of macroscopic properties of the copper crystal studied here. These values were computed by means of the inter-atomic EAM potential.

elastic forces. The main purpose of the coarse domain is to reduce significantly the number of unknowns to handle.

The presented work concentrate on contact mechanics by use of the overlapping scheme presented in [8] with the modification brought by [9]. Here, our main objective is to prove the technical feasibility of using coupled schemes to handle contact problems in 3D [10]. We will show that it is important to replace those missing atoms by (a coarse) continuum domain to account for the elasticity of the full system size. The organization of this paper is as follows. The first section describes in details the implemented model. Then the second section presents a validation example consisting of a hertzian contact case where a spherical rigid body, constituted of real atoms, is indenting a copper crystal. We conclude with a discussion on perspectives.

## 2 Model implementation

The purpose of this section is to describe the molecular dynamics properties relevant to the presented work. All initial states will consider a perfect copper crystal. The considered lattice is cubic face centered (CFC) with dimension $a = 3.615$ Å. The mass of each copper atom is $m = 63.55$ g/mol. The cutoff radius is set to $R_{cut} = 4.95$ Å. Thus, atomic forces are computed only taking into account for the first three neighbors. The interatomic potential used to drive atoms through molecular dynamics is the *Embedded Atom Method* (EAM) [11]. The expression of the associated global potential energy is :

$$E_{tot} = \frac{1}{2} \sum_{i,j \neq i} V(r_{ij}) + \sum_{i} F(\rho_{i}) \quad \text{with} \quad \rho_{i} = \sum_{j \neq i} \phi(r_{ij}) \quad (1)$$

where $V$ is a purely attractive pair potential, $\phi(r)$ is the electronic density approximation function of a neighbor atom with respect to the distance $r$, $\rho_i$ is the electronic density associated with atom $i$, $\mathbf{r}_i$ is the position of atom $i$, $\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i$ and $r_{ij}$ is the norm of the vector $\mathbf{r}_{ij}$.

We use the well known linear momentum balance expression for a continuum with the Lagrangian formulation and no body forces :

$$\nabla \cdot \underline{P} = \rho_0 \dot{\mathbf{v}} \quad (2)$$

where $\underline{P}$ is the first Piola-Kirchhoff stress tensor, $\rho_0$ is the mass per unit of volume at initial time and $\mathbf{v}$ is the velocity field. As a constitutive law for the crystal copper material we use the Cauchy-Born rule (CBR) [12]. The Cauchy-Born rule assumes an homogeneous deformation of the crystal lattice to compute the Piola-Kirchhoff tensor. The derivation of the inter-atomic potential EAM yielded the elastic modulus and Poisson ratio. These elastic properties, as well as wave velocities and mass density are presented in table 1 and they match usual copper constants. Thus, using this method based on the interatomic potential insures a good fit of material properties between the continuum and molecular dynamics.

The coupling method, which was implemented, is based on an overlapping region (Bridging Zone) where the continuum and the atomistic models are both valid. They are combined by using a weighting function for the energy. The two models are glued together by constraining the degrees
of freedom in order to have coherent dynamics. The overlap allows us to change the predominance
of each model so that the constraint can take effect gradually when the overlap is crossed and
consequently minimize the effect of wave reflections. A more detailed description can be found in
[8, 9].

3 Hertz contact simulation

The purpose of this section is to validate the coupling scheme on an hertzian contact simula-
tion. We will introduce a spherical indenter (rigid body) that will be put in contact with a substrate
at the high displacement rate of 10 m/s = 10⁻⁴/fs. The timestep is defined to Δt = 2 fs. To
measure the accuracy of the method all results from this section will be compared to an equiva-
lent simulation performed with molecular dynamics (complete domain is filled with atoms). We
describe the model before the adhesive Hertzian contact. All results presented in this article were
carried out using the open LibMultiScale <fnref n="1" /> parallel coupling framework.

The dimension of the body (substrate) to
be simulated is 90a × 90a × 110a where per-
iodic boundary conditions are used in X and
Y directions. The bottom atoms (in pure mo-
lecular dynamics) or the bottom surface (in
coupled simulation) are held fixed to provide
a boundary condition while the top surface is
leaved as free. When the coupled scheme will
be used, the atomic zone has for dimensions
90a × 90a × 28a - containing 1 014 747 atoms
- with a planar coupling zone as deep as 18a
oriented in the Z direction. Figure 2 presents
the geometry of the substrate. One important
note must be added here : the area involving
uncoupled atoms has a very small number of
atomic layers (10 a). This was done on pur-
pose in order to show the limitation of the cou-
pling method with regard to dislocation pro-
pagation. It is straightforward to enlarge this
area - keeping the surface section constant to
allow contact mechanics simulation to run fur-
ther with no misrepresentation of the disconti-
nuities.

The mesh used here is constituted of te-
trahedral elements with constant sizes. It was

built by decomposition of cubic elements in tetrahedra. Each initial cubic element had the dimensions 3a × 3a × 3a. The mesh contains 324 000 elements and 442 981 nodes, where 64 800 elements and 583 200 atoms are in the Bridging zone. This gives 9 atoms per elements in the coupling area. The bridging size which is set to 18a involves 36 atomic layers. This is a large number in comparison to the overlapping used elsewhere [13]. This size was chosen to obtain a good coupling quality with regards to wave reflection at the interface between molecular dynamics and finite elements (18a is in the low range of bridging zone thickness for minimizing reflections as shown in [9]). We recall again that by having a thin layer of uncoupled atoms close to the surface, we want to demonstrate the limits of the actual method, and especially with regards to the

<figure><img src="image_2.png" /><figcaption>FIG. 2 – Schematic of the coupled domain (substrate) where $a = 3.615$ is the CFC lattice parameter.</figcaption></figure>

<footnote n="1">¹http://libmultiscale.gforge.inria.fr/index.html</footnote>
<figure><img src="image_3.png" /></figure>

FIG. 3 – Kinetic energy measured in the region where contact occurs. The curves represent kinetic energy of that region in the cases where coupling is used, full MD is used and when a reduced MD sample is used.

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

FIG. 4 – Stress component $\sigma_{33}$ for both models. The stress is measured along the line B (resp. A and C) - see figure 2 - on the top (resp. middle and lower) graphic.

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

FIG. 5 – A cut in the middle of the sample along axis Y is presented through the evolution of the simulation. The gray scale presents the magnitude of the displacement field. The left side beam is the solution obtained with the coupled scheme while the right one is the purely molecular dynamics solution. (a) is given for a displacement $\delta$ of 3 Å and (b) for 5 Å.

issue of dislocation propagation. It should also be noted that the subsequent gain brought by the coupling divided the number of manipulated atoms by three.

We bring in contact the substrate body with a rigid body constituted of real copper atoms. This rigid body has been built first by constructing a CFC copper crystal in an identical manner as the substrate. Then the crystal was cut according to a spherical geometry of curvature $R = 50a$. Such a structure is stepped but, as the effect of atomic structure of the indenter was studied by Luan and Robbins [14] and will not be included in this paper.

Now, we present the numerical results obtained when considering that forces between rigid body atoms and substrate atoms are computed by means of the EAM potential described earlier. As our purpose is the validation of the coupling scheme, we need to exhibit measures relevant enough to be compared. We consider for instance the comparison of the residual kinetic energy close to the contact surface. Indeed, when the coupling method is used, well-known spurious wave reflections can impair part of the energy of the molecular domain. Figure 3 shows the kinetic energy of the common molecular zone in three cases. The two first are the coupled and
<figure><img src="image_6.png" /><figcaption>FIG. 6 – Atoms extracted from centrosymmetry criterion with a threshold of 0.1. (a) corresponds to the situation when displacement equals 5 and (b) 6. The lines shows the beginning of the coupling zone.</figcaption></figure>

the fully atomistic simulations. The last one is a simulation for which the MD zone is the same as in the coupled case but the boundary condition has been replaced by fixed displacements on the three bottom layers (instead of a coupling to continuum mechanics). It can be seen that the behavior of the two first simulations are similar until a specific displacement has been reached (indicated with an arrow on the figure). In the third case, the use of a reduced set of atomic layers increased the system rigidity which led to an accumulation of kinetic energy in the contact zone. This demonstrates that the domain size has an important influence on the flow of energy density and therefore on contact mechanics. The use of a coupling method that dissipates and forwards kinetic energy allows to model large system sizes which is paramount to capture the correct system response during contact.

Another important information relative to coupling method quality is the stress continuity. We define the per atom pseudo-stress on atom i :

$$ \underline{\underline{\sigma}}^i = \sum_k \frac{f_{ik}}{r_{ik}} \mathbf{r}_{ik} \otimes \mathbf{r}_{ik} \qquad (3) $$

where the sum over k is over the neighboring atoms, $f_{ik}$ is the force magnitude exerted by the atom k on atom i. This quantity has the dimension of a pressure-distance. To get a fully local expression of the stress, a volume per atom has to be computed. Although this quantity can be approximated by means of Voronoi tessellation, it is more convenient to define a grid to compute the volumic density of the peratom pseudo stress in such a grid. The results for the atomic stress as well as the continuum stress, for displacement $\delta = 1.5$ are presented on figure 4. This figure shows the diagonal component $\sigma_{33}$ of the stress on the three lines A,B and C notified on figure 2. It is clear that continuity of the stress is ensured.

We now return to the divergence between models as observed from the kinetic energy plot (figure 3). To explain this lack of accuracy of the coupling scheme we present a cut into the substrate showing the inner displacement field. This is shown on figure 5 for indenter displacements 1 Å, 3 Å and 5 Å, which correspond to 0.28a, 0.83a and 1.38a respectively. On figure 5(b), it is clear that a discontinuity of the displacement field has been created by the compressive load. The nature of this discontinuity is fully identified as we extracted, on figure 6, the atoms in non symmetrical situation through the centrosymmetry criterion. It originates from a dislocation, which has been stopped by the coupling zone. While the dislocation front propagates deep through the atomic layers in pure MD simulation, it is annihilated in the coupled simulation as observed on figure 6(b).

## 4 conclusion

In the presented work, one of our main aim was to validate the use of our coupling scheme in the context of contact mechanics. The Hertzian contact study showed that the use of the LibMultiScale framework is nicely fitted to contact problems, when the assumption of continuity of the deformation is true in the coupling area. Indeed, coupled simulations compare reasonably well with purely molecular dynamics. The stress continuity as well as other indicators give a good confidence in the coupling scheme for contact mechanics.
We recognize that the coupling model is valid only when some assumptions are verified. The first limitation of the approach stems from dislocation propagations as well as plane slips interactions with the coupling zone. They should be minimized to avoid important inaccuracies of the computations. Secondly, no thermal effects are treated at present time. While molecular dynamics and continuum mechanics can take into account phonon and thermal diffusion, these two phenomena are of very distinct physics and scale. Thus, at present time, one has to limit the coupled simulation at a 0K temperature (or very close to zero).

Nonetheless, and despite of these limitations, in this work we produced a truly 3D parallel coupled simulation of contact problems. This brings us new opportunities to deal with realistic material samples involving up to billions of atoms and coupled in parallel with continuum mechanics.

Références

[1] D. C. Rapaport. *The art of Molecular Dynamics Simulation*. Cambridge University Press, 1995.

[2] H. Yao, M. Ciavarella, and H. Gao. Adhesion maps of spheres corrected for strength limit. *Journal of Colloid and Interface Science*, 315(2):786-790, 2007.

[3] H. J. Kim, W. K. Kim, M. L. Falk, and D. A. Rigney. Md simulation of microstructure evolution during high-velocity sliding between crystalline materials. *Tribology Letters*, 28:299-306, 2007.

[4] O. M. Braun and A. G. Nauvomets. Nanotribology : Microscopic mechanism of friction. *Surface Science Reports*, 60(6-7):79-158, 2006.

[5] H. H. Yu, P. Shrotriya, Y. F. Gao, and K. S. Kim. Micro-plasticity of surface steps under adhesive contact : Part i - surface yielding controlled by single-dislocation nucleation. *Journal of the Mechanics and Physics of Solids*, 55(3):489-516, 2007.

[6] J. Rottler and M. O. Robbins. Macroscopic friction laws and shear yielding of glassy solids. *Computer Physics Communications*, 169(1-3):177-182, 2005.

[7] W. A. Curtin and R. E. Miller. Atomistic/continuum coupling in computational materials science. *Modelling and simulation in materials science and engineering*, 11 :R33–R68, April 2003.

[8] S. P. Xiao and T. Belytschko. A bridging domain method for coupling continua with molecular dynamics. *Computer Methods in Applied Mechanical Engineering*, 193(17-20):1645–1669, 2004.

[9] G. Anciaux, O. Coulaud, J.Roman, and G. Zerah. Ghost force reduction and spectral analysis of the 1d bridging method. submitted to *Journal of Computational Physics*, 2008.

[10] G. Anciaux and J.F. Molinari. Contact mechanics at the nanoscale, a 3d multiscale approach. *International Journal for Numerical Methods in Engineering*, accepted, 2008.

[11] M. S. Daw and M. I. Baskes. Embedded atom method : Derivation and application to impurities, surfaces and other defects in metals. *Physical Review B*, 29:6443-6543, 1984.

[12] E. B. Tadmor, M. Ortiz, and R. Phillips. Quasicontinuum analysis of defects in solids. *Philosophical Magazine*, 73(6):1529-1563, 1996.

[13] B. Q. Luan, S. Hyun, J. F. Molinari, M. O. Robbins, and N. Bernstein. Multiscale modeling of two-dimensional contacts. *Physical Review E*, 74, 2006.

[14] B. Luan and M. O. Robbins. The breakdown of continuum models for mechanical contacts. *Nature*, 435:929-932, 2005.