Computational Materials Science 81 (2014) 98–103
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Multiscale model of elastic nanocontacts Guillaume Michal ⇑, Cheng Lu, A. Kiet Tieu School of Mechanical, Materials and Mechatronic Engineering, University of Wollongong, 2500 NSW, Australia
a r t i c l e
i n f o
Article history: Received 13 February 2013 Received in revised form 4 June 2013 Accepted 25 June 2013 Available online 19 July 2013 Keywords: Nanoscale Adhesion Friction Elasticity Multiscale modelling
a b s t r a c t A multiscale scheme for adhesive elastic contact is proposed to evaluate the interaction mechanisms between two solids. The constitutive law, internal to each body, is described by classical continuum mechanics. The model introduces the atomic structure in the region of contact to evaluate the contact conditions at the coarse scale. Simulations of normal loading and sliding were used to compare quantities such as contact radius, stresses and frictional forces with that from Molecular Dynamics. The radius of contact obtained from the non-adhesive cases show that the model falls in between continuum contact mechanics and molecular dynamics. Stresses resulting from a normal loading are within a few percents of the discrete model for regions further than 15r from the contact interface. The frictional characteristics exhibit essential stick–slip phenomenon and compare favourably with the Molecular Dynamics results. The integration of the non-linearities of the material at the coarse scale and the account of the out of phase motions of atoms within the coarse scale mesh are necessary to improve the solution in regions subject to more severe stresses. Ó 2013 Elsevier B.V. All rights reserved.
1. Introduction and review Conditions such as experienced in Atomic Force Microscopy’s contact mode induce mainly adhesive elastic mechanisms of interaction. The extensive use of rigid tips in most MD simulations of contact AFM shows that limiting the total number of degrees of freedom of the system, even in an important part of the contact, can lead to acceptable results. In this context, is a local but still fully atomistic model required to sensibly represent adhesive elastic contacts? The notion of adhesive contact was ﬁrst synthesised by Bradley [1,2] based on the earlier work of Rayleigh [3,4]. Bradley’s integration showed that contact mechanisms accounting for adhesion are non-local and extend to scales much larger than the atomic spacing. At the same time, the adhesion force can be considered localised to a relatively small part of the asperity. JKR  and DMT  models have been demonstrated to apply at the nanoscale, depending on the properties of the surfaces. However they do not consider the inﬂuence of the atomistic structure and its impact on the spatial repartition of adhesion. Luan and Robbins [7,8], as well as Mo et al. , showed that the use of discrete modelling is more suited. Models such as Prandl–Tomlinson  or Frenkel– Kontorova  are useful to get an insight on the excitation, dispersion and stick–slip mechanisms at the atomistic level but they are not well suited for extensions to complex geometries, complex interaction potentials or lubrication problems. ⇑ Corresponding author. Tel.: +61 620414576157. E-mail address: [email protected]
(G. Michal). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.06.053
Continuum models are limited in their representation of friction phenomena as they do not derive the contact forces ab initio. Instead, they prescribe a relationship between normal and tangent forces such as Amontons’ law. No such deﬁnition is required with an atomistic description. Two atomistically deﬁned surfaces brought together at close range would exhibit at least a non-zero static lateral force in most conditions. Because they allow the integration of friction without the deﬁnition of constitutive laws, they are more ﬂexible and depend solely on the crystal structure and associated pair potentials used in the model. The generality in handling various materials and the possibility of using lubricants with detailed chemical chains and bonds is advantageous. One basis of computational efﬁciency is the acknowledgment that a high degree of detail is not necessary everywhere for the results to be meaningful. Mesh based methods such as Finite Element Analysis make use of varying element sizes to reduce computational costs while providing a ﬁne descriptions in localised regions. Multiscale methods extend this principle further through the connection of theories valid at distinct length and time scales. Most multiscale methods are applied to problems involving locally large stresses and strains, dislocations and crack mechanisms. The mechanical conditions in these problems are considerably more severe than those encountered in elastic contacts. Even though a breakdown of continuum mechanics occurs at the nanoscale, continuum theory still holds relatively well in some cases , illustrating that the limits of applicability of theories in regard to length scales have some ﬂexibilities. Multiscale works on tribological and nano-contact problems have been carried by different groups such as in [12,13]. Luan
G. Michal et al. / Computational Materials Science 81 (2014) 98–103
and Robbins [12,14] based their work on methods that impose displacement boundary conditions at the interface MD/FEM in a way similar to the FEAt method  and CAAD . Contrary to methods which downscale the mesh to atomic separation, the method from Luan et al.  leaves the mesh relatively coarse. Despite not having a mathematical framework as sound as other multiscale methods, especially in the treatment of the energy wave propagation, the method performs rather well for contact problems. In this paper a multiscale model of purely elastic contact is presented to limit the atomistic description to a minimum. The foundation of the description is based on the decomposition of the coarse and ﬁne scales established by Wagner et al. in the Bridging Scale Method (BSM) . Here, continuum mechanics handles the constitutive laws internal to the material, while the discrete modelling is used to concurrently describe the external forces emerging from the contact interactions between each bodies. The presented approach is rather simplistic in regard to available multiscale methods. The important aspect is that the use of atomic interactions is not meant to describe the material itself. It is intended to provide an approximation of the adhesive contact laws and friction laws based on the discrete interactions and the coarse scale position of the atoms present in the contact zone. The model’s description of the mechanics at the interface between bodies falls in between continuum adhesive models and purely discrete models.
At the ﬁrst order the equation of motion for the particle a in body 1 is given by (5), where the internal and external forces are themselves written in terms of coarse and ﬁne scales as shown in Eqs. (6)–(9).
€ a ¼ f a þ fa0 þ ga þ g 0a ma u f a ¼
ua ¼ ua þ u0a
Consider a contact problem involving two interacting bodies noted 1 and 2, composed of n1 and n2 atoms respectively. The equation of motion of an atom a in body 1 is given by (2). The forces acting on the atom a are divided into internal forces f, due to particles b1 in body 1, and external forces g, due to particles b2 in body 2. f and g derive from any pair potential.
€a ¼ ma u
X X f ua ; ub1 þ g ua ; ub2 b1 –a
For small deformations of the bodies and at low temperatures it is assumed that individual atoms oscillate slightly around their ideal equilibrium. On average the internal forces are null when no external interactions exist between the bodies. The internal forces on particle a are written in the form of a Taylor expansion in (3) with ua;b1 ¼ ua ub1 . The separation between the bodies can be large, and the assumption of small displacements between particles of distinct bodies does not hold in general. However, a u b2 þ u0a u0b , and the ﬁne scales u0 are ideally small ua;b2 ¼ u 2 quantities. A Taylor expansion of the pair interaction at 0 0 a;b2 ¼ u a u b2 relative to u0n u a;b2 ¼ ua ub2 is given in (4).
X @fa;b 1 @ n fa;b1 n 1 ua;b1 þ þ u n! @una;b1 a;b1 @ua;b1 b –a
h i X @fa;b f 1 u0a;b1 ¼ ka u01 @ua;b1 b –a
X a ; u b2 g u
g 0a ¼
X @g a;b
h i g g u0a;b2 ¼ ka;1 u0a þ ka;2 u02
n X @g 1 @ g a;b2 0n a;b2 þ a;b2 u0a;b þ þ g u u 2 a;b2 na;b a;b2 n! @ u @u b 2 2
€1 u €2 u
2h i 3 f k ½0 u 1 g 6 1 7 1 h i5 ¼4 þ f 2 2 u g ½0 k2 h i 2h i h i 3 f g g 0 k1 þ k1;1 k1;2 6 7 u1 þ4 h iT h i h i5 0 g f g u2 k2 þ k2;2 k1;2 ð10Þ
h i h i f f k1 and k2 are n1 n1 and n2 n2 internal stiffness matrices at atomic sites due to the small displacements taking place within 1 and g 2 are external force vectors, bodies 1 and 2 respectively. g of length n1 and n2 respectively, taken relative to the coarse displacement of individual atoms. The remaining matrix is the assembly of internal and external stiffnesses for the ﬁne scale motion of h i f the particles. The sub-matrices ki are assumed to be constant h i g for small deformations. However ki;j depends on the effective separation. A FEM mesh is superimposed over each atomistic body. The nodal displacements are represented by a vector d. In the BSM, the projection of the total scale onto the coarse scale is deﬁned as the sum of squared difference for all atoms weighted by their atomic masses . For a body ‘i’, the quadratic form in (11) is obtained by considering the diagonal matrix [Ma,i] and by minimising the functional through the nodal degrees of freedom. The solution to the minimisation is given in (12) where [Mi] is the coarse scale mass matrix obtained from (13).
T min ðui ½Ni di Þ ½M a;i ðui ½Ni di Þ
di ¼ ½M i 1 ½Ni T ½M a;i ui
½Mi ¼ ½N i T ½M a;i ½Ni
Using (12) in (10) leads to the nodal equation of motion for i ¼ ½N i di by conbody ‘i’ in interaction with a body ‘j’ (14). Since u struction, the system is written in terms of nodal displacements in (15).
€ ¼ ½N T ½M i d i i ð4Þ
The approximation leads to the matricial system in (10) for a two body interaction.
h i X @fa;b 1 1 a;b1 ¼ kfa u u @u a;b1 b –a
2. Presentation of the method The Bridging Scale Method (BSM) formulation, as presented by Tang et al. in , is the foundation of the proposed method. In the BSM, the individual displacement ua = u(xa) of atom a, called the a and a ﬁne scale total scale, is decomposed into a coarse scale u ua 0 as shown in (1). The BSM systematically separates the scales through a projection operator. It considers the ﬁne scale only in critical regions to limit computational requirements. The total scale is supposed to be the solution from Molecular Dynamics while the coarse scale is the projection of the total scale on the FEM space of displacements.
" #! h i h h i h i h i i u0 i f f g g i þ g i þ ki þ ki;i ki u ki;j u0j
0 0 € ¼ ½K d þ ½N T g ½Mi d i i i i i þ K i u
G. Michal et al. / Computational Materials Science 81 (2014) 98–103
h i f K i ¼ ½Ni T ki ½Ni
hh i h i h ii 0 f g g K i ¼ ½Ni T ki þ ki;i ki;j
In classical continuum adhesive contact mechanics the effects of internal and external ﬁne scales are discarded altogether, i.e. 0 K i is considered null. In these models the coarse scale interaction , sum of all external atomic interactions, is approximated by an g integral. The internal forces are based on continuum mechanics and K i is approximated by a linear elastic formulation commonly used in FEM. Higher order schemes can be used to include non-linearities of the material. For instance, Luan and Robbins used a quadratic form to describe the coarse internal forces in . In multiscale methods, the ﬁne scale in the internal forces is neglected in the regions away from strong perturbations. Only the coarse scale’s constitutive laws remain in these regions. Looking at the form presented in (17), this simpliﬁcation means that part h i f of ki is set to zero in K 0i . However the ﬁne scale is kept in parts of the domain requiring the use of a discrete description. In the present model the ﬁne scale displacement is discarded entirely as in continuum methods. Only the coarse scale of the external forces is used but it is expressed in terms of the sum of the individual atomic interactions taken at the coarse scale separation. In FEM schemes, contact is often considered with algorithms of interactions along the boundary elements. It implies that the surface topography is always as smooth as the vertices in contact. In the proposed formulation it is the coarse scale motion of the atomic scale topography and its underlying atomic structure that rules the interaction. Since the ﬁne scale is discarded, the equation of motion from MD is not used in practice. The applicability of the model is limited to problems with small deformations where the energy of adhesion between the bodies is lower than the cohesive energy within the bodies themselves. The individual dynamics of atoms must not be a major actor of the mechanics in action, problems at high temperatures cannot be accurately modelled. The proposed approach does not have more restrictions than continuum adhesive contact models. In these models, the Young modulus of the bodies in contact is always fully decoupled from the energy of interaction. Their constitutive laws for internal and external forces are continuum-based. The present model substitutes the integration of the external forces by a discrete sum. The atomic structure plays a role in (a) the force received by an atom due to the presence of the second body and (b) in the force effectively received by an FEM node through the extrapolation of the external atomic forces enclosed in the elements. (Eq. (18)) gives the ﬁnal system of equation of the model.
€ ¼ ½Kd þ ½NT g ½Md
[M] is the assembled consistent mass matrix of all the objects in is interactions. ½K is the assembled coarse scale stiffness matrix. g the external force vector at atomic sites between all bodies. It is easily obtained by using MD schemes that make use of cut-off distances and neighbour lists to improve the efﬁciency. 3. Simulation of elastic contact In this section a two dimensional model of interaction between a tip and a substrate is presented. It follows the work of Luan and Robbins in . X and Z are the in-plane horizontal and vertical axis respectively. The tip and the substrate structures are based on a triangular lattice with an atomic spacing ra/r = 21/6. The tip is a rigid curved plane with a radius of curvature R/
Fig. 1. (a) Atoms of the tip and the substrate for R/r = 100 and position of cuts A, B and C in the substrate. (b) Mesh of the substrate for multiscale simulation.
r = 33,66,100,133 and 166 and a thickness of ﬁve atomic layers (27/6r). The substrate is 400r in width and 200r in depth, see Fig. 1. The substrate’s lower face is ﬁxed in position along X and Z. Its lateral faces are left unconstrained and no periodic boundary conditions are used. Both a full MD and a multiscale model were developed. The substrate was made up of 73083 atoms in the MD case and 6060 atoms with 4369 elements based on 2318 nodes in the multiscale model. A fundamental interest in the original BSM is the possibility of using a reasonably coarser mesh than the one presented in this work. However the coarsening was limited by geometry and resolution considerations. Firstly, the small radius of contact limits the coarsening level to provide a good level of contact conformity between the tip and the substrate. A larger element size reduces the conformity by enforcing a disk/plan tangency condition that rises the inter-atomic forces locally. Secondly, a coarser mesh reduces the resolution of large stress gradients in the region of the contact. Following the work in  the mesh size was set to 4r in the contact region with a gradual coarsening away from the surface of the substrate. In the MD simulation, the internal forces of the substrate are modelled using the LJ12-6 potential with i ¼ 0:15; r ¼ 3Aand a cut off distance Rci/r = 1.5. Subscript ‘i’ denotes a quantity related to the internal forces. In the multiscale model, the substrate is meshed by linear triangular elements (Fig. 1(b)) and internal forces are represented by a classical linear stiffness matrix using a Poisson ratio m = 1/3 and a Young modulus E = 67.04.r2. The contact modulus in this conﬁguration is E⁄ = 75.42r2. To remain within was set such that the elastic domain the maximum loading F max z =RE < 0:025 . Multiscale and MD simulations both use a F max z LJ12-6 potential with e = i/2 and r ¼ 3Ato model external interactions between the atoms of the two bodies. The temperature was set to T = 0.0001/kb and the Langevin thermostat used a damping pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ rate of 0.1s1, where s ¼ mr2 = and m is the atomic mass. The time step was set to DtMD = 0.005s for both the multiscale model and MD. Loading simulations were carried out ﬁrst. The initial contact between the bodies was displacement driven. The tip was displaced by a small amount until the interaction force became repul-
G. Michal et al. / Computational Materials Science 81 (2014) 98–103
sive. Once the contact was initiated, a series of step loadings with DFz/E⁄ = 0.125 was applied until the maximum loading force was reached. For each step loading the system was allowed to relax for 350s to establish the equilibrium. Shearing simulation were subsequently carried out. Allowing for a relaxation period after the initial loading, the ﬁnal load was kept constant while the tip was displaced laterally at a constant velocity of 0.01r/s during 1500s. Forces and stresses were recorded at an interval Dt = s, equivalent to recording the state of the system every 200 steps for a total of 350 recordings per step loading. In the following, the total tangential and normal interaction forces between the bodies are denoted by Ft and Fn respectively. The normal stresses along X and Y are denoted by rxx and ryy respectively. The shear stress is denoted by rxy. 4. Discussion The relationship between loading and contact radius for nonadhesive elastic contacts at the continuum level can be obtained from Hertz’s theory and follows (19) .
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 4F z R=pE
The contact radius of the multiscale model and MD was evaluated by averaging the data over the relaxation period at each step loading. The variation of the contact radius was measured in the adhesive and non-adhesive conditions. Fig. 2(a) shows the results for the non-adhesive cases. Continuous lines are the Hertz’s solution (only for the non-adhesive case), circles are the solution from MD and squares indicate the solutions from the multiscale model. At low loading, the multiscale model follows closely Hertz’s solution as both are based on linear constitutive laws of the material. MD consistently predicts a larger contact radius. As explained by Luan et al. an increase of loading leads to a stiffening of the LJ potential which limits the expansion of the contact radius in MD. This is particularly obvious for R/r = 100 and R/r = 133. Despite having the constitutive laws of the materials
Fig. 3. Stresses for a loading Fz = 3.771 nN and R/r = 100 with adhesion for (a) cut A, (b) cut B and (c) cut C.
Fig. 2. Contact radius with adhesion. Circles and squares represent MD and MU respectively.
based on linear elasticity, the multiscale approach leads to a similar behaviour, albeit the effect is less pronounced than in the MD case. This is due to the non-linearity of the external interactions. Both the non-linearity of the material constitutive laws and the non-linearities of the external interaction play a role in the overall solution. Both multiscale and MD models lead to a stepped curve due to the discreteness of the topography at nanoscale. The smoothness of the Hertz’s solution is due to the assumption of smooth surface inherent to the theory. Fig. 2(b) shows, for the adhesive cases, the same trend due to the differences in the description of the material’s non-linearities. These results emphasise that the presented multiscale scheme sits between Hertz’s theory and a purely atomistic scheme. Normal and shear stresses in the substrate are shown in Fig. 3(a)–(c) for the adhesive case. The stresses are taken at a depth of 30r (cut A), 100r (cut B) and along the depth at the centre line (cut C) for a loading Fz 3.771 nN. The cuts are shown in Fig. 1(a). The equivalence between the loading and the quantity Fz/RE⁄ is given in Table 1. Following the convention used by Luan et al. a minus sign is used to represent the stresses, such that a positive value indicates a compressive stress. All the results are obtained by averaging the 350 recordings of the relaxation period.
G. Michal et al. / Computational Materials Science 81 (2014) 98–103
Table 1 Equivalence between loading Fz 3.771 nN and ratio Fz/RE⁄ and percentage of maximum allowed loading to remain in the purely elastic domain. R/r
% of F max =RE ¼ 0:025 z
33 66 100 133 166
0.01875 0.009375 0.00625 0.0046875 0.00375
75 37.5 25 19.5 15
Table 2 Equivalence between Fz/RE⁄ and the effective loading used for the simulation of friction. R/r
33 66 100 133 166
1.0070 2.0139 3.0209 4.0279 5.0348
2.0139 4.0279 6.0418 8.0557 10.0697
3.0209 6.0418 9.0627 12.0836 15.1045
4.0279 8.0557 12.0836 16.1115 20.1394
5.0348 10.0697 15.1045 20.1394 25.1742
(nN) (nN) (nN) (nN) (nN)
In cut A, the root mean square deviation (rmsd) of each stress curve, relative to their respective absolute maximum stress, was less than 7% in all loading cases tested. Results tend to degrade at the deeper cut B. While the normal stress ryy and the shear stress rxy were a reasonable approximation of MD, the rmsd of the normal rxx was up to more than 20% at this depth. Results for cut C in Fig. 3(c) reveal that the multiscale model can represent the stress levels at the centreline reasonably well for depths larger than approximately 15r. While the rmsd is consistently below 10% along cut C, the multiscale model tends to differ more from MD as higher stresses are considered near the contact region. The inﬂuence of the tip radius on the accuracy of the model at constant loading was not signiﬁcant in the range of radius studied. Considering its simplicity, the multiscale scheme compares favourably with MD to represent stresses but discrepancies exist. The deviations are mostly the result of the lack of non-linear constitutive laws to represent the continuum scale. The frictional forces are now discussed. Loadings of Fz/ RE⁄ = 0.05, 0.01, 0.015, 0.020 and 0.025 were applied for each tip radius before the shearing motion was applied. Table 2 shows the equivalence between these ratios and the loading. Fig. 4 illustrates the lateral force obtained for the adhesive case R/r = 66 and Fz/RE⁄ = 0.015 over 5 periods. MD and multiscale models exhibit a stick–slip behaviour with a periodicity equal to the atom spacing but there is a discrepancy in the amplitude. The appearance of stick–slips in the multiscale model is a result of the existence of the atomic structure in the model. A purely continuum model would not intrinsically exhibit such features. The variations in the frictional forces can be linked to the properties of the model. Firstly, the stiffness of the material in MD is reduced if the material is subject to traction, due to the non-linearity of the potential . As the tip moves in a stick–slip fashion, the region of contact in the substrate alternates locally between tensile and compressive stresses. The non-linear properties of the atomic interactions in MD lead to a reduction of stiffness of the material in the tensile zones and a stiffening in the compression zones. Secondly, since the multiscale model considers only the atomic coarse scale, atoms inside a given element move in phase. The variation of the magnitude of the friction force is due to this in-phase motion. In MD, atoms follow a more random distribution that averages to a lower energy barrier . Reducing the element’s size in the contact zone from 4r to r showed indeed a reduction of the amplitude of the lateral force, but only to a small proportion. In the MD mod-
Fig. 4. Friction force from MD and MU for a tip radius R/r = 66, cut-off Rc/r = 6 and a loading Fz/RE⁄ = 0.015.
Fig. 5. Static (continuous line) and kinetic (dot line) friction forces for several tip radius and loadings with adhesion.
el, because the material becomes stiffer in compression, the break of lateral adhesion occurs at an earlier stage. Static and kinetic friction forces were extracted from the lateral force curve of the multiscale model. The static friction was taken as the average of 15 peak values. Kinetic friction was taken as the average over 14 complete periods. Fig. 5(a) shows the static and kinetic friction forces. The results compare qualitatively well with the results of Luan et al. . A sublinear behaviour of the static friction force is observed as the loading increases while kinetic friction increases exponentially at low loading to reach a steady rate at larger loadings. The section on friction in the paper of Luan et al. involved rough surfaces, a problem that differs noticeably from the single asperity modelled in this work, precluding quantitative comparisons. Most multiscale models that take into account the ﬁne scale motion of atoms demonstrate the production of phonon during stick–slip frictional motion. An important focus of such simulations is to adequately handle the numerical reﬂection at the interfaces between the coarse scale and the small scale regions. Because the presented model discards the ﬁne scale motion entirely, no mechanisms were necessary to limit the reﬂections. Nonetheless, coarse scale stress waves were produced in the contact region at the instant of contact and during subsequent stick–slips in our simulations. These waves propagated away from the contact region in a circular fashion before being reﬂected along the free surfaces of the substrate. The coarsening of the mesh away from the contact region also led to a progressive and partial reﬂection of the waves within the bulk of the material. 5. Conclusion While retaining the simplicity of the continuum model, the proposed multiscale model introduces the description of atomic structure in the region of contact to prescribe the contact conditions
G. Michal et al. / Computational Materials Science 81 (2014) 98–103
without specifying a friction law explicitly. The calculated radius of contact obtained from the non-adhesive cases show that the model falls in between continuum contact mechanics and molecular dynamics. The scheme reproduces qualitatively well results from MD but quantitative differences exist. Stresses resulting from a normal loading are within a few percent of the discrete model for regions further than 15r from the contact interface. A description of the non-linearities of the material at the coarse scale will produce solutions closer to MD for regions subject to more severe stresses. Results for friction exhibit essential stick–slip phenomenon and qualitatively agree with results obtained by Luan et al. . The in-phase motions of atoms inside the elements of the continuum mesh leads to a larger energy barrier to initiate slipping. The ﬁne scale motion of atoms is a key factor in the quantitative description of the friction force, even at low temperature. The introduction of scheme able to represent the ﬁne scale without having to compute it explicitly would be beneﬁcial. A Time History Kernel  for surface atoms or the use of statistical distribution of the ﬁne scale, based on parameters such as local temperature and pressure can provide answers at the cost of a more complex implementation. The velocity interfacial approach proposed by Tang  offers an alternative to the THK with a much reduced complexity.
References                    
R.S. Bradley, Philos. Mag. 11 (1931) 846–849. R.S. Bradley, Philos. Mag. 13 (1932) 853–862. Lord Rayleigh, Philos. Mag. 30 (1890) 285–298. Lord Rayleigh, Philos. Mag. 30 (1890) 456–475. K.L. Johnson, K. Kendall, A.D. Roberts, Proc. Roy. Soc. Lond. A 324 (1971) 301– 313. B.V. Derjaguin, V.M. Muller, Y.P. Toporov, J. Colloid Int. Sci. 53 (1975) 314–325. B. Luan, M.O. Robbins, Nature 435 (2005) 929–932. B. Luan, M.O. Robbins, Phys. Rev. E 74 (2006) 026111. Y. Mo, K.T. Turner, I. Szlufarska, Nature 457 (2009) 07748. G.A. Tomlinson, Philos. Mag. 77 (1929) 905–939. Y.I. Frenkel, T. Kontorova, Zh. Eksp. Teor. Fiz. 8 (1938) 1340. B.Q. Luan, S. Hyun, J.F. Molinari, N. Bernstein, M.O. Robbins, Phys. Rev. E 74 (2006) 046710. C. Yang, U. Tartaglino, B.N.J. Persson, Eur. Phys. J. E 19 (2006) 47–58. B. Luan, M.O. Robbins, Tribol. Lett. 36 (2009) 1–16. P. Kohlhoff, S. amdGumbsch, H.F. Fischmeister, Philos. Mag. A 64 (1991) 851– 878. L.E. Shilkrot, R.E. Miller, W.A. Curtin, J. Mech. Phys. Solids 52 (2004) 755–787. G.J. Wagner, W.K. Liu, J. Comput. Phys. 190 (2003) 249–274. S. Tang, T.Y. Hou, W.K. Liu, Int. J. Numer. Meth. Eng. 65 (2006) 1688–1713. W.K. Liu, H.S. Park, D. Qian, E.G. Karpov, H. Kadowaki, G.J. Wagner, Comput. Methods Appl. Eng. 195 (2006) 1407–1421. S. Tang, J. Comput. Phys. 227 (2008) 4038–4062.