Int. J. Therm. Sci. (1999) 38, 66-74
~) Elsevier, Paris
Molecular dynamics simulation of an evaporating sodium droplet
Anil P. Bhansali ~,
Yildiz Bayazitoglu ~*,
Shigeo Maruyama ba Department o f Me£hanical Engineering and Materials Science, Rice University, 6100 South Main Street, Houston, TX 77005-1892, USA
b Department of Mechanical Engineering, The University of Tokyo, 7-3-I Hongo, Bunkyo-ku, Tokyo 113, Japan (Received 23 January 1998, accepted 3 June 1998)
Abstract--Developments and advancements have recently been made on the nanoscale level, particularly in the area of the thermal sciences. Since continuum mechanics fail in such phenomena, a demand for molecular level analysis has been created.
Molecular dynamics simulation has proven to be a viable means of microscopic analysis, due primarily to the advanced design of high speed computers. A molecular dynamics simulation is performed to analyse the limits of macroscopic behaviour of an isolated evaporating liquid sodium droplet. Lennard-Jones 12-6 potential is used to determine the intermolecular forces. Details of the simulation are presented as well as variations in properties obtained from the simulation. Parameters such as the bulk liquid density, vapour density, vapour pressure, surface tension, and interfacial surface thickness with respect to temperature are determined. Comparisons of the simulation results to the limiting macroscopic properties are made and trends in the data are discussed. ~) Elsevier, Paris
molecular / droplet / evaporation / phase change / simulation
R~sume - - Simulation de r~vaporation d'une gouttelette de sodium par la dynamique moleculaire. Des d~veloppements importants ont 6t~ r~alis~s r~cemment ~. I'~chelle nanom~trique, particuli~rement dans le domaine de la thermique. Puisque la m~canique du continuum fait d~faut pour I'~tude de tels ph~nom~nes, une demande pour une analyse ~ 1'6chelte mol~culaire a vu le jour. La simulation par la dynamique moleculaire s'av~re un moyen efficace pour une analyse microscopique, grgtce principalement aux performances des ordinateurs ~. haute vitesse de calcul. Dans cette ~tude, on effectue une simulation par la dynamique mol~culaire, dans le but de cerner les limites d'un comportement macroscopique d'une gouttelette isol~e de sodium liquide en {vaporation. On utilise le potentiel 12-6 de Lennard-Jones pour d~terminer les forces intermol~culaires. Les d~tails de la simulation, de m~_me que les variations des propri~t~s obtenues de la simulation, sont pr~sent~s. On d~termine, en fonction de la temperature, des param~tres tels que la densit~ du liquide, la densit6 de la vapeur, la tension de vapeur, la tension superficielle ou I'~paisseur de la surface interfaciale. Des comparaisons entre les r~sultats de la simulation et les propri~t~s macroscopiques limitatives sont pr~sent~es et les tendances obtenues sont interpr~t~es. ~) Elsevier, Paris
mol~culaire / gouttelette / ~vaporation / changement de phase / simulation
Nomenclatu re
B(T)
2 nd virial coel~cientds thickness of the interface . . . F i force on atom i N
kB Boltzmann constant 0.381.10 -23 J) m a mass of sodium atom (3.82.10 -26 kg) M molecular weight of sodium (22.99 kg.mo1-1) N number of atoms in the simulation
* Correspondence and reprints.
£
N(~) P(~)
p *
r *
~'c
number of atoms per spherical shell
local pressure within the droplet . . . M P a reduced pressure ( =
Paa/e)
universal gas constant . . . k J - m o l - l . K -1 magnitude of the radial position as
measured from the centre of the drop _~
reduced position as measured from the centre of the drop ( =
r/s)
cut-off radius used in neighbour list
routine . . . A equimolar radius . . . A
rL neighbour list radius in neighbour list
routine . . . /~
rij magnitude of the distance between atom i and atom j
r0 coefficient in the tanh curve fit (equa- tion (7))
r~ radius of tension . . . A T temperature . . . K T* reduced temperature (= kBT/e)
V velocity . . . m.s - 1 v* reduced velocity (= t , / e X / 7 ~ )
v(r) volume of spherical shells used in den- sity and pressure tensor calculations
A.a
V volume of the simulation box . . . A a
Greek Symbols
Tolman's length parameter . . . A A r width of spherical shells . . . A At timestep . . . s At* dimensionless timestep ( A t / a ~
depth of the potential well . . . l
"7 surface tension . . . N.m -~
q~(r),q5 potential function . . . J
p(r) local density . . . kg.m -3 p* reduced density (= Na3/V)
cr first zero crossing of the potential . . . . /~
Subscripts/Superscripts
act actual
i with reference to atom i j with reference to atom j
k with reference to a co-ordinate direction K kinetic term of the normal pressure tensor 1 liquid phase
new scaled velocity component p set point
U configurational term of the normal pressure tensor v vapour phase
assuming infinite radius of curvature
1. I N T R O D U C T I O N
Recent developments in nanoscale phenomena, both spatially a n d in time, have led to a n increasing n u m b e r of studies dealing with small scale behaviour. Nanoscale processes must be understood on an atomic or molecular level a n d t h e n related to the macroscale system of relevance. Advances, particularly in the t h e r m a l science areas, have created this d e m a n d for molecular level analyses. Kaltz et al. [1] have described the need to s t u d y the evaporation of droplets at supercritical conditions as it relates to combustion in cryogenic
rocket motors on a microscale basis. Kotake [2] has d o c u m e n t e d the need for microscale analysis in heat pipe design, laser applications for material processing, a n d t h e r m a l t r e a t m e n t of metal castings, while Matsui and M a t s u m o t o [3] have given much a t t e n t i o n to rarefied gas flows on a nanoscale level. In particular within the metallurgical industry, an emphasis has been placed on the nanoscale behaviour of metals because of their relevance in metal powder p r o d u c t i o n via nucleation from a s u p e r s a t u r a t e d vapour, chemical vapour deposition processes, a n d the crystallisation of various polymers.
Since liquid metals are evident in m a n y diverse nanoscale processes, it has become i m p o r t a n t to u n d e r s t a n d these processes on a nanoscale level. As early as 1952, a t t e n t i o n has been focused on the need for u n d e r s t a n d i n g the kinetics of nucleating vapour [4, 5]. P o u n d [6] observed the nucleation of various types of vapours a n d approximated t h a t the critical nuclei consisted of only between 80 a n d 100 atoms. For drops this small, statistical theory predicts a deviation in the behaviour of such droplets from macroscopic or c o n t i n u u m analyses. Hill et al. [7] confirmed this notion in their s t u d y of nucleating metal vapours by noting t h a t the surface tension could be in s u b s t a n t i a l error since drops of critical size have such small radii of curvature. Nucleation p h e n o m e n a are small in spatial extent in t h a t the critical radii of the nuclei have been found to be of the order 10 -9 m [8]. Furthermore, the time required for formation of critical nuclei has been documented to be of the order of 10 10 s [9].
Molecular dynanfics simulations have become viable as a tool for analysing systems on a nanoscale level. This can be a t t r i b u t e d primarily to the advanced design of high speed computers. Molecular dynamics formulations are deterministic and consist of solving Newton's equations of m o t i o n for each particle to determine properties of the material. Recently, molecular dynamics simulations have been used to investigate the molecular mechanism t h a t governs m a n y heat transfer processes such as the evaporation and condensation at liquid surfaces [10 12], rarefied gas interaction with a solid wall [3], conduction p h e n o m e n a in t h i n fihns [13], and cluster formation a n d growth [12, 14 16]. The m o n a t o m i c molecular dynamics is simulated by Xu, Xiong and Tzou [17] by t r a n s f o r m a t i o n relationships between one dimension to three dimension. The thermal energy transfer near the critical point of Xenon has been analysed by Ishii, Masuda a n d Maekawa [18]
using molecular dynamics simulation. Hydrogen b o n d in water is analysed over a wide range of t e m p e r a t u r e s and densities by Ohara a n d Aihara [19]. In the present investigation, a molecular dynamics simulation is employed to investigate the interracial p h e n o m e n a a n d to determine properties of a n evaporating liquid metal droplet. Details of the simulation will be presented in the following section.
Molecular dynamics simulations consist of four m a j o r phases: (i) the construction of an adequate potential
m
c
I B
I m
L _
O
,ll,
A.P. Bhansali
t h a t governs the intermolecular forces acting between the i n d i v i d u a l particles; (ii) the initialisation of the simulation and run p a r a m e t e r s ; (iii) the calculation of the molecular t r a j e c t o r i e s and velocities of each particle during simulation, and (iv) t h e analysis of the t r a j e c t o r i e s and velocities to d e t e r m i n e the physical p r o p e r t i e s of the s y s t e m (next section).
In simulating the s o d i u m molecules the well-known L e n n a r d - J o n e s 12-6 p o t e n t i a l (figure 1) is used,
¢ ( r i j ) = 4 s [(¢r/rij) 12 -- (a/rij) 6] (1) where r~j is the s e p a r a t i o n distance between a t o m s i and j , a is the first zero of the p o t e n t i a l (3.24 A), mid ¢ is t h e d e p t h of the p o t e n t i a l well (8.27.10 -21 J).
T h e L e n n a r d - J o n e s p o t e n t i a l exhibits a short range repulsion for s e p a r a t i o n distances less t h a n a and long range a t t r a c t i o n for distances greater t h a n a. This p o t e n t i a l was selected because of its ease of applicability, b u t it could possibly have its l i m i t a t i o n s in modelling liquid m e t a l s such as sodium, because of the short-range oscillatory b e h a v i o u r which liquid metals display.
20
15
~- 10
e-
0
-5
-10 i ) , I ,
2 3 4 5 6 7 8
Separation Distance, rij (Angstroms)
Figure 1. Lennard-Jones potential (sodium).
Initialisation of the simulation is d e p e n d e n t on the a c t u a l s y s t e m being considered. Since it was the intent of this s t u d y to investigate t h e r m a l characteristics of an isolated e v a p o r a t i n g liquid s o d i u m d r o p l e t at various t e m p e r a t u r e s , a microcanonical ensemble was employed in which the n u m b e r of particles, t h e volume, and the t o t a l energy of the s y s t e m were the constrained constants [20]. A cubic simulation box served as t h e control volume which contained 864 s o d i u m a t o m s w i t h the walls of the simulation box replaced with cubic
periodic b o u n d a r y conditions (Allen and Tildesley [22]).
The size of t h e simulation box was chosen sufficiently large such t h a t the periodic b o u n d a r i e s did not interfere w i t h the d r o p ( T h o m p s o n et a l . [21]). The length of the box was therefore chosen to be 80 A. T h e simulation was i n i t i a t e d with each a t o m occupying a position of a perfect face-centre cubic (fcc) lattice s t r u c t u r e in a liquid s t a t e at a t e m p e r a t u r e slightly over the fusion t e m p e r a t u r e . W h i l e other initial configurations such as a b o d y - c e n t r e d cubic (bcc) lattice exist, the fcc s t r u c t u r e is quite a d e q u a t e for simulations such as the one in this s t u d y in which a b n o r m a l l y high densities are not considered [20]. Each a t o m was a r b i t r a r i l y assigned a velocity in each co-ordinate direction (v'k) via a r a n d o m n u m b e r g e n e r a t o r and t h e n scaled to the set point t e m p e r a t u r e , :/-~, as
v~k = v ( 2 )
where subscript i corresponds to a p a r t i c u l a r atom, subscript k represents the co-ordinate direction, and T~*t is the a c t u a l t e m p e r a t u r e of the s y s t e m d e t e r m i n e d from
1E z
T2*~t = 5-~ v~k v~k (3)
k
The variables in the simulation were normalised with respect to the s t a n d a r d p a r a m e t e r s used in t y p i c a l soft- sphere models and are d e n o t e d as such with the * s u p e r s c r i p t [21, 22[.
In order to solve for the molecular t r a j e c t o r i e s and velocities, it was necessary to choose a suitable numerical i n t e g r a t o r to solve Newton's equations of motion. Since no dissipative e x t e r n a l forces exist and since the p o t e n t i a l was assumed pairwise additive, the force on each particle could t h e n be related to the p o t e n t i a l in the following manner:
_ ~ - , a 2 ~
F~ = - ~ V ¢ ( r ~ / ) = m ~t 2 (4) i
A 5th order G e a r predictor-corrector a l g o r i t h m was used to solve the equations of motion. This a l g o r i t h m was found to have superior energy conservation charac- teristics to the other algorithms considered for this prob- lem, such as the Verlet or B e e m a n algorithms [23, 24].
T h e simulation was governed by the t e m p e r a t u r e control scheme employed in the simulation as shown in figure 2. A n equilibration p e r i o d was specified at each set point t e m p e r a t u r e (T~) during which the velocity was scaled according to equation (2). The length of the equilibration periods are i n d i c a t e d by the s h a d e d boxes in figure 2.
T h e r m a l equilibration was m o n i t o r e d by using a nearest neighbour routine t h a t tracked the n u m b e r of v a p o u r a t o m s with respect to time. Each a t o m was classified as vapour, liquid, or surface d e p e n d i n g on the n u m b e r of a t o m s t h a t are within a sphere of ra- dius 1.5 a centred at t h a t p a r t i c u l a r atom. A n a t o m was
1.1
I--I!--
1
-it-
t'.O. 9
~
0.70.6
0 . 5 L . . . _ _ ~
0 200
---, ~ 0
400 600 800 I000 0
Dimensionless Time, t*
500 450 400
~
350~
250Z 150 100 50
i i i L i
200 400 600 800 1000
Dimensionless Time
Figure 2. Temperature control strategy. Figure 3. Vapour atom count as a function of time.
m
C
E M
B g
L _
O
considered v a p o u r if it h a d 1-2 neighbours, interfacial if it had 3 - 7 neighbours, and liquid if it had 8 or more neighbour a t o m s [25]. T h e r m a l equilibrium h a d been reached when the n u m b e r of v a p o u r a t o m s h a d become relatively constant, as i n d i c a t e d in figure 3.
T h e velocity scaling was t h e n t e r m i n a t e d and t h e s y s t e m was allowed to proceed in a s t a t e of constant t o t a l energy d u r i n g which the t r a j e c t o r i e s and velocities were a c c u m u l a t e d for the subsequent d e t e r m i n a t i o n of the physical properties. T h e set point t e m p e r a t u r e was t h e n reset and t h e r m a l equilibrium was t h e n achieved at the new desired state. T h e dimensionless t i m e s t e p used in t h e simulation was At* =0.00144 ( A t ~ 10 -1~ s).
The equilibration p e r i o d became increasingly longer as the set point t e m p e r a t u r e of t h e system was increased and it consisted of up to 70 000 timesteps.
To e x p e d i t e the simulation run time, a Verlet neighbour list routine was i m p l e m e n t e d [26-28] to calculate t h e force on each atom. For a s y s t e m of N atoms, there are N ( N - 1)/2 possible force interactions and the calculation of these interactions is e x t r e m e l y time-consuming. Hence, a t r u n c a t e d p o t e n t i a l was used such t h a t the force for s e p a r a t i o n distances greater t h a n a critical cut-off distance, re, equals zero. T h e n for each atom, the routine m a i n t a i n s a list of neighbouring a t o m s t h a t lie within a distance rL of t h a t p a r t i c u l a r a t o m where t y p i c a l l y rL = re + 0.3 a. T h e neighbour list for each particle was a u t o m a t i c a l l y u p d a t e d based on a m a x i m u m particle displacement criterion, as r e p o r t e d by Verlet.
2. PROPERTY EVALUATION
Once t h e r m a l equilibrium had been reached at a desired t e m p e r a t u r e , the density profile t h r o u g h o u t the d r o p could be calculated by d e t e r m i n i n g the n u m b e r of a t o m s in differential spherical shells of equal width, A r = 0.1 a, t h r o u g h the following equation:
p(r) = ( N ( r ) ) / v ( r ) (5) where ( } denotes an ensemble average taken over the d u r a t i o n of t h e simulation for which t h e r m a l equilibrium exists. T h e shell volume is simply represented by:
4 rc A r (3r2 + A r 2 / 4 ) (6)
v ( r ) - a
where r denotes the m i d p o i n t of the shell as measured from the centre of the drop. T h e centre of the d r o p coincided with the center of t h e simulation box because the centre of mass of the d r o p was recalculated at each time step to c o m p e n s a t e for any bulk m o t i o n of the drop. A regression analysis was t h e n performed used [21] to fit the density d a t a with the commonly used hyperbolic t a n g e n t function such t h a t :
p(r) = ~ (pc + pv) -- ~ (Pc - Pv) t a n h (7) where pe and pv are the liquid density at the centre of the d r o p and v a p o u r density at the b o u n d a r y of the container, respectively, a n d r0 is an e s t i m a t e of the d r o p radius. T h e p a r a m e t e r ds is a measure of
A.P. Bhansali
tile thickness of the surface layer a b o u t the equimolar plane at r = r0. It was also necessary to calculate the equinmlar radius, re, as well as the r a d i u s of surface of tension, rs, since b o t h variables were needed in the d e t e r m i n a t i o n of the surface tension. The surface tension can be a p p r o x i m a t e d t h r o u g h the well-known Laplace equation:
(Pc - Pv) r~
(8)~ / - 2
where Pe and Pv represent the n o r m a l components of the liquid pressure within the d r o p and the v a p o u r pressure well displaced from the centre of the drop, respectively.
The equimolar r a d i u s was a p p r o x i m a t e d as:
r e = 4 ~ ( ~ -- p ~ ) / 3 ] (9)
where rn~ is the mass of one s o d i u m a t o m (3.82.10 -26 kg) and V is the volume of the simula- tion box. r~ is defined such t h a t if the limiting densities of two phases were c o n s t a n t up to r = r~ and changed discontinuously at r = r~. Having d e t e r m i n e d r~ via equation (9), r~ could be calculated using equation (10), which is a r e a r r a n g e m e n t of equation (8) and the Tolman equation given by equation (16) [21]:
3")'~ -- [9")/2 -- 4 7 ~ r e ( P g -- P v ) ] 0"5
r~ = Pe - Pv
(10)where 7 ~ is the surface tension in the macroscopic or p l a n a r limit (i.e. infinite curvature).
The liquid and v a p o u r pressures were d e t e r m i n e d via calculation of the I r v i n g - K i r k w o o d pressure tensor [29]. O n l y the n o r m a l c o m p o n e n t of the pressure tensor, PN(r), was d e t e r m i n e d , since the calculation of t h e surface tension requires only this component. A detailed description of the m e t h o d o l o g y is provided by T h o m p s o n et al. [21], so only the m a j o r points will be highlighted here. The n o r m a l c o m p o n e n t can be w r i t t e n as the sum of kinetic, PK (r), and configurational, Pc (r), contributions:
PN(r)
= PK(r) + Pu(r) where the kinetic t e r m is given byPK(r) = kB T~ct
p(r)/m~
(11)
(12) T h e evaluation of the kinetic c o m p o n e n t is straight- forward and can be o b t a i n e d solely from knowledge of the t e m p e r a t u r e and density profile [equation (7)] at a given instant.
The configurational c o n t r i b u t i o n can be w r i t t e n as 1 E 17" >+ijl 1 d ¢ ( r i j ) (13) P u ( r ) - - 4 ~ r 3 I T i j ~ d r , j
k
where
rij
is a vector between any two a t o m s along which the force between t h e two particles acts and ~b is the pairp o t e n t i a l as defined by equation (1). It was possible for this vector to intersect the surface of any of the various spherical shells t h a t were used in the evaluation of the density profile. Any intersection with a shell results in a c o n t r i b u t i o n of the n o r m a l pressure to t h a t p a r t i c u l a r shell, where each shell is defined by r, the radial distance from the centre of t h e drop. The scalar value of each force interaction t h a t was positive was repulsive while any value t h a t was negative was a t t r a c t i v e . A hyperbolic t a n g e n t flmetion, similar to equation (7), was used to fit b o t h the repulsive and a t t r a c t i v e components of the configurational t e r m s separately. These t e r m s will be discussed more in detail in t h e following section. I t should be noted t h a t the t a n h fit for b o t h the kinetic t e r m and the configurational t e r m was a result of time averaged d a t a over 70 000 t i m e s t e p s in the equilibrium state. T h e liquid pressure could then be d e t e r m i n e d in the limit of r as
Pe = PN(r
-+ 0) (14) T h e v a p o u r pressure was d e t e r m i n e d using two different methods. T h e first m e t h o d involved direct calculation of the virial during the simulation. T h e second m e t h o d included the use of the well-known virial equation as derived from kinetic t h e o r yPv-- Pv-RTM [ I + B(T)p~]M
(15) where R is the ideal gas constant, M is the molecular weight of sodium, andB(T)
is the second virial coefficient which is a function of t e m p e r a t u r e only [30].E q u a t i o n (15) incorporates the correction to the ideal gas relation such t h a t the volume of the v a p o u r does not a p p r o a c h zero as the absolute t e m p e r a t u r e approaches zero for finite pressures, since t h e molecules occupy some finite volume regardless of the t e m p e r a t u r e . T h e virial coefficient,
B(T)
is related to the two force functions of the L e n n a r d - J o n e s p o t e n t i a l [31], since the virial coefficients can be expressed in terms of molecular forces from s t a t i s t i c a l mechanics, a n d evaluated u p o n the selection of the p o t e n t i a l function.3. D I S C U S S I O N
It was necessary to validate the accuracy of the molecular dynamics model. This was accomplished by c o m p a r i n g w i t h d a t a in t h e published literature.
Figure ~
represents a plot of the reduced p o t e n t i a l energy and t h e reduced pressure (virial) o b t a i n e d from t h e current model as well as the average values for an identical simulation performed by Halle[23].
T h e p a r t i c u l a r simulation was performed at a reduced density, p* = 0.6, and a reduced t e m p e r a t u r e , T* = 1.54, using a L e n n a r d Jones potential. Excellent agreement between the two simulations is evident as the p r o p e r t i e s o b t a i n e d from t h e current model fluctuate precisely a b o u t the average values given by Halle[23].
:= 0 (7_ 1
~ - 2 t~
-3
Average Values from Haile
Fluctuations: Values from Current Simulation
Potential Energy, U *
i i i i
-50 2000 4000 6000 8000 10000
Timestep Figure 4. Comparison to Haile's data.
T h e focus of this work was to investigate the lim- its of macroscopic b e h a v i o u r of an isolated e v a p o r a t i n g s o d i u m d r o p l e t by using molecular d y n a m i c s simula- tions. T h e r m a l p r o p e r t i e s of b o t h the bulk liquid as well as t h e s u r r o u n d i n g v a p o u r m e d i u m at various t e m p e r a - tures were determined. Figure 5 depicts this e v a p o r a t i o n process. T h e t e m p e r a t u r e s investigated ranged from 380 to 600 K. T h e t h e r m a l p r o p e r t i e s are s u m m a r i s e d in t h e
table.
T h e simulation was i n i t i a t e d from a perfect fcc lattice s t r u c t u r e and given initial r a n d o m velocities. As was previously discussed, the t h e r m a l equilibration for each s t a t e point consisted of anywhere between 30 000- 70 000 t i m e s t e p s (see figure 2) which was followed by 10 000 t i m e s t e p s of r e l a x a t i o n and 70 000 t i m e s t e p s of p r o d u c t i o n simulation during which p r o p e r t i e s for t i m e averaging were accumulated. T h e r m a l equilibrium was m o n i t o r e d by d e t e r m i n i n g the n m n b e r of v a p o u r a t o m s in the s y s t e m at periodic intervals as was previously discussed. The relatively fiat p l a t e a u x seen in figure 3, in which the n u m b e r of v a p o u r a t o m s is relatively constant, represent the equilibrium periods which are associated with a constant t o t a l energy, as shown in figure 6.
The density profile d a t a along with the hyperbolic t a n g e n t fit for each s t a t e point is shown in figure 7. It can readily be observed t h a t the t a n h fit correlates well w i t h the d a t a . However, the d a t a near the centre of the d r o p l e t display a great deal of fluctuation due to t h e fact t h a t the spherical shells used in the density calculations [equations (5) and (6)] become very small as the centre of the d r o p is approached. As expected, the liquid d e n s i t y of the d r o p l e t decreases as the t e m p e r a t u r e increases, while the v a p o u r density exhibits the opposite trend.
Initial Configuration
$$
T = 3 8 6 . 6 K
O
o o
O O (3
O O
o oo ~ °
o 0
T--444.6 K D
o o
og
O o O
O )
T=486.7 K
T=536.3 K
(
Pc
T=582.1 K
Figure 5. Schematic of the evaporation process.
i t can also be observed fl'om figure 7 t h a t the interface between the liquid and the v a p o u r is much more distinct at lower t e m p e r a t u r e s . But, as the t e m p e r a t u r e increases, the interface experiences a g r a d u a l spreading. As a result, the surface thickness exhibits a r a p i d increase. T h o m p s o n et al. [21], however, concluded in their comparison to the work performed by Powles et al. [32], t h a t the length of the cut-off radius could significantly affect the values of the surface thickness values and t h a t special a t t e n t i o n should be given to such values if e x t r e m e l y low cut-off values are used.
T h e macroscopic limiting values for the liquid den- sity corresponding to the values o b t a i n e d in simula- tion shown in the t a b l e are 923.3, 909.6, 899.6, 887.9, 877.1 kg.m -3, respectively, from lowest to highest tem- perature. For lower t e m p e r a t u r e s , the d a t a over pre- dicts the macroscopic limit, b u t as the t e m p e r a t u r e continues to increase (as the d r o p l e t shrinks), the d a t a crosses below the macroscopic limit. This t r e n d was also observed by T h o m p s o n et al. [21]. T h e y at- t r i b u t e d the decrease below the macroscopic limit to
a
as c
# m
~ m
L _
@
ii,!!!N ~ ~::
A.P. Bhansali
-0.5 -1
~-1.5
L~
-2
o
"~ - 2 . 5 - 3
-3.5
- 4
-4.5 i
0 2 0 0
/
i i i T
4 0 0 6 0 0 8 0 0 1 0 0 0
D i m e n s i o n l e s s T i m e Figure 6. Total energy versus time.
a g r a d u a l l y decreasing cohesion within the d r o p as the size decreased. However, this t r e n d was not observed by M a r u y a m a et al. [25] in their s t u d y of an e v a p o r a t i n g argon droplet.
T h e profiles of the configurational t e r m of the n o r m a l pressure tensor, P u ( r ) , defined in equation
(13),
are shown in figure 8 for the two e x t r e m e t e m p e r a t u r e s . It was previously s t a t e d t h a t t h e configurational t e r m was divided into a repulsive c o m p o n e n t for positive force interactions and an a t t r a c t i v e c o m p o n e n t for negative force interactions. T h e profile of the kinetic t e r m of the pressure tensor, PK(r) mirrors the d e n s i t y profile [see equation (12)] and is therefore not reproduced. A hyperbolic t a n g e n t function was again used to fit the d a t a for b o t h repulsive and a t t r a c t i v e c o m p o n e n t s and the limiting values ( r - * 0) were o b t a i n e d by equation (14). T h e repulsive c o m p o n e n t profiles, a t t r i b u t e d to short range interactions, display more u n c e r t a i n t y near the centre of the drop, t h a n its a t t r a c t i v e c o u n t e r p a r t s , a t t r i b u t e d to longer range interactions, due p r i m a r i l y to the fact t h a t fewer short range interactions exist which makes the averaging process less reliable. The v a p o u r pressure d e t e r m i n e d from t h e simulation is c o m p a r e d to the macroscopic bulk values as a function of t e m p e r a t u r e , as shown in figure 9. W h i l e b o t h m e t h o d s for d e t e r m i n i n g the v a p o u r pressure c o m p a r e quite well with each other, t h e y b o t h over predict the macroscopic bulk values at lower t e m p e r a t u r e s but seem to a p p r o a c h the macroscopic limit as t e m p e r a t u r e increases.The limiting macroscopic values of the surface tension d a t a corresponding to the values o b t a i n e d from the simulation via L a p l a c e ' s equation [equation (8)] are 0.062, 0.051, 0.054, 0.032, and 0.025 N.m -1 respectively
from lowest to highest t e m p e r a t u r e . The average of the two values o b t a i n e d for the v a p o u r pressure was used as the value of Pv in equation (8). T h e surface tension d a t a are much lower t h a n the p l a n a r limit which can be e x p e c t e d based on T o l m a n ' s pioneering work [33].
T h o m p s o n et al. [21] also observed such a deviance from the macroscopic limit with the discrepancy being much larger at smaller values of equimolar radius considered in this work. Tolman d o c u m e n t e d an expected decrease in surface tension as the radius of c u r v a t u r e decreased, especially in the range of d r o p sizes considered in this work. Tohnan proposed an equation to account for the variation in the surface tension with drop size
7 26
- - = 1 - - - ( 1 6 )
7~ rs
where 6 is the difference re - rs. Values for 6 are given in the table L
T(K) (r*)
386.6 (0.65) 444.6 (0.74) 486.7 (0.81) 536.3 (0.89) 582.1 (0.97)
TABLE I
Summary of properties from the simulation.
Pen pv re ds rs Pe *Pc 7 5
(kg.m 3) (kg.m 3) (A) (/~) (/~) (MPa (MPa) (Nm) (~) 982.7 2.1 19.8 5.1 14.8 83.4 0.24 0.O625.O5
(0.29)
940.5 6.6 19.6 5.9 14.2 72.6 1 . 0 6 0.051 5.45 (0.99)
907.2 12.8 19.2 6.8 14.2 8 1 . 1 2.10 0.054 4.96
(2.05)
864.6 20.1 18.6 7.8 13.3 57.0 3.84 0.032 5.34 (3.48)
830.6 36.7 16.2 9.4 11.4 5 1 . 1 6.66 0.0254.86 (6.30)
* The above state.
vapour pressure obtained from the virial is listed the vapour pressure obtained from the equation of
The value of the surface tension at T = 486.7 K s t r a y s from the monotonic decline in surface tension with respect to t e m p e r a t u r e . However, u p o n inspection of figure 2 again, it can be seen t h a t the t e m p e r a t u r e experienced a g r a d u a l decrease d u r i n g the p r o d u c t i o n or equilibrium p e r i o d for this s t a t e point. This indicates t h a t while the d a t a are quite good, b o t h longer equilibration times might be needed in order to b e t t e r define an equilibrium s t a t e and longer p r o d u c t i o n periods might be required to make the t i m e averaging more accurate.
As mentioned earlier, the need for molecular d y n a m - ics analyses is quite a p p a r e n t in m e t a l powder produc- tion processes, in p a r t i c u l a r , ones t h a t involve nucle- ation from a s u p e r s a t u r a t e d vapour. If the e v a p o r a t i o n process described here was continued by increasing the set point t e m p e r a t u r e , t h e n the d r o p would eventu- ally e v a p o r a t e and a uniform v a p o u r s t a t e point would exist. Subsequently, if the t e m p e r a t u r e were t h e n sud- denly lowered (to model a s u d d e n expansion process),
1200
1000
800
o. 600
400
200 Increasing Temperature (T* = 0.65, 0.81, 0.97)
I ~ I I I [ -,
2 4 6 8
Reduced Position from Drop Center, r*
10 Figure 7. Density profile at various temperatures.
4 ~ s i v e (T* = 0.65)
~ A
-4
~ Attractive (1"* = 0.65) -6
0 1 2 3 4 5 6 7
Reduced Position from Drop Center, r*
F i g u r e 8. Reduced pressure profiles.
the vapour would exist in the thermodynamically unsta- ble supersaturated state required for the h o m o g e n e o u s nucleation to initiate [34].
4. CONCLUSIONS
A molecular dynamics simulation was employed to investigate the t h e r m a l characteristics of an evaporat-
5
0~:" 4
~d
>
350
O Vixial (from Simulation)
o/
r a //
D Virial Equation ofState (Eq. 15) /
400 450 500 550 600
Temperature, T(K)
Figure 9. Vapour pressure variation with temperature.
ing sodium droplet. Various t h e r m a l properties were o b t a i n e d from the simulation a n d compared to the lim- iting macroscopic values. Some of the trends in the d a t a were also observed in the previously published literature, b u t there are still some questions left un- resolved. While the surface tension values o b t a i n e d by M a r u y a m a et al. [25] accurately predict the macroscopic limit for argon, which is an insulator, the surface tension d a t a obtained for sodium, a liquid metal conductor, do not reveal such a comparison. It has been well docu- m e n t e d t h a t the pair potentials of liquid metals exhibit a n oscillatory behaviour due to ion screening which is not accounted for in the L e n n a r d - J o n e s potential [35].
Since the pioneering work of J o h n s o n et al. [36], there has been a growing interest in developing more suitable pair potentials for liquid metals [3~39]. Such potentials are now being implemented to determine their effects on the t h e r m a l behaviour of liquid metals.
Acknowledgements
This w o r k w a s in part supported by Texas A d v a n c e d Technology P r o g r a m under Grant No. 003604-041 a n d National Science Foundation under Grant No. C T S - 9312379. W e also acknowledge the c o m p u t e r time m a d e available by the Pittsburgh S u p e r c o m p u t i n g Center.
REFERENCES
[1] Kaltz T., Little J., Wong B., Micci M., Long L.N.
Supercritical droplet evaporation modeled using molecular dynamics on parallel processors, in: Euromech Colloquium 324, The Combustion of Drops, Sprays, and Aerosols, Marseilles, France, 1994, pp. 1-9.
W l
O
A.P. Bhansali
[2] Kotake S., Future aspects of molecular heat and mass transfer studies, Thermal Sci Eng. 2 (1) (1994) 12- 20.
[3] Matsui J., Matsumoto Y., Molecular dynamics of gas-surface interaction, ASME/JSME Thermal Engineering Proceedings 4 (1) (1991) 59-164.
[4] Reiss H., Theory of the liquid drop model, Ind. Eng.
Chem. 44 (6) (1952) 1284-1288.
[5] Michaels A.S., in: Proceedings of Nucleation Phe- nomena Conference, A Chemical Study Publications, Wash- ington, DC, 1969.
[6] Pound G.M., Liquid and crystal nucleations, Ind.
Eng. Chem. 44 (6) (1952) 1278-1283.
[7] Hill P.G., Witting H., Demetri E.P., Condensation of metal vapors during rapid expansion, J. Heat Trans.- T. ASME (1963) 303-317.
[8] Lothe J., Pound G.M., Statistical mechanics of nucle- ation, in: Zettlemoyer A.C. (Ed.), Nucleation, Marcel Dekker, Inc, New York, 1969, p. 112.
[9] Mandell M.J., McTague J.P., Rahman A., Crystal nu- cleation in a three-dimensional Lennard-Jones system: a molecular dynamics study, J. Chem. Phys. 64 (9) (1976) 3699-3702.
[10] Yasuoka K., Matsumoto M., Kataoka Y., Evapora- tion and condensation at a liquid surface. I. Argon, J.
Chem. Phys. 101 (9) (1994) 7904-7911.
[11] Matsumoto M., Yasuoka K., Kataoka Y., Micro- scopic Features of Evaporation and Condensation at Liquid Surfaces: Molecular Dynamics Simulation, 1994, Vol. 2, No. 1, pp. 1-6.
[12] Kotake S., Aoki X., Atomic and molecular clusters and their film condensation, in: ICHMT Symposium Proceedings on Molecular and Microscale Heat Transfer in Material Processing and other Applications, 1996, pp. 118- 129.
[13] Wakuri S., Kotake S., Molecular dynamics study of heat conduction in very thin films, ASME/JSME Thermal Engineering Proceedings 4 (1991) t 11-116.
[14] Rey C., Gallego L.J., Iniguez M.P., Alonso J.A., A molecular dynamics study of the evaporation of small argon clusters, Physica B 179 (1992) 273-277.
[15] Nagasaki T., Hijikata K., Miyazaki K., Kaneko H., Itoh H., Fundamental study on formation of metal clusters, in: ICHMT Symposium Proceedings on Molecular and Microscale Heat Transfer in Material Processing and other Applications, 1996, pp. 106-117.
[16] Sugiyama T., Echigo R., Yoshida H., Cluster for- mation around the critical point in two- dimensional Lennard-Jones particle system, in: ICHMT Symposium Pro- ceedings on Molecular and Microscale Heat Transfer in Material Processing and other Applications, 1996, pp. 91- 105.
[17] Xu Y.S., Xiong D.X., Tzou R.D., One dimensional molecular dynamics simulation, in: ICHMT Symposium Proceedings on Molecular and Microscale Heat Transfer in Material Processing and other Applications, 1996, pp. 41- 49.
[18] Ishii K., Masuda S., Maekawa T., Thermofluid dy- namics and molecular dynamics analyses of thermal energy transfer near the critical point, in: ICHMT Symposium Pro- ceedings on Molecular and Microscale Heat Transfer in
Material Processing and other Applications, 1996, pp. 58- 69.
[19] Ohara T., Aihara T., MD studies on dynamic structure of water, in: ICHMT Symposium Proceedings on Molecular and Microscale Heat Transfer in Material Processing and other Applications, 1996, pp. 70-80.
[20] Rowley R.L., Statisical Mechanics for Thermophys- ical Property Calculations, PTR Prentice Hall, Englewood Cliffs, New Jersey, 1994, pp. 41 et 236.
[21] Thompson S.M., Gubbins K.E., Walton J.P.R.B., Chantry R.A.R, Rowlinson J.S., A molecular dynamics study or liquid drops, J. Chem. Phys. 81 (1) (1984) 530-522.
[22] Allen M.P., Tildesley DJ., Computer Simulation of Liquids, Clarendon Press, New York, 1987, pp. 17 et 87.
[23] Halle J.M., Molecular Dynamics Simulation, John Wiley& Sons, inc., New York, 1994, pp. 147-176 & 231.
[24] Amini M., Fincham D., Evaluation of temperature in molecular dynamics simulation, Comput. Phys. Commun.
56 (1990) 313-324.
[25] Maruyama S., Matsumoto S., Ogita A., Surface phenomena of molecular clusters by molecular dynamics method, Thermal Sci. Eng. 2 (1)(1994) 77-84.
[26] Verlet L., Computer experiments on classical fluids. I. Thermodynamic properties of Lennard-Jones Molecules, Phys. Rev. 159 (1967) 98-103.
[27] Arnold A., Mauser N., An efficient method of bookkeeping next neighbours in molecular dynamics simulations, Comput. Phys. Commun. 59 (1990) 267-275.
[28] Chialvo A.A., Debendetti P.B., On the performance of an automated verlet neighbor list algorithm for large systems on a vector processor, Comput. Phys. Commun.
64 (1991) 15-18.
[29] Tsai D.H., The Virial Theorem and Stress Calcu- lation in Molecular Dynamics, J. Chem. Phys. 70 (1979) 1375.
[30] Reed T.M., Gubbins K.E., Applied Statistical Me- chanics, McGraw-Hill, New York, 1973, pp. 202-207.
[31] Hirschfelder J.O., Curtiss C.F., Bird R.B., Molecular Theory of Gases and Liquids, John Wiley & Sons, New York,
1964.
[32] Powles J.G., Fowler R.F., Evans W.A.B., A new method for computing surface tension using a drop of liquid, Chem. Phys. Lett. 98 (1983) 421-426.
[33] Tolman R.C., The effect of droplet size on surface tension, J. Chem. Phys. 17 (1949) 333-337.
[34] Stever H.G., High Speed Aerodynamics and Jet Propulsion, Section F., Princeton University Press, 1954, pp. 532-536.
[35] Ashcroft N.W., Mermin N.D., Solid State Physics, Holt Rinehart and Winston, New York, t976, pp. 337-345.
[36] Johnson M.D., Hutchinson P., March N.H., Ion-ion oscillatory potentials in liquid metals, P. Roy. Soc. Lond.
A 282 (1964) 283-302
[37] Price D.L., Shyu W., Singwi K.S., Tosi M.P., Ar- gonne National Report, ANL 7761, 1970.
[38] March N.H., Structure and Forces in Liquid Metals and Alloys, Can. J. Phys. 65 (1987) 219-240.
[39] Rani M., Pratap A., Saxena N.S., Single Particle and Collective Motions in Liquid Rubidium, Ind. J. Pure Appl.
Phys. 27 (1989) 269-274.
