Nonequilibrium Molecular Dynamics Study of Heat Conduction in Ionic

of ionic systems was investigated using nonequilibrium molecular dynamics simulations. ... The Journal of Physical Chemistry B 2018 122 (5), 1662-...
1 downloads 0 Views 520KB Size
Download PDF
J. Phys. Chem. 1996, 100, 1879-1888

1879

Nonequilibrium Molecular Dynamics Study of Heat Conduction in Ionic Systems§ Fernando Bresme,† Bjørn Hafskjold,*,‡,| and Inge Wold⊥ Departamento de Quı´mica-Fı´sica, Facultad de CC. Quı´micas, UniVersidad Complutense de Madrid, E-28040, Spain, Instituto de Quı´mica-Fı´sica Rocasolano, CSIC, Serrano 119, E-28006 Madrid, Spain, Cluster Science Group, National Institute for AdVanced Interdisciplinary Research, 1-1-4, Higashi, Tsukuba, Ibaraki 305, Japan, and Department of Physical Chemistry, UniVersity of Trondheim-The Norwegian Institute of Technology, N-7034 Trondheim, Norway ReceiVed: May 2, 1995X

The effects of charge strength and potential range on the thermal diffusion factor and thermal conductivity of ionic systems was investigated using nonequilibrium molecular dynamics simulations. The potential model was a Lennard-Jones/spline function with an ionic tail of variable range and magnitude. The range was varied by selecting either a Coulomb or a Yukawa tail, and the magnitude was varied with the ionic charge strength. The study covers the range from a simple binary fluid mixture in the neighborhood of its critical point to a molten salt in the supercritical regime. Our main conclusions are that, under the conditions used here, the ionic charge inhibits the Soret effect in the Coulomb system (except when the valence is less than 0.02) but not in the Yukawa system. The thermal diffusion factor decreases from +1.03 ( 0.13 for zero charge to -0.036 ( 0.007 for univalent ions in the Coulomb case and to +0.098 ( 0.004 in the Yukawa case. (The uncertainties are given as 1 standard deviation of the mean.) The thermal conductivity appears to be insensitive to the potential range for ionic systems. An analysis of the flux contributions reveals that the major contributions to the heat flux are the flux of kinetic energy and the collision energy flux. The potential energy flux contribution is not important, but it increases significantly as the system approaches its critical region.

1. Introduction The transport properties of molten salts are very different from those of nonionic liquids. For instance, the electric conductivity of molten KCl is of the order of 108 times higher than that of pure water (close to the triple points for both liquids), and the thermal conductivity is about 4 times larger than that of argon at similar conditions. Despite the technological importance of molten salts and electrolyte solutions in e.g. metallurgy and batteries, transport properties other than the electric conductivity are still relatively unexplored by experimental and theoretical methods as compared to nonionic liquids. There is, however, an increasing interest in thermal transport data in design methods for electrochemical process units.1 For instance, it was recently pointed out that local heat effects at the electrodes in aluminum electrolysis are substantial2 and may therefore be important for the design and operation of electrolysis cells. Phenomenological descriptions are well developed for dilute electrolytes but not for transport phenomena at higher concentrations and in molten salts. This is mainly due to the existence of cooperative effects that are not contemplated in the theories of electrolyte solutions, based on a single-ion transport process.3 However, recent works on electrolytes near interfaces4-6 and the thermodynamic and structural properties of electrolytes near the critical point7-9 have exposed the different natures of ionic and nonionic systems much more clearly than before. These differences are attributed to the difference in structure, inter§ This paper was originally submitted for the Harold Friedman Festschrift [J. Phys. Chem. 1996, 100 (4)]. * Author to whom correspondence should be addressed at University of Trondheim-The Norwegian Institute of Technology. † Universidad Complutense de Madrid and Instituto de Quı´mica-Fı´sica Rocasolano. ‡ National Institute for Advanced Interdisciplinary Research. | On leave from the University of Trondheim-The Norwegian Institute of Technology. ⊥ University of Trondheim-The Norwegian Institute of Technology. X Abstract published in AdVance ACS Abstracts, December 15, 1995.

0022-3654/96/20100-1879$12.00/0

molecular energies, and range of the potential.10 It is difficult to separate the relative importance of these factors by experimental methods. For instance, the strength and range of the Coulomb potential are not separable, and they dominate the thermodynamic properties of electrolyte solutions and molten salts. The typical molten-salt short-range ordering (large fluctuations in the charge-charge radial distribution function 1/ [g 2 ++(r) - g+-(r)]) is, however, not a consequence of the long range of the ionic potential but rather the fact that like ions repel each other and unlike ions attract each other at short range; an ionic Yukawa model gives essentially the same structure as a Coulomb model at a corresponding state.11 In order to further investigate the different roles of the strength and range of the ionic potential on the system’s properties, it seems worthwhile to study the transport properties by computer simulations. We have chosen to carry out nonequilibrium molecular dynamics (NEMD) simulations of two ionic models, one with the Coulomb potential and the other with the Yukawa potential as the long-range tails, with the objective to search for systematic changes in the thermal conductivity and thermal diffusion factor as a function of the charge strength and range of the interionic potential. Thermal diffusion (the Soret effect) is in itself an interesting phenomenon, and even for simple fluids there is currently no theory that can predict even the sign of the effect.12 That is, given a binary mixture of two simple fluids subject to a temperature gradient and the thermodynamic state, one cannot a priori tell which of the components will migrate to the hot side of the system or how extensive the migration is. A molten salt is in the Gibbs phase rule sense a one-component system because of the electroneutrality condition, in which there should be no Soret effect. The electroneutrality condition is closely related to the long range of the Coulomb potential.13 However, if the charge strength is reduced, the system becomes a twocomponent system in the limit of zero charge. To what extent © 1996 American Chemical Society

1880 J. Phys. Chem., Vol. 100, No. 5, 1996

Bresme et al.

an experimentally feasible temperature gradient is strong enough to separate the charged species and polarize the ionic system is not known. More specifically, the question is: How does the thermal diffusion factor depend on the ionic charge and the range of the ionic potential? A considerable amount of theoretical studies and computer simulations have been done to understand the thermal conductivity and thermal diffusion in simple and molecular liquids.12,14-18 NEMD methods have proved to be a useful tool in such studies, and a number of computer algorithms have been developed for this purpose.19 We have in this work applied the algorithm discussed by Ikeshoji and Hafskjold.17 For the cross-transport coefficients, such as the thermal diffusion factor, NEMD is more efficient than equilibrium MD and use of the Green-Kubo formula.10 Besides, NEMD gives a more direct simulation of experimental situations and, as such, is more instructive in interpretations of experimental data. Application of NEMD techniques to the computation of transport properties in ionic systems is so far limited to some work by Svishchev and Kusalik20 concerned with the computation of electric conductivity of primitive electrolyte solutions and molten salts. As far as we know, no previous studies of heat conduction using NEMD have been done for ionic systems. The remainder of this paper is organized as follows: Section 2 contains the formal relations for the transport processes in ionic systems, which include the definitions of the thermal diffusion factor and thermal conductivity as well as the expressions for the molecular and phenomenological fluxes and their relationships. In section 3 we describe the potential models and the simulation details. Section 4 is divided into three parts, where we first discuss the assumption of local equilibrium in the NEMD system. Next, the features of thermal diffusion and heat conduction in the Yukawa and Coulomb systems are discussed. Finally, we consider in detail the molecular mechanisms for heat conduction.

our binary mixture eq 1 simplifies to

σ)-

σ)-

1

Jq‚∇T -

T2

1

2

∑Jk‚{∇Tµk - Fk}

Tk)1

(1)

where it is assumed that the system is not subject to viscous forces, external magnetic fields, or electronic current and that no chemical reactions take place. In eq 1, ∇T is the temperature gradient, µk is the chemical potential of component k, and the subscript T in the gradient denotes that the temperature dependency of µk is not included in the gradient. The mass flux Jk used here refers to a barycentric frame of reference. The charges, zk, carried by the particles appear in the definition of the external forces Fk

Fk ) zkE

2

I ) ∑zkJk

(4)

k)1

The phenomenological equations for the fluxes and the thermodynamic forces of eq 3 assuming linear relations are given by22

∇T 1 1 Jq ) -Lqq 2 - Lq1 ∇T(µ1 - µ2) + Lqe E T T T

(5)

∇T 1 1 J1 ) -L1q 2 - L11 ∇T(µ1 - µ2) + L1e E T T T

(6)

∇T 1 1 - Le1 ∇T(µ1 - µ2) + Lee E 2 T T T

(7)

I ) -Leq

where the L’s represent the phenomenological coefficients. The NEMD simulations were carried out so that a stationary state state with zero mass flux was reached, which led to zero electric current I as defined by eq 4. Besides, no external electric field was applied to the system, and the third term in eqs 5-7 disappears. Taking this into account we get from eq 6

dµ1 ∇T(µ1 - µ2) -L1q ) ) -sT ∇T L11T d ln w1

(8)

where sT is the Soret coefficient and w1 is the weight fraction of component 1. The thermal diffusion factor is defined by

(

R12 ) sTT ) -

)

∇ ln(w1/w2) ∇ ln T

(

)J1)0

)

∇ ln(x1/x2) ∇ ln T

J1)0

(9)

where x1 and x2 are the molar fractions of species 1 and 2, respectively, R12 ) -R21, and R11 ) R22 ) 0. The thermal diffusion factor can be readily computed in computer simulations, provided that both a temperature gradient and a composition gradient exist. The thermal conductiVity, λ, may be defined in terms of Jq by

Jq ) -λJ1)0∇T

(10)

The thermal conductivity is somewhat more difficult to compute because it requires the heat flux. The computation of Jq is made indirectly via the so-called internal energy flux, JU, which for a local control volume (CV) is given in MD as19

JU ) JKIN. + JPOT. + JCOL.

(11)

where

(2)

where E is the static electric field operating on the system. For

(3)

where we have eliminated the mass flux of species 2 (J2), 2 Jk ) 0. The total electric making use of the relation ∑k)1 current I is given by

2. Formal Expressions for the Transport Processes This paper deals with charged systems subject to a temperature gradient. Each system is composed of an equal number of anions and cations of the same valence. We shall initially treat this as a two-component system and examine the consequences of the electroneutrality condition later on. Three nonconvective transport processes may be observed in this system; heat flux Jq, diffusive mass flux Jk (of a given component k), and electric current I, i.e., charge flux. The expression for the entropy production, σ, is a convenient starting point for describing these processes by irreversible thermodynamics,21

1 1 1 Jq‚∇T - J1‚∇T(µ1 - µ2) + I‚E T T T2

JKIN. )

1

1

∑ mi(vi - v)2(vi - v) Vi∈CV 2

(12)

Heat Conduction in Ionic Systems

JPOT. )

JCOL. ) -

1

1

J. Phys. Chem., Vol. 100, No. 5, 1996 1881

∑ φi(vi - v)

Vi∈CV

(13)

N

∑ ∑((vi - v)‚Fij)rij

2Vi∈CV j)1

(14)

In eqs 12-14, V is the size of the control volume, mi and vi are the mass and velocity, respectively, of particle i, v is the barycentric velocity of the system, φi is the potential energy of particle i in the field of all the other particles, Fij is the force acting on i due to j, rij is the vector from the position of i to that of j, and N is the total number of particles in the system. The internal energy flux and the heat flux are related through

Jq ) JU - (h1 - h2)J1

(15)

where hk is the partial specific enthalpy of component k. When the stationary state has been reached, the mass flux J1 equals zero and the potential contribution to the energy flux (cf. eq 13) becomes independent of the reference state. Under these conditions, Jq ) JU, and we get an expression that relates the measurable, macroscopic heat flux to the microscopic internal energy flux which can be computed in MD simulations. The assumption that a small system adequately represents the thermodynamic limit is implicit in this relation and will be discussed in section 4.1. 3. Nonequilibrium Molecular Dynamics Simulations 3.1. Potential Models. Our objective in this paper is to search for systematic changes in the thermal diffusion factor and thermal conductivity as a function of the charge strength and range of the ionic potential. For this purpose we have considered a potential model that can be divided into two contributions, a short range nonionic contribution modeled by a Lennard-Jones/spline (LJ/s) potential and an ionic long-range contribution,

φij(r,z) ) φij,LJ/s(r) + φij,ionic(r,z)

(16)

The LJ/s contribution is given by23

φij,LJ/s(r) )

{

4ij [(σij/r)12 - (σij/r)6] aij(r - rcij)2 + bij(r - rcij)3 0

if 0 < r < rsij if rsij < r < r cij (17) if rcij < r

where r is the distance between particles i and j, rsij ) (26/7)1/6σij, rcij ) (67/48)rsij, aij ) -(24192/3211)(ij/rsij2), and bij ) -(387072/61009)(ij/rsij3). In this work the two species differ only by their charge and mass. The LJ/s potential becomes that of an isotopic mixture, and the subscript ij may be discarded. This LJ/s potential has been used in previous NEMD works;12,24 therefore it will be a good reference in our forthcoming discussion. Two models were considered for the ionic contribution; one is the Coulomb potential given by

φij,Coul(r,z) )

e2 zizj 4π0 r

Figure 1. Model potential used in this work for z ) 0 and z ) 0.1. For z ) 0.1, only the Coulombic case is shown.

(18)

where e is the electronic charge, 0 the permittivity of vacuum, and zi the valence of ion i. The second is the Yukawa potential,

φij,Yuk(r,z) )

[ (

)]

r e2 zizj -1 exp -κ 4π0 r σca(z)

(19)

where κ controls the screening of the potential. In all cases, we set z1 ) -z2 ) z. To distinguish between the different potentials, we will denote them by LJ/s when z ) 0 and Coulomb potential or Yukawa potential (including the shortrange LJ/s part) for z > 0, depending on which long range model is used for the ionic interactions. The Coulomb potential for z ) 0.1 is shown in Figure 1. The interparticle distance appears in reduced units of the atomic diameter of the uncharged system, i.e., r* ) r/σ. For convenience we have chosen σ ) 3.4 Å, the diameter of argon, which is isoelectronic with K+ and Cl-. In Figure 1, the energy units also refer to the Ar parameters; i.e., φij* ) φij(r,z)/ where  ) 1.7 × 10-21 J. For this choice, the depth of the total potential for cation-anion interactions, *ca(z) ) ca(z)/, ranges from 1 for z ) 0 up to 456.74 for z ) 1.0. From Figure 1 we observe that σ*ca(z) ) σca(z)/σ, defined as the distance at which the cation-anion potential is zero, diminishes as the charge is increased. Because the cationanion distance is the most characteristic length parameter in the system, we will define the reduced number density as F*(z) ) Nσca3(z)/V. In a system with a fixed number of particles, N, and fixed volume, V, this would mean that the reduced density would decrease as the charge increases, which would give an undesired contribution to the changes in the system’s properties as the charge strength varies. We shall instead impose the constraint that F*(z) ) constant. This condition may be fulfilled by changing the volume of the system until the desired value of the reduced density is reached. With the definitions of σca(z) and ca(z), it is possible to construct a set of reduced variables as is usually done for other potentials. Some selected variables are compiled in Table 1. The dependence of σca(z) on the ionic charge is depicted in Figure 2. As in the previous figure, σ*ca(z) and the reduced distance of the potential minimum, r*min ) rmin/σ, refer to the Ar parameter, σ ) 3.4 Å. The main changes in the reduced distances occur at low values of z. An inflection point is observed around z ) 0.15. It is not clear what physical significance these features have. By using σca(z) as a parameter in the Yukawa potential, eq 19, we ensure that both the Coulomb and the Yukawa potentials have zero value at r ) σca(z), and consequently the same ionic diameter. Two additional conditions were imposed on the Yukawa potential in this report. First the potential depth should be equal to the Coulomb potential depth, and the value of the potential at the cutoff (half of the shortest box length in this work) should be smaller than 1.5% of the potential depth. The

1882 J. Phys. Chem., Vol. 100, No. 5, 1996

Bresme et al.

Figure 2. The potential parameters σ*ca and r*min as a function of charge.

TABLE 1: Definitions of Reduced Quantities Used in This Work quantity

definition

density temperature pressure time energy energy flux thermal conductivity

ρ* ) (N/V)σca3(z) T* ) kBT/ca(z) P* ) (P/ca(z))σca3(z) t* ) (t/ca(z))(σca(z)/m1)1/2 U* ) U/ca(z) J*U ) JU(σca3(z)/ca(z))(m1/ca(z))1/2 λ* ) λ(σca2(z)/kB)(m1/ca(z))1/2

TABLE 2: Simulation Conditions and Thermal Properties run J*U 1a 1b 1c 2a 2b 3 4a 4b 5 6 7a 7b 8 9 10 11

0.10 0.10 0.10 0.10 0.12 0.12 0.10 0.15 0.08 0.08 0.10 0.10 0.20 0.20 0.20 0.30

steps/ 106

z

1 1 1 2 2 1 2 2 0.3 0.3 0.3 1 1 1 1 1.2

0.00 0.00 0.00 0.01 0.01 0.02 0.05 0.05 0.25 0.50 1.00 1.00 0.05 0.10 0.50 1.00

κ

λ*

R12

potential

0.7 1.0 1.3 1.6

2.08 ( 0.09 2.05 ( 0.13 2.05 ( 0.11 2.08 2.01 ( 0.08 2.02 ( 0.09 1.86 ( 0.07 1.90 ( 0.05 1.70 ( 0.10 1.76 ( 0.09 1.76 ( 0.02 1.74 ( 0.05 1.87 ( 0.02 1.81 ( 0.02 1.88 ( 0.02 1.87 ( 0.01

1.06 ( 0.20 0.89 ( 0.22 1.14 ( 0.25 0.427 0.35 ( 0.09 0.10 ( 0.07 -0.02 ( 0.03 -0.03 ( 0.02 -0.01 ( 0.04 -0.03 ( 0.04 -0.04 ( 0.01 -0.03 ( 0.01 0.23 ( 0.02 0.13 ( 0.01 0.089 ( 0.008 0.098 ( 0.004

LJ/s LJ/s LJ/s Coulomb Coulomb Coulomb Coulomb Coulomb Coulomb Coulomb Coulomb Coulomb Yukawa Yukawa Yukawa Yukawa

conditions were attained by varying κ and scaling the Yukawa potential (including the LJ/s contribution). In this respect the only difference between the potentials is their range. The resulting κ-values are included in Table 2. 3.2. Heat Exchange Algorithm. The NEMD simulations were made with the heat exchange (HEX) algorithm previously used in studies of Lennard-Jones mixtures.17 The method is described in the following to the extent necessary for the subsequent discussion. The objective is to generate an energy flux, JU, through the system. The MD cell was divided into 32 layers, perpendicular to the direction of the flux. One layer at each end of the box, hereafter referred to as the hot layers, was used for addition of kinetic energy, and two layers in the middle of the box, the cold layers, were used as a sink of kinetic energy (see Figure 3). This method gives a simulation cell that is fully periodic. At each time step a specified amount of kinetic energy, ∆U, was added or subtracted by scaling and shifting the velocities of all the particles in the hot and cold layers, respectively, under the constraint that the momentum for the total system is conserved (and equal to zero by virtue of the initial condition). In this way the barycentric velocity v of the system was also

Figure 3. Layout of the simulation cell. H and C denote hot and cold layers as described in the text. The direction of the heat flux is depicted as well.

zero. The barycentric velocity in each layer was zero at stationary state,18 and the temperature in each layer was computed from the average kinetic energy according to the equipartition principle. Assuming local equilibrium, this quantity was interpreted as a thermodynamic temperature. As a consequence of this procedure, the system developed a gradient of temperature ∇T and an internal energy flux JU. At stationary state the internal energy flux must fulfill the continuity equation,

{

}

∆U JU ) ( ,0,0 2δtS

(20)

where the right-hand side is the source term. The + and signs refer to 0 < x < Lx/2 and Lx/2 < x < Lx, respectively, δt is the time step, and S is the cross-sectional area of the box, i.e., LyLz. Equation 20 served as a test to check the computation of the local fluxes from eqs 11-14. 3.3. Simulation Conditions. The simulations were done at overall temperature T* ) 1.0 and overall number density F* ) 0.4,25 with a mass ratio m1/m2 ) 10. Component 1 was chosen to be the cation. At these conditions the ionic charge was varied from z ) 0 to z ) 1.0, which corresponds to a supercritical molten salt. The number of particles, N, was 512, 256 cations and 256 anions, arranged in a noncubic box, with the proportions Lx ) 2Ly ) 2Lz()L). The runs were started either from a final configuration from a previous NEMD run or from a configuration of an equilibrium simulation. The runs were carried out for (0.3-2) × 106 time steps, allowing a stationary state to be established during the first (0.1-0.4) × 106 time steps. A list of the conditions is given in Table 2. For z > 0.02, we noted that the small Soret effect caused that only 0.3 × 106 steps were enough to obtain reliable results of the thermal conductivity and thermal diffusion factor (as shown by runs 7a and 7b in Table 2). The thermodynamic conditions chosen for the LJ/s system are in the neighborhood of the critical point, which gives large density fluctuations and necessitates longer runs to obtain reliable results. Besides, some care must be taken when setting the energy flux so as to avoid thermodynamic states in the two-phase region.26 The equations of motion were integrated by means of the “velocity Verlet” algorithm27 with time steps δt* ) 0.002 or 0.005, depending on the system. The drift in total energy was about 6% for a run of 2 × 106 time steps, similar to that observed in a previous work.18 The long-range interaction of the Coulomb potential was computed by the Ewald summation method.27 The expressions used to compute the collision energy flux contribution (eq 14) are described in the Appendix. The parameters for the Ewald B|2 < 38, sum were chosen to be B hmax ) {(6,(3,(3}, with |h giving a total of 274 reciprocal vectors. The convergence parameter in the Ewald sum expressions (see Appendix) was R

Heat Conduction in Ionic Systems

J. Phys. Chem., Vol. 100, No. 5, 1996 1883

Figure 5. Pressure profile along the cell for two systems at zero charge. The state at P*1 is located near the critical point, whereas that at P*2 represents a supercritical state.

Figure 4. Temperature, density, and mole fraction of component 1 along the simulation cell as a function of charge strength for the Coulombic potential. The other variables are m1/m2 ) 10.0 and z ) 0.00, 0.01, 0.05 (runs 1b, 2b, and 4a, respectively). The lines in plot c represent the linear regressions used to compute the thermal diffusion factor.

) 4.95/L. With this set of parameter values, the error in the total electrostatic energy of a representative configuration is lower than 0.002%, which we considered acceptable for our purposes. 4. Results 4.1. Gradients, Local Equilibrium, and Equation of State. The stationary state is characterized by a temperature gradient, a density gradient, and in some cases a composition gradient. These gradients are shown for the LJ/s system and two Coulomb systems at low ionic charge, z ) 0.01 and 0.05, in Figure 4. The vicinity of the critical point for z e 0.01 causes the cold layer to be close to a liquid state whereas the hot layer is close to a gas. Given that the liquid is a better heat conductor than the gas, this explains the curvature of the temperature profile. For the same reason, the density profile exhibits a hyperbolic tangentlike shape (cf. Figure 4). The temperature and density profiles are rather insensitive to increments in the charge strength for z < 0.05, but at z ) 0.05, the profiles are almost linear because the system is moved away from the critical region. Although the systems at z ) 0 and 0.01 present very similar temperature and density gradients, their composition profiles show significantly different behaviors (see Figure 4c). This is a clear manifestation that the charge strength inhibits the Soret effect (a further discussion of this is given in the next section). A basic assumption in this work is that local equilibrium prevails in each layer of the MD cell. The analyses by Hafskjold and Ratkje24 for an isotope LJ/s mixture and by Hafskjold et al.28 for a nonideal LJ/s mixture show that this assumption is valid for external forces that are stronger than those considered here. A useful criterion for local equilibrium is that the local thermodynamic properties fulfill the equation of state (EOS) for the system. In the following, we shall test the assumption of local equilibrium for two states of the LJ/s system. For fluid states, the pressure remains uniform along the cell. This means that if there is local equilibrium, the EOS is obtained along an isobar. Figure 5 shows the pressure profile along the left half of the cell for two NEMD simulations, one corresponds to conditions near the critical point (P*1 at T* ≈ 1.0) whereas

Figure 6. Internal energy and density vs temperature for pressures P*1 and P*2 from NEMD simulations (white circles). For P*1, NPT molecular dynamics (black circles) and HMSA integral equation (lines) results are also shown.

the other (P*2 at T* ≈ 1.5) corresponds to a system well into the supercritical region. For each condition, the pressure is found to be constant along the cell, except in the neighborhood of the cold and hot layers, where the EOS is not fulfilled. Parts a and b of Figure 6 show the internal energy and density vs temperature for the two pressures mentioned above. In the same figure, we have included results obtained from an integral equation using the HMSA closure,29 which is very accurate in describing the phase diagram of the LJ system.30 The internal energy obtained by NEMD is in excellent agreement with the HMSA equation at the supercritical pressure, but some disagreement is found at P*1. In order to assign this discrepancy to a lack of local equilibrium or a system-size effect in the NEMD simulation, or a deficiency in the theory, we have performed several NPT molecular dynamics runs using the Evans and Morriss algorithm.31 This ensemble seems the most suitable for our purposes, because just by fixing the pressure and the temperature, one can get the density as a result of the simulation. These plots show that the NEMD method develops density and temperature values that fulfill the equation of state (except near the hot and cold layers). The agreement between NPT and NEMD is very good at high densities (low temperatures), whereas, for low densities (high temperatures), there is some discrepancy between the two sets

1884 J. Phys. Chem., Vol. 100, No. 5, 1996 of MD results. Actually, the agreement between the NPT results and the theory is better than between the two MD results at low densities. The deviation between NPT and NEMD results observed at low densities may be explained by the deviation of the pressure from the average value observed as one approaches the hot layer (cf. Figure 5). This deviation has its origin in the proximity of this thermodynamic state to a phase boundary. It is also reflected in the system’s volume expansivity (β ) -1/F(dF/dT)P), found as the slope of the curves in Figure 6. A visual inspection shows that β increases from P*2 to P*1 as one should expect from the divergence in the isothermal compressibility, κT, at the critical point. Note that β is related with the isothermal compressibility through the thermal pressure coefficient γV, κT ) β/γV. At the critical point γV is positive and finite, whereas κ is positive and infinite. It follows that β is positive infinite at the critical point. In summary, the fulfillment of the EOS warrants the computation of the thermal conductivity along an isobar in just one simulation. Preliminary computations on molten salts show that thermal conductivities obtained in this way behave consistently for pairs of density and temperature.32 Besides, the NEMD method is able to supply important information other than transport coefficients; therefore, it seems worthwhile to investigate the applicability of this NEMD algorithm to the study of the EOS of simple systems. 4.2. Thermal Diffusion Factor and Thermal Conductivity. The thermal diffusion factor R12 is readily obtained from eq 9, provided temperature and composition gradients are present along the simulation cell. Similarly, the thermal conductivity can be obtained from eq 10. Numerical results for both quantities are compiled in Table 2. Because of the proximity to the two-phase region of the phase diagram when z < 0.02, we have to use relatively small internal energy fluxes to avoid a phase separation. Table 2 contains data for three independent simulations with J*U ) 0.10 for z ) 0, runs 1a-c. We observe that the thermal conductivity and the thermal diffusion factor are reproducible, but the uncertainties in the thermal diffusion factor are large because of the large density and composition fluctuations under these thermodynamic conditions. Note that an increase in the charge reduces the uncertainties in R12 for approximately the same J*U. The transport properties were not found to depend significantly on the imposed internal energy flux, as shown by runs 2a-b and 4a-b. This is consistent with previous work,18 where we observed that the transport properties are independent on the imposed flux under similar conditions. Finally, we note that, for large enough charges, one can obtain reliable properties both of the thermal diffusion factor and thermal conductivity using relatively short runs; see runs 7a and 7b. Figure 7 shows the thermal diffusion factor as function of charge, both for the Coulomb and Yukawa potentials. The thermal diffusion factor is positive for neutral and weakly charged systems, which means that the lighter component migrates to the hot region (after the definition adopted in this work for R12; see eq 9), a fact noted in the 19th century in experiments concerned with gas mixtures. The thermal diffusion factor decreases sharply with increasing ionic charge, which inhibits the Soret effect. In a formal sense, an increase in charge strength has the effect of interpolating between a two-component mixture (at zero charge) and a one-component system (at nonzero charge), the latter being a consequence of the electroneutrality condition. Hence, a reduction of the thermal diffusion factor with increasing charge is expected since the Soret effect does not exist for a one-component system. From Figure 7 one can appreciate that the thermal diffusion factor reaches es-

Bresme et al.

Figure 7. Dependence of the thermal diffusion factor on the ionic charge for Coulomb and Yukawa systems at F* ) 0.4, T* ) 1.0, and m1/m2 ) 10.0.

sentially zero value for the Coulomb system at charges as small as 0.05, and only for very small charges, such as 0.01-0.02, is the thermal gradient able to polarize the system. These changes occur in a continuous way; i.e., the thermal diffusion factor is a smooth function of the charge strength. For z > 0.25, the Coulomb system shows slightly negative values for R12, but they are just at the limit of being significantly different from zero (given the statistical uncertainties represented here as 1 standard deviation of the mean). We believe that the correct value is zero and that the negative values may be a consequence of the linear regression used in the computation of R12. Figure 7 also shows the results obtained with the Yukawa potential. As in the Coulomb case, a drop in the thermal diffusion factor is observed, although it now decreases less drastically and never reaches zero. Even for z ) 1.0, the value of R12 is of the order 0.1 for the Yukawa system. There are at least two different interpretations of the behavior observed for the thermal diffusion factor as the ionic charge is increased. The first one is based on a consideration of the energy. Increasing the charge (whether it is for the Coulomb or the Yukawa system) will favor unlike particles as nearest neighbors and reduce the tendency for the two components to separate. The increased repulsion between like particles works in the same direction because each species will tend to spread uniformly over the available volume. The result is that the thermal diffusion factor decreases with increasing charge strength. The second interpretation is based on the electroneutrality condition in the Coulomb case. The observed difference in the thermal diffusion factor between the Coulomb and Yukawa systems must be attributed to the different long range nature of the models. The charge-charge structure factors Szz(k) of the Coulomb and Yukawa systems differ in the limit of zero k values. A Coulomb system must fulfill Szz(0) ) 0, whereas Szz(k) * 0 for the Yukawa system.33 This implies that in a Coulomb system the long-wavelength charge fluctuations, which in our context represents the system’s tendency to polarize, are suppressed, but not in the Yukawa system. This statement can also be represented by the Stillinger-Lovett moment conditions.13 The zeroth moment condition (the local electroneutrality condition) applies to potentials varying asymptotically as r-n, where 1 e n e 3, i.e., not to the Yukawa system. The second moment condition is strictly valid for potentials that vary asymptotically as 1/r, i.e., the Coulomb potential. In a formal sense, the Yukawa system represents a mixture, whereas the Coulomb system represents a one-component system. This explains why the two systems can have different values of R12 for large charge strengths. In addition to the thermal diffusion factor, we report in Table 2 the thermal conductivity, λ*, as function of charge. Figure 8

Heat Conduction in Ionic Systems

J. Phys. Chem., Vol. 100, No. 5, 1996 1885

Figure 8. Dependence of the thermal conductivity on the ionic charge, for Coulomb and Yukawa systems at F* ) 0.4, T* ) 1.0, and m1/m2 ) 10.0.

TABLE 3: Results of Thermal Conductivity Using Different Potential Depths To Get the Reduced Quantities (λ*′ ) λ*(Eca(z)/E)1/2) run

z

*ca(z) ) ca(z)/

λ*

λ*′

1a 1b 1c 2a 2b 3 4a 4b 5 6 7a 7b

0.00 0.00 0.00 0.01 0.01 0.02 0.05 0.05 0.25 0.50 1.00 1.00

1.0 1.0 1.0 1.036 1.036 1.142 1.894 1.894 24.943 104.277 456.741 456.741

2.08 2.05 2.05 2.08 2.01 2.02 1.86 1.90 1.70 1.76 1.76 1.74

2.08 2.05 2.05 2.12 2.05 2.16 2.56 2.61 8.49 17.97 37.61 37.19

shows λ* for the Yukawa and Coulomb systems, as well as the LJ/s system. As pointed out in the Introduction, the thermal conductivity of an ionic system such as KCl is larger than that of a neutral system such as Ar. The decrease in λ* with increasing z at small values of z must be attributed to the different potential depth, ca(z), used to represent the thermal conductivity in reduced units (cf. Table 1). If we use a constant potential depth, , and constant diameter, σ, irrespective of the charge, this effect is eliminated and we get the new values λ*′ reported in Table 3. These values increase monotonically with the ionic charge. For high charges (z g 0.25), the increase is proportional to the potential depth (note the constancy of λ* in Figure 8), which suggests that a corresponding states law would apply for the thermal conductivity with these potentials. As a test of the results for the thermal conductivity, we have compared the data for the LJ/s system with experimental data for Ar in the critical region.34 Our value for the conditions F ≈ 1.14Fc and T/Tc - 1 ≈ 0.1 is λ ) 39.4 mW m-1 K-1 vs λ ≈ 39.5 mW m-1 K-1 for Ar at F ) 1.169Fc and T/Tc - 1 ) 0.1, where we have used the σ and  values of Ar given in section 3.1. The agreement is excellent, although a proper comparison should involve a larger set of data. It would be interesting to extend the computations of λ* to a larger set of densities and temperatures in order to assess to what extent the LJ/s potential is a suitable model to describe simple systems, such as argon, but this will not be pursued here. A remarkable feature of the thermal conductivity and the thermal diffusion factor is the difference in their sensitivity to the range of the potential (note that the thermal conductivity for the Yukawa and Coulomb systems agree to within the combined uncertainties in the entire charge-strength range). Whereas the thermal diffusion factor is strongly influenced by the electroneutrality condition, the heat flux of a dense fluid is dominated by the term JCOL. (see eqs 11 and 14 and section

Figure 9. Collision, kinetic, and potential contributions to the total energy flux as a function of charge strength for the Coulomb system. The symbols have the same meaning as in Figure 4.

4.3). Clearly, this term will in general be different for Yukawa and Coulomb systems, but the results show that, provided the thermodynamic state is suitably chosen, the reduced thermal conductivity will be almost the same for the two systems. This is another manifestation of the law of corresponding states that exists between the Yukawa and Coulomb systems, and which is most clearly seen as a similarity in the short-range structure of the two systems.11 It remains to be seen if such coincidence is observed at other conditions, particularly at densities lower than those studied in this report, where the heat flux is dominated by JKIN.. 4.3. Fluxes and Transport Mechanisms. In this section, all flux contributions are given as reduced (i.e., divided) by the imposed energy flux; therefore, the value of the imposed J*U is ideally 1 in these units. Hereafter, we omit the asterisk that denotes the reduced units. Figure 9 shows the three contributions JKIN., JCOL., JPOT. as functions of ionic charge. The quantities are averages over the two half-cells (making use of the system’s symmetry) for the three systems presented in Figure 4. Irrespective of the charge, the kinetic contribution to the total energy flux decreases from the hot to the cold layer, whereas the collision flux increases in the opposite direction. As observed in a previous paper on LJ systems,18 the ratio JCOL./ JKIN. increases from the hot to the cold layer. These results are consistent with a higher collision rate as the density increases from the hot to the cold layer. The system with z ) 0.05 shows a less pronounced slope in the density gradient (see Figure 4b); therefore, the changes in the collision flux are not so large. In this respect the fluxes seem insensitive to the effect of the ionic charge and depend only on the temperature and density. The contribution JPOT. is relatively small for the LJ/s system, but when it is compared to results at higher temperature at the same density,18 we find that it increases with decreasing temperature. This is compatible with results reported by Wang et al.14 in their study of the thermal conductivity of carbon dioxide, where they observed that the potential contribution is the dominant one near the critical point. At pressure P*1 (run 1b) it accounts for approximately 15% of the total flux as compared with 7% at P*2 and 5% at a reduced pressure of approximately 4 times P*1,18 all at the same average density. A large value in the potential contribution to the energy flux of the LJ/s system has also been observed at a vapor/liquid interface.17 A tentative conclusion is therefore that as the thermodynamic state approaches a phase boundary, the potential energy flux contribution increases relative to the other two contributions. If we carry this over to the Yukawa and Coulomb systems, we must expect that the potential energy flux decreases with increasing charge strength, because the system moves away from the critical region further into the supercritical region. This

1886 J. Phys. Chem., Vol. 100, No. 5, 1996

Bresme et al.

Figure 11. Total energy flux JU and contributions from light (J2) and heavy (J1) species. Black circles correspond to the LJ/s system, whereas white circles represent the Coulomb system at z ) 0.05.

Figure 10. Collision and kinetic contributions from species 1 and 2 to the total energy flux. The symbols have the same meaning as in Figure 4.

is confirmed by Figure 9, which shows that the potential contribution of a system at z ) 0.05 is clearly smaller than that of the neutral system. More details of the transport mechanisms can be obtained by studying the contribution of each species to the fluxes discussed above. Parts a and b of Figure 10 show the kinetic and collision fluxes for each component, i.e., JK,1 and JK,2 such that JK,1 + JK,2 ) JKIN., and similarly for the collision contribution. For z ) 0 and z ) 0.01, the lighter component shows the largest variation in the kinetic contribution along the cell. Besides, the kinetic contributions from light and heavy species converge as one moves from the hot to the cold layer. The picture is somewhat complementary for the collision contributions, but they do not converge from the cold to the hot layer. The changes in JC,1 and JC,2 along the cell are of the same order of magnitude (cf. Figure 10b), and at z ) 0.05 they are essentially the same for the light and heavy components, both in the kinetic and collision contributions. Indeed, for the latter case the ratios JC,2/JC,1 and JK,2/JK,1 are rather constant along the cell. We believe this is due to the absence of the Soret effect in this system, which implies a uniform concentration of both species along the MD cell. The results above agree with previous work,18 where it was concluded that the kinetic energy flux of the lighter species and collision flux of the heavier species play the main role in the mechanism of thermal diffusion. In addition, we find that the collision contribution of the lighter species plays an important role. Figure 11 shows the contributions to the total energy flux from each species for z ) 0 and the Coulomb system with z ) 0.05. The error bars represent 1 standard deviation of the mean, based on a statistical analysis of a time series of the appropriate subaverages, to give an idea of the uncertainties in our results. The difference between the local and imposed energy fluxes (1.0 in this graph) in the x-direction is within the statistical uncertainties in the averages, which is a confirmation of the proper functioning of the algorithm. Conversely to an observation in a previous work,18 there is no sign of convergence of the contributions JU,1 and JU,2 from the hot to the cold layer; indeed, the fluxes are rather constant (within the uncertainties). A similar representation of the contributions from each species to the total energy flux is given in Figure 12. In this case, the Coulomb and the Yukawa systems are compared. Again, JU,1 and JU,2 are found to be constant along the cell, and it is

Figure 12. Influence of potential range on the energy flux. The results correspond to charge z ) 1.0. Black circles refer to the Coulomb system and open circles to the Yukawa system.

remarkable that the contributions are insensitive to the model potential. Hence, it seems that the fluxes depend on changes in density and temperature and, in the case of the potential contribution, also on the proximity to a phase boundary, but these contributions considered as a whole appear to be largely independent on the potential. 5. Conclusions In this work we have analyzed the influence of the potential range and strength on the transport properties of a model consisting of two species and that varies continuously from a Lennard-Jones/spline potential to an ionic potential. The main result of our study is that the Soret effect diminishes as the ionic charge is augmented. The thermal diffusion factor is generally positive, indicating that the lighter species migrates to the hot region of the system. As a function of the ionic strength, the greatest changes in the thermal diffusion factor are observed at low charges. Two different potential models, Yukawa and Coulomb potentials, were used in searching for a possible dependence of transport properties on the range of the potential. Our study reveals that a Coulomb potential inhibits the Soret effect, except for very small charges, whereas the Yukawa potential gives a positive thermal diffusion factor even at the largest charge strength studied in this work (z ) 1.0). The decrease of R12 with ionic charge seems related to the local interactions, because the attraction between unlike particles favors a uniform distribution of the two species. The differences between the Yukawa and the Coulomb systems may be explained by the different behavior of the long-range fluctuations in the charge density. This is expressed by the Stillinger-Lovett sum rules, i.e., the second moment condition and the zeroth (or local electroneutrality) condition. In the context of the Gibbs

Heat Conduction in Ionic Systems phase rule, this implies that the Yukawa system is a twocomponent mixture (except when κ f 0), whereas the Coulomb system is a one-component system (except when z f 0). The question of exactly how and where the transition between a twoand a one-component system occurs, and the somewhat intricate implications this transition has on the Gibbs phase rule and the phase behavior in general, has not been analyzed here. The application of short-ranged potentials to describe molten salts or Coulomb systems is not justified when dealing with properties such as the thermal diffusion factor. Hence, this property is relatively sensitive to the peculiarities of the molecular interactions, as has been noted by other authors.35 Conversely, the thermal conductivity has essentially the same results for the two ionic potentials mentioned above. These conclusions are probably valid for different mass ratios given that the thermal conductivity is quite insensitive to the mass ratio.18 It is not clear to what extent our conclusions will remain valid at lower densities. We have found that, at F* ) 0.4, the main mechanism of the transport of heat is the collision among particles and that a local structure is responsible for the observed thermal conductivity, which explains why the Yukawa and Coulomb systems give the same results. The thermal conductivity seems to follow a simple law of corresponding states. A further analysis of the potential models employed in this report would be very interesting but has not been pursued here. The kinetic energy flux was found to decrease from the hot to the cold layer in the system, whereas the collision energy flux increases in the same direction. These two contributions complement each other along the cell and are by far the dominant contributions to the total energy flux. The potential energy flux is much smaller than the other two, but we found it to increase in the neighborhood of the liquid/vapor phase boundary, in agreement with previous observations. The most significant changes in the kinetic and collision contributions from each species from the hot to the cold layer are a decrease of the kinetic flux for the lighter component and a similar increase of collision flux for both components. In addition, neutral and weakly charged systems are characterized by a decrease of the ratio JK,2/JK,1 from the hot to the cold layer, whereas JC,2/JC,1 increases in the same direction. The absence of the Soret effect, caused by an increase in the charge, tends to make these ratios constant along the cell. The total energy flux contribution from each species is insensitive to the range of the potential, and its behavior depends mainly on the thermal and density gradients and the mass ratio. Finally, we have reported preliminary computations of the equation of state of the LJ/s system using NEMD. Aside from being a further probe that the hypothesis of local equilibrium is justified, the NEMD technique appears to be a powerful tool in obtaining not only transport properties but also the equation of state. Acknowledgment. F.B. acknowledges The Norwegian Research Council for Sciences and the Humanities (NAVF, the cultural exchange program between Spain and Norway) and the Light Metals Centre of Excellence at The Norwegian Institute of Technology for support during his stays at the Department of Physical Chemistry of the University of TrondheimsThe Norwegian Institute of Technology in 1992/3 and 1994, respectively. I.W. was supported by NAVF under Grant No. 100301/431. B.H. thanks The Agency of Industrial Science and Technology (AIST) for financial support to stay at The National Institute for Advanced Interdisciplinary Research, Japan. Some of the computer resources used in this work were provided by TRU at the University of Oslo. Stimulating discussions with

J. Phys. Chem., Vol. 100, No. 5, 1996 1887 Professors John M. Kincaid and Signe Kjelstrup Ratkje are highly appreciated. Appendix: Ewald Sum for the Collision Energy Flux Contributions The Coulomb potential between two ions is

zizje2 zizjc φij ) ) 4π0rij rij

(A1)

where zi is the valence of ion i, e is the electronic charge, 0 is the permittivity of vacuum,  is the dielectric constant of the bj - b ri|, and c ) medium in which ions are immersed, rij ) |r e2/4π0. The total electrostatic energy of a periodic Coulomb system surrounded by a perfect conductor is given by27

φij )

c

N

N

zizj

∑∑′ 2 i)1 j)1 r

erfc(Rrij) +

( )

ij

2cπ

∑

′

1

LxLyLz |hB|