43

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

Special Feature: CAE and Simulation
Research Report

Exact Solution and Simulation Method for Energy Dissipation via Lattice
Vibration at a Friction Interface
Seiji Kajita, Hitoshi Washizu and Toshihide Ohmori
Report received on Oct. 30, 2014

In order to dissipate lattice-vibration energy generated in a molecular dynamics (MD)calculation properly, we developed a numerical method by which an infinitely large lattice system can beconnected to an MD calculation cell. MD calculations of kinetic friction models showed that an enormous
number of solid atoms is required to simulate energy dissipation due to long-wavelength lattice vibrations.
The proposed method derived accurate results without calculating large number of the solid atoms explicitly.
This method would be useful in clarifying non-equilibrium phenomena, such as catalysis and fracture
dynamics, as well as the kinetic friction mechanism.
Atomic Scale Simulation, Friction, Lattice Vibration, Energy Dissipation,
Heat Exchange, Lattice Heat Bath

1. Introduction
Despite industrial demands related to the growing
need for energy conservation, the mechanism involvedin friction remains unclear. One difficulty in clarifying
the friction mechanism stems from the elementary
processes at an atomic scale. For example, if one
supposes a system in which two solids are in sliding
contact, the surface atoms would repeatedly generate
and break chemical bonds at the interface. If the sliding
surfaces are lubricated, the lubricant may react with thesurface and alter the friction properties. These atomic
elementary processes affect friction, but are difficult toobserve the friction interface experimentally.
Simulations are useful for accessing this hidden
friction interface, and molecular dynamics (MD)and electronic state calculations have been used toinvestigate tribology phenomena such as lubricationand tribo-chemical reactions.(1,2)
However, these atomic simulations are not free fromtechnical problems. One of the major problems is theinsufficient treatment of energy dissipation in kinetic
friction. During shear, atoms within the frictional
region absorb the kinetic energy of the sliding material,
and the lost kinetic energy is transferred to otheratomic degrees of freedom, such as lattice vibrations
and electron excitations. In the present study, we
focus on the energy dissipation via lattice vibrations.
In order to simulate a highly non-equilibrium state in© Toyota Central R&D Labs., Inc. 2014

kinetic friction, the heat generated in a calculation cellshould be removed by means of a thermostat, which
exchanges heat with the calculation cell as the rest of
the entire system. Usually, the Langevin equation and
the Nose-Hoover method have been used to express
the thermostat to maintain a specified temperature.(3)Although these methods provide excellent thermostats
for liquid and gas systems, they are not useful in asolid system because the Langevin equation is basedon a Brownian environment and the Nose-Hoover
method cannot retain the ergodic property in harmonic
oscillator systems.(4)We present a method that mimics the heat exchange
via lattice vibrations in a solid friction system. Our
strategy is simple. As shown in Fig. 1(a), in the
MD simulation cell without thermostats, the excited
lattice vibrations by shear propagate to the surface
normal direction and reflect at the boundary. In order
to dissipate the excited lattice vibrations, we derive
the exact solutions for surface atoms in an infinitely
large harmonic-coupled oscillator in Fig. 1(b). Then,as shown in Fig. 1(c), the surface atoms of theinfinite harmonic-coupled oscillator are set at a
calculation-cell boundary to interact with the MD
atoms. Since the derived surface atoms involve infinite
degrees of freedom, the excited lattice vibrations
propagate toward the infinite lattice system and neverreturn to the surface. Therefore, the proposed method
can correctly reproduce the energy dissipation within
http://www.tytlabs.com/review/

44

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

the harmonic approximation.
The remainder of the present paper is organized asfollows. Section 2 provides theoretical and numericaldescriptions for the dynamics of the surface atoms
in the infinitely thick lattice. Section 3 presents thesimulation results of the proposed method with friction
models. Section 4 presents remarks and conclusions of
the present study. Although, for simplicity, we discussa two-dimensional lattice system, this formalism canbe readily extended to a three-dimensional latticesystem.
2. Methods

which N and N' atomic layers are aligned in the surface
normal direction, z, and lateral direction, x, respectively,
as shown in Fig. 2(a). Here, h and x are used to identify
atoms. Periodic boundary conditions are applied in thex direction, and the bottom layer connects to a rigid
body. The bonds between the atoms are approximatedby two-dimensional harmonic springs, where Kx(') andKz(' ) are the spring constants in the x and z directions,
respectively, for a bond extending in the z(x) direction(see Fig. 2(b)). We assume that only the surface atomic
layer is subjected to an external force. Denoting the x
and z displacements of an atom by u, the equation of
motion for the lattice is given as follows:

2. 1 Semi-infinite Dynamic Green’s Function
MethodWe consider a two-dimensional square lattice in(a)Calculation cell
Sliding,(1)

where d is the Kronecker delta, and ∆ refers to the(c)Calculation cell
Sliding(b)

Infinite harmonic oscillators
Infinite harmonic oscillators
Fig. 1

Conceptual schematic diagram of the method for energy dissipation via lattice vibrations. (a) Calculation cell for afriction system by MD simulation without a thermostat. Periodic boundary conditions are applied along the lateralsurface direction. (b) Surface model of a coupled-harmonic oscillator containing an infinite number of atomic layers in
the thickness direction. (c) Method for describing the energy dissipation via lattice vibrations. The surface atoms of theinfinite coupled harmonic oscillators are included and are connected to the bottom atoms in the MD calculation cell.

© Toyota Central R&D Labs., Inc. 2014

http://www.tytlabs.com/review/

45

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

nearest-neighbor atoms connected to an atom labeledby h and x. The quantities m and Fh are the atomic
mass and the external force acting on the surface of
atom h, respectively.
Equation (1) can be solved analytically, and we cantake the limit of an infinite number of atomic layers N
in the solution.(5,6) The solutions for the surface atoms(x = 1) in the infinite N limit are expressed by discrete
Fourier transform as follows:
,

(2)
,

,
,

© Toyota Central R&D Labs., Inc. 2014

2. 2 Friction Decomposition with Lattice Vibration
Modes
By virtue of the analytically exact solutions, the
SI-DLGF formalism provides information about thedegree to which a lattice-vibration mode dissipates the
kinetic energy of a sliding material. Combining Eq. (2)with the thermodynamic relationship between thefriction force and the energy-loss rate, we can derivein thethe equation for the kinetic-friction force(a)Periodic boundaries

=1

=2

...

=N'

=1
=2

...

where l is the wave number index,, and. The term
will be discussed
later herein. The memory kernel function sl includes
the effect of all of the atoms in the infinite lattice.
Hereinafter, we refer to the function sl as thesemi-infinite dynamic Green’s function (SI-DLGF).
Note that, despite calculating only the surface atoms
explicitly described in the Eq. (2), the surface atoms
can move as if they are connected to the infinite lattice,thanks to the SI-DLGF.
is a general solution of Eq. (1) under
The function
the initial conditions for the uniform distortion of the, wherelattice. The explicit form isis the force associated with the uniform distortion. Thissolution is important when applying uniform forces
expresses the normal
through the lattice; namely,is used to balance
force applied to the lattice, andthe shear force in order to achieve a steady state in afriction system. Therefore, the surface atoms describedby the SI-DLGF method act as a heat exchanger but
also a role to apply normal pressure and shear stress.
The second term on the right-hand side of Eq. (2) istroublesome because the computational costs for the
convolution integral are too high. Using a conventional
integral method, the calculation over Nt time steps
requires a large memory space O(Nt2 ) in order tomaintain the force history and O(Nt ) operations.
Fortunately, this difficulty can be overcome by a fast

convolution algorithm that uses the inverse Laplacetransformation to reduce the computational cost
significantly, so that the memory requirement isO(log Nt ) and the operation requirement is O(Nt logNt ).
Details concerning the algorithm can be found in Ref. (6).
For the time integral in Eq. (2), a typical method suchas the Runge-Kutta method and varlet algorithms canbe used.

=Nz

Rigid body
x

(b)

' Kz)
'
(Kx,
(Kx, Kz)
Fig. 2

Schematic diagram of (a) infinite coupled harmonic
oscillators and (b) spring constants for the atomic
bonds.
http://www.tytlabs.com/review/

46

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

stationary state, as follows: (6,7), (3)						,
where w is the frequency. The power spectrum
can be obtained by the history of thesurface forces in the simulations. The analytical formof the real part of the SI-DLGF is given and represents
the coefficient of the energy absorption by a surface
vibration mode. Using Eq. (3), we can determine
the contribution of the lattice vibrational modes to
the kinetic friction by the analysis of a SI-DLGF
simulation. Equation (3) states that the kinetic-friction
force can be decomposed by the contributions of themode indices u, n, and l. Furthermore, the contributionsof the bulk and surface effects to the kinetic friction
are independent of each other because the friction is
the product of the surface term and the bulk term in
Eq. (3), rather than their sum.
3. Results
3. 1 Simple Demonstration: Friction Dependence
on the Number of Atomic Layers
A friction system, in which the upper and lower bodies
have the same lattice constant s and an atomicallyflat surface, as shown in Fig. 3, are used to compare
the MD and SI-DLGF results, i.e., simulations withfinite and infinite atomic layers systems are compared.
The upper layer slides at V = 10 m/s, and the bonds
between the atoms are frozen. The lower layer is a
coupled harmonic oscillator with the spring constantof diamond and the atomic mass of a carbon atom. The
number of the inner atomic layer is N. The interaction
potential between the upper and lower atoms is arepulsive soft core potential. For the calculation of thefinite N system, we introduced a velocity-proportional
damping term that mimics the energy dissipation.
More details of the simulation can be found in Ref. (6).
Figure 4 shows the friction coefficients obtained by
the MD and SI-DLGF simulations. In the case of N = 1,
the friction coefficients are approximately zero. The
friction coefficients increase with N, and the values
converge to those for infinite N, which is obtained© Toyota Central R&D Labs., Inc. 2014

by the SI-DLGF method. Over a thousand atoms areneeded to obtain the converged friction coefficients.
In order to analyze the friction dependence on thenumber of atomic layers more closely, the friction
decomposition in Eq. (3) was used. Owing to thetranslation symmetry of the friction system, thewave number l = 0 is only allowed to be excited.
Figure 5 shows the decomposition of the friction
with the frequency mode, and the lowest frequency
component n = 1 shares a large amount of the total
friction coefficient. The lowest frequency is 2pV/s
and is regulated by the sliding velocity V, which is
generally in the range of 0 to 10 m/s. In the case of
N = 1, the eigenfrequency for the single harmonic
oscillator is on the order of the sound velocity in asolid, which is typically several thousand meters per
second. Therefore, since the eigenfrequency is toohigh to resonate with the sliding velocity, the single
harmonic oscillator cannot dissipate the kinetic energyof the sliding counter object. As N increases, thelattice can be excited at longer wavelengths and lower
frequencies, and energy dissipation eventually occurs
when the lattice vibrations become so slow as toresonate at 2pV/s. These result reveal the importanceof using an extremely large number of solid atoms
when reproducing the energy dissipation via lattice
vibrations.

V

σ

N

Rigid bodyFig. 3

Schematic diagram of a friction system between
flat surfaces. The lower solid is a coupled-oscillator
surface model that consists of N inner atomic
layers. Vertical lines indicate periodic boundaries.
http://www.tytlabs.com/review/

47

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

3. 2 Single Asperity System
The potential effectiveness of the SI-DLGF was
examined using the simple model described in theprevious subsection. Here, we use a more complex but
realistic system involving friction at a single asperity.
The friction system is illustrated in Fig. 6. In the lower
solid, there is a coupled-harmonic oscillator simulated
by MD calculation, and the number of surface atoms
N' and 10 of the inner lattice layers are in alignment.
The atomic mass is the same as that of a siliconatom, and the spring constant is that of bulk silicon.
Atoms simulated by the SI-DLGF with the samebulk parameters as the MD lattice are deployed at thebottom of the MD lattice with spring connections. In

N=infinite
N=2000N=1000Friction coefficient

0.25
0.2
0.15

N=300

0.1
0.05

this manner, the MD lattice can make a conjunctionand exchange heat fluxes with the infinitely large
lattice. In the upper solid, the surface atoms align inthe x direction to form an atomically flat asperity oflength Lc. The upper layer slides at V = 10 m/s, and thebonds between the atoms are frozen. The upper andlower lattices have the same lattice constant s, andthe interaction potential is a repulsive Lennard-Jones
potential. Owing to the periodic boundary conditions,the lower solid has an infinitely large surface area,
whereas the upper solids form asperities at intervalsof N's.
Figure 7 shows the dependence of the kinetic friction
on the number of the surface atoms N'. As N' increases,
the friction coefficient first decreases for N' < 1000.
For an N' of approximately 1000, the friction coefficientincreases and finally converges to a finite value in thelarge N' limit. This strong dependence indicates that alarge number of surface atoms, and not only internal
atoms, is needed for the energy-dissipative surface
lattice vibrations.
In order to clarify the origin of the N' dependence,
Eq. (3) is modified. The surface force term ofis approximated to separate the space

N=100
N=1

0
0

50

100

150

Calculation cell

200

Normal pressure [MPa]V

Fig. 4

Dependence of kinetic friction coefficient on thenumber of atomic layers N. Symbols , , , ,and indicate plots for N=1, 100, 300, 1000, 2000
and infinite, respectively. The infinite data wereobtained by the SI-DLGF method.

Lc
σ

N'

10

Friction coefficient

0.16
0.12
0.08
0.04

Infinite harmonic oscillators
0

Fig. 5

1

2

3

n

4

5

Decomposition of kinetic friction coefficient withfrequency index n. Normal pressure is 100 MPa.

© Toyota Central R&D Labs., Inc. 2014

Fig. 6

Schematic diagram of a friction system for a single
asperity. The yellow atoms are simulated by MD
calculation, whereas the green atoms within thecalculation cell are simulated by the SI-DLGF
method.
http://www.tytlabs.com/review/

48

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

and time variables. (A detailed description is provided
in Ref. (7).) The modified equation is as follows:,(4)

where J0 and kl are a zero-order Bessel function of thefirst kind and the surface wave number, respectively,. The Bessel function termandis derived from the shape of the upper contacting
surface. The term Cu(n) is a quantity related to the
autocorrelation of the temporal fluctuation of thesurface force during shear. Figure 8 shows Cu(n) for
this friction simulation with N' = 128 and 4096. As
in the previous case of the flat contact system, thecontribution of the lowest frequency mode, n = 1, is thelargest. The other n components decay exponentially.
Since the n = 1 component is the largest, we canreasonably consider only the n = 1 term in Eq. (4).
Parameter A in the maximum restriction for |l| in Eq. (4)refers to the ratio of the sliding velocity to the speed ofsound in a solid. The actual value of A in this setup isapproximately 1/1000. Thus, the surface wave number
l modes that satisfy |l| < N'/1000 act asenergy-dissipation channels. When N' is less than 1000,
only the l = 0 mode is excited in order to dissipate theshear energy. In this case, Eq. (4) yields.
As shown in Fig. 8, Cu(n) does not change significantlywith N'. Therefore, the above equation assumes the
kinetic friction to decrease as 1/N'. This assumption
coincides well with the results shown in Fig. 7. When
N' > 1000, |l| ≥ 1 modes begin to contribute to theenergy dissipation. The number of energy dissipativel modes increases by 1 with each increment of N' by
1000, whereas the magnitude of each energy-dissipation
mode decreases as 1/N'. Therefore, the value of thefriction coefficient converges in the large N' limit.
The N' dependency can be intuitively understoodas follows. When N' is small, surface lattice vibration
modes other than l = 0 cannot be excited, because thewavelength of |l | ≥ 1 is too short and the time scale ofthe frequency is too fast for resonance to occur withthe force fluctuation due to the sliding velocity
V. As N' increases, the surface wavelength increases
and the frequency gradually decreases. Finally, the lmodes can resonate to dissipate the kinetic energy ofthe sliding material.
4. Conclusion and Remarks
A method for energy dissipation via lattice
vibrations was proposed herein. We successfully
derived analytical solutions for surface atoms in an

1
10-1

0.15

Cu(n)[ nN 2 ]Friction coefficient

0.2

0.1
0.05
0

10-2
10-3
10-4

100

N'

Fig. 7

1000

10000

Dependence of the friction coefficient on thenumber of the surface atoms N' by SI-DLGF
method. The normal force for an upper surface
		 atom is 0.0473 nN. The number of the uppersurface atoms is 20.
© Toyota Central R&D Labs., Inc. 2014

10-5

1

2

3

4

n
Fig. 8andindicate
Results of Cu(n). Symbols
the results for u = z and u = x for N' = 4096,
respectively. Symbols and indicate the results
for u = z and u = x for N' = 128 respectively.

http://www.tytlabs.com/review/

49

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

infinite harmonic-coupled oscillator and obtained the
SI-DLGF method, which can simulate the dynamics
of the surface atoms on which the existence of an
infinitely large lattice is projected implicitly. Since
the SI-DLGF atoms and MD atoms can coexist in a
calculation cell, the SI-DLGF atoms dissipate thevibrational energy generated by a non-equilibrium
phenomenon occurring in the MD atoms.
We demonstrated the SI-DLGF method in friction
systems, and the strong dependence of friction on thenumber of solid atoms was found. In order to dissipate
the kinetic energy of a sliding counter material, theenergy-dissipative lattice vibrations should have along wavelength and be slow enough to resonate
with the sliding velocity. These results have two
important implications. First, in terms of simulations,the number of solid atoms is a key parameterbecause the energy-dissipative lattice vibrations
are long coherent movements. Second, in terms ofapplications, kinetic friction can be changed by not
only surface modifications but also by modificationsof bulk characteristics, because Eq. (3) implies thatthe contributions of the bulk and surface effects to thekinetic friction are independent.
The SI-DLGF method is not a solid heat bathcompletely, because this method can only treat theinfinitely large lattice at zero Kelvin, but not at a
finite temperature. Introducing a finite temperature to
the Si-DLGF method is possible by extension of theformulation in Eq. (2) to the generalized Langevinframework; this extension will be presented elsewhere.
Another important shortcoming is that the SI-DLGF
method cannot treat the general structure of solids,
such as face-centered cubic and body-centered cubic
structures. No analytical method by which to solve
this problem has yet been established, but numericalmethods for obtaining Green’s functions for general
structures have been proposed.(8) Although there
remain problems to overcome, the presented method
may be a fruitful step for the goal of the proper solid
heat bath, that could contribute to automotive fields,such as tribology, catalysis, and fracture dynamics inatomic simulations.

© Toyota Central R&D Labs., Inc. 2014

References(1) Washizu, H. and Ohmori, T., “Molecular Dynamics
Simulations of Elastohydrodynamic Lubrication OilFilm”, Lubr. Sci., Vol. 22, No. 8 (2010), pp. 323-340.
(2) Zilibotti, G., Corni, S. and Righi, M. C., “Load-induced
Confinement Activates Diamond Lubrication by Water”,
Phys. Rev. Lett., Vol. 111 (2013), 146101.
(3) Hünenberger, P. H., “Thermostat Algorithms for
Molecular Dynamics Simulations”, Adv. Polym. Sci.,
Vol. 173 (2005), pp. 105-149.
(4) Watanabe, H., “Basis of Temperature Control in
Molecular Dynamics Simulations” (in Japanese),<http://apollon.issp.u-tokyo.ac.jp/~watanabe/pdf/
mdthermo.pdf>, (accessed 2014-9-5).
(5) Kajita, S., Washizu, H. and Ohmori, T., “Deep Bulk
Atoms in a Solid Cause Friction”, EPL, Vol. 87(2009), 66002.
(6) Kajita, S., Washizu, H. and Ohmori, T., “Approach of
Semi-infinite Dynamic Lattice Green’s Function andEnergy Dissipation Due to Phonons in Solid Friction
between Commensurate Surfaces”, Phys. Rev. B,
Vol. 82 (2010), 115424.
(7) Kajita, S., Washizu, H. and Ohmori, T., “Simulationof Solid-friction Dependence on Number of Surface
Atoms and Theoretical Approach for Infinite Numberof Atoms”, Phys. Rev. B, Vol. 86 (2012), 075453.
(8) Benassi, A., Vanossi, A., Santoro, G. E. and Tosatti, E.,“Parameter-free Dissipation in Simulated Sliding
Friction”, Phys. Rev. B, Vol. 82 (2010), 081401(R).
Figs. 4 and 5
Reprinted from Phys. Rev. B, Vol. 82 (2010), 115424,
Kajita, S., Washizu, H. and Ohmori, T., Approach of
Semi-infinite Dynamic Lattice Green’s Function andEnergy Dissipation Due to Phonons in Solid Friction
between Commensurate Surfaces, © 2010 American
Physical Society.
Figs. 7 and 8
Reprinted from Phys. Rev. B, Vol. 86 (2012), 075453,
Kajita, S., Washizu, H. and Ohmori, T., Simulation ofSolid-friction Dependence on Number of Surface Atoms
and Theoretical Approach for Infinite Number of Atoms,© 2012 American Physical Society.

http://www.tytlabs.com/review/

50

R&D Review of Toyota CRDL, Vol.45 No.4 (2014) 43-50

Kajita Seiji
Research Fields:
		- Nano Tribology
		 - Atomic Scale Simulations
Academic Degree: Dr.Sci.
Academic Societies:
		 - Japanese Society of Tribologists
		 - The Physical Society of Japan

Hitoshi Washizu
Research Fields:
		 - Molecular Simulation of Interfaces
		- Tribology
Academic Degree: Ph.D.
Academic Societies:
		 - Japanese Society of Tribologists
		 - The Molecular Simulation Society of Japan		 - American Chemical SocietyAward:
		 - Young Scientist Award of Japanese Society of		 Tribologists, 2004

Toshihide Ohmori
Research Fields:
		 - Frontier Research Promotion
		- Tribology
Academic Degree: Dr.Eng.
Academic Societies:
		 - Japanese Society of Tribologists
		 - Society of Automotive Engineers of Japan
Awards:		 - Best Paper Award, Japanese Society of Lubrication,
1991
		 - R&D 100 Award, 1991
		 - Award of Tokai Chemical Industry Association, 1996
		 - Chunichi Industrial Technology Award, 2007
		 - Best Paper Award, Society of Automotive Engineersof Japan, 2012
		 - Technology Award of Japanese Society of Tribologists,
		 2012© Toyota Central R&D Labs., Inc. 2014

http://www.tytlabs.com/review/