Thermal Science & Engineering Vol.15 No. 3 (2007)

Molecular Dynamics Study on Wave Equation of Liquid
Tatiana ZOLOTOUKHINA and Toshihiro IWAKI

Abstract
The equation of wave propagating in fluid is described as a differential equation of the velocity or small element
density from the Navier-Stokes (NS) equation. On the other hand, the molecular dynamics (MD) equation expresses the
motion of a particle constituting the fluid and shows that the particle is always in motion regardless of the existence of
the wave. In this study, we discuss in what way meanings of the velocity and density in the wave equation can be
adopted in the MD system. What are the differentials with respect to time and space in the NS equation for the system
of MD particles? Ordinarily, the physical quantities in the NS equation are obtained as the ensemble and time averages
over the MD system. We investigate the number of particles and duration of time that are sufficient for the averages andsimultaneously confirm whether the averaged values satisfy the differential equation. The two-dimensional MD method
is used for the qualitative understanding. The fluid is assumed to consist of particles connected by the Lennard-Jones
potential. The satisfaction of the differential equation by the MD averaged values is shown by the propagation velocity
in the wave equation. The propagation velocity can be also obtained in another way i.e. by the µobservation¶ of the wave
fronts motion of velocity, density or total energy wave in the fluid. The propagation velocity that has resulted from thewave equation is strongly affected by the ensemble and time averages. On the other hand, when it is obtained from thewave front, it is independent of the ensemble and time averages. We can have the former propagation velocity close to
the latter one when a long time average is used for a large ensemble and a short time average for a small ensemble, i.e.
the product of both averages can be considered changed at same rate as if by a scaling coefficient.
Key Words:

1

Molecular dynamics, Wave equation, Wave front, Liquid

Introduction

Theoretical studies of thermo-fluid phenomena byusing the continuum mechanics have been the subject of
the ongoing research for several decades because of theextended importance in wide range of practical engineering applications. The continuum characteristics ofthe thermo-fluid phenomena is obtained by solving governing equations describing conservations of mass, momentum and energy in an extremely small cubic element
or unit cell (called thereafter the element) under certain
boundary and/or initial conditions. The extremely small
cubic element is considered to be homogeneous and continuous there. It is assumed in the continuum mechanics
that, in the cubic element, the temporal and spatial deviations of the physical quantities such as density orpressure from their average values are very small. This
assumption enables the deviations to be linear. Then we
can ignore the terms higher than second order one of the
Taylor¶s series expansion and represent the deviation by
the first order derivative. Therefore, the governing equations are expressed by differential equations. This dif-̐

ferential equation approach represents the approach tothe physical phenomena that are considered continuous.
The study using the continuum mechanics is the macroscopic one.
However, from the physical viewpoint, the cubic
element consists of many atoms and/or molecules in motion at high velocity. The element is not continuous but
discrete microscopically. Atoms and molecules are assumed to be called particles hereafter. The motion of aparticle is expressed by the MD equation and then the
cubic element is nothing but a limited volume of space.
Particles don¶t (continuously) always stay in a certain
cubic space and may move in or out from the cubic
space to the neighboring elements due to their individual
motion. The cubic element is not homogeneous in relation to the number and type of particles that are moving
inside it.
Due to the duality of approaches, here we have avery important and fundamental problem of connection
between the continuum mechanics and the particle motion outlined by the following questions. How thephysical quantities in the governing equations of theReceived: Feb. 2, 2007, Editor: Taku OHARA
Department of Mechanical and Intellectual System Engineering, Faculty of Engineering, Toyama University (3190 Gofuku,
Toyama 930-8555, JAPAN)-1-- 111

C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No. 3 (2007)

2

In the continuum mechanics, the wave equation for
disturbed velocity (or density waves), u, is expressed by

c2

w 2u
wx 2

i=3

' yi

i=2

x

' xi

'xs

Fig. 1 Particles near the left end of the MD systemand definition of subregion.
controls oftemperatureand pressure

movement of wall end ofon the left side
simulation

t=±5000 t=±4500

t=0 t=500

Fig. 2 Time development of calculation procedure.
time j=1, (t1)
ta

'ts

time j=2, (t2)tb
tn

Fig. 3ttn

Average with respect to time and definition oftime interval.

byc2§ wP ·¨¸© wU ¹

(2)

c2

J

P
U

(3)

or

where P, ȡ and Ȗ are pressure, density and specific heatratio, respectively. There is no wave motion without disturbance, and a region of disturbance propagation shouldhave a coherent change in respective values of the wave.
On the other hand, the particle is always in motion regardless of the wave. The motion of a particle constituting the liquid is expressed by the MD equation

m

w2 x
wt 2

F

(4)

where m and x are the mass and position of the particle,
respectively, and F is the force acting on the particle.

Wave Equation and Particle Motion

w 2u
wt 2

subregion i=1ycontinuum mechanics are derived from characteristics ofthe particle motion? What are definitions of differentials
x and t of the governing equations in the MD system?
Namely, how can the MD method connect particle motion to the continuum mechanics? The answer to the
problem will make clear how small a unit size and how
short a time interval can be applied in the continuum
mechanics in the limit of small scale.
For the problem stated, the governing equations of
the continuum mechanics were derived from the MD
equations [1-4] by Kotake. He concluded that the local
uniformities in the temporal and spatial scales are necessarily demanded in the MD system for the connectionof particle motion to the continuum mechanics. If the
MD system has no such local uniformities, it cannot explain even non-equilibrium and/or unsteady state macroscopic mechanical phenomena. These conclusions lead
us to discuss what the local uniformities are and then
how they are obtained in the real MD simulation.
In order to determine the local uniformity, the planeelastic wave, that is a propagation of disturbed displacement of an extremely small cubic element in anelastic solid continuum, was investigated as one of thetypical non-equilibrium and unsteady state macroscopic
mechanical phenomena [5]. It was found that, if an appropriate average value of particle displacement withrespect to space and time in the MD system was definedas the displacement of a point in the elastic solid continuum, it behaved like in the wave equation of the continuum mechanics.
There are many types of disturbances propagating ina medium besides the displacement wave and are different mediums propagating various types of disturbances.
Therefore, it would be important to compare whether
some or all features and behavior of chosen disturbances
are same as the displacement wave. In the present study,
propagations of disturbed velocity and density waves ina liquid are investigated. They are compared with theplane elastic wave. Although, the wave propagation of
the energy disturbance is not described by the same typeof the wave equation as the motion of the disturbance invelocity or density waves, it is also examined.
We use Lennard-Jones particles and calculate theirmotions by means of the two-dimensional MD method
for understanding of the phenomenon qualitatively.

3(1)

under the conditions of stationary fluid and an isentropicchange, where c is a propagation velocity and is given-2-- 112

Model and Calculation Method
The potential of the Lennard-Jones particl is

ID ,E­§ V ·12 § V ·6 ½
°
°
4H ®¨

‫ޓ‬
¨ r ¸¸ ¨¨ r ¸¸ ¾°¯© D , E ¹ © D , E ¹ °¿C 2007 The Heat Transfer Society of Japan
○

(5)

vx

1
1
kB
Vf
eoe,
ror,
T oT ,P o P,
24H
24H
V
24H
f
24H
m
,
, V UoU
(6)
t
t
v
v
o
o
m
mV 2
24H

0.15

vx

where İ, ı and rĮ,ȕ are the potential depth, the length parameter and the distance between particles Į and ȕ, respectively. We have the following dimensionless energy,pressure, length, temperature, time, velocity and density.

0
0.15

vx

Thermal Science & Engineering Vol.15 No. 3 (2007)0.15

Ǎxi=20, n=99, subregion i=18

Ǎxi=20, n=999, subregion i=18

0

0
vx

Here kB and f are Boltzmann¶s constant and dimension,respectively. We set f=2 for the two-dimensional MD
method.
Figure 1 shows a particle configuration close to and at
the disturbance generation region on the left end of the
MD system. The unfilled circles (ż) show liquid particles, their number being 17360. The number of particlesof two rigid walls marked by filled circle (Ɣ) at the left
and right ends is 280. The periodic boundary condition is
employed in the y direction. The equation (4) of molecular motion is transformed into the difference equationand the time difference ¨t=0.01 is used. Particles are set
in motion under the conditions of temperature
T=1.95×10-2 and pressure P=5.10×10-3. If the particles
are argon atoms, the time difference is equal to 4.4 fs andthe temperature is 55K. As shown in Fig. 2, we put t=0
after a considerably long time elapses. In order to investigate the direction of the wave, the rigid wall on the left
hand side is moved to the distance 0.677 at an angle of60° relative to the x axis during the time duration of
10.38.
The physical quantities in the continuum mechanics
are obtained from the spatial and temporal averages ofthe properties of MD system. Therefore, the point ofmatter is focused on how large space region and how
long time are required for the wave equation (1). For the
spatial average, a series of subregions are defined in theliquid and are named subregion i=1, 2, 3 and so on fromthe left hand side to the right hand side as shown in Fig.
1. Four kinds of subregion width are employed i.e.
¨xi=20, 40, 80 and 160. We always have ¨yi=20. Thenthe average numbers of particles ¢ Ni² are 288, 577,
1153, and 2307, respectively, where Ni is the number ofparticles in a subregion i, the symbols ¢ ² and ¯ mean
the averages with respect to subregion and time, respectively. The distance between the neighboring subregions
is denoted by ¨xs and is used as x in the wave equation(1). We set ¨xs=¨xi. As shown in Fig. 3, the average time
between ta and tb is expressed as tn=n¨t and four kinds ofn =99, 999, 2999, and 9999 are taken. The values averaged with respect to time are defined as those at the
middle time j (tj). The time interval ¨ts between theneighboring average values is set to be equal to n¨t andis used as t in the wave equation (1).
We define the disturbed velocity, density, total energy and pressure in a subregion i at a time j as

Ǎxi=160, n=999, subregion i=30.15

Ǎxi=160, n=9999, subregion i=3

0
0

200

400

600

800

1000
t

Fig. 4 Time development of disturbed velocity vx .

vi , j·1 ( n 1) / 2 § 1 Ni
mvD ¸
¨
¦
¦
n k  ( n 1) / 2 © N i m D 1¹t  k 't(7)

1 ( n 1) / 2 § N i ·¦ ¨ ¸
n k  ( n 1) / 2 © Si ¹t  k 't

(8)jUi , jj

eti , j·1 ( n 1) / 2 § 1 Ni
¨ ¦ etD ¸¦
n k  ( n 1) / 2 © Si D 1 ¹t  k 't

(9)j

Pi , j·1 ( n 1) / 2 § 1 Ni
¨ ¦ pD ¸¦
n k  ( n 1) / 2 © Si D 1 ¹t  k 't
j

pD

Ni

1
1
mD vD2 ¦ rD ,E  FD ,E
2f E 1f(10)

respectively, where vĮ and etĮ are the velocity and totalenergy of particle Į, respectively, and Si is the area of thesubregion i. The derivatives at both sides of Eq. (1) aretransformed into the difference equations. For example,
the left side is

w 2u
wt 2

ui , j 1  2ui , j  ui , j 1
( 'ts ) 2

䇭

(11)

The macroscopic propagation velocity, c, is calculated
from the MD system by use of Eqs. (1) and (11).

-3-- 113

C 2007 The Heat Transfer Society of Japan
○

2

±0.04

Ǎxi=20, n=99, subregion i=18

w vx /wt

etThermal Science & Engineering Vol.15 No. 3 (2007)


2

±0.04

Ǎxi=20, n=999, subregion i=18

w vx /wt

et

±0.05

Ǎxi=160, n=999, subregion i=3×1010

Ǎxi=20, n=999, subregion i=18

0
±5

×10
10

Ǎxi=160, n=999, subregion i=3

±0.04

2



±0.05
Ǎxi=160, n=9999, subregion i=3

w vx /wt

et

0
±4


2

±0.04

w vx /wt

et

±0.05

0
±6

×1010

Ǎxi=160, n=9999, subregion i=3



±0.05
±0.06

Ǎxi=20, n=99, subregion i=18

0.1

0

200

400

600

800

0

1000
t

0

200

400

600

800

1000
t

Fig. 7 Time development of 2vx/t2.

Fig.5 Time development of disturbed total energy et .
2



0

2

×10 5

w vx /wx×10 5

Ǎxi=20, n=99, subregion i=18

00.15

w vx /wx



0

2

×10 5



0

2

Ǎxi=20, n=999, subregion i=18

×10 5
0

Ǎxi=160, n=999, subregion i=30

Ǎxi=160, n=999, subregion i=3

±6

0.15

w vx /wxvy

Ǎxi=20, n=999, subregion i=18

±6

0.15

w vx /wxvy0

Ǎxi=160, n=9999, subregion i=3

0

±5

0

200

400

600

800

1000
t

Ǎxi=160, n=9999, subregion i=3

0

200

400

600

800

1000
t

Fig.8 Time development of 2vx/x2.

Fig.6 Time development of disturbed velocity Ȟy .
3

Ǎxi=20, n=99, subregion i=18

±4



vy

vy

±4

0.15

Results and Discussions

3.1 Averages with respect to particle and time
Figure 4 shows the time development of disturbed
velocity vx and is arranged from the top to the bottom byamounts of average i.e. (¨xi=20, n=99), (¨xi=20, n=999),(¨xi=160, n=999) and (¨xi=160, n=9999). As shown in
the figure, the disturbed velocity increases abruptly att§250 as if an arrival of a wavelike matter because thevalue of the disturbed velocity is almost zero before thearrival. After then, it decreases gradually and is approaching to zero again. Therefore we can consider thatthe abrupt increase is the disturbed velocity wave. As
shown in the top graph of Fig. 4, the value of the disturbed velocity before the arrival of the wave is not ex-

actly zero and this result is different from one in the continuum mechanics. If the amount of average is large,
the disturbed velocity increases gradually as time
elapses as shown in the bottom graph of Fig. 4. The
value of the disturbed velocity before the arrival of the
wave is almost zero as if in the continuum mechanics.
The maximum value of vx decreases and the small fluctuations of high frequency of the curve disappear as theamount of average becomes large. The time of the maximum disturbed velocity comes later for the case of large
subregion (¨xi=160) relative to small subregion (¨xi=20).
It is concluded that the disturbed velocity curve strongly
depends on the amount of average with respect to timeand space. The features of the disturbed velocity for theother cases of the amount of average are-4-- 114

C 2007 The Heat Transfer Society of Japan
○

±1

10

±2

10

±3

10

±4

10

±5

10
4
10

±6

2

w vx /wx ('xi=20)
2
w vx /wx ('xi=160)

0.03

10

2

3

10

0.10
0.08

0.02



2

w vx /wx

0.78

0.80

Ǎxi=20, n=99

tmax ±xmax
tmax ±xt±xmax
t±x

20

10

200

0

200

400

2

c = 1.049

2

c = 1.156

w vx /wt
600

w vx /wx

wave front group velocity

Ǎxi=160, n=9999


400

200

0

8

200

400

600
t

Fig. 11 Position of wave front (for the case of large
amount of average) as changed with time.

200

400

600
x

tmax ±xmax
tmax ±xt±xmaxt±xǍxi=160, n=9999

6

4

2

0

0

0

Fig. 13 Propagation velocity obtained from small amount of
average vs. the wave motion coordinate x.

c

800

0

600
t

Fig. 10 Position of wave front (for the case of small
amount of average) as changed with time.

x

0.76

30

c = 1.138

400

0

40

wave front group velocity

x

600

0.74

U

c
2

0.72

Fig. 12 Relation between density and pressure.

Ǎxi=20, n=99


p=0.9923ȡ-0.7157

0.70

n

w vx /wt

p=1.0180ȡ-0.7280
p=0.9392U±0.6739

0.04

Fig. 9 Maximum values of disturbed velocity vx and its
derivatives 2vx/x2 and 2vx/t2 vs. number of averages.

800

Ǎxi=160, n=9999, subregion i=3

0.06



0.06



2

w vx /wt ('xi=20)v ('x =20)
2w vx /wt ('xi=160) x ivx('xi=160)

0.09

w vx /wt , w vx /wx

vx

disturbed velocity

2



2

10

0.12

pressure p

Thermal Science & Engineering Vol.15 No. 3 (2007)

0

200

400

600
x

Fig. 14 Propagation velocity obtained from large amount of
average relative to the wave motion coordinate x.
-5-- 115

C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No. 3 (2007)

between the top and bottom figures of Fig. 4. We also
investigated the disturbed density, pressure and total energy. The time developments of disturbed density andpressure are very similar to the one of the disturbed velocity but the total energy (Fig. 5) differs from them.
However, a wave-like pattern can be found also in the
total energy curve although it is not as clear as shown inthe bottom two graphs of Fig. 4.
In order to know in which direction the wave propagates, we moved the rigid wall on the left hand side towards the upper right. As shown in Fig. 6, we can findno wave in the y direction but only small fluctuations.
This result shows that the shear wave such as a distortional wave in solid [5] is not produced in a fluid even in
the nano scale system. The wave propagates in the normal direction relative to the rigid wall surface.
Figures 7 and 8 show the second partial derivativesof the disturbed velocity with respect to time and space,respectively. They are also arranged from the top to bottom in the same manner as Figs. 4, 5 and 6. When theamount of average is small, fluctuations of the second
partial derivatives are remarkable and the steep changesof the derivatives due to the arrival of the wave disappear in the fluctuations. On the other hand, if the large
ensemble and long time averages are used, these steep
changes come out clear as a maximum or a minimum
value.
From Figs. 4, 7 and 8, we assembled the maximum
values of the disturbed velocity and its second partial
derivatives with respect to space and time, respectively,
at the time close to arrival of the wave. They are shownin Fig. 9. They strongly depend on the amount of average when n is large. It can be supposed that the magnitude of the disturbance itself and its second partial derivatives with respect to space and time, respectively,
will disappear when a huge amount of average is used.
We also investigated the disturbed density, pressure
and total energy. They have very similar features to the
disturbed velocity.
From these results, it can be concluded that if theamount of average is large, the peak values of physicalquantities decrease and their temporal and spatial derivatives, respectively, also disappear. The degree of the decrease depends on the physical quantity and the amountof average.
The times of the maximum values in Figs. 7 and 8
correspond to the start of the abrupt increase of the disturbed velocity vx as shown in Fig. 4. Therefore, we canfind the wave front and the arrival time from Figs. 7 and8. The times when 2vx/t2 and 2vx/x2 attain their maximums are not always identical. If the amount of averageis large, this time difference will be very small as described later (Figs. 13 and 14).
3.2 Propagation velocities of the wave of dilatationand distortion
The propagation velocity of the wave can be experi-mentally obtained as a numerical value by measuring thearrival times of the wave front at each subregion. Figure
10 shows the relation between the positions of the centerof subregion and the arrival times of the wave front there.
In order to find the arrival time of the wave front, weused two ways, that is, 2vx/t2 (Fig. 7) and 2vx/x2 (Fig.
8) values. The marks ż denoting 2vx/x2 take positions
along a slightly curved line. It can be supposed that they
are located roughly on a straight line. Then, the wave
front velocity c=1.138 is obtained. If the particles areargon atoms, the velocity is 873m/s and is close to a
general sound velocity of 900~1500m/s for liquid.
Therefore, it can be concluded that the value of 1.138 is
the macroscopic propagation velocity. Let us call it apropagation velocity obtained from the wave front hereafter. On the other hand, the marks Ɣ denoting 2vx/t2
are scattered over a large value region. Therefore, thepropagation velocity cannot be obtained from them. If alarge amount of average is used, the marks Ɣ are on astraight line as shown in Fig. 11. Then, we have wave
group velocity c=1.049 and this value is very close to theprevious one of 1.138. Because the velocities obtained
from Figs. 10 and 11 are almost equal, it can be concluded that the propagation velocity obtained from thewave front is independent of the amount of average except for the very small amount of average.
The propagation velocity is also obtained from themechanical property of the liquid, i.e. Eq. (2). Figure 12
shows the relation between density and pressure obtained for the case of ¨xi=160 and n=9999. Although, therelation is not exactly linear, it is approximated as a linear one shown by a straight line. Then we extrapolate
c=1.009 from the estimation of p by the upper linear
dependence. This value is also close to 1.138.
The macroscopic propagation velocity is also calculated in other way i.e. by the wave equation (1) in thecontinuum mechanics. We find first the times when
2vx/t2 and 2vx/x2 attain maximum or peak around
the arrival time of the wave and then obtain the values of2vx/t2 and 2vx/x2 at the times. This makes four valuesi.e.

(w 2vx / wt 2 ) max 䇭䇭: tt(w2 v / wt 2 )(w 2vx / wt 2 ) 䇭䇭䇭䇭: tt(w 2 v / wx2 )xmax(w 2vx / wx 2 ) max 䇭䇭: tt(w 2 v / wx2 )max(w 2 x v / wx 2 ) 䇭䇭䇭䇭: tt(w 2 v / wt 2 )

x

x

x

max

max

(12)
(13)
(14)
(15)

Combining four values, we have four kinds ofpropagation velocities i.e. {Eq. (12)÷Eq. (14)}½, {Eq.
(12)÷Eq. (15)}½, {Eq. (13)÷Eq. (14)}½ and {Eq. (13)÷
Eq. (15)}½ which are shown by marks Ɣ, ¸, ' and ,
respectively, in Figs. 13 and 14. The horizontal lines in
the figures show the values obtained previously from thewave front. For the case of small amount of average, the-6-- 116

C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No. 3 (2007)

wave front group velocity

cvalues shown by marks ¸, ' and  are considerably
scattered but ones shown by marks Ɣ have almost thesame value (Fig. 13). However, the points marked by Ɣ
are remarkably far from the straight line. The propagation velocity obtained by the small amount of average isvery large. It can be concluded that x and t in the wave
equation (1), which have been approximated by the finite
difference elements with gradually decreasing size, are
not respectively identical with an infinitesimally small
length and time used in Taylor¶s series. Therefore, x andt are finite from the physical point of view.
As shown in Fig. 14, if a large amount of average isused, the difference among the values of marks Ɣ, ¸, 'and  decreases and these marks approach to the straightline. It can thus be said that the disturbed velocity defined in the continuum mechanics is averaged with respect to particle in a region and time in the MD system.

20

wave equation
wave front
wave equation ofpressure, 'xi =160

15'xi'xi = 80
= 40

10

5

'xi = 160'xi
= 20
'xi = 20 ~160

0
10

2

10

3

4

10
n

c

20

wave equation
wave front
wave equation ofpressure, 'xi =160

15

wave front group velocity

wave front group velocitycFig. 15 Relation between propagation velocity of waveof disturbed velocity and space and time averages.

Figure 15 shows the effect of the amount of average on
the propagation velocity obtained by the wave equation.
The horizontal broken line indicates the propagation
velocity obtained from the wave front. If a large subregion, that is, a large space average is applied together
with the small time average, the propagation velocity
obtained from the wave equation differs from that one of
the wave front. In addition to this space average, if thetime average increases, then the propagation velocity isclose to one marked by the horizontal broken line. Whena long time average over n=10000 is used, the propagation velocity deviates from the horizontal broken line
again. The reason is not clear at the present time. In order to make a value close to the propagation velocity
obtained from the wave front, it is necessary to select anappropriate time average for a certain space average.
Figure 16 shows the effect of the amount of averageon the propagation velocity of the disturbed density obtained by the wave equation in the same manner as thedisturbed velocity wave. We have the result that is similar to the disturbed velocity wave.
The marks  in Figs. 15 and 16 show the case of the
disturbed pressure obtained from several numerical experiments when ¨xi=160. The results are also very similar to the cases of the disturbed velocity and the disturbed density.
Figure 17 shows the case of the total energy which
can not be expressed by Eq. (1) along with the potential
and kinetic energies obtained from several numericalexperiments when ¨xi=160. Their results are also very
similar to the cases of the disturbed velocity and the disturbed density.
Kotake concluded that the MD system can not explain non-equilibrium and/or unsteady state macroscopic
phenomena if it has no local uniformities as described in
the introductory section. It can be expected that if a large
amount of average is used, the'xi 'xi 'xi = 160= 40 = 80

10

'xi
= 20

5

20

15

'xi = 8010102'xi = 40'xi = 205'xi = 20 ~160

'xi = 20 ~160
0wave equation
wave front
wave equation of
kinetic energy,'xi =160
wave equation ofpotential energy,'xi =160'xi = 160

10

3

0

4

10

10

2

10

3

4

10
n

n

Fig. 16 Relation between propagation velocity of waveof disturbed density and space and time averages.

Fig. 17 Relation between propagation velocity of waveof disturbed energies and space and time averages.

-7-- 117

C 2007 The Heat Transfer Society of Japan○

Thermal Science & Engineering Vol.15 No. 3 (2007)

physical quantities become locally equilibratory and locally uniform. However, Figs. 15, 16 and 17 show thatthe propagation velocity obtained from the wave equation does not approach to one obtained from the wave
front motion only by increase in the amount of average.
The propagation velocity obtained from the wave equation is identical with the one from the wave front when along time average is used for a large ensemble and ashort time average for a small ensemble. It is necessary
to select an appropriate combination of time and space
averages relative to a particle system being studied. Thisconclusion is the same as in the case of propagation of
the disturbed displacement in an elastic solid continuum[5].
The present calculation of the wave front motion inliquid has been done for the case of two-dimensional(2D) sample. In order to validate our conclusions obtained in 2D reduction for the case of three-dimensional
calculations, one has to take into account the force nature
between particles and two-body approximation of LJ
potential. Presence of neighbors in the third dimension
will cause the interaction potential of each particle to beapproximately 20 to 30% deeper compared with 2D one.
The individual particle motions will also change velocity
characteristics. However, change will be localized
mostly in the vibrational motion of each particle andsystem averages used in this study will be applicable andvalid for the 3D system. Wave front propagation in 3D
particle medium as well as in 2D one is a one-dimensional
simultaneous motion of perturbation of some values, spatially and temporally localized within several layers
normal to the wave motion. With deeper potential andshorter bonds at the same values of temperature andpressure in the system as compared with the 2D case, the3D calculation will show a higher rate of perturbationexchange in the unit of time between layers normal todirection of wave front motion. It will result in the increase in the amplitude of perturbation of around 10 to
15%, in our estimation. The group velocity of the wave
as a collective process should be affected less by thedepth of the particle bonds, the shift of the values insome 10% margins can be expected as we have seen in apressure estimation compared for two-body andthree-body potentials. Exact values can be obtained only
by direct 3D calculations, though the dependencies obtained for the wave front motion should remain similarin the 2D and 3D evaluations.

4

Conclusion

In summary, the question whether and at what conditions the wave equation in the continuum liquid can be
precisely described by particle motions in the MD system has been answered positively and conditions havebeen clarified. For a given system of particles of definedsize and population of species, if an appropriate spatial
and temporal average value of particle motion can be
defined, the averaged characteristics of motion undercondition of wave propagation correspond to the wave
equation of the continuum mechanics. Relating to thesize of the MD model, the spatial region average as it
was shown should be directly proportional to the average
size of time interval. Such wave behavior was confirmed
for the propagation of the velocity, density, and pressure
disturbance wave, as well as for the energy disturbancewave. Then, the wave propagation velocity can be extracted from the continuous wave equation. For
two-dimensional Ar model, obtained values of c in interval from 1.138 to 1.009 are close to general soundvelocities in liquid.
Acknowledgements
The authors thank the late Professor S. Kotake for his
useful discussions and comments. They also thank thestudents, D. Nishi, F. Yoshida and A. Itoh for their computational supports.
References[1] Kotake, S., ´Molecular Dynamics Derivation of
Thermo-Fluid-Dynamics Equations 1. Equations of
Motion´, Thermal Science and Engineering, 10-2
(2002) 47-54.
[2] Kotake, S., ´Molecular Dynamics Derivation of
Thermo-Fluid-Dynamics Equations 2. Equation ofEnergy´, Thermal Science and Engineering, 10-3
(2002) 43-50.
[3] Kotake, S., ´Molecular Dynamics Derivation of
Thermo-Fluid-Dynamics Equations 3. Equation ofContinuity´, Thermal Science and Engineering,
10-5 (2002) 1-8.
[4] Kotake, S., ´Author¶s Comments on µMolecular Dynamics Derivation of Thermo-Fluid-Dynamics¶´,
Thermal Science and Engineering, 11-2 (2003)
23-26.
[5] Iwaki, T., ´Molecular Dynamics Study on Wave
Equation of Solid´, Thermal Science and Engineering, 12-2 (2004) 17-24.

-8-- 118

C 2007 The Heat Transfer Society of Japan○