Molecular Dynamics Simulations of the Local Structures and Transport

Article pubs.acs.org/JPCB

Molecular Dynamics Simulations of the Local Structures and Transport Coefficients of Molten Alkali Chlorides Jia Wang,†,‡ Ze Sun,*,‡ Guimin Lu,*,†,‡ and Jianguo Yu‡ †

School of Mechanical and Power Engineering and ‡National Engineering Research Center for Integrated Utilization of Salt Lake Resource, East China University of Science and Technology, Shanghai 200237, China

ABSTRACT: Systematic results from molecular dynamics simulations of molten alkali chlorides (ACl) serials are presented in detail in this paper. The effects of temperature and cationic size on the structures and transport properties of molten salts have been investigated and analyzed. The local structures of molten ACl have been studied via the analysis of radial distribution functions and angular distribution functions. The coordination number of ACl decreases when ACl melts from solid and increases as cationic radius increases. Molten LiCl takes a distorted tetrahedral complex as the microconfiguration, while other melts have the tendency to keep the original local structure of the corresponding crystal. Temperature has no significant effect on the local structures of molten ACls. The results also show that the Tosi−Fumi potential predicts positive temperature dependences for self-diffusion coefficients and ionic conductivity, and negative temperature dependences for both viscosity and thermal conductivity of molten ACls. Ionic diffusivity decreases as cationic radius increases from LiCl to CsCl. The simulation results are in agreement with the experimental data available in the literature.

1. INTRODUCTION Inorganic molten salts are liquid mixtures of ionic species at high temperatures. Due to their chemical and high temperature stability and other excellent performances, inorganic molten salts play an important role in many technological applications, such as molten salt reactors,1 coolants in nuclear reaction processes,2 electrolytes for metal production,3 the pyrochemical treatment of nuclear waste,4 and electricity storage devices.5 Revealing the structures and properties and understanding the relationships between them of high temperature inorganic molten salts have been interests of many researchers due to their wide applications. These are difficult tasks through experimental methods, given the extremely high temperature and corrosive conditions. It is not surprising, therefore, that there are substantial disparities in the reported experimental results for these materials. These disparities are normally attributed to convective and radiative effects in some experimental methods. It is desirable to accurately predict these properties over a wide range of conditions through experiments, while high temperature experiments are laborious and inconvenient. Molecular dynamics (MD) simulation, which is performed on the atomic level, can avoid the experimental problems referred to previously and provide an alternative to obtain the properties of molten salts. The MD method dates © 2014 American Chemical Society

from several decades ago and has successfully predicted the local structures and transport properties (such as diffusion coefficients, viscosity, ionic conductivity and so on) of not only many pure inorganic melts systems but also their mixtures.6−13 Alkali halides have been studied extensively via molecular simulations. Tosi and Fumi14 developed potential parameters for 17 alkali halides long ago by reproducing the properties of crystals. Then the thermodynamic and structural properties of molten alkali halide salts were primarily simulated with these potentials by molecular dynamics simulations. Galamba and coworkers15,16 computed the viscosity and thermal conductivity of molten KCl and NaCl with equilibrium molecular dynamics (EMD) simulations and nonequilibrium molecular dynamics (NEMD) simulations. Nevins and Spera17 examined the calculation conditions to obtain a good compromise between calculation time and viscosity quality when calculating the viscosity of molten NaCl with the EMD method. Ohtori and co-worker18 constructed the ab initio polarizable potentials for the alkali chlorides and then calculated the thermal conductivity of molten LiCl, NaCl, and KCl. The transport properties of Received: May 22, 2014 Revised: August 7, 2014 Published: August 8, 2014 10196

dx.doi.org/10.1021/jp5050332 | J. Phys. Chem. B 2014, 118, 10196−10206

The Journal of Physical Chemistry B

Article

Table 1. Parameters for the Potentials of the Tosi−Fumi Modela σ = σi + σj/Å

a

C/10−79 J·m6

ρ/Å

system

++

+−

−−

LiCl NaCl KCl RbCl CsCl

1.632 2.340 2.926 3.174 3.44

2.401 2.755 3.048 3.172 3.505

3.170 3.170 3.170 3.170 3.170

0.3425 0.317 0.3367 0.3185 0.3338

D/10−99 J·m8

++

+−

−−

++

+−

−−

0.073 1.68 24.3 59.4 152.0

2.0 11.2 48.0 79.0 129.0

111.0 116 124.5 130.0 129.0

0.03 0.8 24.0 82.0 278.0

2.4 13.9 73.0 134.0 250.0

223.0 233.0 250.0 260.0 260.0

b = 0.338 × 10−19 J, A+ + = 1.25, A+ − = 1.00, A− − = 0.75, except for LiCl, where A+ + = 2.00, A+ − = 1.375, A− − = 0.75.

Table 2. Densities for Molten ACl Salts at Different Temperatures (Units: g/cm3) LiCl a

NMT 1100 K 1200 K 1300 K 1400 K 1500 K

NaCl

KCl

RbCl

CsCl

ρsim

ρexpt

ρsim

ρexpt

ρsim

ρexpt

ρsim

ρexpt

ρsim

ρexpt

1.4076 1.3181 1.2793 1.2381 1.1962 1.1431

1.4946 1.4082 1.3650 1.3218 1.2787 1.2355

1.455 1.444 1.393 1.340 1.296 1.244

1.553 1.542 1.4878 1.4335 1.3793 1.325

1.3749 1.3548 1.2988 1.2446 1.1889 1.1345

1.5219 1.4945 1.4362 1.3779 1.3184 1.2609

2.0611 1.9609 1.8873 1.8008 1.7238 1.6485

2.2349 2.1498 2.0646 1.9785 1.8943 1.8092

2.6451 2.4690 2.3796 2.2573 2.1656 2.0466

2.8067 2.6287 2.5239 2.4192 2.3144 2.2097

temp

a

NMT in this paper means Near Melting Temperatures, i.e. 900 K for LiCl, 1080 K for NaCl, 1053 K for KCl, 1000 K for RbCl, and 900 K for CsCl, respectively.

2.2. Simulation Details. The LAMMPS code was used to perform the simulations. All the simulations reported here were performed for a cubic box consisting of 432 ions (216 cations and 216 anions, respectively) for CsCl and 512 ions for LiCl, NaCl, KCl, and RbCl (256 cations and 256 anions, respectively). It has been demonstrated that further increasing the particle numbers has no significant improvement on the results.17,23 The initial configurations were the crystal structure, i.e. fcc for the former four systems and sc for the CsCl system. The periodic boundary conditions method was employed to keep the particle numbers constant and to eliminate the boundary effect. The interionic interactions were truncated at rc = L/2. Here, L meant the ideal box length for the systems after relaxation and was primarily estimated by the experimental density24 at the corresponding temperature. Long-range interactions were handled using the Ewald summation method to eliminate the truncation errors, which means that the interactions within the cutoff distance are computed directly, and interactions outside that distance are computed in reciprocal space. In the Ewald sum, the precision was set equal to 1.0 × 10−5. The initial velocities were randomly assigned and obeyed the Gaussian distribution. The Verlet algorithm25 was used to solve Newton’s equations of motion with a time step Δt = 5 fs. For all the simulations, the systems were first melted at 3000 K and then cooled to the desired temperature; after that the systems were equilibrated in an isothermal−isobaric ensemble (NPT) with a pressure fixed at 0 GPa around the desired temperature using versions of the Nose−Hoover thermostat and barostat.26,27 The production runs were performed in a canonical ensemble (NVT) with the Nose−Hoover thermostat method28,29 at the equilibrated cell volume. The damping parameters which determined how rapidly the temperature and pressure were relaxed were set to be 0.5 ps, being equal to 100 time steps. To control initial oscillations, Nose−Hoover chains were used for the thermostat and barostat, and the chain numbers were 3 in this work.30 The NVT stages lasted up to more than 5 ns for each system in order to ensure good statistics for the computation of the transport coefficients.

molten alkali halide mixtures have also been studied through molecular dynamic simulations, such as those of molten LiF/ KF and LiCl/KCl.19 Although many studies on typical molten alkali halides have been reported, there are few systematic comprehensive investigations about them. In this work, the microstructures and macro transport properties of single molten alkali chlorides ACl (LiCl, NaCl, KCl, RbCl, and CsCl) were studied fully. The effects of cation size on the structure are analyzed, and transport properties are uncovered.

2. METHODOLOGY 2.1. Potential Model. For the interaction between different ions, the following Born−Mayer−Huggins potential was used: Uij(r ) =

⎛ σi + σj − r ⎞ Cij Dij + Aij b exp⎜ ⎟− 6 − 8 ρ r ⎝ ⎠ r r

qiqj

(1)

where the first term describes the electrostatic interactions between ions, qi is the ionic charge, and formal charge (+1 for alkali ions and −1 for chloride) is used; the second term represents Born−Huggins exponential short-range repulsion due to overlap of electron clouds, Aij is the Pauling factor,20 being equal to 1 + (zi/ni) + (zj/nj), in which Z denotes the charge of the ion and n denotes the number of outer electrons, b is a constant, σ is the crystal ionic “radii”, while ρ is the hardness parameter, and the last two terms correspond to the dipole−dipole and dipole−quadrupole dispersion interactions, in which Cij and Dij are the dispersion parameters. The values of short-range repulsive parameters and dispersion parameters used in the present calculation are listed in Table 1. The former four sets are given by Tosi and Fumi,14 and they were determined to give the best agreement with density and compressibility data for the 17 alkali metal halides. The repulsive values for CsCl are given by Dixson and Sangster.21 The dispersion constants Cij and Dij for all the five systems are determined by Mayer.22 10197



Article

Here, the “bond” does not mean a chemical bond, but a near neighbor atom pair, for example, an Li−Cl pair. The angular distribution function can be considered as a three body correlation function and can be calculated according to the position coordinates of atoms stored during the calculation. First, for a centered atom A, its nearest neighbor atoms Cl should be determined through the rmin of the corresponding ACl RDF. Then, taking the centering atom i as vertex and two arbitrary coordinated atoms j and k as end points, the bond angle is calculated through the following equation and finally the angle distribution probability is analyzed to infer the atomic configuration types existing in the system.

During the production stages, the radial distribution functions (RDF) of different ion pairs and angular distribution functions (ADF) were calculated to reveal the structure features of different melts, the mean square displacements (MSD) of different ions were accumulated to compute diffusion coefficients, and the autocorrelation functions (ACF) were built and integrated to give the values of other transport coefficients with 10 million steps. From initial NPT simulations, the liquid densities could be determined by the following equation:

ρ=

NM VENA

(2)

where N is the particle number, M is the molar mass, NA denotes Avogadro’s constant, and VE is the equilibrated volume of the simulation cell at the given temperature in the NPT ensemble simulations. The averaged densities obtained here are listed in Table 2 and compared with experimental values;24 this work underestimated densities by about 5−10%, primarily indicating that the Born−Mayer−Huggins potential could describe molten ACl systems accurately. The densities of LiF−BeF2 obtained by Jabes et al.31 with the transferable rigidion mode potential were approximately 17−20% lower than the experimental values while polarizable ion models (PIM) gave densities within 5% of the experimental value. The induced dipole moments resulting from anion polarization screen the repulsion Coulombic interaction between cations and result in a more compact pack between cations. Further calculations taking polarization effects into account should be able to improve the accuracy of the results. 2.3. Evaluated Properties. Molecular dynamics simulations were used to compute the static and dynamic properties detailed below.25,32 2.3.1. Radial Distribution Functions. The mainly used function for describing a liquid structure is the radial distribution function (RDF), defined as gαβ(r) for an ion pair of species α and β gαβ(r ) =

1 ⎡ dNαβ(r ) ⎤ ⎢ ⎥ 4πρβ r 2 ⎣ dr ⎦

θjik =

∫0

rmin

gαβ(r )r 2dr

cos ⎜ ⎝

+ rik2 − r jk2 ⎞ ⎟ ⎟ 2rijrik ⎠

(5)

In addition to obtaining the structural information, MD is also a powerful tool to determine both individual quantities, such as diffusion coefficients of various species, and collective quantities, such as the thermal conductivity, shear viscosity, ionic conductivity, and so on, for which it is difficult to obtain accurate values from experiments at high temperatures. It is important to run the simulations for enough time to obtain good statistics sufficiently. 2.3.3. Diffusion Coefficient. The slope of the mean square displacement (MSD) versus time in the linear long time limit region is related to the diffusion coefficients of the diffusing atoms. The diffusion coefficients Dα of ions can be calculated from the mean square displacement through the Einstein expression:32 Dα = lim t →∞

|δri(t )|2 6t

(6)

δri(t) is the displacement of a typical α-type ion in time t, and the angular brackets denote ensemble average over all the ions of species α; thus, the diffusion coefficient is an individual property. Then the diffusion coefficients are extracted from the slope of the MSD at long times. 2.3.4. Shear Viscosity. The shear viscosity η is a measure of the resistance of a fluid to being deformed by shear stress. Molecular dynamics can be used to predict viscosity using either equilibrium or nonequilibrium methods. Here, it is calculated from the time integral of the off-diagonal stress tensor autocorrelation function based on the Green−Kubo (GK) formula16,32 in a steady-state equilibrated molecular simulation (EMD):

(3)

where ρβ is the number density of species β and Nαβ(r) is the mean number of β-type ions lying in a sphere of radius r centered on an α-type ion. The coordination number (Nαβ), one of the most important parameters to characterize the local structure, is the mean number of β-type ions lying in a sphere of radius rmin centered on an α-type ion; rmin is the position of the first peak valley of the RDF. Nαβ can be calculated from the RDF through the following equation: Nαβ = 4πρβ

⎛ r2 −1⎜ ij

η=

1 kBTV

∫0

∞

⟨Sxy(0)Sxy(t )⟩dt

(7)

where kB is the Boltzmann constant, T is temperature, V is the simulation cell volume, and Sxy is the xy-component of the stress tensor. Sxy is defined as

(4)

The ionic equilibrium distance corresponds to the position of the first maximum of the RDF. 2.3.2. Angular Distribution Functions. The RDF plays a very important role in the research of melt structure, but it just counts the probability of atoms pairing, containing no orientation information among atoms. To further describe the local structure of molten alkali chloride salts, aside from the RDF, the angular distribution function (ADF)33 of A−Cl−A, which can reflect the bond orientation, has also been analyzed.

Sxy =

N

⎡

i=1

⎣

∑ ⎢⎢mivxivyi +

1 2

⎤

∑ xijf y (rij)⎥⎥ j≠i

⎦

(8)

where mi is the mass of ion i, vxi and vyi are the x-component and y-component of vi (the velocity of ion (i), respectively), xij is the x-component of rij = ri − rj, and f y(rij) is the y-component of the force fij on ion i due to ion j. Each of the three independent off-diagonal components of the stress tensor (i.e., 10198



Article

Figure 1. Radial distribution functions for molten ACl at NMTs.

⎡N 1⎢ 1 JE = ∑ Ei vi + ⎢ V ⎣ i=1 2

Sxy, Sxz, and Syz) provides an independent estimate of shear viscosity, and the averaged value is taken as the final viscosity value. 2.3.5. Thermal Conductivity. Using the GK relationship, the thermal conductivity κ can be calculated from the time integral of the energy flux autocorrelation function through the EMD method:15,32 V κ= kBT 2 V = 3kBT 2

∫0

Ei =

⟨JxE (0)JxE (t)⟩dt

1 [mivi2 + 2

j≠i

∑ Uij(rij)] N

(11)

where mi is the mass of ion i, vi is the velocity of the ion, Uij(rij) is the pair potential between particles i and j, rij is the position vector between particles i and j, and fij is the force on ion i due to ion j. Each of the three independent components (i.e., JxE, JyE, JzE) provides an independent estimate of thermal conductivity, and the averaged value is taken as the final thermal conductivity value. 2.3.6. Ionic Conductivity. Using the GK relationship,32 the ionic conductivity λ can be calculated from the time integral of

∞

⟨JE (0) ·JE (t )⟩dt

(10)

and the energy per particle Ei is defined as

∞

∫0

⎤ ∑ (rijfij)·vi⎥⎥ ⎦ j≠i N

(9)

where kB is the Boltzmann constant, T is temperature, V is the simulation cell volume, JE is the energy flux, and JxE is the xcomponent of the energy flux vector. The energy flux vector JE is defined as 10199



Article

Table 3. Structural Parameters for Molten ACl Systems rmax (Å)a system

ion pair

present

rigid

LiCl

Li−Li Li−Cl Cl−Cl Na−Na Na−Cl Cl−Cl K−K K−Cl Cl−Cl Rb−Rb Rb−Cl Cl−Cl Cs−Cs Cs−Cl Cl−Cl

3.87 2.25 3.76 4.17 2.61 4.05 4.58 2.96 4.58 4.67 3.12 4.74 4.76 3.30 4.90

3.7 2.2 3.7 4.0 2.6 4.0 4.5 2.9 4.4 4.5 3.1 4.7 4.8 3.3 4.8

NaCl

KCl

RbCl

CsCl

b

hmax (Å)a expt

b

2.47c 3.9 2.7 3.9 4.8 3.1 4.8 4.8 3.2 4.8 3.9 3.4 3.9

present

rigid

1.86 3.95 2.13 1.74 3.95 1.83 1.76 3.94 1.72 1.81 4.22 1.76 1.91 4.13 1.75

1.8 3.6 2.0 1.9 4.0 2.0 1.8 3.7 1.8 1.9 4.2 1.8 1.8 4.1 1.8

rmin (Å)a

b

expt

b

1.7 3.3 1.9 2.3 3.0 1.8 1.6 4.4 1.7 1.8 3.3 1.8

present

rigidb

5.28 3.49 5.28 6.09 4.05 6.09 6.86 4.65 6.79 7.22 4.80 7.22 7.55 5.11 7.55

5.4 3.5 5.4 6.0 4.1 6.1 6.8 4.7 6.9 7.1 4.9 7.2 7.5 5.0 7.5

exptb

Nαβa

6.1 3.9 5.8 7.0 4.9 7.0 7.3 5.2 7.3 4.3 5.4 4.3

11.81 4.05 11.86 13.59 4.76 13.54 14.57 5.18 14.28 15.66 5.41 15.71 16.26 5.70 16.34

rmax = the position of the first peak; hmax = the height of the first peak; rmin = the position the first minimum; Nαβ = coordination number. bRigid: Results from rigid-ion simulations in ref 21; expt: experimental data referenced in ref 21. cReference 35. a

the charge flux autocorrelation function through the EMD method: λ=

1 3VkBT

∫0

Table 4. Diffusion Coefficients for Molten ACl Salts at Different Temperatures (in units of 10−9 m2/s)

∞

⟨JZ (0) ·JZ (t )⟩dt

(12)

where kB is the Boltzmann constant, T is temperature, V is the simulation cell volume, and the charge flux vector JZ(t) is defined as

system

ion

NMT

1100 K

1200 K

1300 K

1400 K

1500 K

LiCl

Li Cl Na Cl K Cl Rb Cl Cs Cl

8.39 5.86 8.21 7.39 7.25 7.34 3.96 4.80 2.86 3.31

14.72 10.74 9.36 8.14 8.27 8.49 5.68 5.86 5.13 6.02

17.51 12.93 10.74 10.25 11.34 11.62 7.79 8.73 7.07 7.32

19.48 14.27 12.70 11.31 13.10 14.22 9.87 10.36 8.66 10.55

23.68 19.72 15.33 13.67 14.23 15.20 11.96 12.53 9.80 11.52

27.61 21.91 18.26 17.99 19.83 18.92 14.81 15.85 12.65 15.57

NaCl KCl

n

JZ (t ) =

∑ zie vi(t ) i=1

RbCl

(13) CsCl

where zie and vi are the charge and velocity of atom i, respectively. Each of the three independent components (i.e., JxZ, JyZ, JzZ) of the charge flux vector provides an independent estimate of ionic conductivity, and the averaged value is taken as the final ionic conductivity value.

plotted in Figure 1. All RDFs show the following common features: they all have a relatively high first peak, the amplitude of the fluctuations after the first peak decreases with distance, and they approach unity at larger distances. The RDFs of all pairs approach unity beyond the very weak third peak, i.e. r = 8 Å, indicating the short-range order and long-range disorder of the liquid state. On the other hand, the RDFs of both like and unlike ions still display oscillations beyond 13 Å, which will be analyzed below. In particular, each ACl function is characterized by a sharp first peak followed by a well-defined minimum. For the A−A and Cl−Cl pairs, the first peak is much broader and its minimum is higher, which suggests a stronger coordination structure between unlike ions. All the gAA(r) and gClCl(r) values are antiphase to the corresponding gACl(r); this phenomenon reflects an alternating arrangement of cations and anions, being consistent with the charge ordering arguments.34 As mentioned above, the RDFs of both like and unlike ions still display oscillations beyond 13 Å, larger than the distance at which the all pairs RDF (without distinguishing charge) approaches unity, indicating that charge ordering exists in a longer range than short-range order of nonionic liquids. As shown in Figure 1, due to the existence of short-range repulsion, there are no particles distributed within a region centered at any type of ion. The unlike ions locate at the first shell, the like ions locate at the second shell, and the shapes of

3. RESULTS AND DISCUSSIONS 3.1. Short-Range Structure. 3.1.1. Radial Distribution Functions. To study the local structures of molten alkali chlorides, the simulated radial distribution functions for ACl systems near melting points, corresponding to all pairs, unlike ion pairs (A−Cl), and like ion pairs (A−A and Cl−Cl), are

Figure 2. Angle distribution functions for molten ACl. 10200



Article

Figure 3. Self-diffusion coefficients for molten NaCl, KCl, RbCl, and CsCl (straight line: experimental data, square symbols: simulated results).

However, there are significant differences between the simulated results and the experimental data. Introduction of ionic polarization, which provides a mechanism whereby the Coulombic repulsion between the like-ions is partially screened, may improve the simulated results.21 In consideration of the fact that the coordination number calculation is affected by the upper limit of integral rmin, the coordination numbers of ACl obtained here are in close agreement with those obtained by Baranyai et al.,36 i.e. 4.0, 5.0, 5.8, 5.7, and 5.8 from LiCl to CsCl. All the studies are in agreement with the fact that the coordination number of the first shell for unlike ions decreases from solid to liquid, being smaller than 6, and the Nαβ for like ions increases from solid to liquid, being larger than 12, except for LiCl. When changing the alkali metal ions, the ratio between cationic and anionic radii, σ+/σ−, increases from 0.515 for LiCl to 1.085 for CsCl; the ionic equilibrium distances and coordination numbers for both like ion pairs and unlike ion pairs increase from lithium to cesium, resulting in a gradually increasing coordination structure for ACl serials. The larger is the cation ion, the larger are the coordination numbers of that ion with respect to Cl and to itself and also the corresponding Cl−Cl pair. It is obvious that the local structure is strongly affected by the cationic size, which is mainly reflected in the short-range term. A conclusion can be drawn that the shortrange repulsion resulting from overlap of clouds dominates the microcoordination structure. The RDFs at different temperatures from 1100 to 1500 K with an interval of 100 K have also been calculated to study the temperature dependence of short-range order. The results show that temperature has no significant impact on the local structures of molten ACl serials. As temperature increases, atomic arrangements become somewhat loose and the coordination numbers decrease slightly. 3.1.2. Angular Distribution Functions. Figure 2 shows the angular distribution functions for ACl melts near melting points, all showing a broad peak and a smearing phenomenon

Figure 4. Time evolution of the stress tensor ACF and the running integral for molten KCl at NMT.

all the gAA(r), gClCl(r), and gACl(r) for different systems are similar, respectively, except for a tiny difference in the peak height, suggesting that, for the ACl melts, the Coulomb forces dominate the interionic interactions, and each system has the same microdistribution trend. For each system, there is an overall similarity between gAA(r) and gClCl(r), especially for molten KCl and RbCl, where gAA(r) and gClCl(r) almost fully match at all distances. This phenomenon can be explained in view of the predominance of the repulsive Coulombic interactions in determining the distribution between similarly charged ions. However, for molten LiCl and NaCl, these are still clearly distinct. The height of the gClCl(r) first peak is higher than that for gAA(r), and for molten CsCl, the height of gCsCs(r) is higher, which are ascribed to the ionic size effects. The quantitatively structural information for molten ACl salts and their comparison with previous calculations and experimental data21,35 are summarized in Table 3. The structural results are in good agreement with previous rigid-ion simulations, confirming the validity of the calculations. 10201



Article

Figure 5. Viscosities of molten ACl at different temperatures (line: experimental data, square symbols: simulated results).

Figure 6. Time evolution of the energy flux ACFs for molten ACl at NMTs.

Figure 7. Thermal conductivities of molten NaCl, KCl, and RbCl at different temperatures (line: experimental data, square symbols: simulated results).

which indicate that the coordination bonds orient randomly. Here, “bond” does not mean a real chemical bond but just an atom pair within some distance. The Cl−A−Cl bond angles mainly distribute between 60° and 120°, and the peak values from LiCl to CsCl are 98.18°, 87.27°, 83.64°, 87.27°, and 83.64°, respectively. This shows that, for molten ACl systems, though losing long-range order, at the microlevel there is still some tendency to keep the original crystal local structures

where the bond angle is 90°37 (for the CsCl crystal, it will transform to a NaCl-type crystal structure at 778 K and the bond angle is also 90°). However, the local structures in molten ACl systems have a certain degree of tetrahedron. It must be noted that the Nαβ values calculated in 3.1.1 are statistically averaged values. In fact, for ACl melts, there coexist 3-fold, 410202



Article

= 1.633 for a perfect tetrahedron. It can be rationally speculated that the microconfiguration exists in molten ACl is an anomalous, distorted tetrahedron. The distortion may be because the large Cl ions arrange around small Li ions. The distance between Cl ions is close, and then the interaction is strong. As the σ−/σ+ value decreases from Li to Cs, the cation first coordination shell can accommodate more anions, and the tendency to form an octahedron becomes stronger. However, the complicated interaction between ions leads to a distorted octahedron configuration compared with crystal local 6-fold structures. 3.2. Self-Diffusion Coefficients. The self-diffusion coefficients are calculated from the slope of the mean square displacement in the linear region for each melt. To study the temperature dependence of ionic diffusivity, the self-diffusion coefficients of ACl molten salts at different temperatures are calculated, i.e. NMT, 1100, 1200, 1300, 1400, and 1500 K, and the results are listed in Table 4. In Figure 3, the simulated data are compared with experiment results40,41 for NaCl, KCl, RbCl, and CsCl. The comparisons suggest some difference for the self-diffusion coefficients, especially for molten NaCl and CsCl. The simulation underestimates the self-diffusion coefficients of both ions in molten CsCl and of Na+ diffusivity in molten

Figure 8. Time evolution of the charge flux ACF and the running integral for molten KCl at NMT.

fold, 5-fold, and 6-fold structures and so on, with different proportions, as is the same in the case of some other molten salts.38,39 For molten LiCl the peak of the ADF is 98.18°, combining with the 4.05 coordination number and the fact that the bond angle for a standard tetrahedron is 109.47°. In addition, the ratio between the Cl−Cl and the Li−Cl bond lengths dCl−Cl/ dLi−Cl is calculated to be 1.671 against the ideal value of (8/3)1/2

Figure 9. Ionic conductivities of molten ACl at different temperatures (line: experimental data, square symbols: simulated results). 10203



Article

The shear viscosities of molten ACl salts have been calculated at different temperatures, i.e. near the melting point and from 1100 to 1500 K with an interval of 100 K. The simulated viscosity results for ACl melts at different temperatures are shown in Figure 5, and they are in close agreement with available experimental data,48−51 with the error within 10%, except for molten CsCl, where the simulated results are obviously high. Research about molten RbCl and CsCl is rare from both the experimental and computational respects; therefore, it is difficult to validate the results. However, the simulated results here can be a supplement to the current database providing the variations of basic properties with temperature or be a reference for later research. For all systems, the shear viscosities decrease as the temperature increases, in contrast to the changing rules of self-diffusion coefficients. Among the five alkali chlorides, LiCl melt has the lowest viscosity, and the effects of LiCl on the mobility and fluidity of KCl molten salts have been studied19 extensively due to their wide application. However, the molten NaCl system has a relatively high viscosity, and the addition of LiCl may change the performance of molten NaCl, which needs to be further researched. 3.4. Thermal Conductivity. It is difficult to carry out experiments, and no empirical estimation technique has been proven to be useful to estimate thermal conductivity. The thermal conductivity is the least known property of molten salts, so it is significant to simulate the thermal transport coefficients accurately through molecular dynamic simulations. Figure 6 shows the time evolution of the normalized energy flux ACFs for molten ACl at NMTs. The normalized ACFs of the three components all decay to zero quickly, within 1 ps. Negative correlations within a certain region are observed for molten ACl, except for KCl, and the phenomena are more remarkable for systems in which the differences between the cationic and anionic radii are large. The running integrals reach a plateau after some time, which varies from 1.0 to 8.0 ps for different systems. During thermal conductivity calculations, for energy flux ACFs calculations, a time of 5 ps is adopted for molten LiCl and KCl, while 10 ps is used for molten NaCl, RbCl, and CsCl. The thermal conductivities of molten ACl salts have been calculated at different temperatures, i.e. near the melting point and from 1100 to 1500 K with an interval of 100 K. Figure 7 compares the simulated thermal conductivity values with those available experimental data52 for molten NaCl, KCl, and RbCl. The thermal conductivities calculated here are shown not to agree so well with the experimental data as in shear viscosity, being consistent with those simulated by Ohtori at al.53 For molten NaCl and KCl, the calculated results are somewhat satisfactory, showing errors within 10% compared with experimental results. However, for molten RbCl, the results were about 60−70% larger than experimental data. For systems where the difference between cationic and anionic size is large, such as molten LiCl and molten CsCl, the calculated results of the thermal coefficients are about 2−3 times experimental values. The unsatisfactory results may be due to the omission of polarization and induced-dipole contributions18 and the inevitable temperature fluctuation during MD calculations. However, the results are in agreement with corresponding experimental data in showing negative temperature dependence. Further research including the polarization effect could be carried out to improve the accuracy of the calculated results and

NaCl. However, due to the difficulty and inconvenience of achieving high temperature conditions, experimental measurements of diffusion coefficients may also have large error. In fact, there is a range of values in the experimental self-diffusion coefficients, with different experimental studies showing up to 1 order of magnitude difference. The results of molten LiCl at 900 K herein are consistent with those calculated by Okazaki,42 D(Li+) = 9.39 × 10−9 m2/s, D(Cl−) = 6.52 × 10−9 m2/s at 956 K. For Cl diffusivity in molten NaCl, the result 7.39 × 10−9 m2/s at 1080 K is close to previous calculated values: 7.5 × 10−9 m2/s at 1073 K.43 The results for LiCl melt at 1100 K and KCl melts at 1100 and 1200 K are in close agreement with those obtained by Bengston44 using the FPMD (first-principles molecular dynamics) method. That is, the self-diffusion coefficient of Li+ in molten LiCl at 1096 K is 14.04 × 10−9 m2/s, while the K+ and Cl− diffusivities in molten KCl are 8.39 × 10−9 m2/s and 8.6 × 10−9 m2/s, respectively, at 1096 K and 12.11 × 10−9 m2/s and 11.38 × 10−9 m2/s at 1200 K. All of the above confirms the accuracy of the calculations, indicating that such simulations can fulfill a significant role in augmenting existing experimental work. The ACl melts serials have relatively high diffusion coefficients, so they have excellent ionic conductivity and fluidity. The self-diffusion coefficients for all ions in different systems increase as temperature increases; that is, the ionic diffusivity increases as the temperature increases. In a single melt, the rule is D(Li+) > D(Cl−), D(Na+) > D(Cl−), D(K+) < D(Cl−), D(Rb+) < D(Cl−), D(Cs+) < D(Cl−). At the same temperature, D(Li+) > D(Na+)> D(K+) > D(Rb+) > D(Cs+), following the reverse order of the cationic radii. The mass and size of one ion significantly affect its dynamics.45,46 The ionic size effects on ionic dynamics in molten RbCl reported by Alcaraz and Trullas47 are also present in the other molten alkali chlorides; that is, the smaller the ion is, the larger diffusivity it has. For molten KCl, D(K+) is close to D(Cl−). This may be due to the similar radii and mass they have. Regarding the self-diffusion coefficients of the Rb+ and Cl− ions, the experimental results of Bockris et al.40 predicted the self-diffusion coefficient of the Rb+ ions to be larger than that of the Cl− ions. This is exactly converse in this work, while the relations revealed here should be rational from the point of view that the Rb atom has larger mass and radius. 3.3. Shear Viscosity. The shear viscosity is a measure of the resistance of a fluid to being deformed by shear stress. As mentioned above, in the MD simulations, it can be calculated using the EMD method as the time integral of the shear stress autocorrelation function through the GK formula. The correlation function has been averaged for 10 million MD steps to ensure sufficiently good statistics that repeated runs give the same plateau value. The time to compute the autocorrelation function (ACF) should be long enough to capture the decay of ACF and to ensure that the running integral reaches the plateau. Before calculating the viscosity, the decay and running integral behaviors of the stress tensor ACFs for molten ACl have been tested. Figure 4 shows the time evolution of the normalized stress tensor ACF and the running integral of molten KCl at NMT. The normalized ACFs of the three offdiagonal components all decay to zero after a certain time, which varies from 1.0 to 2.0 ps for different systems. The running integrals reach a plateau after some time, varying from 1.0 to 4.0 ps for different systems. In this work, a time of 5 ps is adopted for stress tensor ACFs calculation. 10204



Article

then complete the current database for molten salts’ physical properties. 3.5. Ionic Conductivity. Figure 8 shows the time evolution of the normalized charge flux ACF and the running integral for molten KCl at NMT. The normalized ACFs of the three components all decay to zero quickly. The running integrals reach a plateau after about 1 ps for all the systems. During ionic conductivity calculations, a time of 5 ps is adopted for charge flux ACFs calculation. The ionic conductivities of molten ACl salts have been calculated at different temperatures, i.e. NMTs and from 1100 to 1500 K with an interval of 100 K. Figure 9 compares the simulated ionic conductivity results of ACl melts with the available experimental data.51 The ionic conductivities obtained here are in good agreement with experiment, especially for molten KCl. In the experimental temperature range, the results of calculation are in agreement with experimental data within the error of 10−15%. Further research including the polarization effect should be able to improve the accuracy of the calculated results.31 Our results show that the ionic conductivities of molten ACl increase with increasing temperature, as can be expected from the increase in the self-diffusion coefficients that accompanies rising temperature.

transport coefficients of molten alkali earth chlorides will be further studied.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected](Z. Sun). *E-mail: [email protected] (G. Lu). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We acknowledge the financial support provided by the National Natural Science Foundation of China (Grant 21206038), the Specialized Research Fund for the Doctoral Program of Higher Education (New Teachers) (Grant 20120074120014), and the Fundamental Research Funds for the Central Universities.

■

REFERENCES

(1) Brun, C. L. Molten Salts and Nuclear Energy Production. J. Nucl. Mater. 2007, 360 (1), 1−5. (2) Waldrop, M. M. Nuclear Energy: Radical Reactors. Nature (London) 2012, 492 (7427), 26−29. (3) Groult, H.; Barhoun, A.; Briot, E.; Lantelme, F.; Julien, C. M. Electrodeposition of Zr on Graphite in Molten Fluorides. J. Fluorine Chem. 2011, 132 (12), 1122−1126. (4) Fukasawa, K.; Uehara, A.; Nagai, T.; Sato, N.; Fujii, T.; Yamana, H. Thermodynamic Properties of Trivalent Lanthanide and Actinide Ions in Molten Mixtures of LiCl and KCl. J. Nucl. Mater. 2012, 424 (1), 17−22. (5) Bradwell, D. J.; Kim, H.; Sirk, A. H.; Sadoway, D. R. MagnesiumAntimony Liquid Metal Battery for Stationary Energy Storage. J. Am. Chem. Soc. 2012, 134 (4), 1895−1897. (6) Trullas, J.; Alcaraz, O.; Gonzalez, L. E.; Silbert, M. Structure and Dynamics of Molten AgCl: The Inclusion of Induced Polarization. J. Phys. Chem. B 2003, 107 (1), 282−290. (7) Urahata, S. M.; Ribeiro, M. C. Molecular Dynamics Simulation of Molten LiNO3 with Flexible and Polarizable Anions. Phys. Chem. Chem. Phys. 2003, 5 (12), 2619−2624. (8) Bitrian, V.; Trullas, J. Molecular Dynamics Study of Polarizable Ion Models for Molten AgBr. J. Phys. Chem. B 2006, 110 (14), 7490− 7499. (9) Okamoto, Y.; Suzuki, S.; Shiwaku, H.; Ikeda-Ohno, A.; Madden, P. A. Local Coordination about La3+ in Molten LaCl3 and Its Mixtures with Alkali Chlorides. J. Phys. Chem. A 2010, 114 (13), 4664−4671. (10) Pauvert, O.; Salanne, M.; Zanghi, D.; Simon, C.; Reguer, S.; Thiaudiere, D.; Bessada, C. Ion Specific Effects on the Structure of Molten AF-ZrF4 Systems (A+ = Li+, Na+, and K+). J. Phys. Chem. B 2011, 115 (29), 9160−9167. (11) Salanne, M.; Simon, C.; Turq, P.; Madden, P. A. Heat-Transport Properties of Molten Fluorides: Determination from First-Principles. J. Fluor. Chem. 2009, 130 (1), 38−44. (12) Levesque, M.; Sarou-Kanian, V.; Salanne, M.; Gobet, M.; Groult, H.; Bessada, C.; Madden, P. A.; Rollet, A. L. Structure and Dynamics in Yttrium-Based Molten Rare Earth Alkali Fluorides. J. Chem. Phys. 2013, 138 (18), 184503. 1−7 (13) Sarou-Kanian, V.; Rollet, A. L.; Salanne, M.; Simon, C.; Bessada, C.; Madden, P. A. Diffusion Coefficients and Local Structure in Basic Molten Fluorides: in Situ NMR Measurements and Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2009, 11 (48), 11501−11506. (14) Fumi, F. G.; Tosi, M. P. Ionic Sizes and Born Repulsive Parameters in the NaCl Type Alkali Halides 1. J. Phys. Chem. Solids 1964, 25 (1), 31−43. Tosi, M. P.; Fumi, F. G. Ionic Sizes and Born Repulsive Parameters in the NaCl Type Alkali Halides 2. J. Phys. Chem. Solids 1964, 25 (1), 45−52.

4. CONCLUSIONS Molecular dynamics simulations of molten ACl (A denotes Li, Na, K, Rb, Cs) serials were performed systematically to complete basic research of simple ionic melts. The radial distribution functions and angular distribution functions were calculated to analyze the local structures of molten ACl. (1) When melting from solid, molten ACl lost long-range order but still kept short-range order and charge order. The local structure was affected by cationic size: as cationic radius increased from Li to Cs, the coordination number of ACl increased from 4.05 to 5.70. (2) Combining the coordination number of Li−Cl and the ADF of molten LiCl, it could be speculated that the main microconfiguration existing in molten LiCl was a distorted tetrahedron. For other molten ACl systems, the melts tended to keep the original crystal local structures, resulting in a distorted octahedron configuration. (3) Below 1500 K, temperature had no significant effect on the local structures of molten ACl. (4) The calculated transport coefficients of molten ACl were in good agreement with experimental values in the literature, generally within 5%−15%. The self-diffusion coefficients of molten ACl serials were relatively high and showed a positive dependence on temperature. At the same temperature, the ionic diffusivity decreased as the cationic size increased. The ionic conductivities of molten ACl also increased as the temperature increased. (5) The shear viscosities of molten ACl decreased as the temperature increased. Though overpredicting thermal conductivities of molten ACl, this work reproduced the negative temperature dependence of this property. Among ACl melts serials, molten LiCl has excellent ionic diffusivity, fluidity, and ionic conductivity, which makes it possible to improve the performance of melt solvents through addition of LiCl. The effect of LiCl on the structures and 10205



Article

(39) Cikit, S.; Akdeniz, Z.; Madden, P. A. Structure and Raman Spectra in Cryolitic Melts: Simulations with an ab Initio Interaction Potential. J. Phys. Chem. B 2014, 118, 1064−1070. (40) Bockris, J. M.; Hooper, G. W. Self-Diffusion in Molten Alkali Halides. Discuss. Faraday Soc. 1961, 32, 218−236. (41) Bockris, J. M.; Richards, S. R.; Nains, L. Self-Diffusion and Structure in Molten Group II Chlorides. J. Phys. Chem. 1965, 69 (5), 1627−1637. (42) Okazaki, S.; Miyamoto, Y.; Okada, I. Molecular Dynamics Study of the Dynamic Structure Factors of Molten LiCl and LiCl-CsCl. Phys. Rev. B 1992, 45 (5), 2055−2062. (43) Lewis, J. W.; Singer, K. Thermodynamic Properties and SelfDiffusion of Molten Sodium Chloride. A Molecular Dynamics Study. J. Chem. Soc., Faraday Trans. 2 1975, 71, 41−53. (44) Bengston, A.; Nam, H. O.; Saha, S.; Sakidja, R.; Morgan, D. First-Principles Molecular Dynamics Modeling of the LiCl-KCl Molten System. Comput. Mater. Sci. 2014, 83, 362−370. (45) Ciccotti, G.; Jacucci, G.; McDonald, I. R. Transport properties of molten alkali halides. Phys. Rev. A 1976, 13, 426−436. (46) Trullas, J.; Padro, J. A. Diffusive transport properties in monovalent and divalent metal-ion halide melts: A computer simulation study. Phys. Rev. B 1997, 55, 12210−12217. (47) Alcaraz, O.; Trullas, J. Mass and size dependence of single ion dynamics in molten monohalides. J. Chem. Phys. 2000, 113 (23), 10635−10641. (48) Wakao, M.; Minami, K.; Nagashima, A. Viscosity Measurements of Molten LiCl in the Temperature Range 886−1275 K. Int. J. Thermophys. 1991, 12 (2), 223−230. (49) Janz, G. J. NIST Properties of Molten Salts Databases, NIST SRD 27; National Bureau of Standards: Boulder, CO, 1992. (50) Janz, G. J. Molten Salts Data as Reference Standards for Density, Surface Tension, Viscosity and Electrical Conductance: KNO3 and NaCl. J. Phys. Chem. Ref. Data 1980, 9 (4), 791−830. (51) Janz, G. J.; Dampier, F. W.; Lakshminarayanan, G. R.; Lorenz, P. K.; Tomkins, R. P. T. Molten Salts: Volume I. Electrical Conductance, Density, and Viscosity Data (NSRDS-NBS-15); National Standard Reference Data Series; National Bureau of Standards: Boulder, CO, 1968. (52) Nagasaka, Y.; Nakazawa, N.; Nagashima, A. Experimental Determination of the Thermal Diffusivity of Molten Alkali Halides by the Forced Rayleigh Scattering Method. I. Molten LiCl, NaCl, KCl, RbCl, CsCl. Int. J. Thermophys. 1992, 13 (4), 555−574. (53) Ohtori, N.; Oono, T.; Takase, K. Thermal Conductivity of Molten Alkali Halides: Temperature and Density Dependence. J. Chem. Phys. 2009, 130 (4), 044505. 1−5.

(15) Galamba, N.; de Castro, C. N.; Ely, J. F. Thermal Conductivity of Molten Alkali Halides from Equilibrium Molecular Dynamics Simulations. J. Chem. Phys. 2004, 120 (18), 8676−8682. (16) Galamba, N.; de Castro, C. N.; Ely, J. F. Shear Viscosity of Molten Alkali Halides from Equilibrium and None Equilibrium Molecular-Dynamics Simulations. J. Chem. Phys. 2005, 120 (22), 224501. 1−9. (17) Nevins, D.; Spera, F. J. Accurate Computation of Shear Viscosity from Equilibrium Molecular Dynamics Simulations. Mol. Simul. 2007, 33 (15), 1261−1266. (18) Ohtori, N.; Salanne, M.; Madden, P. A. Calculations of the Thermal Conductivities of Ionic Materials by Simulation with Polarizable Interaction Potentials. J. Chem. Phys. 2009, 130 (10), 104507. 1−13. (19) Ribeiro, M. C. Chemla Effect in Molten LiCl/KCl and LiF/KF Mixtures. J. Phys. Chem. B 2003, 107 (18), 4392−4402. (20) Pauling, L. The Sizes of Ions and Their Influence on the Properties of Salt-Like Compounds. Z. Kristallogr. 1928, 67, 377−404. Pauling, L. The Influence of Relative Ionic Sizes on the Properties of Ionic Compounds. J. Am. Chem. Soc. 1928, 50 (4), 1036−1045. (21) Dixon, M.; Gillan, M. J. Structure of Molten Alkali Chlorides I. A Molecular Dynamics Study. Philos. Mag. B 1981, 43 (6), 1099− 1112. (22) Mayer, J. E. Dispersion and Polarizability and the van der Waals Potential in the Alkali Halides. J. Chem. Phys. 1933, 1 (4), 270−279. (23) Galamba, N.; Cabral, B. J. C. First principles molecular dynamics of molten NaCl. J. Chem. Phys. 2007, 126, 124502. 1−8. (24) Janz, G. J.; Tomkins, R. P. T.; Allen, C. B.; Downey, J. R., Jr.; Garner, G. L.; Krebs, U.; Singer, S. K. Molten Salts: Chlorides and MixturesElectrical Conductance, Density, Viscosity, and Surface Tension Data. J. Phys. Chem. Ref. Data 1975, 4 (4), 871−1178. (25) Allen, M. P.; Tlldesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (26) Melchionna, S.; Ciccotti, G.; Holian, B. L. Hoover NPT Dynamics for systems varying in size and shape. Mol. Phys. 1993, 78, 533−544. (27) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177−4189. (28) Nose, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (29) Hoover, W. G. Canonical dynamics: Equilibrium phase space distributions. Phys. Rev. A 1985, 31, 1695−1697. (30) Martyna, G.; Klein, M.; Tuckerman, M. Nose-Hoover chains: The canonical ensemble via continuous dynamics. J. Chem. Phys. 1992, 97, 2635−2643. (31) Jabes, B. S.; Agarwal, M.; Chakravarty, C. Structure and Transport Properties of LiF−BeF2 Mixtures: Comparison of Rigid and Polarizable Ion Potentials. J. Chem. Sci. 2012, 124 (1), 261−269. (32) Frenkel, D.; Smit, B. Understanding Molecular Simulation-From Algorithms to Applications; Academic Press: 1996. (33) Owczarek, E.; Hawlicka, E. Molecular Dynamics Study of CaCl2 in Methanol. J. Phys. Chem. B 2006, 110 (45), 22712−22718. (34) Lewis, J. W. E.; Singer, K.; Woodcock, L. V. Thermodynamic and structural properties of liquid ionic salts obtained by Monte Carlo computation. Part 2Eight alkali metal halides. J. Chem. Soc., Faraday Trans. 2 1975, 71, 301−312. (35) Levy, H. A.; Agron, P. A.; Bredig, M. A.; Danford, M. D. X-Ray and Neutron Diffraction Studies of Molten Alkali Halides. Ann. N. Y. Acad. Sci. 1960, 79 (11), 762−780. (36) Baranyai, A.; Ruff, I.; McGreevy, R. L. Monte Carlo simulation of the complete set of molten alkali halides. J. Phys. C: Solid State Phys. 1986, 19, 453−465. (37) Kittle, C. Introduction to Solid State Physics, 8th ed.; John Wiley & Sons, Inc.: 2005. (38) Pauvert, O.; Zanghi, D.; Salanne, M.; Simon, C.; Rakhmatullin, A.; Matsuura, H.; Okamoto, Y.; Vivet, F.; Bessada, C. In Situ Experimental Evidence for a Nonmonotonous Structural Evolution with Composition in the Molten LiF-ZrF4 System. J. Phys. Chem. B 2010, 114, 6472−6479. 10206


Molecular Dynamics Simulations of the Local Structures and Transport

Recommend Documents