Numerical simulations of nanotribology

Guillaume Anciaux, Pedro Romero, Jean-François Molinari

► To cite this version:

Guillaume Anciaux, Pedro Romero, Jean-François Molinari. Numerical simulations of nanotribology. 9e colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01419919

HAL Id: hal-01419919

https://hal.science/hal-01419919

Submitted on 20 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
Numerical simulations of nanotribology

Anciaux Guillaume¹, Romero Pedro¹, Molinari Jean-Francois¹*

¹ 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 : jean-francois.molinari@epfl.ch

**Abstract** — We present a 3D multiscale coupling method to address sliding contact simulations at the nanoscale. The components of the model are a molecular dynamics engine, a finite element program and a coupling scheme. We consider the normal loading of a rough (fractal) surface and discuss the influence of dislocation interactions beneath contacting asperities on the contact response. We then conduct a sliding simulation between two rough surfaces, which yields a direct prediction of the evolution of contact clusters.

**Keywords** — Nanotribology, Concurrent multiscale coupling, Molecular dynamics, Fractal surface, Sliding contact.

# 1 Introduction

Frictional forces are important to many biological functions and to almost all industrial applications. With the recent development of nanotechnologies, which implies that surface mechanisms tend to dominate bulk behavior, contact mechanics, and its by product friction, are taking an ever increasing importance. Perhaps surprisingly, considering that Amonton and Coulomb proposed a model for friction centuries ago, there is still no consensus on the origins of frictional forces. Different factors are thought to play a role: these include asperity locking, plastic activity, third body interactions, lattice dynamics, heating and melting, and possibly others. Numerical simulations may become, in the years ahead, a unique tool to understand the parameters that influence friction.

One can envision a virtual experimental facility in which the role played by surface topography and material parameters on friction and ultimately wear could be investigated. Recent FEM simulations [1] as well as atomistic studies [2] show that contact mechanics is dominated by nanoscale asperities. 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 can resort on Molecular Dynamics (MD). However, MD alone is arguably not an appropriate method as the computational cost is too high to evaluate long range elastic forces and heat flow away from the contact interface. The limitations of purely atomistic or purely continuum simulations, have motivated research in multi-scale simulations that bridge atomistic simulation and continuum modeling [3, 4]. The presented work analyses the problem of contact mechanics using the overlapping scheme presented in [5].

The objective of this paper is to describe the main ingredients of a 3D numerical engine to simulate sliding friction at the nanoscale. The numerical framework is composed of a Finite Element engine, a Molecular Dynamics package, for capturing robustly the atomic interactions at the contacting asperities, and a coupling scheme between the two physical models. We describe the results obtained for normal loading of a fractal rough surface and discuss the importance of
dislocation interactions beneath the surface. In the last section, we show, 3D sliding multiscale simulations results for a sliding test. To the best of our knowledge, the preliminary results presented in this paper are the first published 3D multiscale sliding simulations results.

## 2 Fractal surface model

All initial states were taken as in a perfect copper crystal. The considered lattice is face-centered cubic (FCC) with dimension $a = 3.615$ Å. The mass of each copper atom is $m = 63.55$ g/mol and the cutoff radius is set to $R_{cut} = 4.95$ Å. In order to reduce the molecular dynamics computational cost, we used a coupling method based on an overlapping region where the continuum and the atomistic models are both valid and are glued together by constraining the degrees of freedom in order to have coherent dynamics [6, 5, 7]. The used interatomic potential is the *Embedded Atom Method* (EAM) [8] and thus adhesive forces will be included. 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})$$

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}}$$

where $\underline{P}$ is the first Piola-Kirchhoff stress tensor, $\rho_0$ is the mass per unit 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) [9] and linear elasticity.

In order to model contact mechanics involving a surface with a roughness, we use a cubic-like copper crystal for which the top surface has been cut following a self-affine fractal [10]. The main length $L$ of the sample is $100a$, while the height, thanks to the coupling, is $233a$. The number of handled atoms is $3266563$, and the number of atoms represented by the continuum is $2400000$. These last atoms were modeled with a structured mesh constituted of $17261$ nodes and $12000$ elements. The total area of the top surface is $A_0 = 10000a^2 = 130682^2$. The bridging zone has a size of $30a$, and is positioned at a distance $17a$ from the top surface. The complete model is presented on figure 1. Concerning the fractal surface generation, we used a Voss [11, 12] generation method for which the injected Hurst parameter was $H = 0.7$. After the cutting, to insure equilibrium of the free surface, the complete domain was relaxed.

<figure><img src="image_2.png" /><figcaption>Figure 1: Coupled model designed for the study of rough surface contact.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Set of atoms that are in contact (black points are contacting atoms) with the rigid body according to a distance criterion of $R_{cut} = 4.95$, when displacement is $\delta = 4.9 a$ for (a) and $\delta = 8.9 a$ for (b). Fractional contact area $A/A_0$ is 0.18 for (a) and 0.57 for (b).</figcaption></figure>

## 2.1 Rough Surface under planar normal loading

Starting from this atomistic fractal surface, coupled to a continuum, we built a planar rigid body on top of the domain. Such a plane is to be displaced at the high velocity rate of 10 m/s against the fractal surface. A crucial quantity to be measured is the true contact area, revealing the local asperities that are in contact. To measure such a quantity, we used the same cutoff radius as in the EAM potential ($R_{cut} = 4.95$). Thus, we considered as contacting any atom of the fractal surface currently at a distance to the rigid body plane that is lower than the considered cutoff. Such a criterion allows us to identify the set of atoms, at each integration step, that are in contact. On figure 2, this particular set of atoms is shown for displacements $\delta = 4.9 a$ and $\delta = 8.9 a$. Hence this set can be projected to a 2D grid of constant spacing. The elements of such a grid that own at least one atom are considered as contact surface contributors (for an amount related to the area of the element itself). We could then obtain the fractional contact area during the simulation referred to as $A/A_0$ where $A$ is the true contacting surface. Thus, it is particularly interesting to see on figure 2(b) that the contacting surface is percolating with a contact fraction of 57% which is in good agreement with the theory of percolation [13].

Another interesting measure is the pressure distribution on contacting asperities. Knowing the stress $\sigma$ in all the substrate and especially on the fractal surface, we can compute the force applied to the surface by integration over the surface asperities:

$$W = \int_A \sigma \, dS. \qquad (1)$$

Then the effective pressure applied to the surface is $W/A_0$. Experiments [14], theoretical work [15, 16] and finite element simulations [17, 18] all expect a linear behavior of the contact ratio with respect to load for small applied loads. Figure 3 presents the fractional contact area $A/A_0$ as a function of the normalized load $W/(E'A_0)$ where $E'$ is the effective elastic modulus defined as

$$E' = \frac{E}{2(1 - v^2)} \qquad (2)$$

On figure 3, we can also observe a linear trend represented by the plain top line. The dashed curves present the results issued from finite element simulations [17]. The lower dashed curve stands for the case of elastic constitutive law while the point-dashed curve is obtained with an elasto-perfectly-plastic material. The offset between our results and the finite element ones is
due to adhesion effects since for a compressive load, adhesion forces have already brought the
substrate and the rigid plane into contact. In the presented simulation, the slope of the contact-
load curve is clearly above the elastic behavior and slightly lower than the elasto-perfectly-plastic
one. This shows a hardening of the material as with atomistic models plasticity is explicitely
represented through dislocation interactions. This will be highlighted in the next paragraph which
treats dislocation activity.

In order to evaluate the accuracy pro-
vided by the simulation, we must take into ac-
count the dislocation nucleations. Again, the
centrosymmetry criterion is used to highlight
atoms on dislocation cores. Figure 4 presents
the evolution of such atoms, with a transverse
view, for displacement of 11 Å and 17 Å.
Atoms are colored using the force field. It can
be seen that the free surface itself is naturally
plotted, as surface atoms lack symmetry. The
mesh, represented in wireframe, is drawn to
give a reference to evaluate if dislocation fronts
are close to the coupling zone. On figure 4(b)
it can be observed that two dislocation fronts
propagated through the atomic layers. How-
ever, they still did not reach the coupling zone,
proving that the accuracy of the simulation at
this point is still sufficient. More precisely, the
first dislocation front was nucleated when dis-
placement was δ = 2.2 a, while the second was
nucleated at δ = 3.3 a. This means that we are
no longer in the elastic regime at this point. In

addition, at higher displacements, some disorder is created close to the free surface which creates
dislocation loops that do not propagate, at least initially, deeper in the copper layers (arrows in
figure 4). These dislocation interactions at the contact surface create the hardening behavior that
was mentioned above in conjunction with figure 3. Accuracy is questionable at δ = 5.1 *a* when
the first dislocation is stopped at the coupling interface and at δ = 9 *a* when the dislocation loops
have grown up to the coupling zone and are then blocked (perturbating the hardening measure-
ments). Nonetheless, real contact area predictions show trends in agreement with prior models
and experiments.

<figure><img src="image_4.png" /><figcaption>Figure 3: Fractional contact area $A/A_0$ as a function of the normalized load $W/(E' A_0)$. The top line is the linear adjustment of the points issued from coupled simulation while bottom curves are resulting from finite element simulations for elastic and elasto-perfectly-plastic ($\sigma_y/E' = 0.01$) behavior from [17].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of atoms highlighted by the centrosymmetry criterion. figure (a) is for dislocation $\delta = 3 a$ and (b) for displacement $\delta = 4.9 a$.</figcaption></figure>
## 2.2 Rough on rough sliding simulations

In this section, we describe our recent attempt at simulating sliding between two 3D rough surfaces. We therefore extend the domain presented previously so that it now includes a rough rigid that is to be compressed and slided against the underneath rough and flexible substrate domain. Since we use the same fractal surface generation method we must change the random seed so that we obtain conforming asperities and non matching surfaces. More precisely, the constant section of the model is 64 *a* × 64 *a* and was consequently reduced so as to decrease computational cost. Thus we respectively have 64 *a*, 160 *a* and 20 *a* as the height of the molecular domain (820 200 atoms), of the finite element domain (7681 nodes and 41472 elements) and of the top rigid body (211 806 atoms). The final coupled domain is presented on figure 5.

We sled the top rigid body piece at a velocity of 1000 m/s in the X direction for 13.8 Å (4 lattice distances). We present here preliminary results for the temporal evolution of the contacting clusters with sliding. Figure 6 shows the contacting asperities at initial time (before sliding) in red (grey in black and white manuscript) while the updated clusters (in black) are for a sliding position of 4 lattice distances. The contact area is constant with time and is comprised of small non-connected contact clusters that have fractal properties. We observe that new (few atoms) contact clusters have arisen while some initial clusters have suffered an area reduction. From these simulations, we hope to be able to predict, in the near future, extinction rate of clusters, which is an important parameter in rate and state dependent friction models. Since the sliding velocity is high, and therefore non-physical, the computed frictional forces are not relevant and consequently will not be presented here. The study of frictional forces is the subject of ongoing work and will be published in a following paper.

<figure><img src="image_6.png" /><figcaption>Figure 5: Schematic of the rough on rough domain.</figcaption></figure>

# 3 Conclusion

In the presented work, one of our main aim was to validate the use of our coupling scheme in the context of sliding contact mechanics. We demonstrated that our 3D methodology which couples an MD engine to a FE solver enables the temporal evolution of contact clusters under normal and sliding loading. Several challenges are left to be addressed. Dislocation interactions beneath the surface influence the contact response. Yet, the coupling scheme blocks dislocation propagation at the coupling interface. This drawback in our methodology can be overcome by implementing a dislocation passing algorithm or by absorbing most of the dislocation activity at grain boundaries within a nanostructured material. In addition, sliding contact generates at our chosen velocity, a large amount of heat. We are currently implementing a thermal part into the coupling scheme, which will allow heat flow from the atomic interface to the FE domain.

<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of contacting clusters with respect to sliding.</figcaption></figure>
References

[1] S. Hyun and M.O Robbins. Elastic contact between rough surfaces: Effect of roughness at large and small wavelengths. *Tribology International*, 40(10-12):1413-1422, 2007.

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

[3] R.E. Miller, E. Tadmor, and M. Ortiz. Quasicontinuum simulation of fracture at atomic scale. *Modeling Simul: Mater: Sci: Eng*, 6:607-683, 1998.

[4] V.B. Shenoy, E. Tadmor, R.E. Miller, D. Rodney, R. Phillips, and M. Ortiz. Quasicontinuum models of interfacial structure and deformation. *Physical Review Letters*, 80:742-745, 1998.

[5] 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.

[6] 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.

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

[8] 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.

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

[10] B. N. J Persson, O. Albohr, U. Tartaglino, A. I. Volokitin, and E. Tossati. On the nature of surface roughness with applications to contact mechanics, sealing, rubber friction and adhesion. *Journal of Physics: Condensed Matter*, 17:R1 to R62, 2005.

[11] H. O. Peitgen and D. Saupe. *The science of Fractal image*. Springer-Verlag New York, 1988.

[12] R. F. Voss. *Fundamental algorithmhs in computer graphics*. Springer-Verlag Berlin, 1985.

[13] Dietrich Stauffer and Amnon Aharony. *Introduction to Percolation Theory*. CRC Press, 1994.

[14] J. Y. Kim, A. Baltazar, and S. I. Rokhlin. Ultrasonic assessment of rough surface contact between solids from elastoplastic loading-unloading hysteresis cycle. *Journal of the mechanics and physics of solids*, 52(8):1911-1934, 2004.

[15] A. W. Bush, R. D. Gibson, and T. R. Thomas. Elastic contact of a rough surface. *Wear*, 35(1):87-111, 1975.

[16] B. N. J. Persson. Elastoplastic contact between randomly rough surfaces. *Physical Review Letters*, 87(11), 2001.

[17] L. Pei, S. Hyun, J. F. Molinari, and M. O. Robbins. Finite element modeling of elasto-plastic contact between rough surfaces. *Journal of the Mechanics and Physics of Solids*, 53:2385-2409, 2005.

[18] S. Hyun, L. Pei, J. F. Molinari, and M. O. Robbins. Finite-element analysis of contact between elastic self-affine surfaces. *Physical review E*, 70(2):026117.1, 2004.