Size Effect on Transport Properties of Liquid Argon: A Molecular Dynamics Simulation Study

Green-Kubo Formulas and Molecular Dynamics Simulation

The diffusion coefficient is obtained through two routes: the Green-Kubo formula from velocity auto-correlation (VAC) function:

and the Einstein formula from mean square displacement (MSD):

The shear viscosity is given by a modified Green-Kubo formula for better statistical accuracy:4

where

with αβ = xy, xz, yx, yz, zx, and zy. Similarly the thermal conductivity is given by4

where is the heat flux:

with α = x, y, and z. The energy of molecule i is given by

where ϕ [rij(t)] is the potential energy between particles i and j at time t.

We have chosen argon systems of N = 1728, 6912, 13824, 27648 and 110592 to carry out MD simulations at T = 94.4 and ρ = 1.374 g/cc. The length of cubic simulation box of each system is given in Table 1. The usual Lennard-Jones (LJ) 12-6 potential for the interaction between liquid argon is used with the LJ parameters,5 σ = 0.3405 nm and ε/kB = 119.8 K, where kB is the Boltzmann constant. The interatomic potential is truncated at 1.0 nm, which is the cut-off distance used in many other simulations. Long-range corrections are applied to the energy, pressure, etc. due to the potential truncation.6 The equations of motion were solved a velocity Verlet algorithm7 with a time step of l0−14 second for NVT MD simulations with the determined volumes from the system density. The configurations of argon molecules were stored every time step for further analysis. The systems were fully equilibrated and the equilibrium properties are averaged over 10 blocks of 5,000 time steps (50 ps).

Table 1.aUsing other LJ parameters: σ = 0.3268 nm and ε/kB = 119.8 K. bAt 90 K and 1.374 g/cc.9 cObtained from Lagrange interpolation of experimental results at 94.4 K and saturated vapor pressure.12

 

RESULTS AND DISCUSSION

The mean square displacement (MSD) and the velocity auto-correlation (VAC) function of liquid argon show a normal behavior (not shown) as Rahman depicted them first time in 1964.8 The MSD of liquid argon shows a linear behavior within 3 ps with the initial dent within 0.3 ps,8 and the VAC of liquid argon decays rapidly to 0 within 0.5 ps, has a negative value due to the collision with the neighboring particle, and fluctuates near zero after 1.5 ps.8

Table 1 lists diffusion coefficients of the liquid argon for N = 1728, 6912, 13824, 27648 and 110592 at 94.4 K and 1.374 g/cc obtained from MSD’s using Eq. (1) and VAC’s using Eq. (2) which are in good agreement with the experimental measures (2.43×10−5 cm2/s at 90 K and 1.374 g/cc).9 As the number of argon molecules increases, D obtained from MSD approaches the experimental measure (2.393 2.427 → 2.429 → 2.422 → 2.435) and D obtained from VAC also shows a good accuracy (2.401 → 2.403 → 2.435 → 2.428 → 2.435). The previously reported results for D of liquid argon at the same state were 2.44,10 2.48,11 and 2.42×10−5 cm2/s4 by EMD, and 4.02×10−5 cm2/s by NEMD.10 The results using other LJ parameters [σ = 0.3268 nm and ε/kB = 119.8 K] for N = 13824 are much poorer than those using the original LJ parameters [σ = 0.3405 nm and ε/kB = 119.8K] as shown in Table 1. Obviously, the larger diffusion coefficients (D) are attributed to the smaller size LJ parameter (σ).

Stress auto-correlation (SAC) and heat-flux auto-correlation (HFAC) functions of the liquid argon for N = 1728, 13824, and 110592 at 94.4 K and 1.374 g/cc are plotted in Figs. 1 and 2. Both correlation functions are fluctuating near zero after 5 ps. In the inset of Fig. 1, we plot the detailed SAC functions in the very narrow y-axis around the zero correlation. For N = 1728, the fluctuation of the SAC function is very high, but it lowered with increasing number of argon molecules. The SAC function for N = 13824 is acceptable and that for N = 110592 is more perfect.

Figure 1.Stress auto-correlation functions SAC = of liquid argon at 94.4 K and 1.374 g/cc. The inset shows the detailed behavior of SAC functions.

Figure 2.Heat-flux auto-correlation functions HFAC = of liquid argon at 94.4 K and 1.374 g/cc. The inset shows the detailed behavior of HFAC functions.

Running integrals for η(t) of liquid argon for N = 1728, 13824, and 110592 at 94.4 K and 1.374 g/cc are plotted as a function of time in Fig. 3. All the running integrals for viscosity clearly show plateaus which signify that the corresponding SAC functions have decayed to zero and are fluctuating along the horizontal time axis. As shown in the inset of Fig. 1, all the SAC functions reach zero at about 5 ps and we report the shear viscosities for N = 1728, 6912, 13824, 27648, and 110592 at 94.4 K and 1.374 g/cc in Table 1 by averaging the running integrals for shear viscosity in Fig. 3 for 5~10 ps. As the number of argon molecules increases, η obtained from SAC decreases to the experimental and then increases slightly.

Figure 3.Running integrals for η(t) (mP) of liquid argon at 94.4 K and 1.374 g/cc.

The shear viscosities, η, obtained by MD simulations at 94.4 K and 1.374 g/cc overestimate the experimental measure for all the values of N. η for N = 13824 is closer to the experimental measure12 than those for larger N’s as shown in Table 1 and Fig. 3, but the result for N = 1728 is unreliable due to the high fluctuation of the SAC function as seen in the inset of Fig. 1. Increasing further in N makes the result slightly worse. η = 2.895 mp for N = 13824 is better than the previously reported results for η of liquid argon at the same state, 3.0811 and 3.04 mp4 by EMD, but not for 2.29 mp10 by EMD and 2.13 mp by NEMD.10 η obtained for N = 13824 using other LJ parameters [σ = 0.3268 nm and ε/kB = 119.8 K] is much better than that using the original LJ parameters [σ = 0.3405 nm and ε/kB= 119.8 K] as shown in Table 1.

The situation for HFAC is very similar to that for SAC. In the inset of Fig. 2, the detailed HFAC function in the very narrow y-axis around the zero correlation shows the high fluctuation of the HFAC function for N = 1728, but it lowered with increasing number of argon molecules. The HFAC function for N = 13824 is better than that for N = 1728, and that for N = 110592 is more reliable.

Running integrals for λ(t) of liquid argon for N = 1728, 13824, and 110592 at 94.4 K and 1.374 g/cc are plotted as a function of time in Fig. 4. All the running integrals for thermal conductivity clearly show plateaus which signify that the corresponding HFAC functions have decayed to zero and are fluctuating along the horizontal time axis. As shown in the inset of Fig. 2, all the HFAC functions reach zero at about 5 ps and we report the thermal conductivities for N = 1728, 6912, 13824, 27648, and 110592 at 94.4 K and 1.374 g/cc in Table 1 by averaging the running integrals for thermal conductivity in Fig. 3 for 5~10 ps.

Figure 4.Running integrals for λ(t) (cal/cm·s·K) of liquid argon at 94.4 K and 1.374 g/cc.

The thermal conductivities, λ, obtained by MD simulations at 94.4 K and 1.374 g/cc underestimate the experimental measure12 for all the values of N. η for N = 13824 is the best among all the values of N in Table 1 and Fig. 4. λ for N = 1728 is unreliable due to the high fluctuation of the HFAC function as seen in the inset of Fig. 2. λ for N = 13824 is better than λ for N = 110592 and increasing N also makes the result slightly worse. The previously reported results for λ of liquid argon at the same state were 1.84,10 3.05,11 and 1.55×10−4 cal/cm·s·K4 by EMD, and 2.31×10−4 cal/cm·s·K by NEMD.10 λ obtained for N = 13824 using other LJ parameters [σ = 0.3268 nm and ε/kB = 119.8 K] is much better than that using the original LJ parameters [σ = 0.3405 nm and ε/kB = 119.8 K] as the same for the case of η.

In summary, we have carried out a series of equilibrium molecular dynamics (EMD) simulations of liquid argon at 94.4 K and 1.374 g/cc for the calculation of transport properties as a function of the number of argon molecules (N). While the diffusion coefficients (D) of gaseous argon approach the experimental values with increasing N, the viscosities (η) and thermal conductivities (λ) obtained for N = 1728 are unreliable due to the high fluctuation of the time correlation functions, but those for N = 13824 are rather acceptable. Increasing N further to 110592 brings the EMD results a slightly worse for η and λ. The EMD results for η overestimate and those for λ underestimate the experimental measurements,10 respectively, and it is not expected that further increase in N would give results closer to the experimental measurements. The use of the smaller size LJ parameter (σ) could improve the results for η and λ of liquid argon at 94.4 K and 1.374 g/cc but not for D. Therefore a systematic EMD simulation using various values for LJ parameters (σ and ε) is currently under study.