Formulation of Temperature-Dependent Thermal ... - ACS Publications

Article pubs.acs.org/JPCC

Formulation of Temperature-Dependent Thermal Conductivity of NaF, β‑Na3AlF6, Na5Al3F14, and Molten Na3AlF6 Supported by Equilibrium Molecular Dynamics and Density Functional Theory Aïmen E. Gheribi,*,† Mathieu Salanne,‡ and Patrice Chartrand† †

CRCT - Center for Research in Computational Thermochemistry, Department of Chemical Engineering, Polytechnique Montreal, Box 6079, Station Downtown, Montréal, Québec H3C 3A7, Canada ‡ Sorbonne Universités, UPMC Université Paris 06, CNRS, Laboratoire PHENIX, F-75005 Paris, France ABSTRACT: Numerical modeling of heat transport in the side ledge of a cryolite bath as a function of its thickness requires an accurate knowledge of the thermal transport properties of the cryolite bath as well as the solid phases constituting the side ledge. Unfortunately, due to experimental limitations, there is a lack of experimental data on this subject. The aim of this work is to predict these otherwise unknown properties. The thermal conductivity of major constituents of the side ledge in their solid state, namely: NaF, β− Na3AlF6 and Na5Al3F14 as well as molten Na3AlF6 are formulated as a function of temperature. To achieve this, we have performed a series of equilibrium molecular dynamics simulations (EMD) in both the NPT (isobar-isotherm) and NVT (canonical) statistical ensembles. The proposed approach is purely predictive as the ionic interaction potentials were parametrized on the basis of Density Functional Theory (DFT). The results are then compared to a theoretical model which has been shown to be a very good predictive ability for many ionic and ionocovalent compounds. For solid NaF and liquid Na3AlF6, for which data are available, the predictions were found to be in excellent agreement with experiment.

■

INTRODUCTION The production of aluminum consists basically in the dissolution of alumina in a molten cryolite bath (Na3AlF6− AlF3−CaF2−Al2O3), followed by its electrolysis. The typical working temperature of the electrolysis cells lies between 5 and 10 K above the cryolite bath liquidus temperature (∼1233 K). Consequently, a ledge of frozen electrolyte is formed on the sides of the cells during endothermic events. Controlling the side ledge thickness (a few cm) is essential in order to maintain a reasonable life span of the electrolysis cell, as the ledge acts as a protective layer against chemical attacks from the cryolite bath. The lack of experimental data on thermal transport properties may hinder awareness of the problem in many practical applications. It may limit not only the understanding of the process, but also thereby design improvement and the control of its parameters. In particular, knowledge of thermal transport properties is essential when predicting the temperature dependence of microstructural and thermo-mechanical properties, and also for understanding heat treatment and solidification process. They are directly related to the grain size and to the residual tensile strength within the microstructures. The characterization of thermal transport properties of materials require, strictly speaking, two independent properties among these four: thermal diffusivity (a), thermal effusivity (e), thermal conductivity (K) and heat capacity at constant pressure per unit volume (C̃ P). The latter is equal to the ratio CP,m/Vm, where CP,m and Vm are respectively the molar heat capacity and the molar volume. The molar heat capacity at constant pressure © XXXX American Chemical Society

can be measured accurately by traditional calorimetric methods, such as drop calorimetry, or by Differential Scanning Calorimetry (DSC). The molar volume (density) is usually measured by dilatometry or directly deduced from XRD lattice constants. Although there exist experimental methods to measure directly the thermal conductivity (steady-state method) and effusivity (transient methods), measurement by laser flash is now the standard experimental method to characterize the thermal transport properties of materials. K and e are deduced from the relationship: K = a × ρ × Cp and e = K/√a. As mentioned above, the major limitation in the numerical approach, for predicting the side ledge thickness in aluminum electrolysis cells, is lack of experimental data on thermal transport properties in the side ledge phases and microstructural parameters. The include type of microstructure, grain morphology, grain size, grain orientation and porosity. Theoretical models can then be used to predict the thermal conductivity of actual microstructures from experimental information and analytical representations of the thermal conductivity as a function of temperature for each phase constituting the system. Nowadays, in aluminum electrolysis cells, it is common to introduce in addition to AlF3, other species such as LiF to alter physical and thermal properties (melting point, density, conductivity, etc.) of the electrolyte.1 Received: August 6, 2016 Revised: September 15, 2016

A

DOI: 10.1021/acs.jpcc.6b07959 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

a very good predictive method at temperatures up to the melting point. It is applicable for ionic systems, in particular for halide salts and oxides.8,10 Later, a similar approach10 was developed for molten salts, also showing a good predictive ability for a wide variety of molten salt systems: halides, nitrates, hydroxides, carbonates, nitrites and sulfates.10 Equilibrium Molecular Dynamics can play a major role in both fundamental and applied materials science. This study is an example of this. Here, it is purely predictive matter, given the considerable lack of experimental data. Ab initio molecular dynamics (AIMD) is a purely predictive method, as it requires no experimental information. Note that in AIMD, the atoms are treated classically (Born−Oppenheimer approximation) via Newton’s equations of motion and the electrons are treated using the Density Functional Theory (DFT). AIMD requires a very large computational resources and is time-consuming, especially when thermal transport properties have to be calculated. An alternative to AIMD, while remaining purely predictive, is the so-called hybrid method. The hybrid method, initially launched by Madden et al.,15 consists in the parametrization of the interaction potentials on the basis of atomistic first−principles electronic structure calculations. The dipoles and forces are matched for ions calculated from first− principles for condensed phase ionic configurations. Consequently, this method has the advantage of being purely predictive, like AIMD, with reasonable calculation time requirement. The hybrid EMD method has already shown successful results in predicting the temperature dependentthermal conductivity of ionic compounds and mixtures in both solid and molten salts.11−13,16,17 The main advantage of AIMD compared to the hybrid method is the fact that chemical reactions can be studied. In this paper, we aim to provide a formulation of temperature-dependent thermal transport properties of some important phases constituting the side ledge and molten cryolite (XNaF = 0.75 and XAlF3 = 0.25), These phases are part of the electrolysis. These phases are NaF, β−Na3AlF6, Na5Al3F14 and molten Na3AlF6. The main objective is to equip the aluminum production industry with a knowledge of the thermal conductivity of these pure phases together with complementary models10,18−20 associated with grain size and porosity. The industry will more easily simulate the thermal conductivity of a given side ledge microstructure assemblage and hence better take into account the side ledge thickness variation during the process. To achieve this, a series of EMD simulations, based on hybrid DFT/Classical methods, have been performed in both isobaric−isothermal (NPT) and canonical (NVT) statistical ensembles. A set of polarizable potentials was developed, on the basis of DFT calculations, to describe each ion−ion pair potential within the system. The results of EMD simulations are compared to what is predicted by the mean of selfconsistent thermodynamic models developed in previous work for solids8 and liquids.7 These models were developed to link thermal conductivity to thermodynamic properties and key physical properties of the materials, and it is thus a good basis for discussing the reliability of EMD predictions. Finally, a simple expression for the temperature dependence of thermal conductivity is proposed for each phase studied.

For an electrolyte containing LiF, the phases present or potentially present in the side ledge may be listed as2 • NaF (cubic) • LiF (cubic) • CaF2 (cubic) • α−Na3AlF6(monoclinic) • β−Na3AlF6 (cubic) • Na5Al3F14 (tetragonal) • α−NaCaAlF6(monoclinic) • β−NaCaAlF6 (trigonal) • Na2Ca3Al2F14 (cubic) • Na4Ca4Al7F33 (cubic) • α−Na2LiAlF6(monoclinic) • β−Na2LiAlF6 (cubic) • LiCaAlF6 (hexagonal) • α−Al2O3 (trigonal) • γ−Al2O3(1) (cubic) • γ−Al2O3(2) (tetragonal) • NaAl11O17(hexagonal) • Al4C3 (trigonal) • MgF2 (tetragonal) • Coke (anodic) For these phases, data are missing particularly for thermal diffusivity. For a significant number of phases, especially the most complex ones like Na4Ca4Al7F33, Na2Ca3Al2F14 or Na2Al11O17, only information on structure and density at room temperature is available. In fact, the only experimental data related to the thermal transport properties of the side ledge hitherto are from 1970,3 with the use of a hot wire apparatus.4 The following empirical expression of the thermal conductivity, as a function of temperature and alumina content, was then proposed by Haupin:3 K (T , WAl2O3) = 0.594 + 6.35 × 10−4WAl2O3 + 6.35 × 10−4T[W /(m . K )]

(1)

where WAl2O3 denotes the alumina weight fraction. This formulation is simplistic and takes into account the alumina content (without specifying the phase) and temperature. Hence it cannot really be used as a generic function for the thermal conductivity of the side ledge. Indeed, the thermal conductivity depends strongly on the microstucture and its parameters: its type, grain size, porosity level, grain size distribution and grain orientation. Moreover, in the Haupin3 representation, the thermal conductivity increases with temperature, which is not physically acceptable for electrically insulating materials. The thermal conductivity of insulating materials is influenced by phonon−phonon scattering and should decrease as 1/T.5 This may be due to the experimental technique discussed by DiGuilio6,7 or it could mean that the Haupin formulation is valid only in a narrow range of temperature, near the liquidus temperature. The temperature-dependent thermal conductivity of pure stoichiometric compounds, both in solid and liquid states, can be predicted with appreciable precision using either a theoretical approach7−10 or Equilibrium molecular dynamics(EMD) simulations.7,11−14 For solid phases, a self-consistent thermodynamic method8 has been developed to predict the lattice thermal conductivity above room temperature, provided that reliable experimental data are available for (i) heat capacity, (ii) thermal expansion and (iii) adiabatic bulk modulus in a wide temperature range. This procedure has been shown to be B


Article


■

SIMULATIONS AND THEORETICAL METHODS The molecular dynamics simulations required for predicting the thermodynamic and the thermal transport properties involve: 1. Developing reliable interaction potentials to describe all the ion−ion interactions within the system; 2. Performing a robust parametrization of the interaction potentials, in this case based on DFT calculations ; 3. Developing suitable setup and procedures for simulations; 4. Extracting the information from the generated phase trajectories and their conversion into desired properties by means of suitable statistical relations Interaction potentials. The so-called Dipole Polarizable Ion Model (DIPPIM) potential was used to describe ionic interaction between different species i and j.15 The total potential (Vtot) includes a charge−charge (c−c), dispersion (disp.) overlap repulsion (rep.), self-energy (s-e) and polarization (pol.) contributions. The self-energy term includes a Ewald term (E) plus a self-dipole (s-d) contribution. The polarization contribution itself consists of dipole−dipole (d−d) and charge−dipole (c-d) contributions. Thus, the total interaction potential can be written as Vtot =

with respect to the induced dipoles. The self-dipole contribution is independent of rij and is expressed as Vis − d . =

⎧ V c − d .(r ) = [q μ α g ij(r ) − q μ α g ji(r )] (1)(r ) ij α ij ⎪ ij i j 4 ij j i 4 ij ⎨ ⎪ V ijd − d .(rij) = −μiα μjβ  (2) αβ (rij) ⎩

rα ⎧  (1) α (r) = − 3 ⎪ r ⎪ ⎨ 3rαrβ − r 2δαβ ⎪ (2) ⎪  αβ (r) = ⎩ r5

with δαβ being the Kronecker delta. describes the short-range induction on the dipoles. Similar to the dispersion term, these short-range induction effects can be analytically expressed through the Tang−Toennies damping functions:21

∑ (Vijd − d . + Vijc− d .)

i,j i  

(2)

ij

k=0

=1−e

n

∑ k=0

(3)

(bijrij)k k!

k!

(8)

where reflects the amplitude of this damping at ion j due to the presence of i, and bijD is a range parameter describing the overlap of the charge density upon the induced dipole. DFT parametrization of the interaction potential. According to the DIPPIM formalism as described above, {Bij, αij, Cij6, Cij8, bij, cijn, bijD, CijD, αi} is the set of parameters describing the polarizable interaction potential. The potential parameters describing F−-F−, Na+-F− and Na+-Na+ interactions were already formulated by Salanne et al.23 In the present study, we still need the interaction potential parameters for: Al3+-Al3+ Al3+-Na+ and Al3+-F−. The procedure for determining the potential parameters, based on the force matching technique, was initiated by Madden15 and is extensively described elsewhere.22,24 Here, it is not necessary to mention every detail of the optimization procedure; only the key points are presented. First of all, it is important to highlight that this optimization procedure is based on electronic structure calculations using the Density Functional Theory (DFT). Clearly, the potentials generated are of a predictive nature, since no experimental information was used in the parametrization procedure. The optimization procedure consists first in generating several ionic configurations via EMD simulations, for a small system (100 ions) with periodic boundary conditions, considering arbitrary interaction potential parameters for Al− Al, Na−Al and Al−F. Then, for each ion in each generated structure, the force and the dipole are calculated by DFT in a plane-wave pseudopotential framework. We employ the Perdew−Burke−Ernzerhof (PBE) exchange-correlation functional25 within the generalized gradient approximation (GGA).

where qi is the formal charge on ion i and Bije is the BornMayer-Huggins short-range potential responsible for the overlap repulsion of the electronic clouds at short distances. Cij6 and Cij8 are the dipole−dipole dispersion coefficients and f ijn is the dispersion damping function for short-range correction of interactions between charge, dipole and dispersion interactions. f ijn are described by the Tang−Toennies function:21 −bijrij

(bDijrij)k

cijD

−aijrij

f nij (rij)

4

g4ij(rij) = 1 − cDije−bDrij ∑

In contrast to the polarization term, the first three terms of the potential are purely pairwise additive and depend only on the interionic distance (rij) between the ions i and j: qiqj ⎧ V ijc − c(rij) = ⎪ rij ⎪ ⎪ ⎪ ⎡ C ij ⎤ ij ⎨ disp. ij 6 ij C8 ⎪ V ij (rij) = −⎢f 6 6 (rij) + f 8 8 (rij)⎥ ⎢⎣ rij ⎥⎦ rij ⎪ ⎪ ⎪ Vijrep.(rij) = Bij e−aijrij ⎩

(7)

gijn

s−e

pol.

(6)

In this expression, the Einstein summation convention (i.e summation is carried out over repeated indices) is applied with respect to α and β which denote the Cartesian coordinates (x, (2) y, z).  (1) α and  αβ are respectively the charge−dipole and the dipole−dipole interaction tensors:

i 

+ ∑ Vis − d +

(5)

In the self-dipole contribution, μi and αi are respectively the dipole moment vector and the polarizability of the ion i. Finally, the charge-dipole term is defined as

∑ (Vijc− c + Vijdisp. + Vijrep.) + ∑ (ViE .) i,j>i

1 |μi |2 2αi

(4)

ij

where the coefficient b represents the distance at which the correction begins to be taken into account. For the present type of potential, n is an integer which can assume the values 6 or 8. It should also be noted that the two limit values of f satisfy: f n(0) = 1 and f n(∞) = 0. The expression of the Ewald term is that formulated by Madden et al.22 The dipole moment is obtained by by minimizing the polarization energy, Vpol = Vs−d + Vd−d + Vc−d, C


Article


Table 1. Parameters of the DIPPIM potential describing the interactions within the Al−Na−F system. All the values are in atomic units. The parameters for the F−−F−, F−−Na+ and Na+−Na+ pairs are from Salanne et al.23 ion pair −

−23

F −F F−−Na+23 F−−Al3+ Na+−Al3+ Na+−Na+23 Al3+−Al3+

Bij

aij

Cij6

Cij8

bij

bijD = bjiD

cijD

cjiD

282.3 34.0 52.8 1.0 1.0 1.0

2.444 1.855 1.974 5.0 5.0 5.0

15.0 13.3 0.0 0.0 11.7 0.0

150.0 88.2 0.0 0.0 51.8 0.0

1.90 1.90 1.90 1.9 1.9 1.9

1.00 1.84 1.59 1.9 1.0 1.9

0.0 2.54 1.84 0.0 0.0 0.0

0.0 −0.19 0.0 0.0 0.0 0.0

properties are defined by the couple (K, CP). The other thermal transport properties can be deduced therefrom as explained above. The crystallographic information necessary to build the supercell of initial solid configuration was taken from ref 28. The initial configuration of the molten cryolite has been obtained by the simulation of the mechanical melting of an α− Na3AlF6 (monoclinic, P21/c) supercell followed by a thermal equilibration at 10 K above the liquidus temperature (1283 K). The total number of ions is 512 (256 F−, 256 Na+) for NaF, 640 (324 F−, 162 Na+, 54 Al3+) for β−Na3AlF6, 352 (224 F−, 80 Na+, 48 Al3+) for Na5Al3F14 and 640 (324 F−, 162 Na+, 54 Al3+) for molten Na3AlF6. We have considered a wide range of temperature for the solid phases: 373 K, 573 K, 773 K, 873 K 973 K, 1073 K, 1173 and 1273 K plus four temperatures evenly distributed between the liquidus temperature and about 200 K above it. The simulations were carried out in both the isobaric− isothermal (NPT) and canonical (NVT) statistical ensembles. The first series of simulations in the NPT ensemble was performed to determine the equilibrium volume and to generate thermally equilibrated configurations for the subsequent simulations in the NVT ensemble. The thermodynamic properties Vm, CP, isothermal compressibility (ξT) and volumetric thermal expansion (αV) were also extracted from the simulated phase space trajectories obtained from the simulations in the NPT ensemble. Then, starting with an initial configuration generated in the simulation in the NPT ensemble and thermally equilibrated, a new series of EMD simulations was performed in the NVT ensemble in order to determine the thermal conductivities of mixtures. The volume used in the simulation in the NVT ensemble was fixed equal to the average volume determined by the NPT simulations. The temperature was controlled with a suboptimal Nosé−Hoover thermostat29 for both NPT and NVT simulations while for the NPT simulation, the pressure was controlled by an extension of the Martyna barostat.30 The thermostat and barostat relaxation times were both 0.5 ps. The equations of motion were integrated using the Verlet algorithm with a time step of 1 fs for all simulations. The total simulation time was 1.5 ns for the NPT simulations and 30 ns for NVT simulations. Lastly, all simulations were performed using periodic boundary conditions and the minimum image convention. Derivation of the thermal conductivity and thermophysical properties from the MD simulations. Thermal conductivity is a second-order tensor,  , defined by the Fourier law,

All the DFT calculations were performed with the CPMD code package.26 The parameters of the potential are then determined by minimizing the difference between the predicted and the calculated values on these configurations. The minimization procedure is performed using the Minuit nonlinear minimization routine.27 This minimization procedure allows setting the C6 and C8 dispersion parameters. They are determined ensuring that NPT simulations give the experimental and/or ab initio density at one temperature. In addition, as pointed out by Salanne et al.,23 the cation−cation short-range parameters (CijD) were not fitted because they interact mainly through their formal charges. For more details on the potential parameter optimization procedure, see references.22−24 The values for the different polarizable potential parameters are given in Table 1. The polarizabilities of F− and Na+ are respectively: αF− = 1.17 Å3 and αNa+ = 0.15 Å3 and they were determined by Salanne et al.23 The Al3+ ions are considered to be not polarizable, i.e. αAl3+ = 0. The quality of the force matching technique is represented in Figure 1. For 100 ions (60 F−−10 Al3+−30 Na3+) in the molten

Figure 1. Illustration of the quality of the force matching technique fit on a molten cryolite system thermally equilibrated at T = 1293 K, containing 100 ions (60 F−−10 Al3+−30 Na3+). For each ion, the components of dipole moments (μ) and forces (F) in the x, y and z directions obtained by the present polarizable potential are represented as a function of those predicted by DFT calculations. The solid line indicates the equality of EMD and DFT values. Pearson’s χ2 test describing the overall goodness of the fit is χ2dip. = 0.039 and χ2f. = 0.017 for the dipoles moments and forces, respectively.

cryolite structure, equilibrated at 1293 K, x-,y- and zcomponents of the dipole moments and forces obtained by the present polarizable potential are compared with those predicted by DFT. Uniformly good fits are found for both forces and dipole moments with χ2dip. = 0.039 and χ2f. = 0.017 respectively for dipole moments and forces. This indicates the good quality of the fit of the potential. Details of molecular dynamics simulations. In this work, the series of EMD simulations was primarily performed to determine the thermal transport properties of NaF (cubic, Fm3¯m), β−Na3AlF6 (cubic, Fm3¯m), Na5Al3F14 (tetragonal, P4/mnc) and molten Na3AlF6. Here, the thermal transport

jQ = −∇T

(9)

which relates the heat flux vector in a material, jQ, to the temperature gradient, ∇T. For cubic (and isotropic) and tetragonal structures the thermal conductivity tensor is written as D


Article

The Journal of Physical Chemistry C ⎡1 0 0⎤ ⎢ ⎥ cub(iso.) = K 0⎢ 0 1 0 ⎥  tetra. = ⎣0 0 1⎦

⎡ K⊥ 0 0 ⎤ ⎥ ⎢ ⎢ 0 K⊥ 0 ⎥ ⎥ ⎢ ⎣0 0 K⎦

thermoelectrical effect upon the thermal conductivity. In the case of purely ionic solids, the thermal conductivity is reduced to Lee/T2 since the charge displacement within the solids is negligible. When the charge displacement cannot be neglected, the thermal conductivity is a sum of the contribution of an electrically neutral system and the contribution due to thermoelectrical effects. Theoretical Model for Lattice Thermal Conductivity as a Function of Temperature. In this work, we need first to address the question of the validity and accuracy of the predicted thermal conductivity of compounds formed in the NaF−AlF3 system. For this system, no experimental data are available, except for NaF, for which a few experimental points (at temperatures ≤400 K) can be found. Unfortunately, this is far from the working temperature of the electrolysis cells and cannot be used for comparison. It is therefore not straightforward to discuss the quality of EMD simulation results and to quantify their accuracy. Here, we propose to rely on the theoretical models to analyze and discuss our calculations. One way to overcome this severe lack of experimental data was proposed by Gheribi and Chartrand8 through a selfconsistent thermodynamic method. In summary, this method considers the analytical Slack−Morelli expression5 which links the temperature-dependent lattice thermal conductivity of cubic and isotropic systems to the acoustic Debye temperature (θa) and the Grüeinsen parameter (γ)

(10)

All the components of the thermal conductivity tensors are temperature dependent. For cubic and isotropic systems, including liquids, the thermal conductivity is defined by a single constant K0. For tetragonal systems,  is diagonal but it is non diagonal for monoclinic systems. In the literature, the thermal conductivity at constant temperature is often reported as a single value, even for non cubic and non isotropic systems corresponding to an effective value, ⟨K⟩, for the associated polycrystalline material. For hexagonal, trigonal, tetragonal and orthorhombic symmetry the effective thermal conductivity is given by the average of the trace: ⟨K ⟩ = Tr()/3. The thermal conductivity may be calculated from EMD simulations using the Green−Kubo (G−K) formalism.31 It is based on the fluctuation−dissipation theorem and allows deduction of the thermal conductivity from the energy and charge current fluctuations generated by EMD simulations. Following the G−K method, the diagonal components of the lattice thermal conductivity tensor for a system with two species, NaF in the present work, is written as16 Ks =

ss 2 ⎤ ⎡ 1 ⎢ ss (Lez1) ⎥ − L ; s ∈ {x , y , z} ee L zss1z1 ⎥⎦ T 2 ⎢⎣ −

+

(11) 3+

K (T ) =

For a system with 3 species (F , Na , Al ), the diagonal components of  is expressed as16 ⎡ (Lezss1)2 Lzss2z 2 + (Lezss2)2 Lzss1z1 − 2Lezss1Lezss2Lzss1z 2 ⎤ 1 ⎥; s Ks = 2 ⎢Leess − ⎥⎦ T ⎢⎣ Lzss1z1Lzss1z 2 − (Lzss1z 2)2 ∈ {x , y , z}

×

1 3kBV

∫0

τ

dt

(13)

where n and m represent either e or zi. kB and V are respectively the Boltzmann constant and the equilibrium volume of the system. The two vectors Jse and Jszi are the Cartesian components of the low-wave vector limits of the spatial Fourier transforms of the energy and charge currents of the ith ion. In practice, when calculating the thermal conductivity, the total of the NVT runs is divided into 20 blocks of 0.5 ns each, and the different transport coefficients are calculated independently in each block. The correlation time, τ, is chosen to be large enough to ensure the convergence of the thermal conductivity. τ is then set to 20 ps for all compositions studied. Lastly, the thermal conductivity of a perfectly thermally equilibrated system is defined as an average value over all the different blocks and is expressed as Ks = lim ⟨Ks(τ )⟩all blocks

(14)

τ→∞

vat.(T )1/3 [W/(m K)] T

(15)

All the quantities in this expression are expressed in SI units. ℏ, m¯, and vat. are, respectively, the reduced Plank constant, the average atomic weight, and the average volume per atom. n is the number of atoms per primitive cell; n becomes greater as the structure becomes complex. From a theoretical point of view, the acoustic Debye temperature is obtained by integrating the phonon density of state over the acoustic modes. More practically, θa is deduced from the traditional Debye temperature according to the following relationship:32 θa = θDn−1/3. Therefore, K ∝ n−2/3, and n−2/3 describes the structural dependence of the lattice thermal conductivity. Both γ and θD are also involved in the temperature-dependent theoretical expression of the volumetric thermal expansion (αV), the heat capacity at constant pressure (CP), and the adiabatic bulk modulus BS.8 The methodology proposed by Gheribi and Chartrand8 determines γ and θD by simultaneous optimization of αV, CP, and BS in order to reproduce the available set of reliable experimental data. The thermal conductivity is then predicted using eq 15. This method demonstrated a good predictive capability for different types of materials,8,9,12 including ionic materials and semiconductors. Note that in the case of an extreme lack of reliable experimental data, which is the case for many other fluoride compounds present in the NaF−AlF3−CaF2 system, this method cannot be applied to determine the model parameters. As a compromise strategy, a DFT-based method was developed by Gheribi et al.9,33,34 in order to determine self-consistency of both θD and γ. This method allows the prediction of not only the lattice thermal conductivity but also all the thermodynamic properties as a function of temperature. The method, called “Thermodynami-

(12)

The indices e and z refer respectively to the energy and charge currents. The coefficients L are the transport coefficients. They are defined through the Grenn-Kubo formula:31 sp Lnm =

⎛ kB ⎞ 3 (6.51982 × 10−3)m̅ ⎜ θa⎟ × n1/3 × γ 2 − 0.514γ + 0.228 ⎝ ℏ ⎠

The first term in eq 11 and eq 12, i.e., Lee/T , corresponds to the thermal conductivity of an electrically neutral system, whereas the second term (the rational function) describes the 2

E


Article


be predicted with the TSC method. Here, as an effort to keep this manuscript as short as possible, the formulation of the thermodynamic properties will not be repeated; the reader is directed to refs 8, 9, 33, and 34. Basically, this method considers several contributions to the density of the lattice vibration energy: harmonic, lattice expansion, defects, electronic, magnetic, and anharmonic, and these contributions are described by a few key physical parameters. In the present case, only the harmonic, lattice expansion, and defect contributions are taken into account, and the key parameters are Vm (0 K), θD, γ, and EVac, the energy of formation of the Schottky defect.

cally Self-Consistent (TSC) method”, extends the quasiharmonic approximation (QHA) model35 by incorporating a minimization procedure to ensure the conservation of Maxwell relations. The formalism of this method has been published, and it will not be repeated here. For full details, the reader is directed to references 9, 33, and 34. The Debye temperature is calculated from DFT calculations of the elastic constants, according to the following relation35−37 ⎛ ℏ ⎞⎛ 9n ⎞1/3 θD = ⎜ ⎟⎜ 0 ⎟ ρ−1/2 a0−1/3 ⎝ kB ⎠⎝ 4πν ⎠

(16)

■

where n0, ν, and ρ are, respectively, the number of atoms in a unit cell, the unit cell volume, and the density. a0 is a function of the elastic constants and can be found in ref 38. The procedure for calculating elastic constants is well established in the literature. Basically, it consists in calculating the energy difference between the equilibrium lattice R and a distorted lattice (R′), by applying very small isochoric strains (ϵ), in order to make sure the process stays within the elastic domain. The elastic constants (Cij) are obtained by fitting the energy curve E(ϵ) by quadratic functions, the ϵ2 coefficient being a linear combination of elastic constants. The distortion matrix (D) is a linear application which transforms R into R′, according to R′ = RD. The number of distortion matrices to be considered is equal to the number of elastic constants, i.e., 3 (C11, C12, and C44) and 6 (C11, C12 and C13, C33, C44, and C66), respectively, for cubic and tetragonal symmetries. The distortion matrices can be found in ref 39. For each distortion matrix, 10 different strains (5 positives and 5 negatives) were applied to the equilibrium lattice in order to obtain an accurate E(ϵ) curve. This procedure is then repeated under external pressures (1, 2, 3, and 5 GPa) in order to obtain the Grüneisen parameter. Assuming that the volume dependence of all modes of the phonon frequencies is identical to the volume dependence of the Debye frequency, the Grüneisen parameter is defined as

γ=−

∂ln[θD(V )] ∂ln V

RESULTS AND DISCUSSION For NaF, β-Na3AlF6, Na5Al3F14, and molten cryolite, the thermal conductivity values obtained by the present EMD simulations are presented in Table 2. The physical properties Table 2. Thermal Conductivities, Expressed in W/(m K), Obtained by EMD Simulations, for NaF, β-Na3AlF6, Na5Al3F14, and Molten Cryolite at Temperatures between Room Temperature up to the Melting Point β-Na3AlF6

NaF T

K

373 573 773 873 973 1073 1173 1273

18.93 (0.05) 10.07 (0.07) 6.31 (0.05)  4.64 (0.06) 4.77 (0.05) 4.25 (0.07) 3.86 (0.08)

K

Molten cryolite T 1293 1373 1423 1473

(17)

DFT calculations were performed only for β-Na3AlF6 and Na5Al3F14, for which very few experimental data of heat capacity and molar volume are available. For NaF, there are sufficient experimental thermodynamic data available to determine the theoretical model parameter. Experimental data are not completely unused; in fact, the DFT data are used as “target parameters values” for the theoretical model, and they are then adjusted to reproduce simultaneously the experimental heat capacity and molar volume. DFT calculations are based on the plane-wave basis sets and are performed using the Vienna ab initio simulation package (VASP).40−43 The projected-augmented wave (PAW) approach is employed to represent the core electrons.44,45 The generalized gradient approximation (GGA) parametrized by Perdew, Burke, and Ernzerhof (PBE)46,47 was used as the exchange-correlation functional. A plane-wave kinetic cutoff energy of 520 eV and Monkhorst Pack grid of 8 × 8 × 8 for βNa3AlF6 and 8 × 8 × 4 for Na5Al3F14 with a first-order Methfessel−Paxton smearing parameter σ of 0.02 eV were used to ensure that the force and energy convergence criteria are better than 0.02 eV/Å and 0.01 meV, respectively. The heat capacity, thermal expansion, α = V−1(∂V/∂T)P, and the temperature dependence of the adiabatic bulk modulus can

K(Veq) 0.72 0.74 0.76 0.78

Na5Al3F14

(0.12) (0.13) (0.11) (0.14)

1.85 1.45 1.12 0.97 1.01 0.93  0.64

K⊥

 0.91 (0.13) 0.80 (0.15) 0.75 (0.10) 0.62 (0.15)   (0.1)  Molten cryolite

(0.25) (0.20) (0.15) (0.20) (0.20) (0.15)

K∥  0.67 0.63 0.52 0.48   

(0.09) (0.11) (0.07) (0.11)

(at constant volume) K(V1293K) 0.72 (0.12) 0.71 (0.10) 0.71 (0.11) 0.705 (0.12)

   

   

a

The numerical error in the thermal conductivity calculations is reported in parentheses next to the thermal conductivity value and corresponds to the standard deviation of the mean (SEM). The mean of K was calculated in the region where Ks(τ) has reached a plateau value. For molten cyrolite, in addition to the thermal conductivity at the equilibrium volume K(Veq), the thermal conductivity at constant volume at T = 1293 K (V1293K, is also reported (see discussion on molten cryolite)).

obtained by DFT, allowing the prediction of thermal conductivity for solid compounds, using eq 15, are given in Table 3 (equilibrium molar volume, elastic constant, Debye temperature, and Grüneisen constant). NaF (Cubic, Fm3̅m). In Figure 2, the results of the present EMD simulations are compared with the available experimental data and with the theoretical predictions obtained via eq 15. The parametrization of eq 15 was performed in our prior work8 on the basis of self-consistent thermodynamic optimization of CP(T), αV(T), and BS(T). Reliable experimental data for these are available in the entire range of temperatures, i.e., from 0 K up to the melting temperature (Tm = 1269 K). It was found that the set of parameters {θD = 495 K; γ = 1.56} can, along F


Article


molecular dynamics. When the temperature is below the Debye temperature, the quantum effect upon the vibrational properties of the lattice sets in and are all the more important when the temperature is low.51 Because of their classical nature, EMD simulations give the limit of 3kB for the heat capacity at constant volume and thus overestimate the density of the lattice vibration energy. This results also in the overestimation of the lattice thermal conductivity below the Debye temperature. Haigis et al.52 have perfomed EMD simulations to calculate the thermal conductivity of MgO, which is also an ionic compound, at 300 K. Their calculations overestimated the thermal conductivity by 70% in comparison with the experimental value. However, the Debye temperature of MgO is approximately 725 K, and it is thus natural that qualitatively the overestimation of K increases with the ratio θD/T until the classical limit of 1. There is no clear rule hitherto reported to quantify the overestimation of the lattice thermal conductivity by EMD simulations below the Debye temperature, due to the noninclusion of the quantum effect upon the lattice vibration density energy. For NaF, some important conclusions are (1) EMD simulations confirm the validity of the theoretical prediction at high temperature, (2) EMD simulations at high temperature, i.e., above the Debye temperature, are reliable, and (3) similar to the heat capacity predictions, EMD simulations fail to predict accurately the lattice thermal conductivity at low temperatures, i.e., below the Debye temperature. β-Na3AlF6 (Cubic, Fm3̅m) and Na5Al3F14 (Tetragonal, P4/mnc). NaF is the only compound within the NaF−AlF3 system for which experimental data on thermal transport properties are reported. For β-Na3AlF6 and other systems, our approach is thus purely predictive. It is known that cryolite (Na3AlF6) is found as two polymorphs: a low-temperature phase, the so-called α-cryolite, which crystallizes in a monoclinic structure (P21/n), and the β-cryolite (Elpasolite Fm3̅m) which is stable at high temperatures, more specifically above 836 K.55−58 While α-cryolite is a stoichiometric compound, β-cryolite is more complex. The NaF−AlF3 phase diagram is shown in Figure 3, according to the critical assessment by Chartrand and Pelton, using the CALPHAD method.59 As can be seen in the close-up diagram, β-cryolite is a nonstoichiometric compound. This nonstoichiometry was identified and measured first by Dewing in 197853 and later again confirmed by Dewing in 1997.54 Although the range of compositions where the nonstoichiometry is observed is very small (less than 0.75 mol %) strictly speaking the thermal conductivity should be described as one of a solid solution. It is known that the lattice thermal conductivity of solid solutions could deviate strongly from a simple linear behavior. From a theoretical point of view, the degradation of the lattice thermal conductivity of solid solutions, defined as the relative deviation from linear behavior, is explained by the mass and the strain field fluctuations upon the phonon−phonon scattering cross section.60 The thermal conductivity degradation is more pronounced for semiconductor solid solutions, especially those which are constituted of a single sublattice such as Si− Ge, as opposed to oxide or halide solutions.12,17,60,61 For solid solutions such as Si−Ge, the thermal conductivity curve as a function of composition exhibits a U-shaped curve (concave upward). That is, even for dilute solutions, the thermal conductivity degradation can be large. This typical shape is caused by the conjugation of two factors: (1) the isotope effect, which increases considerably the mass fluctuation contribution

Table 3. Theoretical Model Parameters (Equation 15) Determined by DFT (Target Value) and Optimized to Reproduce Simultaneously the Molar Volume and Heat Capacity, Both As Functions of Temperature (Optimized Value) β−Na3AlF6 Vm(0K) θD γ n

Na5Al3F14

Target

Optimized

Target

Optimized

74.38 508 1.48 40

72.25 515 1.51 

160.62 511 1.53 44

152.37 495 1.58 

a

The number of ions per primitive cell is taken from the Materials Project database.28 The molar volume is expressed in cm3/mol and the Debye temperature in K, the Grüneisen parameter being dimensionless.

Figure 2. EMD simulation results (filled stars) of the thermal conductivity of NaF (cubic, Fm3m ̅ , n = 2) as a function of temperature and up to the melting temperature (Tm = 1269 K) in comparison with experimental data (open triangles;48 open diamonds;49 and open circles50) and theoretical predictions (dashed line) obtained via eq 15, using the set of parameters {θD = 495 K; γ = 1.56} determined by Gheribi and Chartrand8 by a self-consistent optimization of thermodynamic properties. For NaF the SEM was found to be not significant (see Table 2), and it is not reported in this figure.

with an energy of vacancy formation EVac. = 122.8 kJ/mol,8 describe simultaneously CP(T), αV(T), and BS(T) from 0 K up to Tm. As explained previously, experimental data on thermal conductivity and/or thermal diffusivity are available only below 400 K. In this narrow range of temperatures, the theoretical prediction matches remarkably well with the experimental data. Above approximately the Debye temperature, the EMD simulation data are in very good agreement with the theoretical prediction, and thus by extrapolation, it should also be close to the experimental data. At 373 K, the discrepancy between EMD simulations and the experiments or theory are found to be rather large: more than 45%. The excellent predictive capability of EMD simulations at high temperature is not surprising, as the same agreement was found for several ionic compounds.16,17 At low temperature, the overestimation is clearly due to an overestimation of the density of the lattice vibration energy. For NaF, at low temperature, the grain size effect upon the lattice thermal conductivity is negligible unless the grain size is below 50 nm.10 Such an overestimation is intrinsically linked to classical equilibrium G


Article


Figure 3. NaF−AlF3 phase diagram (left) according to the CALPHAD optimization of Chartrand and Pelton.2 The calculated phase diagram (solid line) is compared with various experimental data. The close-up phase diagram (right) shows the thermodynamic modeling of the nonstoichiometry of the β-cryolite in comparison with experimental data of Dewing.53,54

Figure 4. Experimental vs calculated heat capacity (CP) and molar volume (Vm) as a function of temperature for cryolite (Na3AlF6). The calculation of CP(T) and Vm(T) was performed according to the self-consistent thermodynamic method proposed by Gheribi and Chartrand8 using a DFT value of Vm, θD, and γ as an initial guess for the optimization procedure. The thermal expansion predicted using the optimized value of Vm and θD, and γ is shown in the inset plot. For CP, two cases were considered: one with a defect contribution, explaining the CP behavior near the transition temperature, and one with only a vibrational contribution. The experimental data are referenced as follows: Anovitz et al.63 (open squares) and assessed values (Ass.) of Barin64 (open triangles) for CP and Zhou and Kennedy65 (open squares), Dolejs and Baker58 (open circles), and Yang et al.66 (open triangles) for molar volume.

mnc) structures, at 154 K.62 Above this temperature, tetragonal chiloite is stable and purely stoichiometric up to 999 K. It then melts incongruently forming a liquid phase in equilibrium with β-Na3AlF6, as shown in the NaF−AlF3 phase diagram (Figure 3). Therefore, for both chiolite and β-Na3AlF6, strictly speaking, it is impossible to determine the Debye temperature experimentally as these phases are not stable at low temperatures. In this case, DFT takes its full meaning. For both β-cryolite and chiolite, very few reliable experimental data are available for heat capacity and molar volume above the α → β transition temperature (Ttrans.). This lack of experimental data at high temperature is due to experimental difficulties, in addition to safety considerations. The optimization procedure of the theoretical model (eq 15) parameters is as follows: Vm(0 K), θD, and γ are first determined by DFT combined with the TSC method,9,33,34 then adjusted to reproduce simultaneously the available experimental data. The calculated and optimized model parameters are given in

upon the phonon−phonon scattering cross section, and (2) the fact that there is no constraint in the sites where the atom can mix. For oxide and halide solid solutions, the thermal conductivity curve as a function of composition exhibits a more quadratic shape (also concave upward), and in this case, the thermal conductivity degradation of dilute solid solutions is usually small. Moreover, the thermal conductivity degradation tends to be attenuated as the temperature increases.12,17,61 In the present case, given the very small composition range (≤0.75 mol %) of the nonstoichiometry of the β-cryolite and the fact that it is stable at high temperature, it can be reasonably assumed that the thermal conductivity of β-cryolite, as a nonstoichiometric compound, is almost equal to that of stoichiometric β-Na3AlF6. The difference is assumed to be at least within the numerical error. Chiolite (Na5Al3F14) also has two polymorphs. Indeed it was reported by Rocquet et al.62 that chiolite exhibits a structural phase transition, from monoclinic (P21/n) to tetragonal (P4/ H


Article

The Journal of Physical Chemistry C Table 3 for both β-Na3AlF6 and Na5Al3F14 (tetragonal). Figures 4 and 6 show the comparison between the calculated and the experimental heat capacity and molar volume versus temperature. Except for the overestimation of the molar volume, which is inherent to the GGA functionals,67 the DFT method, combined with TSC, shows good predictive capability as θD and γ were slightly adjusted to reproduce the experimental data. For β-Na3AlF6, the heat capacity above the α → β transition temperature is well reproduced. The extrapolation of the heat capacity of β-Na3AlF6 below Ttrans. is also in good agreement with the heat capacity of α-Na3AlF6 at up to about 0.77 Ttrans., where the thermal vacancy defect contribution starts to increase the heat capacity strongly. The energy of formation of thermal vacancy defects was estimated, from heat capacity experimental data, as 168.3 kJ/mol. There is a single experimental datum reported by Dolejs et al.,58 at 1282 K and thus very close to the melting temperature. However, this datum cannot be related to the β-Na3AlF6 as a stoichiometric compound because at this temperature the cryolite is nonstoichiometric (see Figure 3) and contains a significant number of structural defects, increasing the molar volume. Since the heat capacities of the two cryolite allotropes are very similar, it should be the same for the thermal expansion. In the optimization procedure, we have considered only the recent data of Zhou and Kennedy65 for α-Na3AlF6, and we assumed that the molar volume of βNa3AlF6 is

Figure 5. EMD simulation results (filled stars) for the thermal conductivity of β-Na3AlF6 (cubic, Fm3̅m) stoichiometric compound as a function of temperature and up to the melting temperature (Tm = 1282 K) in comparison with the experimental thermal conductivity of the side ledge (open diamond3) and theoretical predictions (dashed line) obtained via eq 15. These were obtained using the set of parameters {θD = 515 K; γ = 1.510} determined by a hybrid DFT method and the self-consistent thermodynamic optimization method8 in order to reproduce simultaneously CP(T) and Vm(T) as presented in Figure 4.

thermal conductivity), chiolite, NaF, fluorite (CaF2), and other compounds within the NaF−AlF3−CaF2 system, such as NaCaAlF6 or Na2Ca3Al2F14. In addition, as no information on the microstructure was reported by Haupin, it is, strictly speaking, not possible to compare our predictions with the experimental data reported by Haupin. The fact that they are in good agreement can however be seen as support of the reliability of our calculations. The same optimization procedure was applied to chiolite (Na5Al3F14). As for to β-cryolite, the molar volume is overestimated by the GGA functional, and both θD and γ are predicted with good accuracy. The comparison between the calculated heat capacity and molar volume and the available experimental data, as a function of temperature, is shown in Figure 6. Above the transition temperature, 154 K, both the heat capacity and the thermal expansion are reproduced accurately. Nevertheless, contrary to the case of cryolite, the model parameter of the tetragonal chiolite fails to predict accurately the heat capacity and thermal expansion of monoclinic chiolite, i.e., Ttrans.. This may be due to the specificity of the monoclinic → tetragonal phase transition of chiolite for which the contribution of the optical mode of phonons is predominant in driving the phase transition, as shown by Rocquet and Couz.71 The two components of the thermal conductivity, tensors K⊥ and K∥, and the effective thermal conductivity of polycrystalline chiloite, ⟨K⟩ = (2K⊥ + K∥)/3, are shown in Figure 7 in comparison with the predicted effective thermal conductivity using eq 15. The configuration is almost identical as for βNa3AlF6: at high temperatures, the theoretical model results are in very good agreement with those of EMD simulations. At lower temperatures, EMD simulations predict a lower value for ⟨K⟩. This may be explained by the same arguments proposed earlier for β-Na3AlF6. In the case of tetragonal chiolite, the underestimation of the thermal conductivity is observed around 573 K, i.e., above the Debye temperature (495 K). Contrary to the case of NaF, the inherent errors in the present EMD simulations of K’s, for both β-Na3AlF6 and

α − Na3AlF6 VMβ − Na3AlF6(Tα → β) = V M (Tα → β) + ΔVα → β

where ΔVα→β is the molar volume change associated with the phase transition. ΔVα→β = (0.89 ± 0.19) cm3/mol was determined by Dolejs et al.58 using the Clausius−Clapeyron relation considering the experimental pressure−temperature unary phase diagram for the α → β transition. At 1283 K, the predicted molar volume is lower than that reported by Dolejs58 by about 1 cm3. This difference can be explained by structural and probably thermal defect contributions, which are not taken into account in our optimization procedure. Then, using optimized values for both θD and γ, the thermal conductivity can be predicted. Theoretical and EMD simulations predictions are compared in Figure 5 with the experimental thermal conductivity of the side ledge at 1250 K.3 The overall agreement between the theoretical model and EMD simulations is very good. At low temperatures, i.e., below the Debye temperature, the agreement is fair, but contrary to NaF, EMD simulations predict a lower value. The contrary may be expected. It would be laborious at the moment to discuss this point and to determine if this is a limitation of classical EMDs or if it corresponds to real physical behavior. Equation 15 was developed based on the assumption that the phonon−phonon scattering results from only anharmonic umklapp processes (Uprocesses). However, even if the normal processes (Nprocesses) conserve phonon momentum and hence do not offer thermal resistance, the populated phonons can participate in U-processes thereby indirectly contributing to the overall thermal conductivity. This could explain the lower thermal conductivity value found at 373 K. However, it remains a hypothesis until further theoretical work and EMD simulations at low temperature are performed. At high temperature, the experimental data reported by Haupin3 are in good agreement with both the theoretical model and EMD simulations. Although the side ledge consists mainly of cryolite (≳80 wt. %) it also contains alumina allotropes (which have a higher I


Article


Figure 6. Experimental vs calculated heat capacity (CP) and molar volume (Vm) as functions of temperature for chiolite (Na5Al3F14). The calculation of CP(T) and Vm(T) was performed according to the self-consistent thermodynamic method proposed by Gheribi and Chartrand8 using DFT values of Vm, θD, and γ as initial guesses for the optimization procedure. The thermal expansion predicted using the optimized value of Vm and θD, and γ is shown in the inset plot together with the molar volume close to the transition temperature. The experimental data are referenced as follows: Stuve et al.68 (open squares), King69 (open circles), and Anovitz et al.63(open triangles), for CP and Gomez−Cuevas70 (open squares) for molar volume.

necessary to combine the present thermal conductivity data of single crystals with reliable models of effective thermal conductivity for an equivalent microstructure. Molten Cryolite. In Figure 8, EMD simulation results are compared with all the available experimental data for molten

Figure 7. EMD simulation results for the two components K⊥ (open circles) and K∥ (open squares) of the thermal conductivity tensor of Na5Al3F14 (tetragonal, P4/mnc) as a function of temperature and up to its incongruent melting point (Ti−m = 999 K). The effective thermal conductivity of polycrystalline chiolite, ⟨K⟩ = (2K⊥ + K∥)/3 (filled stars), is compared with the theoretical prediction (solid line) obtained via eq 15 using the set of parameters {θD = 495 K; γ = 1.58} determined by a hybrid method DFT and the self-consistent thermodynamic optimization method8 to simultaneously reproduce CP(T) and Vm(T), as presented in Figure 6.

Figure 8. EMD simulations results (filled stars) for the thermal conductivity of molten cryolite (Na3AlF6) as a function of temperature and up to 1473 K in comparison with experimental data reported of Polyakov and Mozhaev72 (open circles), Filatov et al.73(open squares) and Khokhlov et al.74 Solid and dashed lines represent the uncertainty. The thermal conductivity at constant volume and 1293 K is also represented (open triangles) and shows that the thermal conductivity of molten cryolite satisfies eq 18.

tetragonal Na5Al3F14, are found to be significant: 10%−20%. This is due to the need to determine six different autocorrelation functions and to the very strong cancellation between the two terms of eq 12. It is interesting to note that the physical parameters describing the thermal conductivity, vat, θD, and γ individually, are close for the corresponding species βNa3AlF6, tetragonal Na5Al3F14, and NaF. Thus, according to eq 15, the difference in thermal conductivity for these compounds is due mainly to n, that is, the difference in the crystallographic structure. The thermal conductivity depends on the microstructure and its parameters: (type, grain size, porosity level, grain size distribution, and grain orientation). The thermal conductivity predicted in this work for NaF, beta-Na3AlF3, and Na5Al3F14 is that of single crystals. In order to predict the thermal conductivity of real systems, i.e., microstructures, it is

cryolite. At first glance, considering the scatter of the data and, even if unreported, the experimental error which usually lies in the range 10%−25% (in some cases more than 50%− 100%),7,75,76 the present EMD simulations predict reliably. In particular, there is appreciable agreement with the recent experimental data of Khokhlov et al.74 For both EMD simulations and experiment, the thermal conductivity increases with temperature. This is not physical behavior for most molten salts.7,75,76 Indeed, after a rigorous analysis of the experimental techniques and the various sources of experimental error, DiGuilio6,76 concluded that for pure molten compounds only J


Article


Figure 9. Entropy of mixing of the NaF−AlF3 molten system at 1293, 1373, and 1473 K according to Chartrand and Pelton.2 The cryolite composition (XAlF3 = 0.25) is marked by a vertical line.

First, in order to confirm or refute the validity of eq 18 for molten cryolite, we have calculated the thermal conductivity following the same procedure as described above, but we fixed the volume to the equilibrium volume at 1293 K. The results are reported in Table 2 and represented in Figure 8 as well. It is

experimental data showing a negative value for (∂K/∂T) should be considered reliable. This criterion was later verified by Gheribi et al.7 for a wide variety of molten salts. From a theoretical point of view, Ohtori et al.11 suggested that the temperature dependence of K is determined mainly by the temperature dependence of the density. This behavior was also confirmed by Takase et al.77 for molten salt mixtures. This empirical observation is mathematically expressed by the following equation

⎛ ∂K ⎞ ⎜ ⎟ = 0 ⎝ ∂T ⎠V

shown that

error, it can reasonably be assumed that the thermal conductivity of molten cryolite obeys eq 18. The increase of the thermal conductivity with temperature is thus an anomaly. In the integration of eq 18 to obtain eq 19, only the vibrational frequency of the liquid quasilattice was taken into account. Previous EMD studies, in particular for halides, also predict a decrease of the thermal conductivity with temperature7,11−13,79 and have shown good agreement with eq 19. For most of the compounds studied by Filatov et al., it was shown that their experimental data are reliable only close to the melting temperature.7,75,76 The experimental data reported by Polyakov and Mozahev72 can be seen reasonably as a cloud of experimental points close to the melting temperature and centered at 0.65 W/(m K). The more recent experimental data set reported by Khokhlov et al.74 is very close to what is predicted by EMD simulations with about the same temperature dependence. It is thus probable that the thermal conductivity of molten cryolite increases slightly with temper-

(18)

The integration of this equation has already been presented in a previous work.7,13 The following expression was found ⎧ ⎡ 1⎤ K (x , T ) = K (x , Tliq)⎨1 − α(x , Tliq) ·⎢γ(x , Tliq) + ⎥ · ⎣ ⎩ 3⎦ ⎫ (T − Tliq)⎬ ⎭ (19)

For liquids, the Grüneisen parameter defines the volume dependence of the vibrational frequency of the liquid quasilattice (ν) and is expressed as γ = −

( ∂∂lnln νv )P .

78

( ∂∂KT )V is slightly negative. Considering the inherent

According to

this last equation, it is clear that the thermal conductivity should decrease with temperature, no matter the composition, since α and γ are both strictly positive parameters. It is not the case for molten cryolite, according to both: experiment and EMD simulations.

ature, even if

( ∂∂KT )V ≈ 0. From a practical point of view, EMD

simulations in the NPT ensemble give, for molten cryolite: α = 3.50 × 10−4 K−1 and γ = 0.45. Therefore, according to eq 19, if only the vibrational frequency of the liquid quasilattice is taken K


The Journal of Physical Chemistry C into account,

vib.

( ∂∂KT )P

■

≃ −1.975 × 10−4 W/(m K2); the

( )

*E-mail: [email protected].

≃ +3.335 × 10−4 W/(m K2). However, an addi-

Notes

The authors declare no competing financial interest.

tional “excess” physical contribution of vibrational frequency of the liquid quasilattice has to be taken into account to explain both the recent experimental work of Khokhlov et al.74 and the present EMD simulations. This excess contribution is large, of a magnitude of about twice the contribution induced by vibrational frequency:

excess

( ∂∂KT )P

AUTHOR INFORMATION

Corresponding Author

present EMD simulations show a value of ∂K EMD ∂T P

Article

■

ACKNOWLEDGMENTS This research was supported by funds from the Natural Sciences and Engineering Research Council of Canada (NSERC) and Rio Tinto Alcan. Computations were made on the supercomputer Briaré at the Université de Montréal, managed by Calcul-Québec and Compute Canada. The operations of this supercomputer are funded by the Canada Foundation for Innovation (CFI), NanoQuébec, RMGA, and the Fonds de recherche du Québec - Nature et technologies (FRQ-NT).

≃ +5.310 × 10−4 W/(m K2).

One question remains: what is the physical origin of this excess contribution? An initial answer to the question can be found in the large temperature dependence of the entropy of mixing close to the composition of the cryolite in the NaF−AlF3 system presented in Figure 9 of Chartrand and Pelton.2 Contrary to systems for which eq 19 has shown a good predictive ability,7,13 for the NaF−AlF3 system, it shows (1) a strong deviation from ideality in terms of entropy and enthalpy of mixing and (2) close to the cryolite composition, the enthalpy and entropy of mixing depend strongly on temperature. The temperature dependence of the entropy and enthalpy of mixing is the signature of a change in the vibrational state of the system as a function of temperature and can be induced by several factors. The development of theoretical models for describing the excess contribution of the temperature dependence of thermal conductivity is beyond the scope of the present paper and will be treated in a future contribution. It requires a deep analysis of the thermal conductivity of higher order local structures and its relationship with viscosity, for key systems for which the thermal conductivity increases with temperature.

■

REFERENCES

(1) Totten, G.; MacKenzie, D. Handbook of Aluminum: Alloy Production and Materials Manufacturing; Handbook of Aluminum; CRC Press, 2003; Vol. 2. (2) Chartrand, P.; Pelton, A. D. A predictive thermodynamic model for the Al-NaF-AlF3-CaF2-Al2O3 system. Light Met. (Warrendale, PA, U. S.) 2002, 245−252. CAPLUS AN 2002:208681(Journal). (3) Haupin, W. E. Calculating thickness of containing walls frozen from melt. JOM 1971, 23, 41−44. (4) Haupin, W. E. Hot wire method for rapid determination of thermal conductivity. Am. Ceram. Soc. Bull. 1960, 39, 139−41. (5) Morelli, D. T.; Slack, G. A. In High Thermal Conductivity Materials; Shindé, S. L., Goela, J. S., Eds.; Springer New York: New York, NY, 2006; Chapter High Lattice Thermal Conductivity Solids, pp 37−68. (6) DiGuilio, R.; Teja, A. A rough hard-sphere model for the thermal conductivity of molten salts. Int. J. Thermophys. 1992, 13, 855−871. (7) Gheribi, A. E.; Torres, J. A.; 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. (8) Gheribi, A. E.; Chartrand, P. Application of the {CALPHAD} method to predict the thermal conductivity in dielectric and semiconductor crystals. CALPHAD: Comput. Coupling Phase Diagrams Thermochem. 2012, 39, 70−79. (9) Gheribi, A. E.; Seifitokaldani, A.; Wu, P.; Chartrand, P. An ab initio method for the prediction of the lattice thermal transport properties of oxide systems: Case study of Li2O and K2O. J. Appl. Phys. 2015, 118, 145101. (10) Gheribi, A. E.; Chartrand, P. Effect of Grain Boundaries on the Lattice Thermal Transport Properties of Insulating Materials: A Predictive Model. J. Am. Ceram. Soc. 2015, 98, 888−897. (11) Ohtori, N.; Oono, T.; Takase, K. Thermal conductivity of molten alkali halides: Temperature and density dependence. J. Chem. Phys. 2009, 130, 044505. (12) Gheribi, A. E.; Poncsák, S.; St-Pierre, R.; Kiss, L. I.; Chartrand, P. Thermal conductivity of halide solid solutions: Measurement and prediction. J. Chem. Phys. 2014, 141, 104508. (13) Gheribi, A. E.; Chartrand, P. Thermal conductivity of molten salt mixtures: Theoretical model supported by equilibrium molecular dynamics simulations. J. Chem. Phys. 2016, 144, 084506. (14) Gheribi, A.; Corradini, D.; Dewan, L.; Chartrand, P.; Simon, C.; Madden, P.; Salanne, M. Prediction of the thermophysical properties of molten salt fast reactor fuel from first-principles. Mol. Phys. 2014, 112, 1305−1312. (15) Madden, P. A.; Heaton, R.; Aguado, A.; Jahn, S. From firstprinciples to material properties. J. Mol. Struct.: THEOCHEM 2006, 771, 9−18. Modelling Structure and Reactivity: the seventh triennial conference of the World Association of Theoritical and Computational Chemists (WATOC 2005).

■

CONCLUSION This study is part of the theoretical modeling of the thermal conductivity of the side ledge and molten electrolyte in aluminum electrolysis cells. Among the number of phases constituting the side ledge, only a few reliable data of thermal transport properties are available. Interionic potentials were developed based on the density functional theory (DFT) and a series of equilibrium molecular dynamic simulations carried out in both the NPT and NVT ensembles. The thermal conductivities of NaF (cubic, Fm3̅m), β-Na3AlF6 (cubic, Fm3m ̅ ), Na5Al3F14 (tetragonal, P4/mnc), and molten cryolite were determined as a function of temperature. For NaF and molten cryolite, for which data were available, the thermal conductivity obtained by EMD simulation was found to be close matched to experiment. For β-Na3AlF6 and Na5Al3F14, our calculations are purely predictive as no experimental data are reported in the literature. However, to analyze the reliability of the results generated for these two compounds, the thermal conductivity obtained by the present EMD simulations was compared to that obtained by a theoretical model which has shown a good predictive ability for many halides, oxides, and semiconductors. The theoretical model was also parametrized, based on a hybrid method; the parameters were determined on a DFT basis, then slightly adjusted to reproduce simultaneously the available experimental data for heat capacity and molar volume as functions of temperature. The agreement between EMD simulations and theoretical model results was found to be excellent. L


Article


properties of orthorhombic crystals: Application to TiSi2. J. Appl. Phys. 1998, 84, 4891. (40) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561. (41) Kresse, G.; Hafner, J. Ab initio molecular-dynamics simulation of ̆ the liquid-metalamorphous-semiconductor transition in germanium. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 49, 14251−14269. (42) Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Comput. Mater. Sci. 1996, 6, 15−50. (43) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 11169−11186. (44) Blöchl, P. E. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953−17979. (45) Kresse, G.; Joubert, D. From ultrasoft pseudopotentials to the projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1999, 59, 1758−1775. (46) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (47) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396−1396. (48) Petrov, A. V.; Tsypkina, N. S.; Logachev, Y. A. Temperature dependence of the thermal conductivity of alkali metal halides at elevated temperatures. Fiz. Tverd. Tela (Leningrad) 1974, 16, 65−70. (49) Smirnov, I. A. Thermal conductivity of sodium fluoride single crystals with potassium, lithium, and chlorine impurities. Fiz. Tverd. Tela (Leningrad) 1967, 9, 1845−7. (50) Haakansson, B.; Ross, R. G. Thermal conductivity of solid sodium fluoride under high pressure. Int. J. Thermophys. 1985, 6, 353− 65. (51) Fultz, B. Vibrational thermodynamics of materials. Prog. Mater. Sci. 2010, 55, 247−352. (52) Haigis, V.; Salanne, M.; Jahn, S. Thermal conductivity of MgO, MgSiO3 perovskite and post-perovskite in the Earth’s deep mantle. Earth Planet. Sci. Lett. 2012, 355−356, 102−108. (53) Dewing, E. W. Thermodynamics of the System NaF-AlF3: Part V. Solid Solution in Cryolite. Metall. Trans. B 1978, 9, 687−690. (54) Dewing, E. W. Thermodynamics of the system NaF-AlF3: Part VII. Non-stoichiometric solid cryolite. Metall. Mater. Trans. B 1997, 28, 1095−1097. (55) Khaïroun, S.; Tressaud, A.; Grannec, J.; Dance, J. M.; Yacoubi, A. Structural phase transitions in elpasolite- and cryolite-type fluorides. Phase Transitions 1988, 13, 157−163. (56) Spearing, D. R.; Stebbins, J. F.; Farnan, I. Diffusion and the dynamics of displacive phase transitions in cryolite (Na3AlF6) and chiolite (Na5Al3F14): Multi-nuclear NMR studies. Phys. Chem. Miner. 1994, 21, 373−386. (57) Bruno, M.; Herstad, O.; Holm, J. L. Calorimetric Investigations of Some Phases in the System Sodium Fluoride - Aluminium Fluoride. J. Therm. Anal. Calorim. 1999, 56, 51−57. (58) Dolejs, D.; Baker, D. R. Phase transitions and volumetric properties of cryolite, Na3AlF6: Differential thermal analysis to 100 MPa. Am. Mineral. 2006, 91, 97−103. (59) Saunders, N.; Miodownik, A. CALPHAD (Calculation of Phase Diagrams): A Comprehensive Guide; Pergamon Materials Series; Elsevier Science, 1998. (60) Callaway, J.; von Baeyer, H. C. Effect of Point Imperfections on Lattice Thermal Conductivity. Phys. Rev. 1960, 120, 1149−1154. (61) MURABAYASHI, M. Thermal Conductivity of Ceramic Solid Solutions. J. Nucl. Sci. Technol. 1970, 7, 559−563. (62) Rocquet, P.; Couzi, M.; Tressaud, A.; Chaminade, J. P.; Hauw, C. Structural phase transition in chiolite Na 5 Al 3 F 14. I. Raman scattering and X-ray diffraction study. J. Phys. C: Solid State Phys. 1985, 18, 6555. (63) Anovitz, L. M.; Hemingway, B. S.; Westrum, E. F.; Metz, G. W.; Essene, E. J. Heat capacity measurements for cryolite (Na3AlF6) and

(16) Salanne, M.; Marrocchelli, D.; Merlet, C.; Ohtori, N.; Madden, P. A. Thermal conductivity of ionic systems from equilibrium molecular dynamics. J. Phys.: Condens. Matter 2011, 23, 102101. (17) Gheribi, A. E.; Salanne, M.; Chartrand, P. Thermal transport properties of halide solid solutions: Experiments vs equilibrium molecular dynamics. J. Chem. Phys. 2015, 142, 124109. (18) Gheribi, A. E.; Gardarein, J.-L.; Rigollet, F.; Chartrand, P. Evidence of second order transition induced by the porosity in the thermal conductivity of sintered metals. APL Mater. 2014, 2, 076105. (19) Gheribi, A. E.; Gardarein, J.-L.; Autissier, E.; Rigollet, F.; Richou, M.; Chartrand, P. Experimental study of the thermal conductivity of sintered tungsten: Evidence of a critical behaviour with porosity. Appl. Phys. Lett. 2015, 107, 094102. (20) Gheribi, A. E.; Autissier, E.; Gardarein, J.-L.; Richou, M. Thermal transport properties of multiphase sintered metals microstructures. The copper-tungsten system: Experiments and modeling. J. Appl. Phys. 2016, 119, 145104. (21) Tang, K. T.; Toennies, J. P. An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients. J. Chem. Phys. 1984, 80, 3726. (22) Madden, P. A.; Heaton, R.; Aguado, A.; Jahn, S. From firstprinciples to material properties. J. Mol. Struct.: THEOCHEM 2006, 771, 9−18. Modelling Structure and Reactivity: the seventh triennial conference of the World Association of Theoritical and Computational Chemists (WATOC 2005). (23) Salanne, M.; Simon, C.; Turq, P.; Madden, P. A. Heat-transport properties of molten fluorides: Determination from first-principles. J. Fluorine Chem. 2009, 130, 38−44. (24) Heaton, R. J.; Brookes, R.; Madden, P. A.; Salanne, M.; Simon, C.; Turq, P. A First-Principles Description of Liquid BeF2 and Its Mixtures with LiF: 1. Potential Development and Pure BeF2. J. Phys. Chem. B 2006, 110, 11454−11460. (25) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (26) CPMD, CPMD, external link: http://www.cpmd.org/. Copyright IBM Corp. 1990−2008, Copyright MPI fuer Festkoerperforschung Stuttgart 1997−2001. http://www.cpmd.org/. (27) Minuit user’s guide. URL http://seal.cern.ch/documents/ minuit/mnusersguide.pdf (accessed 2004). (28) Materials project.2016: https://materialsproject.org/. (29) Nose, S. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (30) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177−4189. (31) Sindzingre, P.; Gillan, M. J. A computer simulation study of transport coefficients in alkali halides. J. Phys.: Condens. Matter 1990, 2, 7033. (32) Morelli, D.; Slack, G. In High Thermal Conductivity Materials; Shindé, S., Goela, J., Eds.; Springer: New York, 2006; pp 37−68. (33) Seifitokaldani, A.; Gheribi, A. E. Thermodynamically selfconsistent method to predict thermophysical properties of ionic oxides. Comput. Mater. Sci. 2015, 108, 17−26. (34) Seifitokaldani, A.; Gheribi, A. E.; Dollé, M.; Chartrand, P. Thermophysical properties of titanium and vanadium nitrides: Thermodynamically self-consistent approach coupled with density functional theory. J. Alloys Compd. 2016, 662, 240−251. (35) Blanco, M.; Francisco, E.; Luaña, V. GIBBS: isothermal-isobaric thermodynamics of solids from energy curves using a quasi-harmonic Debye model. Comput. Phys. Commun. 2004, 158, 57−72. (36) Jasiukiewicz, C.; Karpus, V. Debye temperature of cubic crystals. Solid State Commun. 2003, 128, 167−169. (37) Ghosh, G. A first-principles study of cementite (Fe3C) and its alloyed counterparts: Elastic constants, elastic anisotropies, and isotropic elastic moduli. AIP Adv. 2015, 5, 087102. (38) Joardar, P.; Chatterjee, S.; Chakraborty, S. Indian J. Phy. 1980, A54, 433. (39) Ravindran, P.; Fast, L.; Korzhavyi, P. A.; Johansson, B.; Wills, J.; Eriksson, O. Density functional theory for calculation of elastic M


Article

The Journal of Physical Chemistry C reactions in the system NaFeAlSiOF. Geochim. Cosmochim. Acta 1987, 51, 3087−3103. (64) Barin, I. Thermochemical Data of Pure Substances; Wiley-VCH Verlag GmbH, 2008; pp 1080−1237. (65) Zhou, Q.; Kennedy, B. J. High-temperature powder synchrotron diffraction studies of synthetic cryolite Na3AlF6. J. Solid State Chem. 2004, 177, 654−659. (66) Yang, H.; Ghose, S.; Hatch, D. M. Ferroelastic phase transition in cryolite, Na3AlF6, a mixed fluoride perovskite: High temperature single crystal X-ray diffraction study and symmetry analysis of the transition mechanism. Phys. Chem. Miner. 1993, 19, 528−544. (67) He, L.; Liu, F.; Hautier, G.; Oliveira, M. J. T.; Marques, M. A. L.; Vila, F. D.; Rehr, J. J.; Rignanese, G.-M.; Zhou, A. Accuracy of generalized gradient approximation functionals for density-functional perturbation theory calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 064305. (68) Stuve, J. M.; Ferrante, M. J. Low-temperature heat capacities and high-temperature enthalpies of chiolite (Na5Al3F14). 1980, 11. (69) King, E. G. Low Temperature Heat Capacities and Entropies at 298.15° K. of Cryolite, Anhydrous Aluminum Fluoride and Sodium Fluoride. J. Am. Chem. Soc. 1957, 79, 2056−2057. (70) Gomez-Cuevas, A.; Perez-Mato, J. M.; Bocanegra, E. H.; Couzi, M.; Chaminade, J. P. The structural phase transition in chiolite Na 5 Al 3 F 14. III. An ultrasonic study. J. Phys. C: Solid State Phys. 1988, 21, 3641. (71) Rocquet, P.; Couzi, M. Structural phase transition in chiolite Na 5 Al 3 F 14. II. A model for the pseudo-proper ferroelastic transition. J. Phys. C: Solid State Phys. 1985, 18, 6571. (72) Polyakov, P. V.; Mozhaev, V. M. Thermal conductivity of melts in the sodium fluoride-aluminum fluoride system. Teplofiz. Vys. Temp. 1975, 13, 661−3. (73) Filatov, E. S.; Khokhlov, V. A.; Kodintseva, A. O. Rasplavy 2006, 29−33. (74) Khokhlov, V. A.; Filatov, E. S.; Solheim, A.; Thonstad, J. Thermal conductivity in cryolitic melts - new data and its influence on heat transfer in aluminum cells. Light Met. (Warrendale, Pa.) 1998, 501−506. (75) DiGuilio, R.; Teja, A. A rough hard-sphere model for the thermal conductivity of molten salts. Int. J. Thermophys. 1992, 13, 855−871. (76) DiGuilio, R. M. The Thermal Conductivity of Molten Salts and Concentrated Aqueous Salt Solutions. Ph.D. thesis, Georgia Institute of Technology, Atlanta, 1991. (77) Takase, K.; Matsumoto, Y.; Sato, K.; Ohtori, N. Thermal conductivity in molten alkali halides: composition dependence in mixtures of (NaK)Cl. Mol. Simul. 2012, 38, 432−436. (78) Kamal, I.; McLaughlin, E. Pressure and volume dependence of the thermal conductivity of liquids. Trans. Faraday Soc. 1964, 60, 809− 816. (79) Takahashi, T.; Kikuchi, T. Porosity dependence on thermal diffusivity and thermal conductivity of lithium oxide Li2O from 200 to 900°C. J. Nucl. Mater. 1980, 91, 93−102.

N


Formulation of Temperature-Dependent Thermal ... - ACS Publications

Recommend Documents