- Email: [email protected]
Thermal performance of a binary carbonate molten eutectic salt for high-temperature energy storage applications
Applied Energy 262 (2020) 114418
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Thermal performance of a binary carbonate molten eutectic salt for hightemperature energy storage applications
T
Gechuanqi Pana,b,c,e, Xiaolan Weid, Chao Yua,b,c, Yutong Lua,b,c,e, Jiang Lie, Jing Dinga,b,c, Weilong Wanga,b,c, Jinyue Yanf,g a
School of Data and Computer Science, Sun Yat-Sen University, Guangzhou, PR China School of Materials Science and Engineering, Sun Yat-Sen University, Guangzhou, PR China c School of Intelligent Systems Engineering, Sun Yat-Sen University, Guangzhou, PR China d School of Chemistry and Chemical Engineering, South China University of Technology, Guangzhou, PR China e National Supercomputer Center in Guangzhou, Guangzhou, PR China f Division of Energy Processes, Royal Institute of Technology, Stockholm, Sweden g School of Business, Society and Energy, Mälardalen University, Västerås, Sweden b
HIGHLIGHTS
GRAPHICAL ABSTRACT
simulation is used to com• Molecular pute structures-properties of Na CO 2
• •
3
K2CO3. Temperature and component dependence are studied from microscopic view. The temperature-thermophysical properties-composition correlation are obtained.
ARTICLE INFO
ABSTRACT
Keywords: Energy storage Molten carbonates Thermophysical properties Microstructures Concentrating solar power
Molten carbonate eutectic salts are promising thermal storage and heat transfer fluid materials in solar thermal power plant with the feature of large specific heat capacity, wide operating temperature range and little corrosive. The high-temperature properties of molten carbonates should be determined accurately over the entire operating temperature for energy system design. In this paper, molecular dynamic simulation is used to study temperature and component dependence of microstructures and thermophysical properties of the binary carbonate molten salt. Negative linear temperature dependence of densities and thermal conductivities of binary mixtures of different components is confirmed with respect to the distances of ion clusters. Besides, positive linear temperature dependence of self-diffusion coefficient is also obtained. When temperature is constant, densities and thermal conductivities of binary mixtures are linearly related with components. Self-diffusion coefficients of CO32− firstly increase and then decrease with increasing mole fraction of Na2CO3. The temperature-thermophysical properties-composition correlation formulas are obtained, and the database of thermophysical properties of molten carbonate salts over the entire operating temperature is complemented, which will provide the essential data for heat transfer and storage system design, operation, and optimization in CSP.
E-mail addresses: [email protected] (J. Ding), [email protected] (W. Wang). https://doi.org/10.1016/j.apenergy.2019.114418 Received 24 October 2019; Received in revised form 16 December 2019; Accepted 17 December 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.
Applied Energy 262 (2020) 114418
G. Pan, et al.
Nomenclature
z length dimension in the z direction Kb, Kθ, Kφ, r0, θ0, φ0 bond, angle, improper parameters length dimension in the x and y direction Lx , Ly v velocity m atomic mass λ thermal conductivity V equilibrium volume NA Avogadro’s constant
born-type pair potential Uij qi, qj charge of ion distance between particle i and j rij Aij, ρ, σij born repulsive parameters Cij, Dij Van der Waals parameters Ubond, Uangle, Uimproper bond, angle, improper potential t simulation time N particle number M molar mass i particle i j particle J k slab k x x dimension hot hot slab T temperature
Subscripts b θ ψ y cold
bond term angle term improper term y dimension cold slab
1. Introduction Molten carbonates, which have large specific heat capacity, wide operating temperature range and little corrosive, have been considered as very promising thermal storage and heat transfer fluid materials [1,2] in solar thermal power plant [3,4], fuel cell [5,6], nuclear fuel reprocessing etc. [7,8]. For efficient utilization of high-temperature heat, thermodynamic properties and transport properties of molten carbonates should be determined accurately over the entire operating temperature extent for device and system design in relevant industry. Some properties are hard or impossible to be obtained by means of experiments because of high-temperature and corrosion problems. Large discrepancies among the experimental value of thermal conductivity have been observed. Moreover, components of carbonates are varying due to evaporation, thermal decomposition, etc.. Thus, it is desired to study the temperature and component dependence of microstructures and thermophysical properties. MD has been demonstrated to be a powerful tool for the determination of the microscopic and macroscopic properties of molten salts in recent years [9,10]. Tissen and Janssen [11,12] determined the parameters of pair interaction of carbonates including electrostatic interaction, repulsion interaction and induction energy at the Hartree-Fock SCF level and the fitting parameters represented structures and diffusion coefficients well although dispersion interaction was ignored. They observed significant changes of structures from molten Li2CO3 to K2CO3 [13]. Archer [14] created an empirical model for CaCO3 which was constructed by fitting to structures and properties of aragonite and calcite. The new empirical model that the instability of the model of Fisler [15] was removed was stable for a wide range of carbonate structures, and reproduced experimental results with a reasonable accuracy. Costa [16] used the fluctuating charge model to performed simulations Li2CO3 and K2CO3, whose parameters were fitted by ab initio
Fig. 1. RNEMD method schematic view.
Table 2 Numbers of ions in simulated Na2CO3-K2CO3 of different components. Mole fraction of Na2CO3
S-1 S-2 S-3 S-4 S-5 S-6 S-7
0.87 0.75 0.62 0.50 0.38 0.25 0.13
Numbers of ions total
Na+
K+
C4+
O2−
9216
2673 2304 1905 1536 1167 768 399
399 768 1167 1536 1905 2304 2673
1536
4608
Table 1 Parameters of intermolecular and intramolecular interaction for Na2CO3-K2CO3. Atom
qi (e)
Aii (Kcal·mol−1)
ρii (Å)
σ (Å)
Cii (Å6·Kcal·mol−1)
Dii (Å8·Kcal·mol−1)
Na K C O
1.00 1.00 1.54 −1.18
6.0791 6.0791 10.9523 3.2648
0.2899
1.14 1.39 1.10 1.33
0.00
0.00
Atoms CeO OeCeO OeCeOeO
Kb (Kcal·mol−1 ·Å−2) = 730.641 Kθ (Kcal·mol−1 ·rad−2) = 20.607 Kφ (Kcal·mol−1 ·rad−2) = 6.959
2
r0(Å) = 1.16 θo (deg) = 120.0 φo (deg) = 0.00
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 3 RDFs of molten Na2CO3-K2CO3 of different components.
calculations of an isolated carbonate anion. Ottochian used Tissen and Janssen’s model (JT model) to study electric conductivity and diffusion coefficients of Li2CO3-Na2CO3 mixtures and found diffusion of lithium ions depended on its concentration level with respect to sodium ions in
the mixtures [17]. Corradini [18] parameterized a non-polarizable force field for Li2CO3-K2CO3 mixture by fitting densities and radial distribution functions to experiment and first principle molecular dynamics and presented an analysis of the temperature and pressure 3
Applied Energy 262 (2020) 114418
G. Pan, et al.
dependence of transport properties of the Li2CO3-Na2CO3 mixtures including self-diffusion coefficients, viscosity, and ionic conductivity. Recently, we have used the JT model to simulate local structures and thermophysical properties of pure alkali carbonate and obtained the variation law of microstructures and thermophysical properties with temperature [19,20]. However, to our best knowledge, thermal properties of molten carbonate salts are quite limited, and especially there is no comprehensive database for the entire operation temperature for heat transfer system design, operation and optimization in a concentrating solar power (CSP) station. Since it is hard to obtain all the properties due to the temperature and acuuracy limitation of measuring devices, it is necessary to obtain the effects of temperature and component on microstructures and thermophysical properties of molten salts under high-temperature conditions. This paper aims to gain temperature-thermophysical propertiescomposition correlation and complement the database of thermophysical properties of molten carbonate salts. Reverse nonequilibrium molecular dynamics (RNEMD) method is implemented for determining thermal conductivity. Moreover, the temperature and component dependence is investigated by correlating transport properties with local structures.
3.1.2. Coordination number Coordination number is the average number of β-type ions locating 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 by integrating the RDF as follow:
N
jik
+ Aij exp
rij
ij
rij
Cij
Dij
rij6
rij8
Uimproper = K (
1
using
the
following 1
Cij = (Cii Cjj ) 2 , Dij = (Dii Djj ) 2 and
ij
mixing
= 2(
ii
1
+
rules: jj
2 r jk
2rij rik
(7)
i Na, K , C , O
(8)
VE NA
3.2.2. Specific heat capacity The specific heat capacity can be determined simultaneously with density and computed from enthalpy using Eq. (9):
Cp = (
The JT model has been demonstrated to be able to reproduce the properties of the molten alkali carbonates [21]. The parameters of intermolecular interaction used in this work are the same to Tissen and Janssen [11,12], while the bond parameters are referred to Corradini [18], and angle and improper parameters are referred to Archer [14]. The detailed values are shown in Table 1. For the Born-type pair potential (Eq. (1)), the cross terms were
computed
rij2 + rik2
1
where Ni and Mi stand for the particle number and the molar mass of species i, while NA denotes Avogadro’s constant, VE corresponds to the average equilibrium volume.
(4)
2
cos
Ni Mi
(3) 0)
=
=
(2)
2 0)
Uangle = K (
(6)
3.2.1. Density The density of molten carbonate mixture can be determined simply by implementing an isothermal-isobaric ensemble (NPT) running and calculated as follow:
(1)
r0 ) 2
Ubond = Kb (r
g (r ) r 2dr
3.2. Thermophysical properties
Tissen and Janssen’s non-polarizable model is chosen for simulations of Na2CO3-K2CO3 mixture. The JT model is represented by an intermolecular interaction (born-type, Eq. (1)) plus a flexible carbonate intramolecular interaction (Eqs. (2)–(4)). The first two terms of right of Eq. (1) stand for columbic potential and overlapping repulsion while the last two terms correspond to dispersion energies.
qi qj
0
3.1.3. Angular distribution function In order to obtain the image of local structure of molten carbonates, angular distribution function (ADF), which can reflect the bond orientation, also should be obtained. ADF depends only on the coordinates of atoms exported during the simulations. For a centered atom i, and its two nearest neighbor atoms j and k, the bond angle can be calculated as follow:
2. Forcefield
Uij =
rmin
=4
H )P T
(
H )P = T
(U + PV ) T P
(9)
where U is the total energy of the system comprising kinetic energy and potential energy. Thus, specific heat capacity is the slope of enthalpy versus temperature.
1
Aij = (Aii Ajj ) 2 ,
) 1.
1
3. Evaluated properties and methods 3.1. Structure information 3.1.1. Partial radial distribution function Partial radial distribution function (RDF) is often used to characterize the local structure of molten salts, which is defined as Eq. (5):
g (r ) =
dN (r )
1 4
r2
dr
(5)
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 position of first peak of RDF is the average ionic equilibrium distance rmax between species α and β.
Fig. 2. RDF sample of molten Na2CO3-K2CO3 mixtures of different components at 1200 K. 4
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 4 Coordination number curves of molten Na2CO3-K2CO3 of different components.
3.2.3. Self-diffusion coefficient In order to study the ion transport properties of molten Na2CO3K2CO3 mixture, the ion self-diffusion coefficients are calculated. The ion self-diffusion coefficients can be calculated by either the time-
integration of the velocity autocorrelation function (ACF) or the Einstein equation based time-dependent mean-squared displacement (MSD). The latter method is employed here. The corresponding calculation formula is expressed as Eq. (10): 5
Applied Energy 262 (2020) 114418
G. Pan, et al.
D=
1 d |ri (t ) 6 dt
ri (0)|2
revealing the characteristic of short-range order and long-range disorder. For the unlike ion pair, the first peak of RDF is sharper indicating more stable coordination structure. When temperature rises, the height and width of first peak of all RDFs decreases and broadens, corresponding to worse thermal stabilities. Average distances between C4+ and O2− hardly change because of strong covalent bond. Meanwhile, the bond angle of OeCeO has no change indicating that spatial tetrahedral structure of CO32− is quite stable. When temperature rises, the average ionic equilibrium distance rmax of M+-C4+ shifts to the left while the average ionic equilibrium distance rmax of C4+-C4+ shifts to the right. It can be concluded that CO32− group centered ion clusters are formed in molten Na2CO3K2CO3. The distances of ions decrease within ion clusters, while the distances between ion clusters increase with rising temperature. A RDF sample at 1200 K of different components is plotted in Fig. 2. With decreasing mole fraction of Na2CO3, the average ionic equilibrium distance rmax of C4+-C4+ shifts to the right while no significant changes are found for other ion pairs. With decreasing mole fraction of Na2CO3, the height of first peak of C4+-C4+ decreases slightly while the heights of first peak of other ion pairs increase; the height of first peak valley of C4+-C4+ increases slightly while the heights of first peak valley of other ion pairs decrease. For the sake of further studying effect of component on microstructure of Na2CO3-K2CO3, coordination number curves will be discussed in the next section with RDFs.
(10)
where the angle brackets 〈···〉 is the ensemble average and ri is the position of ion i. 3.2.4. Thermal conductivity Thermal conductivity is determined by RNEMD method based on linear response theory in this work:
jz =
T z
(11)
The simulation box is artificially equally divided into N slabs (Fig. 1) along the direction of interest (z dimension). At the start, a Nosé-Hoover thermostat is applied and the systems are equilibrated at target temperature for 1 ns. Then the thermostat is canceled and kinetic energies are switched unphysically between the bottom (1th) and middle (N/2 + 1th) slab in the direction of interest at a certain frequency in the microcanonical ensemble (NVE). The transferred kinetic energy is accumulated over at least 5 ns. Once the system achieves dynamic equilibrium (i.e. the temperature profile is linear), the temperature gradient in the direction of interest is monitored and averaged. Thermal conductivity can be derived as follow:
=
transfer
m 2 (vhot 2
2tL x L y
2 vcold ) T
z
(12)
5.1.2. Coordination number It is observed from Table 4 that all coordination numbers decrease with rising temperature. Coordination number curves of Na+-C4+ and K+-C4+ have oblivious inflection point compared to other ion pairs, which indicate more stable coordination structures. A sample of coordination number curve at 1200 K of different
where the denominator denotes total transferred kinetic energy during time t and the angle brackets 〈···〉 denotes the ensemble average. 4. Simulation details All simulations are performed by use of the open source molecular simulation package LAMMPS [22]. All simulated systems consist of 9216 ions (see Table 2). The crystal models of binary Na2CO3-K2CO3 mixtures of different components are built by replacing cations of pure Na2CO3. The periodic boundary conditions method (PBC) is employed and the interionic interactions are truncated at 20 Å which is slightly smaller than half of side length. Long-range interactions are dealt in reciprocal space with pppm method and the precision is set to 1.0 × 10−6. The initial velocities are assigned obeying the Gaussian distribution. The Verlet algorithm is used to solve Newton’s equations of motion and time step is set to 1.0 fs. For molten carbonates of different components, the systems are firstly heated to a very high temperature (e.g. 2500 K). Then the system is cooled slowly to the target temperatures follow by an equilibrium process in NPT ensemble under atmospheric pressure by use of NoséHoover thermostat and barostat. The damping parameters for thermostat and barostat are set to 0.1 ps and 0.5 ps, respectively. The equilibrium process lasts at least 1 ns for each system for the sake of good statistics.
Fig. 3. Coordination number curve of molten Na2CO3-K2CO3 mixtures of different components at 1200 K.
5. Results and discussions 5.1. Local structures
Table 5 Coordination numbers of molten Na2CO3-K2CO3 at 1200 K.
RDFs, coordination number curves and ADFs of molten Na2CO3K2CO3 mixtures of different components are computed in order to study temperature and component dependence of microstructures of molten binary carbonate salts. The results are showed and analyzed in the following section.
S-1 S-2 S-3 S-4 S-5 S-6 S-7
5.1.1. Partial radial distribution function It is observed from Table 3 that all RDFs of molten Na2CO3-K2CO3 mixtures have a high and narrow first peak. Then the vibration amplitude after the first peak decreases gradually with increasing distance 6
Na-Na
Na-K
Na-C
K-K
K-C
C-C
6.21 5.64 4.28 3.37 2.55 1.59 0.70
1.04 1.95 3.15 4.19 5.06 6.08 7.01
3.39 3.34 3.21 3.05 2.98 2.74 2.67
1.08 1.90 3.02 3.89 4.75 5.96 6.72
5.01 4.60 4.44 4.37 4.12 3.97 3.82
13.19 13.70 13.48 13.44 13.71 13.31 13.61
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 6 Angular distribution functions of molten Na2CO3-K2CO3 mixtures.
between Na+-C4+ and K+-C4+ and K+-C4+ is easier to be formed. Fig. 2 showcases that the heights of first peak of Na+-C4+ and K+-C4+ decrease while the heights of first peak valley of Na+-C4+ and K+-C4+ increase with increasing mole fraction of Na2CO3. The comprehensive consequence is that NNa-C and NK-C increase slightly.
components is plotted in Fig. 3. Corresponding coordination numbers of molten Na2CO3-K2CO3 at 1200 K are summarized in Table 5. With increasing mole fraction of Na2CO3, coordination numbers of Na+-K+, K+-K+ decrease while that of Na+-Na+, Na+-C4+ and K+-C4+ increase. It's worth noting that NK-C also increases with increasing mole fraction of Na2CO3, which indicates that there is competitive binding 7
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 7 Heat capacities of molten Na2CO3-K2CO3 with different component. Salt
Cp (simulation)/J·mol−1·K−1
Cp (experiment)/J·mol−1·K−1
Pure Na2CO3 S-1 S-2 S-3 S-4 S-5 S-6 S-7 Pure K2CO3
178.7 179.3 180.2 180.3 180.9 181.0 181.5 180.1 177.9
187.1 — — — — — — — 199.2
Fig. 4. ADF sample at 1200 K for molten Na2CO3-K2CO3 mixtures of different components.
5.1.3. Angular distribution function Table 6 is angular distribution functions of molten Na2CO3-K2CO3 of different components. The bond angles of O-C-O are distributed between 30 and 150 degree revealing that C atom and O atom are not in the same plane and CO32− is spatial tetrahedral structure. Significant changes are hardly found when temperature increases continually. An ADF sample at 1200 K of different components is plotted in Fig. 4 and the result showcases that no significant changes happen with the different components. Fig. 6. The variation of specific heat capacity with temperature.
5.2. Thermophysical properties
linear temperature dependences for densities of molten Na2CO3-K2CO3 are observed and this can be explained in conjunction with variation of microstructure: with rising temperature, distances of ion clusters increase leading to bigger volume. When temperature is constant, densities are linearly related with components. Three dimension cloud map of the temperature-density-composition correlation is plotted as Fig. 5b. The fit of density with temperature and component is shown as Eq. (13),
5.2.1. Densities For the computation of density, a wide operating temperature range is taken into consideration. The simulated and experimental [23] densities of liquid phase are compared in Fig. 5. The simulated densities of molten Na2CO3-K2CO3 of different components are slightly underestimated within 11.2% of experimental value. The results coincide with pure alkali carbonates of our previous works [19,20]. The errors arise from the potential model of selection since the model takes no account of dispersion forces (Cij and Dij) which have characteristics of attraction. On account of acceptable discrepancies, the potential model is chosen to calculate the other thermophysical properties. Negative
mixture
( a )
= 2.1435
3.9963 × 10 4T + 5.7870 × 10 2xNa2 CO3
( b)
Fig. 5. Temperature and component of density dependence of molten Na2CO3-K2CO3: (a) Temperature (b) Composition and Temperature. 8
(13)
Applied Energy 262 (2020) 114418
G. Pan, et al.
Fig. 7. MSDs of molten Na2CO3-K2CO3 (58–42 mol%) (a) of Na+. (b) K+. (c) CO32−.
Fig. 8. Self-diffusion coefficients of molten Na2CO3-K2CO3 (58–42 mol%) (a) of Na+. (b) K+. (c) CO32−.
5.2.2. Specific heat capacity Specific heat capacities are calculated through the procedure described in Section 3.2.2. A series of simulations in NPT ensembles over a range of temperatures are carried out. Table 7 showcases the specific heat capacities of molten Na2CO3-K2CO3 of different components, plus simulation and experimental value of two pure alkali carbonates [24]. No obvious change is found with respect to component. The variation of specific heat capacity with temperature are computed and shown as Fig. 6. Weak Negative linear temperature dependence of specific heat capacity is observed. The slope is so small that we can suggest specific heat capacity of Na2CO3-K2CO3 be a constant value regardless of temperature.
5.2.3. Self-diffusion coefficient The ion self-diffusion coefficients of molten Na2CO3-K2CO3 of different components are calculated in order to study the ion transport properties by means of Einstein equation based on time-dependent MSD mentioned above in Section 3.2.3. Fig. 7 exhibits the MSDs of Na+, K+ and CO32− of molten Na2CO3-K2CO3 (58–42 mol%) from 1000 K to 1400 K with an interval of 50 K. The MSD curves are nonlinear at the beginning and become linear after certain time, i.e. the characteristic time, which are around 0.5 and 1.0 ps for cations and anions, respectively. Fig. 8 exhibits the corresponding self-diffusion coefficients. Selfdiffusion coefficients increase linearly with rising temperature are observed, which is in agreement with Spedding and Mills’s experiment [25]. Besides, the self-diffusion coefficients of cations are larger than 9
Applied Energy 262 (2020) 114418
G. Pan, et al.
Fig. 9. MSDs of molten Na2CO3-K2CO3 of different components at 1250 K (a) of Na+. (b) K+. (c) CO32−.
CO32− since Na+ and K+ are single atom entity with small mass while CO32− is multi-atom entity with large mass. The movement of CO32− is strongly limited by intramolecular interaction. In order to study component dependence of self-diffusion coefficient of molten Na2CO3-K2CO3 mixtures, MSDs and self-diffusion coefficients of different components at 1250 K are plotted as Figs. 9 and 10. It is observed from Fig. 10 that self-diffusion coefficients of all ions depend on composition and self-diffusion coefficients of Na+ and K+ increase linearly while self-diffusion coefficient of CO32− firstly increases and achieves maximum at the eutectic composition ( Na2 CO3 = 0.5), then decreases with increasing mole fraction of Na2CO3. The fit of self-diffusion coefficients with components at 1250 K are shown as Eqs. (14)–(16),
Fig. 10. Self-diffusion coefficients of molten Na2CO3-K2CO3 of different components at 1250 K (a) of Na+. (b) K+. (c) CO32−.
DNa+ = 2.94457 + 0.02782XNa2 CO3
(14)
D K + = 6.69685 + 0.01424XNa2 CO3
(15)
DCO32 = 2.54323 + 0.02720XNa2 CO3
2 2.51763 × 10 4XNa 2 CO3
(16)
5.2.4. Thermal conductivity Thermal conductivities are determined by RNEMD method mentioned above in Section 3.2.4. For molten carbonates, kinetic energy swap rate between 40 and 200 time steps yields fairly linear 10
Applied Energy 262 (2020) 114418
G. Pan, et al.
(a) S-2
(b) S-4
(c) S-6
Fig. 11. Temperature profiles of molten Na2CO3-K2CO3 of different components with kinetic energy swap rates equal to 200 time steps.
Fig. 12. Thermal conductivities of molten Na2CO3-K2CO3 mixture vs. (a) temperature and (b) component.
temperature profile, which obeys linear response theory and can obtain accurate result. The resulting temperature profiles of molten Na2CO3K2CO3 of different components with kinetic energy swap rates equal to 200 time steps are plotted in Fig. 11. Since there is no experimental data of thermal conductivity for molten Na2CO3-K2CO3 mixtures, simulated thermal conductivities for different components are compared between simulated and experimental value [26] of pure alkali carbonates (Fig. 12a). Negative linear temperature dependences for thermal conductivities are observed and this can be explained in conjunction with variation of microstructure: when temperature increases, the distances of ion clusters become larger leading to obstruction of heat transfer through vibration between ions. When temperature is constant, thermal conductivities are linearly related with components. Three dimension cloud map is plotted as Fig. 12b standing for the temperature-thermal conductivities-composition correlation. The fit of thermal conductivity with temperature and component is shown as Eq. (17),
= 0.4752
6.9533 × 10 5T + 1.9548 × 10 1xNa2 CO3
consisting of harmonic bond, angle, improper terms. The properties of molten Na2CO3-K2CO3 are assessed including density, specific heat capacity, self-diffusion coefficient and thermal conductivity. The results show that the selected potential model and parameters are suitable to reproduce the interactions of alkali carbonates. Temperature and component dependence of microstructures and thermophysical properties of molten Na2CO3-K2CO3 mixtures are studied. The calculated densities are slightly underestimated within 11.2% of the experimental value which coincide with pure alkali carbonates of our previous works. Negative linear temperature dependence of densities and thermal conductivities of binary mixtures of different components is confirmed since the distances of ion clusters increase when temperature increases, while positive linear temperature dependences for self-diffusion coefficient are obtained. The results show that specific heat capacities keep constant with temperature. Densities and thermal conductivities of molten Na2CO3-K2CO3 mixtures are linearly related with components. Self-diffusion coefficients of Na+ and K+ increase linearly while self-diffusion coefficients of CO32− firstly increase and then decrease with increasing mole fraction of Na2CO3. Reverse nonequilibrium molecular dynamics method is suitable for determining thermal conductivity of molten Na2CO3-K2CO3. The temperature-thermophysical properties-composition correlation formulas are obtained, and the database of thermophysical properties of molten carbonate salts over the entire operating temperature extent is complemented, which will provide the essential data for system design, operation, and optimization in CSP.
(17)
6. Conclusions In this paper, molecular dynamics simulations of molten Na2CO3K2CO3 mixtures of different components is conducted by using Borntype pair potential plus flexible carbonate intramolecular interaction 11
Applied Energy 262 (2020) 114418
G. Pan, et al.
CRediT authorship contribution statement Gechuanqi Pan: Methodology, Software, Data curation, Writing original draft. Xiaolan Wei: Formal analysis. Chao Yu: Writing - review & editing. Yutong Lu: Supervision. Jiang Li: Software. Jing Ding: Conceptualization, Supervision. Weilong Wang: Writing - review & editing. Jinyue Yan: Writing - review & editing.
[8] [9] [10]
Declaration of Competing Interest
[11]
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
[12] [13]
Acknowledgements
[14]
This work is supported by the funding of Nature Science Foundation of China (U1707603), and Nature Science Foundation of China (51436009), Science and Technology Planning Project of Guangdong Province (2015A010106006), and Nature Science Foundation of Guangdong (2016A030313362).
[15]
References
[19]
[16] [17] [18]
[1] Navarrete N, Mondragón R, Wen D, Navarro ME, Ding Y, Juliá JE. Thermal energy storage of molten salt–based nanofluid containing nano-encapsulated metal alloy phase change materials. Energy 2019;167:912–20. [2] Olivares RI, Chen C, Wright S. The thermal stability of molten lithium–sodium–potassium carbonate and the influence of additives on the melting point. J Sol Energ-T ASME 2012;134:041002. [3] Hao W, He X, Mi Y. Achieving high performance in intermediate temperature direct carbon fuel cells with renewable carbon as a fuel source. Appl Energ 2014;135:174–81. [4] Sang L, Cai M, Ren N, Wu Y, Burda C, Ma C. Improving the thermal properties of ternary carbonates for concentrating solar power through simple chemical modifications by adding sodium hydroxide and nitrate. Sol Energy Mater Sol Cells 2014;124:61–6. [5] Yang Z, Liao T, Zhou Y, Lin G, Chen J. Performance evaluation and parametric optimum design of a molten carbonate fuel cell-thermophotovoltaic cell hybrid system. Energy Convers Manage 2016;128:28–33. [6] Brun CL. Molten salts and nuclear energy production. J Nucl Mater 2007;360:1–5. [7] Griffiths TR, Volkovich VA, Yakimov SM, May I, Sharrad CA, Charnock JM.
[20] [21] [22] [23] [24] [25] [26]
12
Reprocessing spent nuclear fuel using molten carbonates and subsequent precipitation of rare earth fission products using phosphate. J Alloys Compd 2006;418:116–21. Griffiths TR, Volkovich VA. A new technology for the nuclear industry for the complete and continuous pyrochemical reprocessing of spent nuclear fuel: Catalyst enhanced molten salt oxidation. Nucl Technol 2008;163:382–400. Galamba N, Ca NDC, Ely JF. Equilibrium and nonequilibrium molecular dynamics simulations of the thermal conductivity of molten alkali halides. J Chem Phys 2007;126:1879. Kowsari MH, Alavi S, Ashrafizaadeh M, Najafi B. Molecular dynamics simulation of imidazolium-based ionic liquids. II. Transport coefficients. J Chem Phys 2009;130:147. Janssen GJM, Tissen JTWM. Pair potentials from ab initio calculations for use in MD simulations of molten alkali carbonates. Mol Simulat 1990;5:83–98. Tissen JTWM, Janssen GJM. Molecular-dynamics simulation of molten alkali carbonates. Mol Phys 1990;71:413–26. Tissen JTWM, Janssen GJM, Dereerden PV. Molecular dynamics simulation of binary mixtures of molten alkali carbonates. Mol Phys 1990;82:101–11. Archer TD, Birse SEA, Dove MT, Redfern SAT, Gale JD, Cygan RT. An interatomic potential model for carbonates allowing for polarization effects. Phys Chem Miner 2003;30:416–24. Fisler DK, Gale JD, Cygan RT. A shell model for the simulation of rhombohedral carbonate minerals and their point defects. Am Mineral 2000;85:217–24. Costa MF. Molecular dynamics of molten Li2CO3–K2CO3. J Mol Liq 2008;138:61–8. Ottochian A, Ricca C, Labat F, Adamo C. Molecular dynamics simulations of a lithium/sodium carbonate mixture. J Mol Model 2016;22:61. Corradini D, Coudert F, Vuilleumier R. Insight into the Li2CO3–K2CO3 eutectic mixture from classical molecular dynamics: thermodynamics, structure, and dynamics. J Chem Phys 2016;144:1363–5. Ding J, Du L, Pan G, Lu J, Wei X, Li J, et al. Molecular dynamics simulations of the local structures and thermodynamic properties on molten alkali carbonate K2CO3. Appl Energy 2018;220:536–44. Ding J, Pan G, Du L, Lu J, Wang W, Wei X, et al. Molecular dynamics simulations of the local structures and transport properties of Na2CO3 and K2CO3. Appl Energy 2018;227:555–63. Koishi T, Si Kawase, Tamaki S, Ebisuzaki T. Computer simulation of molten Li2CO3K2CO3 mixtures. J Phys Soc Jpn 2000;69:3291–6. Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117:1–19. Spedding PL. Densities and molar volumes of molten alkali carbonate binary mixtures. J Electrochem Soc 1970;117. Janz GJ, Neuenschwander E, Kelly FJ. High-temperature heat content and related properties for Li2CO3, Na2CO3, K2CO3, and the ternary eutectic mixture. T Faraday Soc 1963;59:841–5. Spedding PL, Mills R. Tracer diffusion measurements in mixtures of molten alkali carbonates. J Electrochem Soc 1966;113. Gheribi AE, Torres JA, Chartrand P. Recommended values for the thermal conductivity of molten salts between the melting and boiling points. Sol Energy Mater Sol Cells 2014;126:11–25.
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Thermal performance of a binary carbonate molten eutectic salt for hightemperature energy storage applications
T
Gechuanqi Pana,b,c,e, Xiaolan Weid, Chao Yua,b,c, Yutong Lua,b,c,e, Jiang Lie, Jing Dinga,b,c, Weilong Wanga,b,c, Jinyue Yanf,g a
School of Data and Computer Science, Sun Yat-Sen University, Guangzhou, PR China School of Materials Science and Engineering, Sun Yat-Sen University, Guangzhou, PR China c School of Intelligent Systems Engineering, Sun Yat-Sen University, Guangzhou, PR China d School of Chemistry and Chemical Engineering, South China University of Technology, Guangzhou, PR China e National Supercomputer Center in Guangzhou, Guangzhou, PR China f Division of Energy Processes, Royal Institute of Technology, Stockholm, Sweden g School of Business, Society and Energy, Mälardalen University, Västerås, Sweden b
HIGHLIGHTS
GRAPHICAL ABSTRACT
simulation is used to com• Molecular pute structures-properties of Na CO 2
• •
3
K2CO3. Temperature and component dependence are studied from microscopic view. The temperature-thermophysical properties-composition correlation are obtained.
ARTICLE INFO
ABSTRACT
Keywords: Energy storage Molten carbonates Thermophysical properties Microstructures Concentrating solar power
Molten carbonate eutectic salts are promising thermal storage and heat transfer fluid materials in solar thermal power plant with the feature of large specific heat capacity, wide operating temperature range and little corrosive. The high-temperature properties of molten carbonates should be determined accurately over the entire operating temperature for energy system design. In this paper, molecular dynamic simulation is used to study temperature and component dependence of microstructures and thermophysical properties of the binary carbonate molten salt. Negative linear temperature dependence of densities and thermal conductivities of binary mixtures of different components is confirmed with respect to the distances of ion clusters. Besides, positive linear temperature dependence of self-diffusion coefficient is also obtained. When temperature is constant, densities and thermal conductivities of binary mixtures are linearly related with components. Self-diffusion coefficients of CO32− firstly increase and then decrease with increasing mole fraction of Na2CO3. The temperature-thermophysical properties-composition correlation formulas are obtained, and the database of thermophysical properties of molten carbonate salts over the entire operating temperature is complemented, which will provide the essential data for heat transfer and storage system design, operation, and optimization in CSP.
E-mail addresses: [email protected] (J. Ding), [email protected] (W. Wang). https://doi.org/10.1016/j.apenergy.2019.114418 Received 24 October 2019; Received in revised form 16 December 2019; Accepted 17 December 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.
Applied Energy 262 (2020) 114418
G. Pan, et al.
Nomenclature
z length dimension in the z direction Kb, Kθ, Kφ, r0, θ0, φ0 bond, angle, improper parameters length dimension in the x and y direction Lx , Ly v velocity m atomic mass λ thermal conductivity V equilibrium volume NA Avogadro’s constant
born-type pair potential Uij qi, qj charge of ion distance between particle i and j rij Aij, ρ, σij born repulsive parameters Cij, Dij Van der Waals parameters Ubond, Uangle, Uimproper bond, angle, improper potential t simulation time N particle number M molar mass i particle i j particle J k slab k x x dimension hot hot slab T temperature
Subscripts b θ ψ y cold
bond term angle term improper term y dimension cold slab
1. Introduction Molten carbonates, which have large specific heat capacity, wide operating temperature range and little corrosive, have been considered as very promising thermal storage and heat transfer fluid materials [1,2] in solar thermal power plant [3,4], fuel cell [5,6], nuclear fuel reprocessing etc. [7,8]. For efficient utilization of high-temperature heat, thermodynamic properties and transport properties of molten carbonates should be determined accurately over the entire operating temperature extent for device and system design in relevant industry. Some properties are hard or impossible to be obtained by means of experiments because of high-temperature and corrosion problems. Large discrepancies among the experimental value of thermal conductivity have been observed. Moreover, components of carbonates are varying due to evaporation, thermal decomposition, etc.. Thus, it is desired to study the temperature and component dependence of microstructures and thermophysical properties. MD has been demonstrated to be a powerful tool for the determination of the microscopic and macroscopic properties of molten salts in recent years [9,10]. Tissen and Janssen [11,12] determined the parameters of pair interaction of carbonates including electrostatic interaction, repulsion interaction and induction energy at the Hartree-Fock SCF level and the fitting parameters represented structures and diffusion coefficients well although dispersion interaction was ignored. They observed significant changes of structures from molten Li2CO3 to K2CO3 [13]. Archer [14] created an empirical model for CaCO3 which was constructed by fitting to structures and properties of aragonite and calcite. The new empirical model that the instability of the model of Fisler [15] was removed was stable for a wide range of carbonate structures, and reproduced experimental results with a reasonable accuracy. Costa [16] used the fluctuating charge model to performed simulations Li2CO3 and K2CO3, whose parameters were fitted by ab initio
Fig. 1. RNEMD method schematic view.
Table 2 Numbers of ions in simulated Na2CO3-K2CO3 of different components. Mole fraction of Na2CO3
S-1 S-2 S-3 S-4 S-5 S-6 S-7
0.87 0.75 0.62 0.50 0.38 0.25 0.13
Numbers of ions total
Na+
K+
C4+
O2−
9216
2673 2304 1905 1536 1167 768 399
399 768 1167 1536 1905 2304 2673
1536
4608
Table 1 Parameters of intermolecular and intramolecular interaction for Na2CO3-K2CO3. Atom
qi (e)
Aii (Kcal·mol−1)
ρii (Å)
σ (Å)
Cii (Å6·Kcal·mol−1)
Dii (Å8·Kcal·mol−1)
Na K C O
1.00 1.00 1.54 −1.18
6.0791 6.0791 10.9523 3.2648
0.2899
1.14 1.39 1.10 1.33
0.00
0.00
Atoms CeO OeCeO OeCeOeO
Kb (Kcal·mol−1 ·Å−2) = 730.641 Kθ (Kcal·mol−1 ·rad−2) = 20.607 Kφ (Kcal·mol−1 ·rad−2) = 6.959
2
r0(Å) = 1.16 θo (deg) = 120.0 φo (deg) = 0.00
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 3 RDFs of molten Na2CO3-K2CO3 of different components.
calculations of an isolated carbonate anion. Ottochian used Tissen and Janssen’s model (JT model) to study electric conductivity and diffusion coefficients of Li2CO3-Na2CO3 mixtures and found diffusion of lithium ions depended on its concentration level with respect to sodium ions in
the mixtures [17]. Corradini [18] parameterized a non-polarizable force field for Li2CO3-K2CO3 mixture by fitting densities and radial distribution functions to experiment and first principle molecular dynamics and presented an analysis of the temperature and pressure 3
Applied Energy 262 (2020) 114418
G. Pan, et al.
dependence of transport properties of the Li2CO3-Na2CO3 mixtures including self-diffusion coefficients, viscosity, and ionic conductivity. Recently, we have used the JT model to simulate local structures and thermophysical properties of pure alkali carbonate and obtained the variation law of microstructures and thermophysical properties with temperature [19,20]. However, to our best knowledge, thermal properties of molten carbonate salts are quite limited, and especially there is no comprehensive database for the entire operation temperature for heat transfer system design, operation and optimization in a concentrating solar power (CSP) station. Since it is hard to obtain all the properties due to the temperature and acuuracy limitation of measuring devices, it is necessary to obtain the effects of temperature and component on microstructures and thermophysical properties of molten salts under high-temperature conditions. This paper aims to gain temperature-thermophysical propertiescomposition correlation and complement the database of thermophysical properties of molten carbonate salts. Reverse nonequilibrium molecular dynamics (RNEMD) method is implemented for determining thermal conductivity. Moreover, the temperature and component dependence is investigated by correlating transport properties with local structures.
3.1.2. Coordination number Coordination number is the average number of β-type ions locating 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 by integrating the RDF as follow:
N
jik
+ Aij exp
rij
ij
rij
Cij
Dij
rij6
rij8
Uimproper = K (
1
using
the
following 1
Cij = (Cii Cjj ) 2 , Dij = (Dii Djj ) 2 and
ij
mixing
= 2(
ii
1
+
rules: jj
2 r jk
2rij rik
(7)
i Na, K , C , O
(8)
VE NA
3.2.2. Specific heat capacity The specific heat capacity can be determined simultaneously with density and computed from enthalpy using Eq. (9):
Cp = (
The JT model has been demonstrated to be able to reproduce the properties of the molten alkali carbonates [21]. The parameters of intermolecular interaction used in this work are the same to Tissen and Janssen [11,12], while the bond parameters are referred to Corradini [18], and angle and improper parameters are referred to Archer [14]. The detailed values are shown in Table 1. For the Born-type pair potential (Eq. (1)), the cross terms were
computed
rij2 + rik2
1
where Ni and Mi stand for the particle number and the molar mass of species i, while NA denotes Avogadro’s constant, VE corresponds to the average equilibrium volume.
(4)
2
cos
Ni Mi
(3) 0)
=
=
(2)
2 0)
Uangle = K (
(6)
3.2.1. Density The density of molten carbonate mixture can be determined simply by implementing an isothermal-isobaric ensemble (NPT) running and calculated as follow:
(1)
r0 ) 2
Ubond = Kb (r
g (r ) r 2dr
3.2. Thermophysical properties
Tissen and Janssen’s non-polarizable model is chosen for simulations of Na2CO3-K2CO3 mixture. The JT model is represented by an intermolecular interaction (born-type, Eq. (1)) plus a flexible carbonate intramolecular interaction (Eqs. (2)–(4)). The first two terms of right of Eq. (1) stand for columbic potential and overlapping repulsion while the last two terms correspond to dispersion energies.
qi qj
0
3.1.3. Angular distribution function In order to obtain the image of local structure of molten carbonates, angular distribution function (ADF), which can reflect the bond orientation, also should be obtained. ADF depends only on the coordinates of atoms exported during the simulations. For a centered atom i, and its two nearest neighbor atoms j and k, the bond angle can be calculated as follow:
2. Forcefield
Uij =
rmin
=4
H )P T
(
H )P = T
(U + PV ) T P
(9)
where U is the total energy of the system comprising kinetic energy and potential energy. Thus, specific heat capacity is the slope of enthalpy versus temperature.
1
Aij = (Aii Ajj ) 2 ,
) 1.
1
3. Evaluated properties and methods 3.1. Structure information 3.1.1. Partial radial distribution function Partial radial distribution function (RDF) is often used to characterize the local structure of molten salts, which is defined as Eq. (5):
g (r ) =
dN (r )
1 4
r2
dr
(5)
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 position of first peak of RDF is the average ionic equilibrium distance rmax between species α and β.
Fig. 2. RDF sample of molten Na2CO3-K2CO3 mixtures of different components at 1200 K. 4
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 4 Coordination number curves of molten Na2CO3-K2CO3 of different components.
3.2.3. Self-diffusion coefficient In order to study the ion transport properties of molten Na2CO3K2CO3 mixture, the ion self-diffusion coefficients are calculated. The ion self-diffusion coefficients can be calculated by either the time-
integration of the velocity autocorrelation function (ACF) or the Einstein equation based time-dependent mean-squared displacement (MSD). The latter method is employed here. The corresponding calculation formula is expressed as Eq. (10): 5
Applied Energy 262 (2020) 114418
G. Pan, et al.
D=
1 d |ri (t ) 6 dt
ri (0)|2
revealing the characteristic of short-range order and long-range disorder. For the unlike ion pair, the first peak of RDF is sharper indicating more stable coordination structure. When temperature rises, the height and width of first peak of all RDFs decreases and broadens, corresponding to worse thermal stabilities. Average distances between C4+ and O2− hardly change because of strong covalent bond. Meanwhile, the bond angle of OeCeO has no change indicating that spatial tetrahedral structure of CO32− is quite stable. When temperature rises, the average ionic equilibrium distance rmax of M+-C4+ shifts to the left while the average ionic equilibrium distance rmax of C4+-C4+ shifts to the right. It can be concluded that CO32− group centered ion clusters are formed in molten Na2CO3K2CO3. The distances of ions decrease within ion clusters, while the distances between ion clusters increase with rising temperature. A RDF sample at 1200 K of different components is plotted in Fig. 2. With decreasing mole fraction of Na2CO3, the average ionic equilibrium distance rmax of C4+-C4+ shifts to the right while no significant changes are found for other ion pairs. With decreasing mole fraction of Na2CO3, the height of first peak of C4+-C4+ decreases slightly while the heights of first peak of other ion pairs increase; the height of first peak valley of C4+-C4+ increases slightly while the heights of first peak valley of other ion pairs decrease. For the sake of further studying effect of component on microstructure of Na2CO3-K2CO3, coordination number curves will be discussed in the next section with RDFs.
(10)
where the angle brackets 〈···〉 is the ensemble average and ri is the position of ion i. 3.2.4. Thermal conductivity Thermal conductivity is determined by RNEMD method based on linear response theory in this work:
jz =
T z
(11)
The simulation box is artificially equally divided into N slabs (Fig. 1) along the direction of interest (z dimension). At the start, a Nosé-Hoover thermostat is applied and the systems are equilibrated at target temperature for 1 ns. Then the thermostat is canceled and kinetic energies are switched unphysically between the bottom (1th) and middle (N/2 + 1th) slab in the direction of interest at a certain frequency in the microcanonical ensemble (NVE). The transferred kinetic energy is accumulated over at least 5 ns. Once the system achieves dynamic equilibrium (i.e. the temperature profile is linear), the temperature gradient in the direction of interest is monitored and averaged. Thermal conductivity can be derived as follow:
=
transfer
m 2 (vhot 2
2tL x L y
2 vcold ) T
z
(12)
5.1.2. Coordination number It is observed from Table 4 that all coordination numbers decrease with rising temperature. Coordination number curves of Na+-C4+ and K+-C4+ have oblivious inflection point compared to other ion pairs, which indicate more stable coordination structures. A sample of coordination number curve at 1200 K of different
where the denominator denotes total transferred kinetic energy during time t and the angle brackets 〈···〉 denotes the ensemble average. 4. Simulation details All simulations are performed by use of the open source molecular simulation package LAMMPS [22]. All simulated systems consist of 9216 ions (see Table 2). The crystal models of binary Na2CO3-K2CO3 mixtures of different components are built by replacing cations of pure Na2CO3. The periodic boundary conditions method (PBC) is employed and the interionic interactions are truncated at 20 Å which is slightly smaller than half of side length. Long-range interactions are dealt in reciprocal space with pppm method and the precision is set to 1.0 × 10−6. The initial velocities are assigned obeying the Gaussian distribution. The Verlet algorithm is used to solve Newton’s equations of motion and time step is set to 1.0 fs. For molten carbonates of different components, the systems are firstly heated to a very high temperature (e.g. 2500 K). Then the system is cooled slowly to the target temperatures follow by an equilibrium process in NPT ensemble under atmospheric pressure by use of NoséHoover thermostat and barostat. The damping parameters for thermostat and barostat are set to 0.1 ps and 0.5 ps, respectively. The equilibrium process lasts at least 1 ns for each system for the sake of good statistics.
Fig. 3. Coordination number curve of molten Na2CO3-K2CO3 mixtures of different components at 1200 K.
5. Results and discussions 5.1. Local structures
Table 5 Coordination numbers of molten Na2CO3-K2CO3 at 1200 K.
RDFs, coordination number curves and ADFs of molten Na2CO3K2CO3 mixtures of different components are computed in order to study temperature and component dependence of microstructures of molten binary carbonate salts. The results are showed and analyzed in the following section.
S-1 S-2 S-3 S-4 S-5 S-6 S-7
5.1.1. Partial radial distribution function It is observed from Table 3 that all RDFs of molten Na2CO3-K2CO3 mixtures have a high and narrow first peak. Then the vibration amplitude after the first peak decreases gradually with increasing distance 6
Na-Na
Na-K
Na-C
K-K
K-C
C-C
6.21 5.64 4.28 3.37 2.55 1.59 0.70
1.04 1.95 3.15 4.19 5.06 6.08 7.01
3.39 3.34 3.21 3.05 2.98 2.74 2.67
1.08 1.90 3.02 3.89 4.75 5.96 6.72
5.01 4.60 4.44 4.37 4.12 3.97 3.82
13.19 13.70 13.48 13.44 13.71 13.31 13.61
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 6 Angular distribution functions of molten Na2CO3-K2CO3 mixtures.
between Na+-C4+ and K+-C4+ and K+-C4+ is easier to be formed. Fig. 2 showcases that the heights of first peak of Na+-C4+ and K+-C4+ decrease while the heights of first peak valley of Na+-C4+ and K+-C4+ increase with increasing mole fraction of Na2CO3. The comprehensive consequence is that NNa-C and NK-C increase slightly.
components is plotted in Fig. 3. Corresponding coordination numbers of molten Na2CO3-K2CO3 at 1200 K are summarized in Table 5. With increasing mole fraction of Na2CO3, coordination numbers of Na+-K+, K+-K+ decrease while that of Na+-Na+, Na+-C4+ and K+-C4+ increase. It's worth noting that NK-C also increases with increasing mole fraction of Na2CO3, which indicates that there is competitive binding 7
Applied Energy 262 (2020) 114418
G. Pan, et al.
Table 7 Heat capacities of molten Na2CO3-K2CO3 with different component. Salt
Cp (simulation)/J·mol−1·K−1
Cp (experiment)/J·mol−1·K−1
Pure Na2CO3 S-1 S-2 S-3 S-4 S-5 S-6 S-7 Pure K2CO3
178.7 179.3 180.2 180.3 180.9 181.0 181.5 180.1 177.9
187.1 — — — — — — — 199.2
Fig. 4. ADF sample at 1200 K for molten Na2CO3-K2CO3 mixtures of different components.
5.1.3. Angular distribution function Table 6 is angular distribution functions of molten Na2CO3-K2CO3 of different components. The bond angles of O-C-O are distributed between 30 and 150 degree revealing that C atom and O atom are not in the same plane and CO32− is spatial tetrahedral structure. Significant changes are hardly found when temperature increases continually. An ADF sample at 1200 K of different components is plotted in Fig. 4 and the result showcases that no significant changes happen with the different components. Fig. 6. The variation of specific heat capacity with temperature.
5.2. Thermophysical properties
linear temperature dependences for densities of molten Na2CO3-K2CO3 are observed and this can be explained in conjunction with variation of microstructure: with rising temperature, distances of ion clusters increase leading to bigger volume. When temperature is constant, densities are linearly related with components. Three dimension cloud map of the temperature-density-composition correlation is plotted as Fig. 5b. The fit of density with temperature and component is shown as Eq. (13),
5.2.1. Densities For the computation of density, a wide operating temperature range is taken into consideration. The simulated and experimental [23] densities of liquid phase are compared in Fig. 5. The simulated densities of molten Na2CO3-K2CO3 of different components are slightly underestimated within 11.2% of experimental value. The results coincide with pure alkali carbonates of our previous works [19,20]. The errors arise from the potential model of selection since the model takes no account of dispersion forces (Cij and Dij) which have characteristics of attraction. On account of acceptable discrepancies, the potential model is chosen to calculate the other thermophysical properties. Negative
mixture
( a )
= 2.1435
3.9963 × 10 4T + 5.7870 × 10 2xNa2 CO3
( b)
Fig. 5. Temperature and component of density dependence of molten Na2CO3-K2CO3: (a) Temperature (b) Composition and Temperature. 8
(13)
Applied Energy 262 (2020) 114418
G. Pan, et al.
Fig. 7. MSDs of molten Na2CO3-K2CO3 (58–42 mol%) (a) of Na+. (b) K+. (c) CO32−.
Fig. 8. Self-diffusion coefficients of molten Na2CO3-K2CO3 (58–42 mol%) (a) of Na+. (b) K+. (c) CO32−.
5.2.2. Specific heat capacity Specific heat capacities are calculated through the procedure described in Section 3.2.2. A series of simulations in NPT ensembles over a range of temperatures are carried out. Table 7 showcases the specific heat capacities of molten Na2CO3-K2CO3 of different components, plus simulation and experimental value of two pure alkali carbonates [24]. No obvious change is found with respect to component. The variation of specific heat capacity with temperature are computed and shown as Fig. 6. Weak Negative linear temperature dependence of specific heat capacity is observed. The slope is so small that we can suggest specific heat capacity of Na2CO3-K2CO3 be a constant value regardless of temperature.
5.2.3. Self-diffusion coefficient The ion self-diffusion coefficients of molten Na2CO3-K2CO3 of different components are calculated in order to study the ion transport properties by means of Einstein equation based on time-dependent MSD mentioned above in Section 3.2.3. Fig. 7 exhibits the MSDs of Na+, K+ and CO32− of molten Na2CO3-K2CO3 (58–42 mol%) from 1000 K to 1400 K with an interval of 50 K. The MSD curves are nonlinear at the beginning and become linear after certain time, i.e. the characteristic time, which are around 0.5 and 1.0 ps for cations and anions, respectively. Fig. 8 exhibits the corresponding self-diffusion coefficients. Selfdiffusion coefficients increase linearly with rising temperature are observed, which is in agreement with Spedding and Mills’s experiment [25]. Besides, the self-diffusion coefficients of cations are larger than 9
Applied Energy 262 (2020) 114418
G. Pan, et al.
Fig. 9. MSDs of molten Na2CO3-K2CO3 of different components at 1250 K (a) of Na+. (b) K+. (c) CO32−.
CO32− since Na+ and K+ are single atom entity with small mass while CO32− is multi-atom entity with large mass. The movement of CO32− is strongly limited by intramolecular interaction. In order to study component dependence of self-diffusion coefficient of molten Na2CO3-K2CO3 mixtures, MSDs and self-diffusion coefficients of different components at 1250 K are plotted as Figs. 9 and 10. It is observed from Fig. 10 that self-diffusion coefficients of all ions depend on composition and self-diffusion coefficients of Na+ and K+ increase linearly while self-diffusion coefficient of CO32− firstly increases and achieves maximum at the eutectic composition ( Na2 CO3 = 0.5), then decreases with increasing mole fraction of Na2CO3. The fit of self-diffusion coefficients with components at 1250 K are shown as Eqs. (14)–(16),
Fig. 10. Self-diffusion coefficients of molten Na2CO3-K2CO3 of different components at 1250 K (a) of Na+. (b) K+. (c) CO32−.
DNa+ = 2.94457 + 0.02782XNa2 CO3
(14)
D K + = 6.69685 + 0.01424XNa2 CO3
(15)
DCO32 = 2.54323 + 0.02720XNa2 CO3
2 2.51763 × 10 4XNa 2 CO3
(16)
5.2.4. Thermal conductivity Thermal conductivities are determined by RNEMD method mentioned above in Section 3.2.4. For molten carbonates, kinetic energy swap rate between 40 and 200 time steps yields fairly linear 10
Applied Energy 262 (2020) 114418
G. Pan, et al.
(a) S-2
(b) S-4
(c) S-6
Fig. 11. Temperature profiles of molten Na2CO3-K2CO3 of different components with kinetic energy swap rates equal to 200 time steps.
Fig. 12. Thermal conductivities of molten Na2CO3-K2CO3 mixture vs. (a) temperature and (b) component.
temperature profile, which obeys linear response theory and can obtain accurate result. The resulting temperature profiles of molten Na2CO3K2CO3 of different components with kinetic energy swap rates equal to 200 time steps are plotted in Fig. 11. Since there is no experimental data of thermal conductivity for molten Na2CO3-K2CO3 mixtures, simulated thermal conductivities for different components are compared between simulated and experimental value [26] of pure alkali carbonates (Fig. 12a). Negative linear temperature dependences for thermal conductivities are observed and this can be explained in conjunction with variation of microstructure: when temperature increases, the distances of ion clusters become larger leading to obstruction of heat transfer through vibration between ions. When temperature is constant, thermal conductivities are linearly related with components. Three dimension cloud map is plotted as Fig. 12b standing for the temperature-thermal conductivities-composition correlation. The fit of thermal conductivity with temperature and component is shown as Eq. (17),
= 0.4752
6.9533 × 10 5T + 1.9548 × 10 1xNa2 CO3
consisting of harmonic bond, angle, improper terms. The properties of molten Na2CO3-K2CO3 are assessed including density, specific heat capacity, self-diffusion coefficient and thermal conductivity. The results show that the selected potential model and parameters are suitable to reproduce the interactions of alkali carbonates. Temperature and component dependence of microstructures and thermophysical properties of molten Na2CO3-K2CO3 mixtures are studied. The calculated densities are slightly underestimated within 11.2% of the experimental value which coincide with pure alkali carbonates of our previous works. Negative linear temperature dependence of densities and thermal conductivities of binary mixtures of different components is confirmed since the distances of ion clusters increase when temperature increases, while positive linear temperature dependences for self-diffusion coefficient are obtained. The results show that specific heat capacities keep constant with temperature. Densities and thermal conductivities of molten Na2CO3-K2CO3 mixtures are linearly related with components. Self-diffusion coefficients of Na+ and K+ increase linearly while self-diffusion coefficients of CO32− firstly increase and then decrease with increasing mole fraction of Na2CO3. Reverse nonequilibrium molecular dynamics method is suitable for determining thermal conductivity of molten Na2CO3-K2CO3. The temperature-thermophysical properties-composition correlation formulas are obtained, and the database of thermophysical properties of molten carbonate salts over the entire operating temperature extent is complemented, which will provide the essential data for system design, operation, and optimization in CSP.
(17)
6. Conclusions In this paper, molecular dynamics simulations of molten Na2CO3K2CO3 mixtures of different components is conducted by using Borntype pair potential plus flexible carbonate intramolecular interaction 11
Applied Energy 262 (2020) 114418
G. Pan, et al.
CRediT authorship contribution statement Gechuanqi Pan: Methodology, Software, Data curation, Writing original draft. Xiaolan Wei: Formal analysis. Chao Yu: Writing - review & editing. Yutong Lu: Supervision. Jiang Li: Software. Jing Ding: Conceptualization, Supervision. Weilong Wang: Writing - review & editing. Jinyue Yan: Writing - review & editing.
[8] [9] [10]
Declaration of Competing Interest
[11]
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
[12] [13]
Acknowledgements
[14]
This work is supported by the funding of Nature Science Foundation of China (U1707603), and Nature Science Foundation of China (51436009), Science and Technology Planning Project of Guangdong Province (2015A010106006), and Nature Science Foundation of Guangdong (2016A030313362).
[15]
References
[19]
[16] [17] [18]
[1] Navarrete N, Mondragón R, Wen D, Navarro ME, Ding Y, Juliá JE. Thermal energy storage of molten salt–based nanofluid containing nano-encapsulated metal alloy phase change materials. Energy 2019;167:912–20. [2] Olivares RI, Chen C, Wright S. The thermal stability of molten lithium–sodium–potassium carbonate and the influence of additives on the melting point. J Sol Energ-T ASME 2012;134:041002. [3] Hao W, He X, Mi Y. Achieving high performance in intermediate temperature direct carbon fuel cells with renewable carbon as a fuel source. Appl Energ 2014;135:174–81. [4] Sang L, Cai M, Ren N, Wu Y, Burda C, Ma C. Improving the thermal properties of ternary carbonates for concentrating solar power through simple chemical modifications by adding sodium hydroxide and nitrate. Sol Energy Mater Sol Cells 2014;124:61–6. [5] Yang Z, Liao T, Zhou Y, Lin G, Chen J. Performance evaluation and parametric optimum design of a molten carbonate fuel cell-thermophotovoltaic cell hybrid system. Energy Convers Manage 2016;128:28–33. [6] Brun CL. Molten salts and nuclear energy production. J Nucl Mater 2007;360:1–5. [7] Griffiths TR, Volkovich VA, Yakimov SM, May I, Sharrad CA, Charnock JM.
[20] [21] [22] [23] [24] [25] [26]
12
Reprocessing spent nuclear fuel using molten carbonates and subsequent precipitation of rare earth fission products using phosphate. J Alloys Compd 2006;418:116–21. Griffiths TR, Volkovich VA. A new technology for the nuclear industry for the complete and continuous pyrochemical reprocessing of spent nuclear fuel: Catalyst enhanced molten salt oxidation. Nucl Technol 2008;163:382–400. Galamba N, Ca NDC, Ely JF. Equilibrium and nonequilibrium molecular dynamics simulations of the thermal conductivity of molten alkali halides. J Chem Phys 2007;126:1879. Kowsari MH, Alavi S, Ashrafizaadeh M, Najafi B. Molecular dynamics simulation of imidazolium-based ionic liquids. II. Transport coefficients. J Chem Phys 2009;130:147. Janssen GJM, Tissen JTWM. Pair potentials from ab initio calculations for use in MD simulations of molten alkali carbonates. Mol Simulat 1990;5:83–98. Tissen JTWM, Janssen GJM. Molecular-dynamics simulation of molten alkali carbonates. Mol Phys 1990;71:413–26. Tissen JTWM, Janssen GJM, Dereerden PV. Molecular dynamics simulation of binary mixtures of molten alkali carbonates. Mol Phys 1990;82:101–11. Archer TD, Birse SEA, Dove MT, Redfern SAT, Gale JD, Cygan RT. An interatomic potential model for carbonates allowing for polarization effects. Phys Chem Miner 2003;30:416–24. Fisler DK, Gale JD, Cygan RT. A shell model for the simulation of rhombohedral carbonate minerals and their point defects. Am Mineral 2000;85:217–24. Costa MF. Molecular dynamics of molten Li2CO3–K2CO3. J Mol Liq 2008;138:61–8. Ottochian A, Ricca C, Labat F, Adamo C. Molecular dynamics simulations of a lithium/sodium carbonate mixture. J Mol Model 2016;22:61. Corradini D, Coudert F, Vuilleumier R. Insight into the Li2CO3–K2CO3 eutectic mixture from classical molecular dynamics: thermodynamics, structure, and dynamics. J Chem Phys 2016;144:1363–5. Ding J, Du L, Pan G, Lu J, Wei X, Li J, et al. Molecular dynamics simulations of the local structures and thermodynamic properties on molten alkali carbonate K2CO3. Appl Energy 2018;220:536–44. Ding J, Pan G, Du L, Lu J, Wang W, Wei X, et al. Molecular dynamics simulations of the local structures and transport properties of Na2CO3 and K2CO3. Appl Energy 2018;227:555–63. Koishi T, Si Kawase, Tamaki S, Ebisuzaki T. Computer simulation of molten Li2CO3K2CO3 mixtures. J Phys Soc Jpn 2000;69:3291–6. Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117:1–19. Spedding PL. Densities and molar volumes of molten alkali carbonate binary mixtures. J Electrochem Soc 1970;117. Janz GJ, Neuenschwander E, Kelly FJ. High-temperature heat content and related properties for Li2CO3, Na2CO3, K2CO3, and the ternary eutectic mixture. T Faraday Soc 1963;59:841–5. Spedding PL, Mills R. Tracer diffusion measurements in mixtures of molten alkali carbonates. J Electrochem Soc 1966;113. Gheribi AE, Torres JA, Chartrand P. Recommended values for the thermal conductivity of molten salts between the melting and boiling points. Sol Energy Mater Sol Cells 2014;126:11–25.