Thermal and Transport Properties of Molten Chloride Salts with Polarization Effect on Microstructure

: Molten chloride salt is recognized as a promising heat transfer and storage medium in concentrating solar power in recent years, but there is a serious lack for thermal property data of molten chloride salts. In this work, local structures and thermal properties for molten chloride salt—including NaCl, MgCl 2 , and ZnCl 2 —were precisely simulated by Born–Mayer–Huggins (BMH) potential in a rigid ion model (RIM) and a polarizable ion model (PIM). Compared with experimental data, distances between cations, densities, and heat capacities of molten chloride slats calculated from PIM agree remarkably better than those from RIM. The polarization effect brings an extra contribution to screen large repulsive Coulombic interaction of cation–cation, and then it makes shorter distance between cations, larger density and lower heat capacity. For NaCl, MgCl 2 , and ZnCl 2 , PIM simulation deviations of distances between cations are respectively 3.8%, 3.7%, and 0.3%. The deviations of density and heat capacity for NaCl between PIM simulation and experiments are only 0.6% and 2.2%, and those for MgCl 2 and ZnCl 2 are 0.7–10.7%. As the temperature rises, the distance between cations increases and the structure turns into loose state, so the density and thermal conductivity decrease, while the ionic self-diffusion coefﬁcient increases, which also agree well with the experimental results.


Introduction
Concentrating solar thermal power (CSP) is a promising technique for high temperature solar energy utilization [1]. Molten salts have the advantages of wide working temperature range, high heat capacity, low viscosity, low vapor pressure, and low material prices [2]; consequently, they are widely used as heat transfer and storage materials in CSP. Molten salts include alkali or alkali earth halides (chloride salts, fluoride salts, and so on), and oxyacid salts (carbonates, nitrates, and so on). In recent years, molten chloride salts including LiCl, NaCl, KCl, MgCl 2 , and ZnCl 2 have been researched as promising heat transfer materials with high operating temperature [3]. Another unique advantage is the tremendous nature resources worldwide, and it is reported that there are 38.81 billion tons magnesium chlorides in Qinghai salt lake of China [4].
Many researchers experimentally reported the microstructure and properties of molten chloride salts. McGreevy [5] investigated the experimental structure of molten alkali halides by the technique of neutron scattering and energy-dispersive X-ray diffraction (EDXRD). Biggin and Enderby [6] researched the structural characteristic of ZnCl 2 at melting point by neutron diffraction technique. Besides, high-energy X-ray diffraction method was taken to obtain the structure of ZnCl 2 by Zeidler et al. [7], showing the result of cornersharing and edge-sharing tetrahedral unit. Pavlatou et al. [8] studied the vibrational spectra of network-forming ionic melts (LaCl 3 and ZnCl 2 ) by the Raman and infrared spectra.
where the first term is electrostatic potential, the second term is repulsion potential and the last two terms represent dipole dispersion energies. The potential parameters in RIM of NaCl and ZnCl 2 are reported in the research [31,32]. A new series of parameters in RIM of MgCl 2 used in this paper are put forward. The interactions between Cl-Cl pairs are adopted from Seo et al. [32] and the rest pair parameters are adopted from Yuen et al. [33], which are listed in Table 1. Table 1. Potential parameters of RIM for molten chloride salts.

Model Pair
A ij (kJ·mol −1 ) B ij (Å −1 ) σ ij (Å) C ij (Å 6 ·kJ·mol −1 ) D ij (Å 8  In order to improve the force field, the polarization ion model (PIM) is researched by Madden and Salanne [34,35]. The PIM is shown by Wilson and Madden [28] with the shell model. In this model, the polarizability of ion is combined by representing an ion as the two-particle system with a charged core and a charged shell connected by an atomic spring. The induction effect is included in the polarization term. The pair potentials are described as where the α i is the polarizability of the i ion and the f ij n (r ij ) is the Tang-Toennies dispersion damping functions [36]. In Equation (3), the first term is the energy cost for the deformation of electron cloud and the last two terms represent the charge-dipole and dipole-dipole energies.
The polarization effect is the deformation of the outer electron cloud caused by the influence of an external electric field [37]. The deformation of the electron cloud and the induced dipole of anions result in the structural changes of ions. Thus, the potential parameters in PIM of molten chloride slats should be restudy. Salanne et al. [35] used the first-principles calculation method to obtain the PIM parameters of NaCl; also, Sharma and Wilson [30] provided the PIM parameters of MgCl 2 and ZnCl 2 , which are listed in Table 2. Table 2. Potential parameters of PIM for molten chloride salts.

Model
Pair

Simulation Details and Conditions
The RIM simulations were carried out by LAMMPS code [38] and the PIM simulations were performed by CP2K package [39]. All the simulations were carried out under the Periodic Boundary Condition (PBC), removing the side effects. The PPPM solver was used to conduct long-range corrections and the precise was set equal to 1.0 × 10 −6 . The initial velocity was set up to obey Gaussian distribution. Newtonian equation of motion was solved by Verlet computing method [40]. The PIM simulations include the iterations for induced dipoles in CP2K, so that the computational complexity of PIM is higher than that of RIM. The number of ions, space group and the details of unit cell are listed in Table 3. The sizes of initial simulation boxes in RIM are 45 distance was set to 20 Å and 10 Å for RIM and PIM, which is slightly shorter than the half of edge length. Although the ions in the PIM simulations are less than the RIM simulations, the PIM simulated results have higher accuracy, which will be illustrated in Section 3. The experimental melting points of NaCl, MgCl 2 , and ZnCl 2 are 1073 K, 987 K, and 591 K [14]. For all the MD simulations, the system was firstly heated to 1700 K and then the system was cooled down to different temperature for simulations in NPT ensemble with the Nosé-Hoover thermostat and barostat at a constant pressure of 0.1 MP. The damping parameters were 100 fs and 500 fs for thermostat and barostat, respectively. After that, th system was simulated on the equilibrated cell volume in NVT ensemble. The simulated time step was set to 1 fs. In order to get precise calculation results, the equilibrium time lasted 1 ns for each simulation. Partial radial distribution function (RDF) is the typical local structure characteristic of liquid molten salts, which defines the distribution probability of j-type ions centered on i-type ion within a radius r where ρ j denotes the number density of species j and n ij denotes the number of j-type ions locating in a globular of radius r around the center of an i-type ion.

Coordination Number
The quantity of j-type ions locating in a globular of radius r min around the center of an i-type ion is the definition of coordination number (CN), and r min is the first peak valley position of RDF. N ij can be calculated as

Density
The densities of molten salts could be computed from the formula where M is the molar mass, N is the number of particles, N a is the Avogadro's constant, and V denotes the equilibrated volume of the simulation system at the given temperature in NPT ensemble.

Specific Heat Capacity
Specific heat capacities were calculated from enthalpies with the formula where U is the total internal energy of the system consisting of kinetic energy and potential energy. Thus, the value of heat capacity can be calculated from the slope of enthalpy variation with temperature.

Thermal Conductivity
The reverse non-equilibrium molecular dynamics (RNEMD) method [22][23][24] can be applied to forecast thermal conductivities of molten salts. RNEMD method executes kinetic energy changes between the center and the bottom slabs and measures the inducing temperature gradient with linear response theory. At first, the system is equilibrated at object temperature for 1 ns. After that, the system is simulated in NVE ensemble to swap kinetic energy between different slabs for 5 ns. Thermal conductivity is calculated from the total transformation of kinetic energy and the ultimate temperature gradient [22] as where t is the swap time of kinetic energy, L x and L y denote the lengths of simulation box in x and y directions, respectively.

Self-Diffusion Coefficient
Ion self-diffusion coefficient could reveal the ion transport characters of molten salts. In the cubic cell, the time-dependent mean square displacement (MSD) method [41] was carried out to compute the ion self-diffusion coefficient. In the recent researches [42,43], this method has been carried out to obtain the self-diffusion coefficients of molten salts, which agree well with the experiment results. MSD and ion self-diffusion coefficient are calculated as where the brackets <···> denote the time averages and r i is the position of ion i.

Local Structures in PIM
To study the relation between microstructure and thermodynamics, the simulated partial radial distribution function (RDF) and coordination number (CN) are illustrated as Figures 1 and 2 and Table 4. It can be seen that all the RDFs have comparatively high first peaks and the height of other peaks decrease gradually, which demonstrates the disorder microstructures of molten chloride salts in long distance. For the ion pair of cation-anion, the first peak is sharper and the first valley is lower than that of cation-cation, indicating that cation-anion pair has a stronger interaction. Figure 1 and Table 4 show the curves and characteristic values of RDFs by PIM simulation. Compared with the experiments by neutron diffraction technique [5], the first peak position of Na + -Na + from PIM is slightly overestimated by 3.8% at 1100 K. Near the melting point, RDFs of MgCl 2 in PIM also have a good agreement with the experimental values [25]. For ZnCl 2 , the first peak position of Zn 2+ -Zn 2+ at 600 K is 3.80 Å [6], and the PIM simulation result is 3.79 Å. Thus, the PIM can provide the reasonable force field between particles and obtain the accurate structures of molten chloride salts. When the temperature increases, the first peak positions of RDFs for cation-anion decrease while those for cation-cation increase, showing that the distances between cations and anions become closer while the distances of cations become larger.

Polarization Effect on Structures
To research polarization effect on structures of molten chloride salts, RIM and PIM simulations are performed with the same conditions in Section 3.2. The simulation results of RDF and CN in RIM and PIM are displayed in Figure 3 and Table 5. The first peak position of Na + -Na + is 4.05 Å from PIM and 4.10 Å from RIM, which are slightly different from experiment value of 3.90 Å. However, there is a large difference between PIM simulated distances of MgCl2/ZnCl2 and those in RIM. When polarization effect is included, the first peak position of Mg 2+ -Mg 2+ changed from 4.38 Å in RIM to 3.95 Å in PIM, comparing with the experiment value of 3.81 Å. In Figure 3b, the first peak position of Zn 2+ -Zn 2+ in PIM is 3.79 Å and fits very well with the experiment value of 3.80 Å, while that in RIM is 4.30 Å.  The coordination number curves of NaCl, MgCl 2 and ZnCl 2 in diverse working temperature condition are displayed in Figure 2. It can be observed that the coordination number diminishes from low temperature to high temperature in all simulations. For NaCl, the coordination number of Na + -Cl − from PIM is 5.06 at 1100 K, which agrees well with the experimental value of 4.83 [5]. As the temperature rises, the coordination number of Mg 2+ -Cl − changes from 5.21 to 4.72, while that of Zn 2+ -Cl − changes from 5.36 to 5.02. According to the analysis of RDF and CN, the structure of molten chloride salts becomes looser under high temperature condition.

Polarization Effect on Structures
To research polarization effect on structures of molten chloride salts, RIM and PIM simulations are performed with the same conditions in Section 3.2. The simulation results of RDF and CN in RIM and PIM are displayed in Figure 3 and Table 5. The first peak position of Na + -Na + is 4.05 Å from PIM and 4.10 Å from RIM, which are slightly different from experiment value of 3.90 Å. However, there is a large difference between PIM simulated distances of MgCl 2 /ZnCl 2 and those in RIM. When polarization effect is included, the first peak position of Mg 2+ -Mg 2+ changed from 4.38 Å in RIM to 3.95 Å in PIM, comparing with the experiment value of 3.81 Å. In Figure 3b, the first peak position of Zn 2+ -Zn 2+ in PIM is 3.79 Å and fits very well with the experiment value of 3.80 Å, while that in RIM is 4.30 Å.   In addition, the coordination number of Zn 2+ -Zn 2+ in PIM is larger than that in RIM, indicating that the cation clusters in PIM are more compact than the cation clusters in RIM. That is to say, the large repulsive coulombic interactions between divalent cations lead to the larger distance between ions in RIM. As soon as the induced dipoles of anions are considered, the large Coulombic force of cation-cation can be screened by the polarization effect, which brings an extra contribution to take the shorter distance of cation-cation. The higher polarizability of Zn 2+ leads to the stronger polarizable force between ions, so that the PIM performs better to derive the first peak position of Zn 2+ -Zn 2+ . In Figure 3, the second peaks of Mg 2+ -Mg 2+ and Zn 2+ -Zn 2+ in PIM show characteristic of the intermediateranged order, which reveals the network structure of MgCl2 and ZnCl2. Besides, the snapshots of network structures for molten salts are plotted in Figure 4. Thus, the PIM simulations are more reasonable and correspond better to the physical laws.  In addition, the coordination number of Zn 2+ -Zn 2+ in PIM is larger than that in RIM, indicating that the cation clusters in PIM are more compact than the cation clusters in RIM. That is to say, the large repulsive coulombic interactions between divalent cations lead to the larger distance between ions in RIM. As soon as the induced dipoles of anions are considered, the large Coulombic force of cation-cation can be screened by the polarization effect, which brings an extra contribution to take the shorter distance of cation-cation. The higher polarizability of Zn 2+ leads to the stronger polarizable force between ions, so that the PIM performs better to derive the first peak position of Zn 2+ -Zn 2+ . In Figure 3, the second peaks of Mg 2+ -Mg 2+ and Zn 2+ -Zn 2+ in PIM show characteristic of the intermediate-ranged order, which reveals the network structure of MgCl 2 and ZnCl 2 . Besides, the snapshots of network structures for molten salts are plotted in Figure 4. Thus, the PIM simulations are more reasonable and correspond better to the physical laws.

Density
Density is the critical thermodynamics property of molten chloride salts, which reflects the change of the system volume. The experimental melting points of NaCl, MgCl2, and ZnCl2 are 1073 K, 987 K, and 591 K [14]. The systems were heated to molten state at

Density
Density is the critical thermodynamics property of molten chloride salts, which reflects the change of the system volume. The experimental melting points of NaCl, MgCl 2 , and ZnCl 2 are 1073 K, 987 K, and 591 K [14]. The systems were heated to molten state at 1700 K in NPT ensemble. After that, the systems in both models were cooled down to the target temperature and equilibrated in the NPT ensemble for 1 ns. Specific heat capacity was also computed with the similar method.
In Figure 5, the simulated densities of molten chloride salts are compared with the experimental results by Yaffe et al. [44], Janz et al. [14] and Angell et al. [45]. It is obvious that the density is in linear negative correlation with temperature. For molten ZnCl 2 , the PIM simulated densities decrease from 2.50 g/cm 3 at 600 K to 2.28 g/cm 3 at 1000 K, which is slightly smaller than the experiment values by 1.1%. When the temperature is over 773 K, the RIM simulated densities of ZnCl 2 become into a gross difference from experiments. Besides, the slope of density in RIM is different from that in PIM for MgCl 2 and ZnCl 2 . Combined with the analysis of polarization effect on local structure in Section 3.2, the large Coulombic interactions of ions induce the longer distance between cations in RIM for lack of the polarization effect. In addition, the looser coordination structures and longer distance between cations in RIM cause the greatly larger volume of system at high temperature, so the densities from RIM are much smaller than real values. As listed in Table 6, the fit formulas of PIM simulation results agree well from the experimental data, with the mean percentage errors of 0.3%, 9.5%, and 1.1% for NaCl, MgCl 2 , and ZnCl 2 respectively. The error of MgCl 2 is larger than others, which may be caused by the error of the fitting potential parameters searched for MgCl 2 . The polarization effect included in the potential of MgCl 2 can bring extra contribution to screen repulsive Coulombic force, but the repulsive force is still slightly larger and induced larger volume for MgCl 2 .

Specific Heat Capacity
As shown in Figure 6, the enthalpy variations of molten chloride salts are in linear relation with temperature, which indicates the specific heat capacity is constant from Equation (8). Listed in Table 7, the specific heat capacities of NaCl, MgCl 2 and ZnCl 2 in PIM are respectively 65.5 J/(mol·K), 99.7 J/(mol·K) and 111.7 J/(mol·K) and their errors are 2.2%, 7.8% and 10.7%. Simulated heat capacities of NaCl in RIM and PIM are both in good agreement with experiment [24], while the results of MgCl 2 and ZnCl 2 in RIM are 38.8% and 32.9% larger than the experimental data [14,46]. In conclusion, the polarization effect can reflect the accurate ability of thermal storage for molten divalent chloride (MgCl 2 and ZnCl 2 ).

Thermal Conductivity
The RNEMD method was used to compute the thermal conductivities. To input a heat flux and obtain a temperature distribution, the simulation box was evenly split into 20 bins along the 'z' axis. The 1st bin was defined as the 'cool' bin, and the 11th bin was defined as the 'hot' bin. The heat flux was generated by swapping the kinetic energy of an atom in the cool bin and one in the hot bin, which results in the increasing temperature in the hot bin and decreasing temperature in the cool bin [22]. After equilibrating for 1 ns, the simulations transferred the kinetic energy over 5 ns in NVE ensemble. Different swap rates were tried for the simulations, and it was found that kinetic energy swap rates between 100 and 400 time steps brought about the linear temperature gradients. The total average kinetic energy exchange as well as the slopes of linear temperature gradient in Figure 7 were substituted into Equation (10) to calculate the thermal conductivities. The swapping rates are same for all molten chloride salts at 400 fs.

Thermal Conductivity
The RNEMD method was used to compute the thermal conductivities. To input a heat flux and obtain a temperature distribution, the simulation box was evenly split into 20 bins along the 'z' axis. The 1st bin was defined as the 'cool' bin, and the 11th bin was defined as the 'hot' bin. The heat flux was generated by swapping the kinetic energy of an atom in the cool bin and one in the hot bin, which results in the increasing temperature in the hot bin and decreasing temperature in the cool bin [22]. After equilibrating for 1 ns, the simulations transferred the kinetic energy over 5 ns in NVE ensemble. Different swap rates were tried for the simulations, and it was found that kinetic energy swap rates between 100 and 400 time steps brought about the linear temperature gradients. The total average kinetic energy exchange as well as the slopes of linear temperature gradient in Figure 7 were substituted into Equation (10) to calculate the thermal conductivities. The swapping rates are same for all molten chloride salts at 400 fs. In order to approximate the finite-size effect on thermal conductivity of molten chloride salt by RNEMD method, a series of different size simulations are carried out, and the lengths in x and y dimension are the same while the length of z dimension is changed. The size-dependent thermal conductivities of molten chloride salt are displayed in Figure 8. It is obvious that the simulated thermal conductivity increases non-linearly with increasing length of z dimension and approaches a stable value when the length of z dimension is over 100 Å. In order to approximate the finite-size effect on thermal conductivity of molten chloride salt by RNEMD method, a series of different size simulations are carried out, and the lengths in x and y dimension are the same while the length of z dimension is changed. The size-dependent thermal conductivities of molten chloride salt are displayed in Figure 8. It is obvious that the simulated thermal conductivity increases non-linearly with increasing length of z dimension and approaches a stable value when the length of z dimension is over 100 Å.
As shown in Figure 9, the thermal conductivities of molten chloride salts are in negative linear correlation with temperature. For molten NaCl, the simulated results of RNEMD decrease from 0.571 W/(m·K) to 0.424 W/(m·K). There is a mean percentage error of 8.1% between simulated results and experiments from Nagasaka et al. [15], so it proves that RNEMD method is suitable for calculating the thermal conductivities. Compared with the experimental values from Singh et al. [47] and Cornwell [48], the RNEMD simulated results of MgCl 2 and ZnCl 2 agree well and are completely predicted over the whole operating temperature. In general, the simulated thermal conductivities of molten chloride salts in different temperature are in the range of 0.213~0.571 W/(m·K). As shown in Figure 9, the thermal conductivities of molten chloride salts are in negative linear correlation with temperature. For molten NaCl, the simulated results of RNEMD decrease from 0.571 W/(m·K) to 0.424 W/(m·K). There is a mean percentage error of 8.1% between simulated results and experiments from Nagasaka et al. [15], so it proves that RNEMD method is suitable for calculating the thermal conductivities. Compared with the experimental values from Singh et al. [47] and Cornwell [48], the RNEMD simulated results of MgCl2 and ZnCl2 agree well and are completely predicted over the whole operating temperature. In general, the simulated thermal conductivities of molten chloride salts in different temperature are in the range of 0.213~0.571 W/(m·K).   As shown in Figure 9, the thermal conductivities of molten chloride salts are in negative linear correlation with temperature. For molten NaCl, the simulated results of RNEMD decrease from 0.571 W/(m·K) to 0.424 W/(m·K). There is a mean percentage error of 8.1% between simulated results and experiments from Nagasaka et al. [15], so it proves that RNEMD method is suitable for calculating the thermal conductivities. Compared with the experimental values from Singh et al. [47] and Cornwell [48], the RNEMD simulated results of MgCl2 and ZnCl2 agree well and are completely predicted over the whole operating temperature. In general, the simulated thermal conductivities of molten chloride salts in different temperature are in the range of 0.213~0.571 W/(m·K).

Ionic Self-Diffusion Coefficient
In order to research the transport property of molten chloride salt ions, the ionic self-diffusion coefficients were calculated for 1 ns in NVT ensemble with the mean square displacement (MSD) method. It is evident that the MSD presents positive linear temperature dependence in Figure 10. Compared with experiments from Janz et al. [14], there are mean percentage errors of 5.9% and 5.7% between simulated results and experiments of Na + and Cl − from 1100 K to 1500 K in Figure 11a, respectively. In Figure 11b, both Zn 2+ and Cl − are difficult to move near melting point with quite small self-diffusion coefficients, because of the glassy state with the network structure [49], while Na + and Cl − can move more easily at melting point without the network structure. When the temperature increases, the ions become easier to move for the break of network structure. It can be seen that self-diffusion coefficients of all cations are larger than Cl − in Figure 11, resulting from the limited mobility of Cl − by the larger mass and ionic radius. Na and Cl from 1100 K to 1500 K in Figure 11a, respectively. In Figure 11b, both Zn and Cl − are difficult to move near melting point with quite small self-diffusion coefficients, because of the glassy state with the network structure [49], while Na + and Cl − can move more easily at melting point without the network structure. When the temperature increases, the ions become easier to move for the break of network structure. It can be seen that self-diffusion coefficients of all cations are larger than Cl − in Figure 11, resulting from the limited mobility of Cl − by the larger mass and ionic radius.

Conclusions
Microstructure and thermopyhsical properties of molten chloride salts (NaCl, MgCl2, and ZnCl2) were researched by molecular dynamics simulations with RIM and PIM. Compared with results from RIM, PIM simulations reveal that the polarization effect brings an extra contribution to screen the large Coulombic interactions, which makes the shorter distance between cations. For divalent chloride, the high charge density of divalent cation leads to more deformation of the electron cloud, so that the polarizable force between anion and cation is stronger. The RIM simulations are unable to precisely predict the densities and heat capacities of divalent chloride (MgCl2 and ZnCl2) for lack of considering

Conclusions
Microstructure and thermopyhsical properties of molten chloride salts (NaCl, MgCl 2 , and ZnCl 2 ) were researched by molecular dynamics simulations with RIM and PIM. Compared with results from RIM, PIM simulations reveal that the polarization effect brings an extra contribution to screen the large Coulombic interactions, which makes the shorter distance between cations. For divalent chloride, the high charge density of divalent cation leads to more deformation of the electron cloud, so that the polarizable force between anion and cation is stronger. The RIM simulations are unable to precisely predict the densities and heat capacities of divalent chloride (MgCl 2 and ZnCl 2 ) for lack of considering polarization effect, while PIM simulations improve the accuracy of simulated densities and heat capacities. Thus, the PIM performs better for molten divalent chloride. There are 0.3%, 1.1%, and 9.5% mean errors between PIM simulated densities and experiments for NaCl, ZnCl 2 , and MgCl 2 respectively, and those of heat capacities are 2.2-10.7%. As the temperature rises, the distance between cations increases and the structure turns into loose state, so the density and thermal conductivity decrease, while the ionic self-diffusion coefficient increases, which agree well with the experimental results. This work provides an alternative way to build the property database of molten chloride salts for the precise design and regulation of high-temperature thermal storage system.

Conflicts of Interest:
The authors declare no conflict of interest.

A
BMH repulsive parameter B BMH repulsive parameter C Van der Waals parameter D Van der Waals parameter L length dimension n number of particles q charge of ion r distance between ions T temperature of the system t simulation time U Born-Mayer-Huggins-Tosi-Fumi potential V volume of the system v velocity Greek symbols α polarizability of the ion ρ number density σ BMH repulsive parameters Subscripts i ion i j ion j k slab k x direction x y direction y