Computational Materials Science 40 (2007) 382–388 www.elsevier.com/locate/commatsci
Conceptual study on the thermal equilibrium state of nanoscale systems Hyun-Kue Choi a, Soon-Ho Choi b
b,*
a Korea Coast Guard, Songdo-dong, Yeonsu-ku, Incheon 573-701, South Korea Korea Maritime University, M258, Marine Engineering Department, #1, Dongsam-dong, Youngdo-ku, Pusan 606-791, South Korea
Received 9 July 2006; received in revised form 18 January 2007; accepted 24 January 2007 Available online 21 March 2007
Abstract According to classical thermodynamics, a state of equilibrium signifies that in a system the thermophysical properties must be equal regardless where the measurement is taken. However, some recent studies have reported that there is a temperature difference over a liquid–vapor interface even though the system is in a saturated equilibrium state when reduced to nanoscale size. Those conclusions lead to questions as to whether the definition of a thermal equilibrium state should be altered for nanoscale systems. In this study, we carefully investigated the temperature profile over the liquid–vapor interface using the molecular dynamics (MD) simulation and found that the temperature discontinuities between a liquid phase and a vapor phase, over an interface, resulted from the methodological difference of the temperature calculation. In other words, the results of this study indicate that there are no sufficient proofs that the classical definition of a thermal equilibrium state should be changed even if the system of interest would be reduced to nanoscale size. 2007 Elsevier B.V. All rights reserved. PACS: 68.10; 68.35.M; 71.15.D Keywords: Thermal equilibrium state; Liquid–vapor interface; Molecular dynamics; Nanoscale system
1. Introduction As well known, the performances of computers have been dramatically improved in the last two decades and one cannot predict the limit of its improvement. With the improved performance and popularization of computers, the studies that use computer simulations have exerted a deep influence on various research fields. More importantly, it is not rare that the computer simulations have played an essential role in investigating and evaluating the characteristics of nanoscale phenomena in scientific and engineering applications [1,2]. Recently, the studies to investigate the characteristics of a liquid–vapor interface have been performed since it is closely related with mass *
Corresponding author. Present address: The College of Marine Science, Gyeongsang National University, 445, Inpyeong-dong, Tongyeoung-si, Gyeongsangnamdo 650-160, South Korea. Tel.: +82 51 555 1801. E-mail address:
[email protected] (S.-H. Choi). 0927-0256/$ - see front matter 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2007.01.011
and heat transfer phenomena, in which we took notice of the specific studies conducted on the temperature boundary layer at a liquid–vapor interface. These studies have reported that, when a system is reduced to a nanoscale size, there is a temperature difference over the interface even though the system is in a thermal equilibrium state [3–6]. The reported temperature discontinuities over a liquid– vapor interface arouse highly intriguing questions as to whether in case of nanoscale systems the definition of a thermal equilibrium state should be changed. According to classical thermodynamics, if a system is in a state of equilibrium, the properties of interest should display the same value no matter where a measurement is taken [7]. For example, in saturated liquid–vapor mixture, the temperatures of a liquid region and a vapor region should be identical. Therefore, if the reported temperature difference over an interface is proven, it is necessary to adopt a new criterion of defining a thermal equilibrium state in nanoscale systems. Although there are many molecular dynamics (MD) studies that surveyed the features
H.-K. Choi, S.-H. Choi / Computational Materials Science 40 (2007) 382–388
occurred over the interface of a liquid–vapor [8–11], a solid– liquid [12–15], and a solid–solid [16–19], there is no study within the authors’ reviews on the related literatures that ascertains the temperature discontinuity over the interface when a system is in a thermal equilibrium state except for the previously mentioned studies [3–6]. Regarding the temperature discontinuity over a liquid– vapor interface in an equilibrium state that is the main focus of the present study, Rao and Levesque first reported that the temperature of a liquid–vapor interface was higher than that of a liquid region [20]. As to non-equilibrium state, Maruyama and Kimura [11,12] observed a temperature jump at the solid–liquid interface in their MD simulation for analyzing the processes of evaporation and condensation. They estimated the quantitative temperature jump at the solid–liquid interface that is often called Kapitza resistance [21] and calculated the equivalent liquid thickness corresponding to the thermal resistance at an interface. Even though the system studied by Maruyama et al. was non-equilibrium state, since a high temperature region for evaporation and a low temperature region for condensation coexist in the system, any remarkable temperature discontinuity over the liquid–vapor interface was not seen. At nearly the same time with Maruyama et al., Fang and Ward also performed the experimental study for investigating the characteristics of a liquid–vapor interface to be evaporated in a vacuum state [21]. In contrast to the results of Maruyama et al., they reported the noticeable temperature discontinuity over the interface, in which the temperature of a vapor region was higher than that of a liquid region. However, their results could not be directly compared with the ones by Maruyama et al. because their experimental system was physically quite different from the simulation system investigated by Maruyama et al. Furthermore, the experimental system by Fang and Ward cannot be regarded in an equilibrium state due to the evaporating interface. Therefore, there is a strong need to clear the ambiguity on the temperature discontinuity over a liquid–vapor interface in a saturated state, i.e., in a thermal equilibrium state. In the present study, we carefully investigated the temperature profiles over the liquid–vapor interface using the molecular dynamics (MD) simulation and argue that the reported temperature difference between a liquid phase and a vapor phase in the previous studies [3– 5] might have resulted from the methodological difference in calculating the temperature profiles of a system. In addition, if a system is at least in the state of a thermal equilibrium and the temperature is properly calculated, there are no sufficient proofs that the temperature discontinuity over a liquid–vapor interface appears when a system size reduces to nanoscale. 2. Simulation system and method For the simulation system to confirm a temperature discontinuity over a liquid–vapor interface, argon molecules were selected as a target material since its behavior is well
383
described by the simple Lennard-Jones (L-J) potential of Eq. (1). r 12 r6 uðrÞ ¼ 4 e ; ð1Þ r r ˚ , and r is the intermolecwhere e = 1.67 · 1021 J, r = 3.4 A ular distance. By differentiating Eq. (1) with respect to r, the intermolecular force is given as r 12 1 r6 1 : ð2Þ F ðrÞ ¼ 48 e r 2 r r In the MD simulation, the time step (Dt) for iteration is very important because an inappropriate selection will not produce an intended state. According to our simulation experiences [17,18,22–27] and other studies [3–15], one femto second (1 fs) is reasonable for a solid state, 5 fs for a liquid–vapor mixture or a liquid state, and 10 fs for a vapor state if the temperature and the pressure of a simulation system are moderate far from a critical state (Dt = 5 fs in this study). However, we do not repeat the detailed methods for the MD simulation to calculate an intermolecular potential, an intermolecular force, and the time evolution of a system since there are too many references and textbooks [28–30] to get information on them. Fig. 1 shows the initial arrangement and the time evolution of the simulation system with the 3-dimensional periodic boundary conditions (3 PBC). Since we aimed at producing the state of a saturated equilibrium state, the liquid region was placed in the middle of the system. In the simulation of a liquid film contacted with vapor molecules, however, it should be paid attention to the selections of a liquid film thickness and an intermolecular length for the formation of a stable liquid film. According to Weng et al. the liquid film could not be maintained if the initial liquid layers were arranged less than 10r [10]. Choi et al. reported that, for the stable suspension of a liquid film as shown in Fig. 1, the intermolecular length, l, calculated from the experimental values of argon provides the criterion of a minimum value [24]. That is, the liquid film will be broken if the initial intermolecular length is selected shorter than the value calculated from the experimental density. Since the fluid molecules have random motions independent of an initial arrangement [22–25], the 8000 molecules (20 · 20 · 20) of argon were arranged in the middle region as the simple cubic structure (SCS) unlike the case of a solid state that is generally arranged with the face centered cubic (FCC). The intermolecular length can be easily calculated from the experimental data since there is only one molecule in a cubic if the SCS is assumed. We intended to simulate the state of the 5% quality (v = 0.05), and hence the inter˚ for the molecular length was selected as l = 3.698 A ˚ saturated state of 100 K, and l = 3.768 A for the saturated state of 110 K [31]. The corresponding system sizes ˚ · 73.97 A ˚ · 357.74 A ˚ and were, respectively, 73.97 A ˚ ˚ ˚ 75.37 A · 75.37 A · 211.88 A in x, y, and z directions. The simulation system was first temperature-controlled through the velocity scaling method performed every 50
384
H.-K. Choi, S.-H. Choi / Computational Materials Science 40 (2007) 382–388
Fig. 1. Snap shots to show the time evolution of a simulation system for the saturated state of 100 K. Total number of the molecules are 8000 and the size ˚ · 73.97 A ˚ · 357.74 A ˚ in x, y, and z of the system was selected so as to have the volume equal to the quality of 0.05, which has the dimensions of 73.97 A ˚ · 75.37 A ˚ · 211.88 A ˚. directions. For the saturated state of 110 K, the dimension is 75.37 A
iterations during the initial run (100 000 steps), which corresponds to 500 ps. After being confirmed that the system was settled down to the desired temperature, the system was relaxed during another 500 ps for the elimination of any undesirable effect that might be resulted from the artificial temperature control in the initial simulation period. Then, the third run was also relaxed for the data production during the same period without any manipulation. For the comparison with the third run, the fourth run was performed in the case of 100 K; however, no meaningful difference was found. Therefore, all data shown in the present study were taken from only the third run. 3. Temperature calculation and simulation results As described previously, Guo et al. [3], Chen et al. [4], Wang et al. [5], and Rao and Levesque [20] reported that the temperature of a vapor region was different from that of a liquid region in the system if the system would be reduced to the nanoscale even in saturated state. These observations are considerably interesting and controversial because the definition of an equilibrium state from classical thermodynamics is no more applicable to nanoscale systems. If those observations were true, it is newly required what a criterion ensures a thermal equilibrium state of nanoscale liquid–vapor system. However, in contrast to the above studies [3–5,20], Choi et al. [24] recently reported that the observed temperature discontinuity over a liquid– vapor interface in the MD simulations had been to blame for the calculation method of a vapor region temperature. In the MD simulations, the temperature and the number
density of a system are evaluated from the time-averaged kinetic energy and the number of molecules, which are given as [28–30]: ! N sys n mv2j 1X 1 X T sys ¼ ; ð3Þ n i¼1 N sys j¼1 fk B qsys ¼
N sys ; V sys
ð4Þ
where kB is the Boltzmann constant (1.38 · 1023 J K1); m is the mass of a argon molecule (6.634 · 1026 kg); v is the velocity of a molecule; n is the iteration number of one run; Nsys is the total number of the molecules in a system; f is the number of degree of freedom including a transnational, a rotational and a vibrational motion, which would be three in this study since argon is a monatomic molecule [28–30]; Vsys is the system volume. However, if we try to evaluate the local profiles of the temperature and the number density, the system should be divided into the slabs with the same thickness, and then the temperatures and the number densities of each slab are calculated from ! N slab n mv2j 1X 1 X T slab ¼ ; ð5Þ n i¼1 N slab j¼1 3k B ! n 1 1X N slab ; ð6Þ qslab ¼ V slab n i¼1 where Nslab is the molecule’s number contained in a slab; Vslab is the volume of a sliced slab. It should be noticed
H.-K. Choi, S.-H. Choi / Computational Materials Science 40 (2007) 382–388
385
that, for keeping the slab temperature from diverging, Tslab is set to zero in the simulation if there is no molecule in the slab. According to Choi et al.’s study performed with 1000 molecules [24], they indicated that there were so many instants to exist nothing within the slabs sliced in a low-density region during a simulation because the MD simulations are typically performed with hundreds of or thousands of molecules; therefore, the temperature of the sliced slab in a vapor region has to be calculated by considering the instants not to exist a molecule. That is, it is necessary in the calculation of a local temperature that the time to exist no molecule should be excluded from the overall simulation time as T slab
! N slab n X mv2j 1 1 X ¼ ; n nempty i¼1 N slab j¼1 3k B
ð7Þ
where nempty is the instants to exist nothing in the slab of a vapor region. However, the calculated temperatures from Eqs. (5) and (7) will have the same results in a liquid region because a molecule certainly exists at any time in a highdensity region. Fig. 2 depicts the temperature profiles of the system with the 8000 molecules as shown in Fig. 1. From this figure, it is apparent that there is no temperature discontinuity over an interface even in a nanoscale system if Eq. (7) is used for the temperature calculation of a vapor region. The difference between the temperatures calculated from Eqs. (5) and (7) was fairly large. Although the suggestion by Choi et al. which was our previous study, is convincible, they did not attempt more additional analysis for verifying their suggestion. Actually, if Eq. (7) should be used to calculate the temperature of a vapor region, the density of it also has to be calculated by ( ) n X 1 1 qslab ¼ N slab : ð8Þ V slab n nempty i¼1 The number density calculated by Eq. (8) for a vapor region will bring higher value than that by Eq. (6) and Choi et al. [24] should have provided the detailed interpretation on both the number densities and the temperature profiles in the vapor region. In the present study, for justifying our previous suggestion, we performed the simulations with a larger system compared with the previous one and tried to analyze the simulation results in more detail. As a matter of fact, if the temperature profile calculated by Eq. (5) in Fig. 2 was real, it means that the high energetic molecules emerged from a liquid region does not interact with the low energetic molecules in a vapor region in spite of intermolecular collisions, and vice versa. This unrealistic situation cannot be imagined at all compared with our physical intuition or experimental observations and it is caused by including the instants not to exist anything when calculat-
Fig. 2. Temperature and density profiles of the system with the 8000 molecules as shown in Fig. 1. The upper one (a) is for Tsat = 100 K with v = 0.05; the lower one (b) for Tsat = 110 K with v = 0.05. The dashed lines are the temperature calculated by using Eq. (5); the solid lines by using Eq. (7); the smooth lines in the middle are the number density profiles. The temperature profiles were calculated based on the slab ˚ . The two horizontal dotted lines on each graph are the thickness of 1 A experimental number density evaluated from Ref. [31]; the upper ones are the value of a saturated liquid and the lower ones are that of a saturated vapor.
ing the local temperatures of each slab in the vapor region. We believe that it is meaningless to measure any physical quantity in an empty space because ‘‘There is nothing to exist’’ means ‘‘There is nothing to measure’’. However, in any real experiment for a macroscale system that has the order of 1023 molecules, there cannot be any gap between two calculated temperatures by Eq. (5) and by Eq. (7) since any slab even with an angstrom thickness in a vapor region certainly has some molecules. We also show the number density profiles in Fig. 2 for the saturate temperatures of 100 K and 110 K when the ˚ . From the comparison of slab thickness was sliced as 1 A Fig. 2 (8000 molecules) and our previous results (1000 molecules) [20], we confirmed that the fluctuation of a temperature in a vapor region was considerably decreased and the gap between two temperature profiles calculated from Eqs. (5) and (7), which will be simply described as DT
386
H.-K. Choi, S.-H. Choi / Computational Materials Science 40 (2007) 382–388
hereinafter, was also reduced if we perform the simulations with a larger system. Furthermore, it was observed that the fluctuations of the temperature profiles and DT in a vapor region were significantly reduced as the slab thickness increases. Fig. 3 shows the temperature profiles for the same saturated temperatures but the slabs were sliced with the thick˚ , in which the temperature profiles marked by ness of 2 A Eq. (5) are duplicated from Fig. 2 for the clear comparison. Although Fig. 3b seemingly appears to be the temperature gradient between two phases, it was intentionally selected from Fig. 2b and the temperature gradient between two phases is only from the effect of the fluctuations as explained in below. Comparing Fig. 3 with Fig. 2, it is clearly seen that DT, which the temperature difference between the temperature profiles evaluated from Eqs. (5) and (7), is fairly dimin-
Fig. 3. Temperature profiles in a vapor region when the slabs were sliced ˚ . The blank circles are calculated by using Eq. (5); the solid squares as 2 A ˚ , it can be seen by Eq. (7). Comparing with Fig. 2 that was sliced as 1 A that the gap of two temperature profiles calculated from Eqs. (5) and (7) is reduced from DT = 20–25 K to DT = 5–10 K in the case of Tsat = 100 K and is nearly negligible DT = Max. about 0.2 K in the case of Tsat = 110 K, which results from the increased molecules due to the thicker slab thickness in a vapor region.
ished. Furthermore, we found that DT is completely disappeared when the slab thickness was selected as more ˚ for Tsat = 100 K and 3 A ˚ for Tsat = 110 K. In than 5 A other words, if we select the slab thickness more than two times of L-J molecule’s diameter for the temperature profiles, the temperature difference between two phases is completely disappeared except for some fluctuation. However, the reduced DT is quite natural because the chance to exist nothing in a vapor region would be appreciably reduced if the slab thickness gets increased. Therefore, it is reasonable to exclude the instants to exist nothing within the sliced slab when calculating its time-averaged temperature. Fig. 4 depicts the enlarged density profiles of the section ‘‘A’’ and ‘‘B’’ indicated in Fig. 2, which were evalu-
Fig. 4. Number density profile in a vapor region when the slabs were ˚ . The blank circles are calculated from sliced with the thickness of as 2 A Eq. (6); the solid squares from Eq. (8). Comparing with Fig. 3, the density profile adversely behaviors, i.e., the evaluated number density value goes far from the experimental value when excluding the instants to exist nothing in the slab for the local density calculation. The difference between the two density profiles from Eqs. (6) and (8) is negligibly small, which results from the increased molecules due to the thicker slab in a vapor region. The difference of the two densities was not distinguishable for the ˚ in the case of Tsat = 100 K, and 3 A ˚ in the case of slab thickness of 5 A Tsat = 110 K. The two horizontal dotted lines on each graph are the experimental number density evaluated from Ref. [31].
H.-K. Choi, S.-H. Choi / Computational Materials Science 40 (2007) 382–388
ated from the same concept as the temperature calculation ˚ and 1 A ˚ . As seen in Fig. 4, the for the slabs sliced as 2 A number density profiles calculated from Eq. (8) for the ˚ goes up from the value calculated slab thickness of 1 A from Eq. (6) and the experimental number density and its trend is completely opposed to the temperature profile. ˚ , the number density from When the slab thickness is 2 A Eq. (8) is slightly larger than that from Eq. (6), however, the difference between two number densities is negligible. Moreover, the number density from Eq. (8) for the slab ˚ is nearly same as the number density thickness of 2 A ˚ , which is also from Eq. (6) for the slab thickness of 1 A the opposite trend in comparison with a temperature profile. This observation is quite confusing and interesting if considering that the temperature and the number density are both intensive property. Although we cannot surely ascertain why this contradictory observation is obtained in a temperature and a density, it might be resulted from the intrinsic attribution of each property. The temperature is calculated from the kinetic energy of molecules, and hence the existence of a molecule is necessarily prerequisite, which means that the temperature is dependent on the existence of matter but not on a space. Therefore, a pure empty space has no physical meaning for the temperature evaluation. However, in the case of a density, it can be seen from Eq. (6) or Eq. (8) that even empty space has a physical meaning because the density is simultaneously dependent with a space and the molecules contained in it. Therefore, on the contrary to the case of a temperature, the instants to exist nothing in the slab should not be excluded when calculating the local density evaluation. The following thought experiment would be helpful to understand that our statement is sound and rational. Consider the system as shown in Fig. 5. Although only one molecule exists in the space between the walls in Fig. 5a, the temperature and the density of the space can be measured because there is a physical entity with a kinetic energy. On the other hand, Fig. 5b shows the purely empty space. In this case, the density would be well defined as zero and no one brings this definition into question. However, if trying to define the temperature of an empty space, an intractable ambiguity arises because there is no molecule for the calculation of a temperature. Then, is the temperature of the empty space absolutely zero? Otherwise, is it meaningless to define the temperature of the empty space? This ambiguity can be cleared by assuming that we insert a thermometer into the empty space. Even though the initial temperature of a thermometer itself is different from the wall temperature, the heat transfer by radiation will be occurred between the thermometer and the walls. Due to this heat transfer, the thermometer eventually reaches to the wall temperature and displays the value of Tw. Therefore, we suggest that, in the MD simulations, the empty instants of no molecule in a slab should be excluded for the temperature calculation while they should be included for the density calculation.
387
Fig. 5. The systems for a though experiment. (a) Not empty space. Although there is only one molecule, the temperature and the density of the space between walls will be defined because matter is inside. (b) Purely empty space. Nevertheless, the density of the space would be well defined as zero. Then, how should one define the temperature of the space since no molecule exists in it? Is it absolutely zero or Tw? However, if one inserts a thermometer in the empty space, it would display the value of Tw due to the thermal equilibrium between the walls and the inserted thermometer after some time elapse.
4. Conclusion According to the recent MD studies, it was reported that the temperatures of each phase in a liquid–vapor mixture were differently observed even though the system would be thermally equilibrated when the system reduces to a nanoscale size. These reports raise a contradictory question on the definition of a state of equilibrium from the viewpoint of classical thermodynamics. If it is true, a new criterion on a thermal equilibrium state for nanoscale systems is require, which would be a highly challengeable research subject. However, this conceptual study suggests that the observed temperature discontinuity over a liquid–vapor interface in some earlier MD studies might have resulted from the methodological difference to calculate the temperature in a vapor region. In briefly, the instants to exist no molecules should be excluded in the time-averaged temperature of a local space while those instants should be included in the time-averaged density of a local space. Therefore, the time-averaged properties must be carefully taken after convincing the detailed characteristics of them in the MD simulations with only hundreds or thousands of molecules. Furthermore, if adopting our suggestion proposed in the present study, the thermal equilibrium state defined by classical thermodynamics can be still applicable even to nanoscale systems.
388
H.-K. Choi, S.-H. Choi / Computational Materials Science 40 (2007) 382–388
Nevertheless, we do not deny the possibility to be observed the temperature profiles indicated by Eq. (5) in Fig. 2. So far, although there is no developed device to measure the temperature of a nanoscale vapor space, the temperature discontinuity will be observed over a liquid– vapor region if the device has the responsiveness of fs and a time-averaged display function. Consequently, this study still leaves an intellectually stimulating question about whether a real experiment obtains a temperature discontinuity over a liquid–vapor interface if a measurement device, which has the time-averaged characteristic and can measure the temperature of a nanoscale vapor region, will be used in coming future. Actually, if there comes such a day, what criterion defines a thermal equilibrium state of nanoscale systems? Otherwise, we cannot define a thermal equilibrium state of nanoscale systems? Acknowledgements This work was supported by Korea Research Foundation Grant (KRF-2004- 037-D00003). The author, SoonHo Choi, deeply thank for KRF. References [1] F.C. Chou, J.R. Lukes, S.G. Liang, K. Takahashi, C.L. Tien, Annu. Rev. Heat Transfer 10 (1995) 141. [2] D.G. Cahill, W.K. Ford, K.E. Goodson, G.D. Mahan, A. Majumdar, H.J. Maris, R. Merlin, S.R. Phillpot, Appl. Phys. Rev. 93 (2003). [3] Z.Y. Guo, D.X. Xiong, C. Yang, M. Chen, Z.X. Li, Int. J. Therm. Sci. 39 (2000) 481. [4] M. Chen, Z.Y. Guo, X.G. Liang, Micro. Thermophys. Eng. 5 (2001) 1. [5] Z.J. Wang, M. Chen, Z.Y. Guo, C. Yang, Fluid Phase Equilibria 183–184 (2001) 321.
[6] G.J. Van Wylen, R.E. Sonntag, Fundamentals of Classical Thermodynamics, John Wiley and Sons, Hoboken, 1991. [7] M.J.P. Nijmeijer, A.F. Bakker, C. Bruin, J.H. Sikkenk, J. Chem. Phys. 89 (1988) 3789. [8] S. Maruyama, S. Matsumoto, A. Ogita, Therm. Sci. Eng. 2 (1999) 77. [9] J.E. Lee, S.H. Park, O.M. Kwon, KSME Int. J. 16 (2002) 1477. [10] J.G. Weng, S.H. Park, J.R. Lukes, C.L. Tien, J. Chem. Phys. 113 (2000) 5917. [11] S. Maruyama, T. Kimura, Therm. Sci. Eng. 7 (1999) 63. [12] T. Kimura, S. Maruyama, Therm. Sci. Eng. 8 (2000) 7. [13] T. Ohara, T. Yatsunami, Micro. Thermophys. Eng. 7 (2003) 1. [14] T. Ohara, D. Suzuki, Micro. Thermophys. Eng. 4 (2000) 189. [15] A.R. Abramson, C.L. Tien, A. Majumdar, ASME J. Heat Transfer 124 (2002) 963. [16] M. Matumoto, H. Yakayabashi, T. Makino, Trans. JSME B 68 (2002) 87. [17] S.-H. Choi, S. Maruyama, K.-K. Kim, J.-H. Lee, J. Kor. Phys. Soc. 44 (2004) 317. [18] S.-H. Choi, S. Maruyama, Int. J. Therm. Sci. 44 (2005) 547. [19] S.-P. Jang, S.U.S. Choi, Appl. Phys. Lett. 84 (2004) 4316. [20] M. Rao, B. Levesque, J. Chem. Phys. 65 (1976) 3233. [21] G. Fang, C.A. Ward, Phys. Rev. E 59 (1999) 417. [22] H.-M. Kim, K.-H. Park, H.-K. Choi, S.-H. Choi, J. KOSME 29 (2005) 34. [23] H.-K. Choi, H.-M. Kim, S.-Y. Choi, K.-K. Kim, S.-H. Choi, J. KOSME 29 (2005) 60. [24] H.-K. Choi, C.-S. Song, H.-M. Kim, J.-H. Lee, S.-H. Choi, Trans. KSME B 29 (2005) 159. [25] S.-H. Choi, J. KOSME 28 (2004) 441. [26] S.-H. Choi, S. Maruyam, J. Kor. Phys. Soc. 45 (2004) 897. [27] S.-H. Choi, S. Maruyama, J.H. Lee, K.K. Kim, J. Kor. Phys. Soc. 43 (2003) 747. [28] D. Frenkel, B. Smit, Understanding Molecular Dynamics Simulation-From Algorithms to Applications, Academic Press, San Diego, 1996. [29] J.M. Haile, Molecular Dynamics Simulation-Elementary Methods, John Wiley and Sons, New York, 1992. [30] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford Univ, Press, New York, 2004. [31] JSME, Thermophysical Properties of Fluids, Tokyo, JSME, 1983.