Nanotribology: A Model for Molecular Dynamics


The microscopic model goto top

Geometry: flat substrates and periodic boundary conditions parallel to the surfaces
draw08.gif

or curved substrate(s) to model an "open" contact (Persson 1989), when the z coordinate of the bottom rigid substrate is defined by

z = – 1

2
hx(dn)rsl
1–cos 2πx

Lx

1

2
hy(dn)rsl
1–cos 2πy

Ly

,
 
where Lx,y is the size of the system in the x or y direction and hx,y are the corresponding curvature parameters (the z coordinate of the top rigid substrate is defined in a similar way)

figmodel.gif

 

Why we must use Langevin equations? goto top

To have thermal equilibrium in a 3D model using Newton equations, one has to solve ~106 or more equations. Therefore, realistic simulation times are ~10τ0 or less, that is too small. Also, anyway the approach with solely Newton equations cannot include e-h damping (excitation of electron-hole pairs in the metal or semiconductor substrate).

 

Due to the driving force, energy is pumped into the system. This energy must finally be removed from the driven system. Indeed, the friction is due to energy losses. These energy losses are produced at the interface (but where within the lubricant? in the utmost substrate layers? – see below simulation results), and then the energy must go out from the interface into the substrates, where it will be absorbed and transformed into excitation of internal degrees of freedom of the substrates (phonons, e-h pairs) and finally to heat which has to be removed from the system. Thus, we cannot use Newton equations the external driving  f  will lead to T → ∞  finally and we have to use Langevin equations with an external damping (but ηext= ?, and where does it operate? – see below).

A standard approach is to use Langevin equations for only few layers far from the interface. However, in simulation there always exists a competition "large system ↔ long times". Besides, the most important is a detailed modeling of the interface itself, so there are no reasons to include too many substrate layers far from the interface. Thus, it is reasonable to use Langevin equations for all lubricant and substrate atoms.

 

External (viscous) damping in Langevin equations goto top

But what is ηext?  If only the thermally equilibrium state is of interest, the value of  ηext  is irrelevant (although the rate of approach to the equilibrium is maximal when ηext~ ω0). In 1D and 2D models, one typically uses  ηext= const, that is reasonable. But the results show that the value of  ηext  itself (i.e., is the system overdamped or underdamped) is very important. Robbins et al used Langevin damping applied only to the degrees of freedom that are perpendicular to the sliding direction. However, for the nonequilibrium state of a 3D model such an approach may lead to wrong results. The external damping must depend on  z (the distances from the substrates) and on the atomic velocity  v (the probability of excitation of phonons in the substrates must depend on the relative velocity of the lubricant atom with respect to the substrate).

We use

ηext= η1(z) [ηph(v) + ηeh]

 

where  η1(z)  is trivial,  η1(z) = 1– tanh [(zz*)/z*],  where  z* is a parameter of the model (z*=2.12 in the simulation presented below) so that for the atoms in the s-layer, where  z z*, we have  η11,  while for the atoms in the first (closest to the substrate) lubricant layer we obtain  η1~ 0.1

The phonon damping  ηph(v): when an atom adsorbed on the crystal surface, vibrates with a frequency  ω  near its equilibrium position, the oscillation decays due to creation of phonons in the substrate with the rate (for references go here)

η(ω) =  

π


2

 

mα


mS

ω2ρ(ω)

 

where the surface local density of phonon states may be described approximately by the function

ρ(ω) =

32


π

 

ω2(ωm2 ω2)3/2


ωm6

 

and we took  ωm=15  for the maximum (Debye) phonon frequency of the solid substrate.

Then we use  ηph(v) = η(ωwash),  where ωwash=2πv/a  is the "washboard" frequency, and  a=as  for the motion along the substrate and  a=z*  for the motion in the z direction.

For the damping due to creation of electron-hole pairs in the metal or semiconductor substrate, we used ηeh= 0.01ωs= 0.049

Dependence of the external damping on
the coordinate  z the frequency  ω = ωwash= 2πv/a
draw09.gif draw10.gif

Stochastic equations for velocity-dependent damping goto top

Langevin equations with a velocity-dependent damping coefficient take a special form as described below. First, recall that a general form of the stochastic equation for a measurable variable  q ≡ {ql} is the following,

dql(t) = Kl(q) dt +
m 
Glm(q) dwm(t)    
(A1)
where   dwm(t) = 0  and   dwl(tdwm(t) = δlmdt
 
This set of equations is equivalent to the Fokker-Planck equation
 
f(q, t)

t
= –

l 
 

ql
  [Kl(q) f(q, t)] + 1

2


kl 
2

qkql
æ
è


m 
Gkm(qGlm(q) f (q, t) ö
ø
 
 
for the distribution function  f(q, t).  If we put in Eq.(A1)  q1= x  and  q2= vK1(x,v) = v,  K2(x,v) = – η(x,v) vU(x)/m,  G11= G12= G21= 0, then the unknown function G22(x,v) is coupled with the random force in Langevin equation by the relationship  δF(x,v,t)/m = G22(x,v) dw2(t)/dt,  or  δF(x,v,t′) δF(x,v,t) = m2G222(x,vδ(tt′). Then the FP equation takes the form
  f

t
+ v f

x
U(x)

m
f

v
=

v
é
ë
η(x,v) æ
è
v + 1

2η(x,v)

v
G222(x,v) ö
ø
 f(x,v,t ù
û
(A2)
Because the Maxwell-Boltzmann distribution  f(x,v) ~ exp{-[(1/2)mv2 + U(x)]/kBT}  must satisfy Eq.(A2), we obtain the equation on G22(x,v),
  G222(x,v)

v
 –  mv

kBT
G222(x,v) + 2η(x,v) v = 0.
 
This equation has the solution G222(x,v) = 2ηR(x,v,T) kBT/m,  where
ηR(z,v,T) =



0 
 eε 
    ~              
η(z,v(ε))    and  
 
~
v
2
 
(ε)  =  v2 + 2kBT

m
ε
 

In the case when the external damping does not depend on the velocity, we have ηR(x) = η(x),  i.e., the standard expression. A deviation of  ηR(ω)  from  η(ω)  is essential at small frequencies/velocities for high temperatures only  (T>0.5  for the parameters used in simulations).

anju-add.gif Figure: the coefficient ηR as a function of the frequency ω = 2πv/as at different temperatures

Motion equations goto top

We use Langevin equations for all "mobile" atoms,
 mα

..

riα

= fiα(int) + 2

S=1 
fiα, S
(M1)
where α=s  or  α=l  for "substrate" or "lubricant" atoms respectively,
fiα(int) = –

riα
all

iα′ 
Vαα(riαriα )
(M2)
The last term in Eq. (M1) describes the interaction of a "mobile" s- or l-atom with the bottom (S=1) and top (S=2) substrates,
fiα, S = fiα, S(int) + fiα, S(fric) + fiα, S(ran)
(M3)
fiα, S(int) = –

riα
  NS

i=1 
Vsα(RiSriα )
(M4)
f, S(fric) = – mαη(...) (  

 .

riαRS
)   
(M5)
where η(...) is the "external" damping coefficient,
η(...) = η( zrel, vrel ),      zrel = (–1)(S–1)(ziαZS),      vrel =

 .

riα

 .

RS
(M6)

RS ≡ { XS,YS,ZS } is the center of mass coordinate of the S-th substrate (for the bottom substrate we took R1 ≡ 0).

fiα, S(ran) describes the random force acting on the i-th atom from the S-th substrate:

fiα, S(ran)(t) fiα′, S(ran)(t′) = 2ηR(...)mαkBT δiiδαα′δSSδ(tt′)
(M7)
Recall that the function ηR(zrel) in Eq. (M7) coincides with the external damping coefficient η(zrel)  if the latter does not depend on the velocity (but may depend on the coordinate). Otherwise, if the external damping depends on  vrel,  the two coefficients are coupled by the relationship
ηR(z, v, T) =



0 
 eε η(z,
~
v
(ε)),      
~
v
2
 
(ε) = v2 + 2kBT

mα
ε
(M8)
For the top rigid substrate we use Newton equation of motion,
MS
..
R2
= NS fext + FS
(M9)

where  MS=NSmS  is the mass of the rigid substrate,  fext={f, 0, fload} is the external force applied to it, and  FS = – ∑all f, S=2  according to third Newton law (conservation of the total momentum of the system).

 

Parameters of the model goto top

All atoms interact via a 6-12 Lennard-Jones pairwise potential
V(r) = Vαα é
ë
æ
è
rαα

r
ö
ø
12

 
– 2 æ
è
rαα

r
ö
ø
6

 
ù
û
 

but the parameters are different for different kinds of atoms:

for the s-s interaction we took Vss= 3 = fixed  (the maximal energy parameter, fixes the base energy parameter) and rss= 3 = fixed  (the base length parameter),

for the l-l interaction  rll= 4.14 (an "incommensurate" case),

for the s-l interaction,  Vls= 1/3 (a "weak" substrate-lubricant interaction) and  rls= 0.5(rss+ rll) = 3.57

Important: VslVll !

either Vsl << Vll → solid lubricant → minimal  fs  and  ηeff  (e.g.,  Vll=1  for the "hard" lubricant),
or Vsl >> Vll → commensurate → melting/freezing mechanism of stick-slip (e.g.,  Vll= 1/9 for the "soft" lubricant).
Masses: ml=ms=mS=1 → a typical frequency is  ωs= [Vss′′(rss)/ms]1/2 = 4.9, the typical period is τs= 2π/ωs= 1.28, and the typical lubricant frequency is ωl= 2.05.

Typically the load was  fload= –0.1

 

Units goto top

We use dimensionless unit ("natural units", n.u.). Connection with the SI system of units:

Let a real system is characterized by

  • the amplitude of interaction in the substrates Vss measured in eV,
  • the substrate lattice constant  as measured in Å, and
  • the mass of lubricant atoms ml measured in proton masses mp.

    Introduce the following coefficients:

  • νe = Vss /Vss = Vss /3,
  • νr = as /as = as /3, and
  • νm = ml /(100ml) = ml /100.

    Then for a typical system we have νe ~ νr ~ νm ~ 1,  and

    For example, the load  fload= –0.1 n.u. typically used in our simulation corresponds to the pressure  P = –fload /as2 = 1.11·10–2  n.u. = 1.78·109  N/m2. To compare, the maximum pressure above which the plastic deformation begins, is

    Velocities: a typical value when the transition from a stick-slip motion to smooth sliding is observed  experimentally, is vc ~ 1 μm/s = 10–9 n.u.

     

    Next: Simulation results goto top

    goto main Back to main page or tribology page


    Last updated on April 1, 2004 by O.Braun.  Translated from LATEX by TTH