2000WJM_EMS_microHeatTransf

W. J. Minkowycz and E. M. Sparrow (Eds), Advances in Numerical Heat Transfer, vol. 2, Chap. 6, pp. 189-226, Taylor & Francis, New York, 2000. CHAPTER SIX MOLECULAR DYNAMICS METHOD FOR MICROSCALE HEAT TRANSFER Shigeo Maruyama 1 INTRODUCTION Molecular level understandings and treatments have been recognized to be more and more important in heat and mass transfer research. A new field, “Molecular Thermophysical Engineering,” has a variety of applications in further development of macroscopic heat transfer theory and in handling the extreme heat transfer situations related to advanced technologies. For example, studies of basic mechanisms of heat transfer such as in phase change heat transfer demand the microscopic understanding of liquid-solid contact phenomena. The nucleation theory of liquid droplet in vapor or of vapor bubble in liquid sometimes needs to take account of nuclei in size of molecular clusters. The efficient heat transfer in three-phase interface (evaporation and condensation of liquid on the solid surface) becomes the singular problem in the macroscopic treatment. Some modeling of the heat transfer based on the correct understandings of molecular level phenomena seems to be necessary. The effect of the surfactant on the heat and mass transfer through liquid-vapor interface is also an example of the direct effect of molecular scale phenomena on the macroscopic problem. The surface treatment of the solid surface has a similar effect. Even though there has been much effort of extending our macroscopic analysis to extremely microscopic conditions in space (micrometer scale and nanometer scale system), time (microsecond, nanosecond and picosecond technology), and rate (extremely high heat flux), there is a certain limitation in the extrapolations. Here, the development of the molecular dynamics (MD) computer simulation technique has shown the possibility of taking care of such microscale phenomena from the other direction. The MD methods have long been used and are well developed as a tool in statistical mechanics and chemistry. However, it is a new challenge to extend the method to the spatial and temporal scale of 2 macroscopic heat transfer phenomena. On the other hand, by developments of high energy-flux devices such as laser beam and electron beam, more physically reasonable treatment of heat transfer is being required. The thin film technology developed in the semiconductor industry demands the prediction of heat transfer characteristics of nanometer scale materials. In this chapter, one of the promising numerical techniques, the classical molecular dynamics method, is first overviewed with a special emphasis on applications to heat transfer problems in section 2 in order to give the minimum knowledge of the method to a reader not familiar with it. The van der Waals interaction potential for rare gas, effective pair potential for water and many-body potential for silicon and carbon are discussed in detail. Then, the molecular scale representation of the liquid-vapor interface is discussed in section 3. The surface tension, Young-Laplace equation, and condensation coefficient are discussed from the viewpoint of molecular scale phenomena. Section 4 deals with the solid-liquid-vapor interactions. MD simulations of liquid droplet in contact with solid surface and a vapor bubble on solid surface are introduced. The validity of Young’s equation of contact angle is also discussed. Then, demonstrations of real heat transfer phenomena are discussed in section 4. Since heat transfer is intrinsically a non-equilibrium phenomenon, the non-equilibrium MD simulations for constant heat flux system and the homogeneous nucleation of liquid droplet in supersaturated vapor and nucleation of vapor bubble in liquid are discussed. Then, the heterogeneous nucleation of vapor bubble on the surface is also discussed. Some interesting non-equilibrium MD simulations dealing with the formation of molecular structures are introduced in section 5.4. Finally, in section 6, future developments of molecular scale heat transfer are discussed. 2 MOLECULAR DYNAMICS METHOD Knowledge of statistical mechanical gas dynamics has been helpful to understand the relationship between molecular motion and macroscopic gas dynamics phenomena [1]. Recently, a direct simulation method using the Monte Carlo technique (DSMC) developed by Bird [2] has been widely used for the practical simulations of rarefied gas dynamics. In the other extreme, statistical mechanical treatment of solid-state matters has been well developed as solid state physics [e.g. 3]. For example, the direct simulation of the Boltzmann equation of phonon is being developed and applied to the heat conduction analysis of thin film [4] for example. However, when we need to take care of liquid or inter-phase phenomenon, which is inevitable for phase-change heat transfer, the statistical mechanics approach is not as much developed as for the gas-dynamics statistics and the solid-state statistics. The most powerful tool for the investigation of the microscopic phenomena in heat transfer is the MD method [e.g. 5]. In principal, the MD method can be applied to all phases of gas, liquid and solid and to interfaces of these three phases. 2.1 Equation of Motion and Potential Function 3 In the MD method, the classical equations of motion (Newton's equations) are solved for atoms and molecules as i i i i dt d m r F r ∂ Φ∂ −==2 2 , (1) where mi, ri, Fi are mass, position vector, force vector of molecule i, respectively, and Φ is the potential of the system. This classical form of equation of motion is known to be a good approximation of the Schrödinger equation when the mass of atom is not too small and the system temperature is not too low. Equation (1) itself should be questioned when applied to light molecules such as hydrogen and helium and/or at very low temperature. Once the potential of a system is obtained, it is straightforward to numerically solve Eq. (1). In principal, any of gas, liquid, solid states, and inter-phase phenomena can be solved without the knowledge of "thermo-physical properties" such as thermal conductivity, viscosity, latent heat, saturation temperature and surface tension. The potential of a system ),...,( N21 rrrΦ can often be reasonably assumed to be the sum of the effective pair potential φ(rij) as )( ij�� > = i ij rφΦ , (2) where rij is the distance between molecules i and j. It should be noted that the assumption of Eq. (2) is often employed for simplicity even though the validity is questionable. The covalent system such as carbon and silicon cannot accept the pair-potential approximation. 2.2 Examples of Potential Forms In order to simulate practical molecules, the determination of the suitable potential function is very important. Here, the well-known Lennard-Jones potential for inert gas and for a statistical mechanical model system is introduced; also introduced are potential forms for water and many-body potential for silicon and carbon. The interaction potential forms between metal atoms are intentionally excluded because the luck of the effective technique of handling free electron for heat conduction prevents from the reasonable treatment of heat conduction through solid metal. 2.2.1 Lennard-Jones potential. An example of the pair potential is the well-known Lennard-Jones (12-6) potential function expressed as � � � � � � � � � � � � −� � � � = 612 4)( rr r σσεφ , (3) 4 where ε and σ are energy and length scales, respectively, and r is the intermolecular distance as shown in Fig. 1. The intermolecular potential of inert monatomic molecules such as Ne, Ar, Kr and Xe is known to be reasonably expressed by this function. Typical values of σ and ε for each molecule are listed in Table 1. Moreover, many computational and statistical mechanical studies have been performed with this potential as the model pair potential. Here, the equation of motion can be non-dimensionalized by choosing σ, ε and m as length, energy and mass scale, respectively. The reduced formulas for typical physical properties are listed in Table 2. When a simulation system consists of only Lennard-Jones molecules, the non-dimensional analysis has an advantage in order not to repeat practically the same simulation. Then, molecules are called Lennard-Jones molecules, and argon parameters σ = 0.34 nm, ε = 1.67×10-21 J, and τ = 2.2 ×10-12 s are used to describe dimensional values in order to illustrate the physical meaning. The phase-diagram of Lennard-Jones system [6] is useful for a design of a simulation. The critical and triplet temperatures are Tc* = 1.35 and Tt* = 0.68, or Tc = 163 K and Tt = 82 K with argon property [7]. For the efficient calculation of potential, which is the most CPU demanding, Lennard-Jones function in Eq. (3) is often cutoff at the intermolecular distance rC = 2.5 σ to 5.5 σ. However, for pressure or stress calculations, the contribution to potential from far-away molecules can result in a considerable error as demonstrated for surface tension [8]. In order to reduce this discrepancy, σ 1.5σ 2σ –ε 0 ε 2ε Intermolecular Distance, r Po te nt ia l E ne rg y, φ( r) r σ6 2 Figure 1 Lennard-Jones (12-6) potential. Table 1 Parameters for Lennard-Jones potential for inert molecules. σ [nm] ε [J] ε/kB [K] Ne 0.274 0.50×10-21 36.2 Ar 0.340 1.67×10-21 121 Kr 0.365 2.25×10-21 163 Xe 0.398 3.20×10-21 232 5 several forms of smooth connection of cutoff have been proposed such as in Eq. (4) by Stoddard & Ford [9]. ( ) ( ) ( ) � � � � � � � � −−�� � � � −+ � � � � � � � � � � −� � � � = −−−− 6* C 12* C 2 C 6* C 12* C 612 47364 rr r rrr rr r σσεφ (4) 2.2.2 Effective pair potential for water. The effective pair potential form for liquid water has been intensively studied. The classical ST2 potential proposed in 1974 by Stillinger and Rahman [10] based on BNS model [11] was widely used in the 1980s. The rigid water molecule was modeled as Fig. 2a, with the distance of OH just 0.1 nm and the angle of HOH the tetrahedral angle θt = ( )3/1cos2 1− ≅ 109.47°. Point charges at four sites shown in Fig. 2a were assumed: positive charge of 0.235 7 e each on hydrogen sites and two negative charges at positions of lone electron pairs (tetrahedral directions). They modeled the potential function as the summation of Coulomb potential between charges and the Lennard-Jones potential between oxygen atoms. Hence, the effective pair potential of molecules at R1 and R2 are expressed as Table 2 Reduced properties for Lennard-Jones system. Property Reduced Form Length r* = r/σ Time t* = t/τ = t(ε/mσ2)1/2 Temperature T* = kBT/ε Force f* = fσ/ε Energy φ* = φ/ε Pressure P* = Pσ3/ε Number density N* = Nσ3 Density ρ* = σ3ρ/m Surface tension γ* = γσ2/ε rOH rOM O +qH +qH -qM -qM ∠HOH rOM O +qH +qH ∠HOH -qM rOH r26 O6 H1 H2 H3 H4 -qM -qM +qH +qH +qH +qH O5 r14 r56 r16 (a) (b) (c) Figure 2 Water potential structures for (a) 5 sites model, ST2, (b) 4 sites and 3 sites models, TIP4P, CC, SPC, SPC/E, (c) definition of interatomic length of MCY and CC potential. 6 ��+ � � � � � � � � �� � � � −�� � � � = i j ij ji r qq RS RR 0 12 6 12 OO 12 12 OO OO2112 4 )(4),( πε σσ εφ RR , (5) where R12 represents the distance of oxygen atoms, and σOO and εOO are Lennard-Jones parameters. The Coulombic interaction is the sum of 16 pairs of point charges. S(R12) is the modulation function to reduce the Coulombic force when two molecules are very close. Later, much simpler forms of SPC (Simple Point Charge) [12] and SPC/E (Extended SPC) [13] potentials were introduced by Berendsen et al. SPC/E potential employed the configuration in Fig. 2b, with charges on oxygen and hydrogen equal to –0.8476 and +0.4238 e, respectively. Lennard-Jones function of oxygen-oxygen interaction was used as ST2 as in Eq. (5) but without the modulation function S(R12). TIP4P potential proposed by Jorgensen et al. [14] employed the structure of water molecule as rOH = 0.09572 nm and ∠HOH = 104.52° based on the experimentally assigned value for the isolated molecule. The positive point charges q were on hydrogen atoms, and the negative charge –2q was set at rOM from the oxygen atom on the bisector of the HOH angle, as in Fig. 2b. The function can be written as Eq. (5) without S(R12) function. The parameters listed in Table 3 were optimized for thermodynamic data such as density, potential energy, specific heat, evaporation energy, self-diffusion coefficient and thermal conductivity, and structure data such as the radial distribution function and neutron diffraction results at 25 °C and 1atm. This potential is regarded as one of the OPLS (optimized potential for liquid simulations) set covering liquid alcohols and other molecules with hydroxyl groups developed by Jorgensen [15]. MYC potential [16] and CC potential [17] were based on ab initio quantum molecular calculations of water dimer with the elaborate treatment of electron correlation energy. The assumed structure and the distribution of charges are the same as TIP4P as shown in Fig. 2b with a different length rOM and amount of charge as in Table 3. For CC potential, the interaction of molecules is parameterized as follows. Table 3 Potential parameters for water. ST2 SPC/E TIP4P CC rOH [nm] 0.100 0.100 0.095 72 0.095 72 ∠HOH [°] 109.47 109.47 104.52 104.52 σOO [nm] 0.310 0.316 6 0.315 4 N/A εOO ×10-21 [J] 0.526 05 1.079 7 1.077 2 N/A rOM qHa [nm] [C] 0.08 0.235 7 e 0 0.423 8 e 0.015 0.52 e 0.024 994 0.185 59 e qM [C] -0.235 7 e -0.847 6 e -1.04 e -0.371 18 e aCharge of electron e = 1.60219×10-19 C 7 [ ] [ ] [ ])exp()exp()exp()exp( )exp()exp()exp()exp( )exp()exp()exp()exp( )exp( 4 ),( 4543542641644 4533532631633 2422321421322 5611 0 2112 rbrbrbrba rbrbrbrba rbrbrbrba rba r qq i j ij ji −+−+−+−− −+−+−+−+ −+−+−+−+ −+=�� πε φ RR (6) a1 = 315.708 ×10-17 [J], b1 = 47.555 [1/nm], a2 = 2.4873 ×10-17 [J], b2 = 38.446 [1/nm], a3 = 1.4694 ×10-17 [J], b3 = 31.763 [1/nm], a4 = 0.3181 ×10-17 [J], b4 = 24.806 [1/nm]. Among these rigid water models, SPC/E, TIP4P and CC potentials are well accepted in recent simulations of liquid water such as the demonstration of the excellent agreement of surface tension with experimental results using SPC/E potential [18]. Because all of these rigid water models are “effective” pair potential optimized for liquid water, it must be always questioned if these are applicable to small clusters, wider range of thermodynamics condition, or liquid-vapor interface. Even though the experimental permanent dipole moment of isolated water is 1.85 D1, most rigid models employ higher value such as 2.351 D for SPC/E to effectively model the induced dipole moment at liquid phase. The direct inclusion of the polarizability to the water models results in the many-body potential, which requires the iterative calculation of polarization depending on surrounding molecules. The polarizable potential based on TIP4P [19], MCY [20] and SPC [21] are used to simulate the structure of small clusters and transition of monomer to bulk properties. On the other hand, flexible water models with spring [22] or Morse type [23] intramolecular potential are examined seeking for the demonstration of vibrational spectrum shift and for the reasonable prediction of dielectric constant. 2.2.3 Many-body potential for carbon and silicon. The approximation of pair potential cannot be applied for atoms with covalent chemical bond such as silicon and carbon. SW potential for silicon proposed by Stillinger and Weber in 1985 [24] was made of two-body term and three-body term that stabilize the diamond structure of silicon. Tersoff [25, 26] proposed a many-body potential function for silicon, carbon, germanium and combinations of these atoms. For simulations of solid silicon, this potential [26] is widely used. Brenner modified the Tersoff potential for carbon and extended it for a hydrocarbon system [28]. A simplified form of Brenner potential removing rather complicated ‘conjugate terms’ is widely used for studies of fullerene [29, 30] and carbon-nanotube. Both Tersoff potential and the simplified Brenner potential can be expressed as following in a unified form. The total potential energy of a system is expressed as the sum of every chemical bond as 1 1 D = 3.3357×10-30 Cm in SI unit. 8 { }� � < −= i jij ijijijij rVbrVrf )( A * RC )()()(Φ , (7) where the summation is for every chemical bond. VR(r) and VA(r) are repulsive and attractive parts of the Morse type potential, respectively. ( ){ }eeCR 2exp1)()( RrSSDrfrV −−−= β (8) ( ){ }eeCA /2exp1)()( RrSS SDrfrV −−−= β (9) The cutoff function fC(r) is a simple decaying function centered at r = R with the half width of D. ( ) ( ) ( )��� �� � � +> +<<−�� � � − π − −< = DRr DRrDRDRr DRr rf 0 /)( 2 sin 2 1 2 1 1 )(C (10) Finally, b*ij term expresses the modification of the attractive force VA(r) depending on θijk, the bond angle between bonds i-j and i-k. 2 * jiij ij bb b + = , δ θ − ≠ � � � � � � � � � � � � += � n ijk jik ikC n ij grfab )()(1 ),( (11) 0.15 0.2 0.25 0.3 –4 –2 0 2 4 Distance rij [nm] Po te nt ia l E ne rg y φ ij [e V] (Re, De) θijk=45° θijk=90° θijk=126.7° and 2–body θijk=180° R rij θijk i j k Figure 3 Many-body characteristics of Tersoff potential for silicon. 9 22 2 2 2 )cos( 1)( θ θ −+ −+= hd c d cg (12) Parameter constants for Tersoff potential for silicon (improved elastic properties) [26] and carbon and Brenner potential for carbon are listed in Table 4. In order to illustrate the characteristic of Tersoff and Brenner potential function, a potential energy contribution from a bond is expressed in Fig. 3. The Tersoff parameters for silicon are assumed and the energy of i-j bond under the influence of the third atom k, { })()()(' ARC ijijijij rVbrVrf −=φ is drawn. The effect of the third atom k is negligible only when the angle θijk is 126.7°. 2.3 Integration of the Newtonian Equation The integration of the equation of motion is straightforward. Unlike the simulation of fluid dynamics, simpler integration scheme is usually preferred [5]. Verlet’s integration scheme, as follows, can be simply derived by the Taylor series expansion of the equation of motion. ( ) ( ) ( ) ( ) iiiii mttttttt )(2 2 Frrr ∆+∆−−=∆+ (13) ( ) ( ) ( ){ } tttttt iii ∆∆−−∆+= 2rrv (14) where ∆t is the time step. A bit modified leap-frog method, as follows, is widely used in practical simulations [5]. After the velocity of each molecule is calculated Table 4 Parameters for Tersoff potential and Brenner potential. Tersoff (Si) Tersoff (C) Brenner (C) De [eV] 2.6660 5.1644 6.325 Re [nm] 0.2295 0.1447 0.1315 S 1.4316 1.5769 1.29 β [nm-1] 14.656 19.640 1.5 A 1.1000×10-6 1.5724×10-7 1.1304×10-2 N 7.8734×10-1 7.2751×10-1 1 δ 1/(2n) 1/(2n) 0.80469 C 1.0039×105 3.8049×104 19 D 1.6217×101 4.384 2.5 H -5.9825×10-1 -5.7058×10-1 -1 R [nm] 0.285 0.195 0.185 D [nm] 0.015 0.015 0.015 10 as Eq. (15), the position is calculated as Eq. (16). ( ) i i ii m tttttt Fvv ∆+� � � � � � ∆ −=� � � � � � ∆ + 22 (15) ( ) ( ) � � � � � � ++=+ 2 tttttt iii ∆∆∆ vrr (16) Typical time step ∆t is about 0.005 τ or 10 fs with argon property