# Noncontinuum Approach for Interfaces and Free Surfaces

More detail about what is actually happening at an interface requires a different approach, because the continuum assumption breaks down. For these situations, a molecular dynamic simulation is needed. Research at the molecular level on vapor formation, condensation, and other interface phenomena is fairly scarce due to the limitations of experimental instruments, because the usual macroscopic continuum numerical models do not describe the phenomena at such a small scale. A promising approach to analyzing phenomena at the molecular level is molecular dynamic simulation (MDS). By solving Newton’s equation of motion for every molecule in the system numerically, information on the macroscopic system can be found. The basic idea behind MDS is that as every substance is made of particles (atoms or molecules), if the dynamics of these particles, such as velocity, position, and interaction force, can be determined, then physical macroscopic properties such as temperature, pressure, volume, etc. may also be found.

The scale treated by MDS is only about 10 nm in space and 10 ns in time with the latest computers. The system needs to be decomposed into elements that can be examined using MDS. Since MDS gives only the transient position and velocity of the molecules in the system, macroscopic quantities such as temperature, internal energy, pressure, heat flux, and momentum flux are obtained by averaging the MDS system to interpret the phenomena.

In MDS, each atom is treated as a point mass and the interactions between atoms are described by Newton’s second law of motion. The equation of motion for each atom or molecule under consideration may be expressed as ${{\mathbf{F}}_{i}}={{m}_{i}}{{\mathbf{a}}_{i}}={{m}_{i}}\frac{{{d}^{2}}{{\mathbf{r}}_{i}}}{d{{t}^{2}}}=\sum\limits_{j\ne i}^{N}{-\frac{d{{\phi }_{ij}}(r)}{dr}\frac{{{\mathbf{r}}_{ij}}}{{{r}_{ij}}}}$ (1)

where Fi is the total force vector acting on the ith molecule; N is the number of molecules; is the pair potential between the ith molecule and the jth molecule; rij represents the position vector and rij is the distance between the ith molecule and the jth molecule; mj is the mass of the ith molecule; and ai is the acceleration vector of the ith molecule. The pair potential approximation does not give the correct potential energy between two molecules when they are isolated in space. Therefore, more accurate models are required to give the correct overall potential of the system. The most basic model is the Lennard-Jones potential model for spherical molecules, which is presented in Chapter 1, eq. (1.28).

MDS has been used to determine the evaporation and condensation coefficients of fluids, such as water and argon, because experimental determination of these values is very difficult. This is due to the large temperature increase in the interface region, which is made up of only a few molecule layers and cannot be directly measured. Some examples of MDS in recent research are: Yasuoka and Matsumoto (1995) used MDS to obtain the condensation coefficient of argon, methanol, and water at certain temperatures; Matsumoto (1998) used MDS to explore evaporation-condensation dynamics of pure fluids under both equilibrium and non-equilibrium conditions; Wu and Pan (2003) investigated the evaporation of a thin argon layer into a vacuum using MDS with the Lennard-Jones potential; Yang and Pan (2005) examined the evaporation of a thin layer of water and determined an evaporation coefficient.

Yang and Pan (2005) used the TIP4P model as the effective pair potential between water molecules. The model consisted of four sites: two positively charged hydrogen atoms, one negatively charged hydrogen atom, and one Lennard-Jones site for oxygen. A total of 1000 water molecules were used for the simulation with dimensions of 25 X 25 Å. A periodic boundary condition was used for all six boundaries. The simulation shows that the hydrogen bond in water has a significant effect on the molecules at the interface and may reduce the evaporation coefficient. This study showed that a combination of MDS and the classic Schrage model might provide a methodology to find the evaporation coefficient at different temperatures and pressures.

Figure 1(a) shows the density distribution for several different temperatures in the normal direction to the interface at equilibrium. The figure shows three zones: liquid, interface, and vapor. The density profile may be fitted to a hyperbolic-tangent function, as shown in Fig. 1 (b): $\rho \left( z \right)=\frac{1}{2}\left( {{\rho }_{\ell }}+{{\rho }_{v}} \right)-\frac{1}{2}\left( {{\rho }_{\ell }}-{{\rho }_{v}} \right)\times \tanh \left[ \frac{z-{{z}_{0}}}{d} \right]$ (2)

where ${{\rho }_{\ell }}$ and ρv are the liquid and vapor densities, respectively; z0 is the position of the Gibbs dividing surface d is a measure of the interfacial thickness and is related to the interfacial thickness, δI, by

 δI = 2.1792d (3)

The interfacial thickness ranges from 4 to 8 Å and increases with temperature by

 δI = − 10.95 + 0.04068T (4)

where δI and T are in Å and K, respectively. The interfacial thickness for water is much smaller than that for argon, which is ranged from 16 to 20Å.

Figure 1 (a) Density profile at various different temperatures, (b) density profile at 423.15K and its fitted curve of hyperbolic tangent form (Yang and Pan, 2005; Reproduced by permission of Routledge/Taylor & Francis Group, LLC).

Molecular dynamics simulation can apply not only to liquid-vapor interfaces, but also to solid-liquid interfaces, to provide detailed information at the microscopic level. For example, conventional disjoining pressure theory, presented in Section 2.6.4, provides a base for determining the effects of wall-fluid attractive forces but ignores many details of the variation in thin liquid film. Molecular dynamics simulation was used by Carey and Wemhoff (2005) to examine the detailed features of the structure and thermophysics of thin liquid films on solid surfaces. The film features and equilibrium vapor pressure were determined for films of varying thickness and were compared to the conventional disjoining pressure results.

The molecular dynamics simulation used combines traditional NVE-type MD simulation with a stochastic boundary condition to represent the equilibrium of pressure. This method allows monatomic fluids that interact via the Lennard-Jones potential ${{\phi }_{LJ}}\left( r \right)=4{{\varepsilon }_{LJ}}\left[ {{\left( \frac{{{\sigma }_{LJ}}}{r} \right)}^{12}}-{{\left( \frac{{{\sigma }_{LJ}}}{r} \right)}^{6}} \right]$ (5)

where σLJ and εLJ are the Lennard-Jones length and energy parameters and r is the distance between two molecules. The system obtains equilibrium by raising the temperature from an initially low value to the desired equilibrium temperature by velocity rescaling, so that the system kinetic energy follows the relation $\frac{1}{2}m\sum\limits_{i=1}^{N}{c_{i}^{2}=\frac{3}{2}N{{k}_{B}}T}$ (6)

where ci is the speed of molecule i and N is the number of molecules in the system. A flux boundary at z = L is used to adjust the pressure to an equilibrium value. For ideal gases, the vapor boundary molecular flux J is related to pressure p and temperature T by kinetic theory relation $J=\sqrt{\frac{{{k}_{B}}T}{2\pi m}}\left( \frac{p}{{{k}_{B}}T} \right)$ (7)

where m is the molecular mass of the injected species. Figure 2 shows the variation of pv,δ/psat with film thickness, which was predicted by the MDS for thickness less than 6 nm. Also shown is the variation predicted by the conventional disjoining pressure theory, presented in Section 2.6.4. It can be seen that the MDS results generally agree with tabulated bulk saturation pressure within about 10% for film thickness greater than 1.5 nm, with the exception of 4.9 nm.

The thermodynamic model only offers a prediction of the relationship between vapor pressure and film thickness at a given temperature. The molecular dynamics simulation, however, provides detailed information on the density and pressure variation in the thin liquid film and the variation of the equilibrium vapor pressure with temperature and thickness. Shi et al. (2005) used molecular dynamic simulations to analyze the contact angles of water droplets on a platinum surface over a range of temperature. To obtain high accuracy for the calculation of long-range forces, the P3M and Ewald summation methods were used. They found that the contact angle decays as the system temperature or the solid wall potential increases. In high temperature regions the contact angle decays very rapidly. The trends of the contact angle for water-platinum, argon-virtual solid wall and water-aluminum cases are similar.

## References

Carey, V.P., and Wemhoff, A. P., 2005, “Disjoining Pressure Effects in Ultra-Thin Liquid Films in Micropassages – Comparison of Thermodynamic Theory with Predictions of Molecular Dynamics Simulations,” IMECE2005-80234, Proceedings of 2005 ASME International Mechanical Engineering Congress and Exposition, Orlando, FL.

Faghri, A., and Zhang, Y., 2006, Transport Phenomena in Multiphase Systems, Elsevier, Burlington, MA.

Matsumoto, M., 1998, “Molecular Dynamics of Fluid Phase Change,” Fluid Phase Equilibria, Vol. 144, pp. 307-314.

Yasuoka, K., and Matsumoto, M., 1995, “Molecular simulation of evaporation and condensation. I. Self condensation and molecular exchange,” Proceedings of ASME/JSME Thermal Engineering Conference, pp. 459-464.

Shi, B., Sinha, S., Dhir, V.K., 2005, “Molecular Simulation of the Contact Angle of Water Droplet on a Platinum Surface,” IMECE2005-80136, Proceedings of 2005 ASME International Mechanical Engineering Congress and Exposition, Orlando, FL.

Wu, Y. W., and Pan, C., 2003, “A Molecular Dynamics Simulation of Bubble Nucleation in Homogenous Liquid under Heating with Constant Mean Negative Pressure,” Microscale Thermophysical Engineering, Vol. 7, pp. 137-151.

Yang, T.H., and Pan, C., 2005, “Molecular Dynamics Simulation of a Thin Water Layer Evaporation and Evaporation Coefficient,” International Journal of Heat and Mass Transfer, Vol. 48, pp. 3516-3526.