Dissertation StefanieHerzog

Modeling and Molecular Simulation of Electrolyte Systems

Modellierung und Molekularsimulation von Elektrolytsystemen

Der Technischen Fakultt der

Universitt Erlangen-Nrnberg

zur Erlangung des Grades

DOKTOR-INGENIEUR

vorgelegt von

Stefanie Elke Herzog

Erlangen 2008

Als Dissertation genehmigt

von der Technischen Fakultt der

Universitt Erlangen-Nrnberg

Tag der Einreichung: 09.05.2008

Tag der Promotion: 26.09.2008

Dekan: Prof. Dr.-Ing. habil. Johannes Huber

Berichterstatter: Prof. Dr.-Ing. Wolfgang Arlt

Prof. Dr.-Ing. Joachim Gro

III

Danksagung Diese Arbeit entstand whrend meiner Ttigkeit als wissenschaftliche Mitarbeiterin an der Technischen Universitt Berlin von Oktober 2003 bis Mrz 2004 und an der Friedrich-Alexander-Universitt Erlangen-Nrnberg von April 2004 bis Juni 2007 unter Leitung von Herrn Prof. Wolfgang Arlt. Ich danke meinem Doktorvater Herrn Prof. Wolfgang Arlt fr die Mglichkeit, eigenverantwortlich in die Tiefen der molekularen Wechselwirkungen zwischen Ionen und dipolaren Moleklen einzusteigen, sowie fr die Mglichkeit, Verant-wortung in der Lehre am Lehrstuhl zu bernehmen. Meine Zeit am Lehrstuhl hat mich auch persnlich reifen lassen. Auch dafr danke ich meinem Doktorvater ganz herzlich.

Herrn Prof. Joachim Gro von der TU Delft danke ich nicht nur fr die bernahme des zweiten Gutachtens, sondern insbesondere fr seine seit Beginn meiner Diplomarbeit kontinuierliche Bereitschaft zu wissenschaftlichen Diskussionen ber Fluidtheorien, Mole-kularsimulationen und vieles mehr. Darber hinaus danke ich Herrn Prof. Joachim Gro fr sein Engagement bei meinem Einstieg in die Promotion, die Bereitschaft zu langen Diskussionen am Telefon, die Mglichkeiten zu Aufenthalten in Delft, die Vermittlung wissenschaftlicher Kontakte auf Tagungen sowie die Beratungen zu meinem Berufseinstieg.

Bei Herrn Prof. Rainer Buchholz vom Lehrstuhl fr Bioverfahrenstechnik und Herrn Prof. Donat-Peter Hder vom Lehrstuhl fr kophysiologie der Pflanzen bedanke ich mich fr die Beteiligung als Prfungsvorsitzender bzw. als weiteres prfungsberechtigtes Mitglied am Prfungsausschuss.

Herrn Prof. Richard Sadus von der Swinburne University of Technology, Melbourne danke ich fr wissenschaftliche Diskussionen bei Besuchen am Lehrstuhl sowie auf Tagungen. Einen ganz herzlichen Dank spreche ich auch an Herrn Prof. Athanassios Panagiotopoulos und Herrn Dr. Philip Lenart von der Princeton University fr die Einfhrung in die grokanonischen Monte Carlo-Simulationen mit histogram reweighting und mixed-field finite-size scaling und die berlassung des entsprechenden Fortran-Codes aus. Ich danke meinen Studenten sowie ehemaligen Berliner und Erlanger Kollegen. Ganz besonders mchte ich mich bei Frau Prof. Irina Smirnova und bei Frau Dr. Steffi Hiller bedanken. Die beiden haben mir schon whrend meiner Zeit als Thermo-Tutorin in Berlin gezeigt, wie man am Lehrstuhl lehrt, forscht und Spa hat. Vielen Dank dafr! Ganz herzlichen Dank auch an meine ehemaligen Brokollegen Supakij Suttiruengwong und Jrn Rolker. Durch ihre ausgeglichene Art habe ich Stress als weniger belastend empfunden. Vielen Dank an Oliver Spuhl, Jrn Rolker und Dirk-Uwe Astrath fr die Leidensgenossenschaft beim Umzug nach Erlangen. Dirk-Uwe, Dir auch ganz lieben Dank fr wissenschaftliche und andere Diskussionen sowie fr den Draht zur Brcke. Alexander, Elisabeth, Edel und Petra danke ich ganz besonders fr die liebe Untersttzung am Tag meiner Promotion und natrlich auch schon davor.

IV

Tobias, Andrea, Thomas, Jozo, Marta, Hella, Sebo, Irina, Christa, Bernd, Opa, Alexandra und Peter, ich habe mich sehr gefreut, dass Ihr bei meiner Prfung wart und mit mir gefeiert habt. Elmar, lieben Dank fr den Platz auf Deinem Sofa und die Untersttzung beim Vortragsfeinschliff.

Ein ganz besonderer Dank gebhrt meinen Eltern, Monika und Wolfgang Herzog. Ihr seid eine groe Sttze fr mich.

Dir lieber Uwe danke ich fr Deine Liebe, Deine Untersttzung in allen Fragen und Lagen sowie ungezhlte Fahrten an meine wechselnden Wohnorte. Durch Dich erfahre ich unerschpflichen Rckhalt. Dafr danke ich Dir!

V

Table of Contents

DANKSAGUNG......................................................................................................................III

NOMENCLATURE ............................................................................................................... VII

KURZFASSUNG..................................................................................................................XIII

ABSTRACT ........................................................................................................................ XIV

1 INTRODUCTION AND OBJECTIVES.................................................................................. 1

1.1 INTRODUCTION ............................................................................................................. 1 1.2 OBJECTIVES .................................................................................................................. 1

2 THERMODYNAMIC BASIS ............................................................................................... 3

2.1 THERMAL EQUATIONS OF STATE AND CALCULATION OF THERMODYNAMIC QUANTITIES.................................................................................................................. 3

2.2 PHASE EQUILIBRIUM CONDITIONS AND SOLUTION ....................................................... 4 2.3 THERMODYNAMICS OF ELECTROLYTE SYSTEMS .......................................................... 6

2.3.1 Electrolytes in a Solvent .................................................................................... 6 2.3.2 Thermodynamic Properties of Electrolyte Systems ........................................... 7

3 STATISTICAL MECHANICS ........................................................................................... 11

3.1 INTRODUCTION ........................................................................................................... 11 3.1.1 A Statistical Mechanical System...................................................................... 11 3.1.2 Statistical Mechanical Ensembles ................................................................... 13

3.2 STATISTICAL MECHANICS FOR DENSE SYSTEMS ........................................................ 14 3.2.1 Distribution Functions..................................................................................... 15 3.2.2 Integral Equation Theories.............................................................................. 16 3.2.3 Microthermodynamics ..................................................................................... 17 3.2.4 Perturbation Theories ..................................................................................... 19

3.3 THEORIES FOR THE DESCRIPTION OF ELECTROSTATIC INTERACTIONS ........................ 20 3.3.1 Electrostatic Pair Potentials ........................................................................... 20 3.3.2 Historical Theories .......................................................................................... 22 3.3.3 Perturbation Theories ..................................................................................... 23 3.3.4 Integral Equation Theories.............................................................................. 24

3.4 MOLECULAR SIMULATIONS ........................................................................................ 26 3.4.1 Introduction ..................................................................................................... 27 3.4.2 The Concept of Monte Carlo Simulations ....................................................... 29 3.4.3 Monte Carlo Simulations in the Grand Canonical Ensemble ......................... 32 3.4.4 Histogram Reweighting for a Pure Component System .................................. 35 3.4.5 Histogram Reweighting for Binary Mixtures .................................................. 37

VI

3.4.6 Determination of the Critical Point from Finite-Size GCMC Simulations ......38 3.4.7 Ewald Summation for Long-Range Interactions..............................................40

3.5 LITERATURE OVERVIEW ON COMPUTER SIMULATION APPROACHES FOR SIMPLE MODEL ELECTROLYTE SYSTEMS.................................................................................45

3.6 LITERATURE OVERVIEW ON THERMODYNAMIC MODELS FOR ELECTROLYTE SYSTEMS .....................................................................................................................47

4 MODELING OF AQUEOUS ELECTROLYTE SYSTEMS USING THE NP MSA....................51

4.1 MODEL AND PROCEDURE ............................................................................................51 4.1.1 PC-SAFT EOS for Non-Electrostatic Interactions...........................................52 4.1.2 Semirestricted Non-Primitive Mean Spherical Approximation for Electrostatic Interactions ..............................................................................................52 4.1.3 Calculation of Phase Equilibria for Electrolyte Systems.................................57 4.1.4 Parameters for the np MSA + PC-SAFT EOS .................................................58

4.2 RESULTS AND DISCUSSION ..........................................................................................59 4.2.1 Pure Component Parameters for Water...........................................................59 4.2.2 Application of the np MSA + PC-SAFT EOS to Aqueous Electrolyte Systems at Infinite Dilution............................................................................................60 4.2.3 Application of the np MSA + PC-SAFT EOS to Aqueous Electrolyte Systems at Finite Salt Concentration.............................................................................63 4.2.4 Limitations of the semirestricted np MSA + PC-SAFT EOS............................73

5 MONTE CARLO SIMULATIONS FOR MODEL ELECTROLYTE SYSTEMS .......................75

5.1 MODEL FLUID AND METHOD.......................................................................................75 5.1.1 Investigated Fluids ...........................................................................................75 5.1.2 Simulation Method ...........................................................................................76

5.2 RESULTS AND DISCUSSION ..........................................................................................80 5.2.1 Pure Stockmayer Fluid.....................................................................................80 5.2.2 Pure Lennard-Jones Salt ..................................................................................84 5.2.3 Stockmayer/Lennard-Jones Mixture ................................................................88 5.2.4 Mixture of a Stockmayer and a Charged Lennard-Jones Fluid.......................89

6 SUMMARY AND PERSPECTIVES .....................................................................................95

APPENDIX A REFERENCE STATE FOR ION ACTIVITY COEFFICIENTS.........................99

APPENDIX B REDUCED QUANTITIES FOR MOLECULAR SIMULATIONS ....................102

APPENDIX C EWALD SUMMATION FOR IDEAL DIPOLES ...........................................103

APPENDIX D MONTE CARLO SIMULATION RESULTS ................................................104

BIBLIOGRAPHY ..................................................................................................................125

VII

Nomenclature Capital Letters

A [J] Helmholtz energy

Ci [-] weight, iterated quantity in histogram reweighting, connected with the partition function of the grand canonical ensemble

F [-] number of degrees of freedom FH [-] quantity for determination of starting values for the solution of

the semirestricted np MSA

G [J] Gibbs energy

Gis [J] Gibbs energy of solvation H [J] enthalpy

K [-] number of components

K [-] chemical equilibrium constant for the dissociation of a salt in a solvent

K [-] total number of observations in a grand canonical Monte Carlo simulation, K = NUf(N,U)

L* [-] reduced simulation box length, L* = L/ M [kg kmol-1] molar mass

M [-] ordering parameter for mixed-field finite-size scaling

Ni [-] number of molecules/particles Q(N,V,T) [-] partition function for the canonical ensemble R [J mol-1 K-1] ideal gas constant, R = 8.314 Jmol-1K-1

S [J K-1] entropy T [K] temperature

T* [-] reduced temperature in molecular simulations, T* = kT/

U [J] potential energy U [J] total energy or Hamiltonian of a system

V [m] volume

Z [-] compressibility factor, Z = pv/kT

VIII

Small Letters

a1, a2 [-] auxiliary parameters for the solution of the semirestricted np MSA

b0, b1, b2 [-] energy parameters for the semirestricted np MSA (solution of the semirestricted np MSA)

c(r) [-] direct correlation function ci [mol salt

l solvent-1] salt molarity

d02, d22 [-] charge density of the ions, and dipole density of the solvent in the semirestricted np MSA (These two parameters describe the properties of the hard-sphere ion-dipole system.)

e [esu] electron/elementary charge fi [Pa] fugacity of i f(N,U) [-] number of observations of a certain number of particles N

associated with a certain potential energy U in a grand canonical Monte Carlo simulation

g2(r), g(r) [-] radial pair distribution function h(r) [-] total correlation function k [J K-1] Boltzmann constant, k = 1.3806610-23 J K-1

kmax [-] Ewald summation parameter k10, k11 [-] auxiliary parameters for the solution of the semirestricted np

MSA

mi [mol salt kg solvent-1]

salt molality

ni [mol] number of moles

ne,ion [-] number of electrons of an ion

n [-] description of the arrangement of periodic boxes around the central one with ( )0,0,0=n in Ewald summation

p [Pa] pressure

ip [N s] linear momentum of particle i

IX

p() [-] probability for the density p(N,U;,) [-] probability distribution function in the grand canonical

ensemble

qi [esu] charge, qi = zie

q* [-] reduced charge for the Lennard-Jones 1-1 salt in Monte Carlo

simulations,

qq =*

r [] radial distance between molecules or molecular sites s [-] field mixing parameter for mixed-field finite-size scaling

u [J] pair potential

v [m mol-1] molar volume xi, yi [-] mole fraction of component i

y1 [-] auxiliary parameters for the solution of the semirestricted np MSA

zi [-] valence/charge number of ions

Greek Letters

ion [3] polarizability of an ion [1/J] inverse temperature, = 1/kT

i [-] activity coefficient of component i

* (m) [-] molality based mean ionic activity coefficient [-1] screening parameter in the primitive mean spherical

approximation

(N,p,T) [-] partition function of the isobaric-isothermal ensemble [J] Lennard-Jones energy parameter [J] dispersion parameter for the PC-SAFT EOS

[-] dielectric constant/relative static permittivity

0 [C2 N-1 m-2] permittivity of free space 0 = 8.8541910-12 C2 N-1 m-2 AiBi [J] pure component parameter for the PC-SAFT EOS (potential

depth for association)

X

AiBj association energy between cations and water to describe the additional solvation in the semirestricted np MSA + PC-SAFT EOS

[-] Ewald summation parameter (screening width for the point charges and point dipole moments)

AiBi [-] pure component parameter for the PC-SAFT EOS (effective association volume)

i [J] chemical potential of component i

D,i [D] dipole moment of component i D* [-] reduced dipole moment in Monte Carlo simulations,

3

* DD =

i [-] number of moles of ion i from 1 mol of salt

[-] number of moles of anion and cation from 1 mol of salt,

= + + -

[] correlation length (goes to infinity at the critical point)

3 [-] packing fraction for spherical molecules, = 33 6 iix pi

(,V,T) [-] grand canonical partition function [kg m-3] density

[-3] number density, = N/V e [C m-3] electrical charge density (Debye-Hckel theory) [] molecular/ionic/hard-sphere diameter [] Lennard-Jones size parameter [] size parameter in the PC-SAFT equation of state standard deviation

ion [] ion diameter for the semirestricted np MSA term (adjustable parameter in the semirestricted np MSA + PC-SAFT EOS)

i [-] fugacity coefficient of component i [J C-1] electrical potential

(m)

[-] molality based osmotic coefficient

XI

i [-] Euler angle between molecular sites

i [-] perturbation term of order i (N,V,U) [-] microcanonical partition function

Indices superscript

assoc association

c based on compressibility equation, see equation (3-19) CC charge-charge

CD charge-dipole

DD dipole-dipole

disp dispersion

hc hard-chain

id ideal gas

LJ Lennard-Jones

m molality based

(np) MSA refers to the (semirestricted non-primitive) mean spherical approximation

res residual

v based on virial equation, see equation (3-18) ' phase 1

'' phase 2

phase

Indices subscript

c critical

(c) molarity based i component i

(m) molality based new state before a Monte Carlo move (position/energy)

XII

old state after a Monte Carlo move (position/energy) real real space term in Ewald summation

reciprocal reciprocal space (Fourier) term in Ewald summation self self space term in Ewald summation SM Stockmayer

soft contrary to repulsive solv solvent

surf surface term in Ewald summation (not necessary if tinfoil boundary conditions are used)

(x) rational (mole fraction based) + positively charged (ion) - negatively charged (ion)

Abbreviations

AAD average absolute deviation for calculated quantities i,calc respective to the corresponding experimental quantities i,exp,

1001%1 exp,

,exp,

= =

N

i i

calcii

NAAD

CGS Gaussian system of units (centimeter, gram, second) EOS equation of state

GCMC grand canonical Monte Carlo simulation

MKS International System of Units (SI), MKS denote meter, kilogram and second

np MSA non-primitive mean spherical approximation

NRTL non-random two liquids (activity coefficient model) p MSA primitive mean spherical approximation

PC-SAFT Perturbed-Chain Statistical Association Fluid Theory

SAFT Statistical Association Fluid Theory

SAFT-VR(E) (electrolyte version) of the statistical associating fluid theory for fluids with attractive potentials of variable range

XIII

Kurzfassung Wssrige Elektrolytsysteme sind fr viele industrielle Prozesse von Bedeutung. Das Problem an den bisher fr diese Systeme verwendeten Aktivittskoeffizientenmodellen und Zustandsgleichungen ist, dass sie die elektrostatischen Wechselwirkungen nur ungengend beschreiben. Der Hauptgrund dafr ist die implizite Bercksichtigung des Lsungsmittels durch seine Dielektrizittskonstante. Daher widmet sich diese Arbeit der fluid-theoreti-schen Modellierung und molekularen Simulation von Elektrolytsystemen unter expliziter Bercksichtigung des Lsungsmittels als molekulare Komponente.

Der Modellierungsansatz kombiniert die perturbed-chain statistical associating fluid theory Zustandsgleichung (PC-SAFT-Zustandsgleichung) mit der semirestricted non-primitive mean spherical approximation (np MSA). Fr Wasser werden die Reinstoff-parameter der semirestricted np MSA + PC-SAFT-Zustandsgleichung an Flssigdichten und Dampfdruckdaten angepasst. Die Alkalihalogenid-Salze werden mithilfe ionenspezi-fischer Parameter aus der Literatur (Kristalldurchmesser von Pauling und Dispersions-parameter aus einer Korrelation von Mavroyannis et al.) beschrieben. Zur Validierung des Modells und zur Bestimmung der anzupassenden Parameter wird die semirestricted np MSA + PC-SAFT-Zustandsgleichung auf unendlich verdnnte wssrige Elektrolytsysteme angewendet. Hierbei zeigt sich, dass der salz-spezifische Ionendurchmesser fr den semirestricted np MSA-Term ion und zustzlich fr 11 der 19 untersuchten Systeme die Kationen-Wasser-Assoziationsenergie AiBj/k angepasst werden mssen. Diese Gren werden an mittlere Ionenaktivittskoeffizienten und osmotische Koeffizienten bei T = 25C ber den gesamten Lslichkeitsbereich des Salzes angepasst. Das Modell kann zur Beschreibung von Systemdrcken und Flssigdichten bis zur Sttigungsgrenze ohne temperaturabhngige Parameter bis 100C extrapoliert werden.

Der Molekularsimulationsansatz nutzt die Monte Carlo-Methode im grokanonischen Ensemble (GCMC) mit histogram reweighting und mixed-field finite-size scaling, um Phasengleichgewichtsdaten fr Modell-Elektrolytsysteme zu erzeugen. Das Modell-elektrolytsystem besteht aus Stockmayer-Moleklen, die das dipolare Lsungsmittel reprsentieren, sowie aus Lennard-Jones-1-1-Salzen, die aus zwei Wechselwirkungszentren gleicher absoluter Ladung bestehen, welche das Anion und das Kation verkrpern. Die Lennard-Jones-Parameter und sind fr alle Molekle und Ionen gleich. Die Simulationsmethode wird auf reine Lennard-Jones-1-1-Salze mit reduzierten Ladungen von |q*| = 2,5 5,0 und 7,0 angewendet. Phasengleichgewichte werden ausgehend vom kritischen Punkt ermittelt. Fr das Stockmayer/Lennard-Jones-1-1-Salz-System mit einem reduzierten Dipolmoment fr das Stockmayer-Fluid von D* = 2,0 werden Phasengleich-gewichte fr T* = 1,6 und q* = 5,0 sowie T* = 1,8 und q* = 5,0 und 7,0 bestimmt. Die Simulationsdaten erffnen die Mglichkeit eine verbesserte statistisch-mechanische Theorie fr Ionen-Dipol-Wechselwirkungen zu entwickeln.

XIV

Abstract Aqueous electrolyte systems appear in many industrial processes. Activity coefficient models and equations of state for these systems suffer from an insufficient description of the electrostatic interactions. The main deficiency is that the solvent is regarded implicitly by its dielectric constant. Therefore, this work is devoted to fluid-theoretical modeling and molecular simulation of electrolyte systems by accounting for the solvent explicitly, as a molecular species.

The modeling approach combines the perturbed chain statistical associating fluid theory equation of state (PC-SAFT EOS) with the semirestricted non-primitive mean spherical approximation (np MSA). Parameters for water for the semirestricted np MSA + PC-SAFT EOS are adjusted to liquid density and vapor pressure data. The alkali halide salts are described using ion-specific parameters (Pauling crystal diameters and dispersion parameters from a correlation of Mavroyannis et al.). The semirestricted np MSA + PC-SAFT EOS is applied to aqueous electrolyte systems at infinite dilution of the salts, in order to validate the model and identify the parameters which need to be adjusted. Thereby it is found that the salt-specific ion diameter for the semirestrictted np MSA term ion and additionally for 11 out of the 19 systems the cation-water solvation energy AiBj/k have to be adjusted. These parameters are adjusted to mean ionic activity coefficient data and to osmotic coefficient data at T = 25C over the full solubility range of the salt. In order to model system pressures and liquid densities up to the solubility limit of the salt, the model can be extrapolated up to 100C without any temperature-dependent parameters.

The molecular simulation approach uses Monte Carlo simulations in the grand canonical ensemble (GCMC) together with histogram reweighting and mixed-field finite-size scaling in order to generate phase equilibrium data for a model electrolyte system. The model electrolyte system consists of Stockmayer molecules, representing the dipolar solvent, and Lennard-Jones 1-1 salt molecules consisting of two sites representing the anion and the cation, where both sites have the same absolute charge. The Lennard-Jones parameters and are equal for all molecules and ions in all systems. The simulation method is applied to pure Lennard-Jones salts with reduced charges of |q*| = 2.5, 5.0 and 7.0. Phase equilibria are determined starting from the critical point. For the Stockmayer/Lennard-Jones 1-1 salt system with a reduced dipole moment for the Stockmayer fluid of D* = 2.0, phase coexistence properties for T* = 1.6 and q* = 5.0 and T* = 1.8 and q* = 5.0 and 7.0 are determined. The simulation data offer the possibility to develop an improved statistical mechanical theory for ion-dipole interactions.

1

1 Introduction and Objectives 1.1 Introduction

Electrolyte systems are ubiquitous in nature and important for many industrial processes. One example for an industrial process, which involves electrolytes, is sour gas treatment. The chemical absorption of sour gases in aqueous amine solutions leads to the formation of ions, which can significantly change the thermodynamic properties of the system. In principle, for all processes which involve organic and inorganic salts as well as acids, which dissociate and form ions, the thermodynamic model has to account for the electrostatic interactions present between the ions and the solvent molecules. The parameters for the thermodynamic model need to describe the system for the range of concentration of all components, pressure, and temperature, which are relevant for operation.

Industrially used thermodynamic models for electrolyte systems, e.g. the Pitzer equation [1,2] and the electrolyte NRTL model [3] still rest on the pioneering work of Debye and Hckel [4] from the twenties of the last century. This theory describes the long-range electrostatic interactions between ions in an electrolyte system. The solvent is implicitly represented by its dielectric constant, and the ions are considered to be point charges. On account of the assumptions of Debye and Hckel this model can describe the mean ionic activity coefficient only up to a salt concentration of 0.01 mol salt/kg solvent [5]. Therefore, the significance of the Debye-Hckel theory is rather the correct limiting behavior at infinite dilution than an accurate description of electrolyte systems at finite salt concentration. The inherent problem of models which consider the solvent (mixture) implicitly by its dielectric constant is that effects like ion solvation cannot be described. Furthermore, the dielectric constant is a bulk property which depends on the density, temperature and composition of the solvent system, what is not accounted for in the models.

To describe all electrostatic interactions, which are present in electrolyte systems, the solvent(s) should be explicitly considered as a molecular species. An electrolyte theory which considers the solvent explicitly is referred to as a model on the Born-Oppenheimer level. Statistical mechanics offers concepts for the description of electrostatic interactions between ions and solvent molecules, namely perturbation theories and integral equation theories. These statistical mechanical theories need to be validated with molecular simulation data. Unfortunately, few molecular simulation data for model electrolyte systems on the Born-Oppenheimer level are available. Therefore, these statistical mechanical theories are seldom validated and consequently not used in engineering thermodynamics models, and rarely further developed.

1.2 Objectives This work focuses on introducing statistical mechanical concepts for electrostatic interactions on the Born-Oppenheimer level into engineering equations of state in order to develop a model which describes electrolyte systems. This work restricts to aqueous electrolyte systems

2 1 Introduction and Objectives

because a vast resource of experimental data is available. The extension of a thermodynamic model in order to account for interactions due to point charges seems to be straightforward, but as Prausnitz et al. mention in [5], p. 508: Extension of a thermodynamic framework for nonelectrolytes to contain also electrolytes is not a trivial task. It is a common misconception that such extension is only a small detail, a little perturbation, like adding a short tail to a big dog. Not so. Extension to include also electrolytes requires concepts and constraints (e.g. electroneutrality) that can be mastered only with patient and devoted study. Besides these conceptual difficulties for extending thermodynamic models to account for electrolyte systems, the electrostatic interactions between ions and also between solvent molecules and ions are long-ranged and strong at short range [6]. This causes problems in molecular simulations [6], and slow convergence of the terms in electrolyte perturbation theories [7].

In principle, one can think of two different approaches for introducing an electrostatic term into an engineering equation of state. The first is to use a statistical mechanical theory for electrostatic interactions from literature together with an engineering equation of state for nonelectrolytes for modeling of electrolyte systems. If there is no suitable statistical mecha-nical model for the description of electrostatic interactions between ions and dipolar solvent molecules, such a model needs to be developed based on molecular simulations. For this purpose molecular simulations need to be carried out. This work covers modeling and molecular simulation of electrolyte systems.

For the first approach the semirestricted non-primitive mean spherical approximation (semirestricted np MSA) [8,9] is used to extend the perturbed-chain SAFT equation of state (PC-SAFT EOS) [10-12] in order to model aqueous alkali halide systems over the full solubility range. The semirestricted np MSA is an integral equation theory which describes excess thermodynamic quantities over a hard-sphere reference system. It considers ion-ion, ion-dipole and dipole-dipole interactions, assuming that these are the main interactions between charged ions and dipolar solvent molecules. This work focuses on validating the model at infinite and finite concentration of the salt, on the adjustable salt-specific parameters, and investigates the temperature extrapolation capability of the model.

Due to limitations of the semirestricted np MSA in the description of aqueous electrolyte systems stemming from the models closures and the neglect of the orientation dependence between the molecules, molecular simulations for model electrolyte systems are conducted. The model electrolyte comprises soft repulsion and dispersive interactions between ions and solvent molecules, as well as electrostatic interactions between the point charges of the ions and the point dipole moments of the solvent. This work uses the Monte Carlo method in the grand canonical ensemble together with histogram reweighting. The simulation results offer the possibility to develop a statistical mechanical theory for ion-ion, ion-dipole and dipole-dipole interactions based on a reference fluid which is more realistic compared to the hard-sphere reference fluid which was used for developing the np MSA.

3

2 Thermodynamic Basis This chapter provides a brief overview of essential concepts of thermodynamics. In section 2.1, the idea of thermal equations of state is introduced. Furthermore, it is shown how thermo-dynamic quantities can be calculated from the Helmholtz energy. In section 2.2, the phase equilibrium problem and a few aspects of its solution are elucidated. In section 2.3, some important aspects of the thermodynamics of electrolyte systems are discussed. This includes the dissolution of a salt in a solvent and the definition of thermodynamic quantities which are peculiar for electrolyte systems.

2.1 Thermal Equations of State and Calculation of Thermodynamic Quan-tities

Thermal equations of state relate the pressure p, the temperature T, the molar volume v and the mole fractions xi of a fluid. In the following sections the word thermal is left out for bre-vity. It is advantageous to use equations of state (EOS) for phase equilibrium calculations as opposed to Gibbs excess energy models because they can be used over wide ranges of temperature and pressure and can be applied to components ranging from light gases to heavy liquids. EOS are frequently used to calculate all kinds of fluid-phase equilibria without conceptual difficulties [13,14]. For phase equilibrium calculations only the residual Helmholtz energy1 Ares is needed which is the total Helmholtz energy minus the ideal gas contribution. The ideal gas contribution is given by the ideal gas equation of state. When working with residual quantities, one needs to pay attention to the reference state, which means, whether the ideal gas boundary condition v or p0 is met. The conversion for all thermodynamic potentials between the two ideal gas reference states is explained by Gross [10]. The equations of state which are considered in this work are based on statistical mechanics. Because the ideal gas part is known for the compressibility factor Z, and is not needed for the calculation of the fugacity coefficient i, an EOS for phase equilibrium calculations only needs to deliver the residual part of the Helmholtz energy Ares.

( ) ( ) ( )NkT

xvTANkT

xvTANkT

xvTA iid

iires

,,,,,,

= (2-1)

Developing an equation of state in terms of the residual Helmholtz energy Ares is convenient because all thermodynamic quantities can be determined using partial derivatives with respect to the canonical variables of A. In this work, the partial derivatives of the Helmholtz energy A with respect to its canonical variables are briefly presented, in order to show, how an EOS in terms of A can be used, for phase equilibrium calculations. A detailed derivation of the partial derivatives of A can be found in the work of Gross [10].

4 2 Thermodynamic Basis

The compressibility factor Z is determined by partial derivation with respect to the number density . From the compressibility factor Z, the system pressure p can be calculated using the definition of the compressibility factor Z = pv/kT.

( )ixT

res

NkTA

Z

,

1

+=

(2-2)

The chemical potential i is defined by the partial derivative of the Helmholtz energy with respect to the number of moles of component i.

ijnvTii

n

A

=

,,

(2-3)

For phase equilibrium calculations usually the mole fractions xi instead of the number of moles ni or the number of molecules Ni are used. The differentiation of the Helmholtz energy can be carried out with respect to xi, under implicit inclusion of the summation relation [15].

( ) ( ) ( )

=

++=

N

kxvT

k

res

k

xvTi

resresres

i

kjij

x

NkTA

xx

NkTA

ZNkTA

kTvT

1,,,,

1, (2-4)

From the residual chemical potential ires the fugacity coefficient i can be calculated.

( )Z

kTvTresi

i ln,ln = (2-5)

For phase equilibrium calculations only the pressure p and the fugacity coefficients i for every component in each phase are needed. The equations for determination caloric properties from the Helmholtz energy are provided in the thesis of Gross [10].

2.2 Phase Equilibrium Conditions and Solution

A closed thermodynamic system is in equilibrium, if the extensive thermodynamic potential functions, which are the internal energy U(S,V,Ni), the enthalpy H(S,P,Ni), the Gibbs energy G(T,P,Ni) and the Helmholtz energy A(T,V,Ni) are in minimum. These minimum potential conditions can be replaced by more useful criteria in terms of the intensive quantities, which are the temperature T, the pressure p, and the chemical potential i [5]. These quantities need to be equal in all phases. This means that all phases are in the thermal and mechanical and chemical equilibrium

1 The Helmholtz energy A is one of the four thermodynamic potentials. Its canonical variables

are the volume V, the temperature T, and the numbers of particles Ni of component i.

2 Thermodynamic Basis 5

piTTT === ''' (2-6) pippp === ''' (2-7)

pi iii === ''' for i = 1K (2-8)

If a phase equilibrium problem is to be solved, F variables need to be specified. F is the num-ber of degrees of freedom and depends on the number of independent intensive variables, which is (K+1) and on the number of equilibrium relations given by equations (2-6) to (2-8), which is (-1)(K+2). The degree of freedom F is the difference of the number of independent intensive variables and the number of equilibrium conditions, and equals

2+= piKF (2-9) Equation (2-9) is known as Gibbs phase rule. The equality of the chemical potentials i can be replaced by the equality of the fugacities fi without loss of generality [5]. Equation (2-8) can be replaced by

piiii fff === ''' for i = 1K (2-10)

If a phase equilibrium problem is solved, the isofugacity conditions (2-10) as well as the mechanical (2-7) and thermal (2-6) equilibrium need to be solved. The fugacity fi can be expressed by the fugacity coefficient i by

pxf iii = (2-11)

If there are only two phases, with the mechanical equilibrium imposed, the isofugacity rela-tion can be written as

'''''' iiii xx = for i = 1K (2-12)

This set of equations needs to be solved iteratively for every component K. The fugacity coefficient i can thereby be determined from an equation of state, see equations (2-4) and (2-5). Besides the phase equilibrium condition in equation (2-12), the summation relation needs to be fulfilled in every phase

=

=

K

iix

11 pi (2-13)

For equations of state, which are not explicit in density, the density needs to be determined iteratively. For numerical robustness, the phase equilibrium iteration which is based on the solution of (2-12) and (2-13), is done in two loops, an outer isofugacity iteration under consid-eration of the summation relation and the inner density iteration. More details about phase equilibrium calculation using equations of state (including flow charts) can be found in Dohrn [16].


2.3 Thermodynamics of Electrolyte Systems

Electrolytes have some peculiarities that make their thermodynamic modeling different from that of classical substances in chemical engineering, e.g. hydrocarbons. These peculiarities result from the fact, that electrolytes dissolve partly or completely in a solvent and form indi-vidual ions, whereas the whole solution is electroneutral. Electrolytes have a dual nature. On the one hand, they form individual ions, which can be detected in an electrolyte solution. On the other hand, the macroscopic charge separation leads to such high counteracting forces, that effectively electroneutrality prevails. This at the same time implies that the effect of an individual ion species on thermophysical properties of the whole solution cannot trivially be singled out.

In this section, the solution behavior of electrolytes is explained and the concentration scales for electrolytes are introduced. In addition, electrolyte solutions are classified according to their solution behavior. Thereafter, thermodynamic properties for electrolyte solutions are introduced.

2.3.1 Electrolytes in a Solvent

If a salt is dissolved, it can dissociate into positively-charged cations and negatively-charged anions. The degree of dissociation depends on the salt and on the type of solvent, see for a solvent classification Barthel et al. [17], p. 3. If a salt completely dissociates, the solution contains only the solvent and the individual ions. However, in case of incomplete dissolution, the molecular form of the electrolyte is also present. This indicates that ion mole fractions xion for electrolyte solutions are ambiguous because the degree of dissociation must be known. Furthermore, even for complete dissolution the type of electrolyte determines how many moles of ions are formed from 1 mol of salt. Therefore, the dissociation behavior of a salt in a solvent needs to be studied first.

The dissociation reaction for a salt M+X- can be written as

+++ +

zzsolvent XMXM (2-14)

where the variables + and - stand for the number of moles of cation and anion, which are formed from 1 mol of salt respectively. Electrolytes are often classified according to these i, e.g. for a 1-1 electrolyte + = - = 1, which means that for complete dissociation of 1 mol of salt, 1 mol of cations and 1 mol of anions are formed. z+ and z- are the valence of the ions. Electrolyte systems are in all practical terms electroneutral, which means that no net charge exists, according to

0=+++ zz (2-15)

However, not all salts dissociate completely. A classification for the dissociation is provided by Robinson and Stokes [18], p. 49 et seqq. For a number of associated electrolytes, Robinson and Stokes use the term weak electrolyte to emphasize that these electrolytes can exist as undissociated or covalent molecules as well as ions. All acids belong to this group, because


they vaporize in their molecular form. For salts which do not completely dissociate, the chemical equilibrium constant K describes the ratio of undissociated electrolyte and dissoci-ated ions. In some salt/solvent systems, ion-pairing occurs as a result of the electrostatic inter-actions. Robinson and Stokes group almost all salts in non-aqueous solvents into this ion-pairing group.

In some salt/solvent systems the salt is believed to exist only in the form of cations and an-ions. Aqueous alkali halide systems, aqueous alkaline-earth halides, aqueous perchlorates, and some aqueous transition-metal halides and perchlorates belong to this group. These systems are of great importance for testing thermodynamic models, because the concentration of the ions in the electrolyte solution is known. Even for these systems, complete dissociation is only fulfilled up to a certain concentration limit, which is temperature dependent.

The concentration of electrolytes is usually given in molality mi as the number of moles of salt per kg of solvent. Another electrolyte concentration measure is molarity ci which is the number of moles of salt per volume of solvent. Molarity is used less often because the density of the solvent is needed for conversion to another concentration scale. The advantage of mo-lality mi and molarity ci is that these scales are independent of the degree of dissociation. Their disadvantage is that they go to infinity, if the salt mole fraction xsalt tends to 1.

For this work, the conversion between the mole fraction xi and the molality mi is used fre-quently, because in phase equilibrium calculation mole fractions xi are used but the experi-mental data is given in terms of molality mi. Assuming that the salt is fully dissociated the conversion is

iisolvi

i

mMkg

gx

1000

1

+

=

(2-16)

with Msolv denoting the molar mass of the solvent in g/mol, and is the total number of moles of ions from 1 mol of salt ( = i). 2.3.2 Thermodynamic Properties of Electrolyte Systems

Because individual ions cannot be isolated in a solution, their influence on the thermophysical properties is not measurable2. This means that only quantities are experimentally accessible which describe the overall effect of positively and negatively charged ions in the solution, e.g. the mean ionic activity coefficient which is defined as

2 This is generally accepted, except by the group of Wilczek-Vera, Rodil and Vera, who deter-

mined individual ion activity coefficients [19-21] from electromotive force data. However, their work was strongly offended by Malatesta [22,23]. Because individual ion activity coefficients are only published by one group and their work is in dispute, we constrain to the classical view which accepts mean ionic activity coefficients.


( ) 1

++ = (2-17)

while + and - are the individual ion activity coefficients. The individual ion activity coeffi-cients can be defined in the mole fraction scale, in the molality scale and in the molarity scale. They describe the deviation of the individual ion chemical potential from a standard state, according to

( ) ( )( ) += xiixii xRT ln

( ) ( )

+= mi

kgmol

imii

mRT

1ln

( ) ( )

+= ci

lmoli

ciic

RT 1

ln (2-18)

For the standard state of the chemical potential in equation (2-18) the hypothetical ideal solu-tion at unit concentration can be used. Unit concentration means that xi = 1, mi = 1 mol/kg, and ci = 1 mol/l. However, electrolytes are not ideal at unit concentration, e.g. at mi = 1 mol/kg and they are nonexistent as a pure liquid (xi = 1) at T and p. Therefore, the unsymmetrical reference is used, where the activity coefficients are defined as 1 for the infinite dilute solution, which means i(x)*1 as xi0, i(m)*1 as mi0 and i(c)*1 as ci0 in equation (2-18). Mean ionic activity coefficients in experimental data collections [18,24] are provided in the molality scale.

The rational ion activity coefficient i(x)* is defined as the ratio of two fugacity coefficients. The nominator is given as the fugacity coefficient for the ion at T and p and in the denomi-nator is the fugacity coefficient at same T and p but at reference concentration, i.e. at infinite dilution. For the unsymmetrical convention according to equation (2-18) the activity coefficient i(x)* reads

( )( )0,,

,,

*)(=

=

ii

iixi

xpTxpT

(2-19)

Using equations (2-17), (2-19) and the conversion from mole fraction to molality, the mean ionic activity coefficient in the molality scale * (m) can be calculated by [25]

( )( )

( )( )

+

=

=

+=

++

++ 0,,

,,

0,,,,

11*

)(xpT

xpTxpT

xpTMm Solvent

m (2-20)

where the indices + and stand for the cation and anion respectively. In Appendix A the conversion of activity coefficients using different concentration scales and reference states is presented for a better understanding of equation (2-20).


The second important thermodynamic quantity is the osmotic coefficient in the molality scale

(m) which represents the activity coefficient of the solvent. It can be calculated using the

fugacity coefficients of the solvent solv.

( ) ( )( )

( )solv

solvsolv

solvsolvsolv

im

MmxpT

xpTx

nPT

=

=

1,,,,ln

,,

(2-21)

In contrast to the activity coefficient of the solvent solv, the osmotic coefficient (m) is more sensitive to the presence of the electrolyte.

In conclusion it can be said, that thermodynamic properties for electrolyte systems can be calculated using common thermodynamic concepts for non-electrolytes. However, attention has to be paid to the different concentration scales and the unsymmetrical convention of the ion activity coefficients. When calculating the mean ionic activity coefficient * (m) (2-20), each ion is treated as a separate species. The concentration of each ion is known for salt/solvent systems that approach complete dissociation. Thermodynamic models for electro-lyte systems need to take into account the strong electrostatic interactions of the ions and the solvent. Statistical mechanics delivers some concepts for modeling of these interactions. They are described in the next chapter.

11

3 Statistical Mechanics Statistical mechanics relates microscopic properties of individual atoms and molecules to macroscopic observable quantities. Quantum and classical mechanics deals with forces and motions of electrons, atoms and molecules, and statistics is used to produce macroscopic aver-ages. The purpose of this chapter is to introduce the concepts of statistical mechanics, which are used to study dense systems with electrostatic interactions. Two important methods are used to study such systems, one are statistical mechanical theories, e.g. integral equation theo-ries and perturbation theories and the other ones are molecular simulations. Both methods are used in this work and explained in this chapter.

In section 3.1 some important terms of statistical mechanics are introduced. Section 3.2 focuses on models for the description of dense systems. In section 3.3 statistical mechanical models for dense systems with electrostatic interactions are presented. In section 3.4 the Monte Carlo simulation method is explained, which is a molecular simulation method based on random moves of particles in a simulation box. The main focus of this section is directed to the simulation of phase coexistence points using Monte Carlo simulations in the grand canonical ensemble and histogram reweighting. The chapter closes with a literature survey on molecular simulation approaches for electrolyte systems, and on thermodynamic models for electrolyte systems.

3.1 Introduction

The objective of this section is to introduce some important terms of statistical mechanics, which are necessary to understand the methods which are explained later in this work. These are terms like effective pair potential, pairwise additivity assumption and statistical mecha-nical ensemble.

3.1.1 A Statistical Mechanical System

A system in the statistical mechanical point of view consists of N material particles, which could be molecules, charged ions, or colloidal particles. Such a system is called an N-body system and is described by the linear momenta ip and the center-of-mass positions ir of all N particles. The total energy or the Hamiltonian of such a system is the sum of the kinetic and the potential energy contribution. The kinetic energy is associated with the motion of the particles, which can be divided into translation, rotation and vibration. It is the main form of energy in low-density systems, e.g. gases. At increasing density the potential energy resulting from the interactions between the particles becomes dominant. The strength of the potential energy between particles is determined by the distance, and orientation between them; they thus depend on their center-of-mass positions ir and their orientations, as expressed by the

Euler angle i .

The total potential energy of a system can be determined by evaluating the interactions bet-ween all molecules or their interaction sites respectively

12 3 Statistical Mechanics

( ) ( )( ) ++

+=

> >>

>

i ij ijkkkjjii

i ijjjii

iii

NN

rrru

rrururU

,,,,,

,,,,,

3

21

(3-1)

The symbols >i ij

and > >>i ij ijk

denote the summation over distinct pairs or triples of i, j

and k without counting any pair or triple twice. The first term in equation (3-1) represents the effect of an external field. The double sum in equation (3-1) represents the contribution to the total potential energy which results from the pair interactions. Pair interactions depend on the distance and the orientations of each pair of interaction sites in the system. The third part is the potential energy due to three body interactions, which depend on the distances and orientations of each triple of molecular sites. Usually, three body and higher body interactions deliver only a small amount to the total potential energy and are therefore neglected in this work. Their contribution to the total potential energy can also be considered by so called effective pair potentials. These effective pair potentials include higher body effects and other possible deficiencies of the molecular model by adjusting their parameters to macroscopic fluid properties.

In the absence of an external field and neglecting the higher order terms equation (3-1) can be simplified to

( ) ( )>

=

i ijijij rurU 2 (3-2)

This is known as pairwise additivity assumption, and is used in this work. In equation (3-2) the pair potential u2(rij) is abbreviated assuming that it is orientation independent. In the following sections the lower case index 2 is omitted for brevity.

One example for an effective pair potential is the 12-6 Lennard-Jones potential which is used frequently to describe the interactions of molecular sites.

( )

=

612

4ijij

ijLJ

rrru

(3-3)

This pair potential is not orientation dependent. Two parameters are needed: the size parameter and the energy parameter . The first term in equation (3-3) describes the soft repulsion of molecular sites. The second term describes the dispersion.

The simplest pair potential to describe dense systems is the hard-sphere potential. This pair potential was an impetus for the development of statistical mechanical models for liquid systems, because these systems are (in their microstructure) dominated by repulsive inter-actions. It is described by

3 Statistical Mechanics 13

( ) +=ru for r ( ) 0=ru for >r (3-4)

Only the hard-sphere diameter is needed for the pair potential in equation (3-4). 3.1.2 Statistical Mechanical Ensembles

A statistical mechanical ensemble is a virtual collection of an infinite number of physical systems with a set of defined properties, say, defined number of particles N, system volume V, and system energy U, but different microstates. This concept was introduced by Gibbs in 1902 as the conceptual foundation of statistical mechanics [26], p. 21. Each statistical mechanical ensemble is characterized by a number of constant state properties. In the canonical ensemble the temperature T, the volume V, and the number of particles N in the system are fixed. Even if these state properties are fixed, there is an infinite number of possibilities for the center-of mass positions and the momenta of each of the N particles in the system. Constant number of particles N, constant volume V and constant energy U give the microcanonical ensemble. If the number of particles N, the temperature T and the pressure p are constant the ensemble is called isobaric-isothermal. The grand canonical ensemble is an open ensemble, where the number of particles N is not fixed, but the chemical potential i, the volume V and the temperature T.

For every one of these ensembles, statistical mechanics provides a link of microscopic states of the ensemble with macroscopic properties. In the canonical ensemble (with fixed N,V,T) which corresponds to a closed system with energy exchange over the system boundaries, the required information on the microscopic states is collected in the partition function of the canonical ensemble, defined as

( ) ( )

=

j

j

kTVNU

TVNQ ,exp,, (3-5)

where Uj denotes the total energy of the systems when it is in state j. In actual fact, in classical statistical mechanics, the partition function cannot be correctly represented as a sum of discrete terms. Because the position and momentum variables can vary continuously, the number of microstates for a given macrostate is uncountable. Therefore, the sum in equation (3-5) can be better represented by an integral. If the Hamiltonian of an N-particle classical system UN comprises only potential energy and translation, the canonical partition function can be written as

( ) ( ) NNNNNN rddrdpdpkTrpU

NTVNQ

= 113

,

exp!

1,, (3-6)

with denoting the de Broglie wave length.

The partition function Q(N,V,T) is in turn related to the Helmholtz energy according to


( )TVNQkTA ,,ln= (3-7) The knowledge of the partition function of the canonical ensemble is thus sufficient to calculate any macroscopic thermodynamic quantity because A(N,V,T) is one of the fundamental potentials, see section 2.1.

The microcanonical ensemble represents an isolated system with no exchange of matter and energy, and the partition function (N,V,U) is here connected to the entropy S.

( )UVNkS ,,ln = (3-8) For the grand canonical ensemble representing an open system the partition function (,V,T) is connected to the product of pV.

( )TVkTpV ,,ln = (3-9) In the isobaric-isothermal ensemble the partition function (N,p,T) is connected to the Gibbs energy G.

( )TpNkTG ,,ln = (3-10) The equations for the partition functions of the microcanonical Q(N,V,T), the grand canonical (,V,T) and the isobaric-isothermal ensemble (N,p,T) can be found in statistical mechanics textbooks, e.g. Lee [26] and McQuarrie [27]. The equations (3-7)-(3-10) that link the partition functions for the different ensembles to thermodynamic quantities are needed to understand the histogram reweighting method which is introduced in section 3.4.

In principle, the canonical partition function would be sufficient to calculate all thermodyna-mic quantities. However, this is only possible for very simple systems, e.g. ideal-gas systems where the molecules have no potential energy. An alternative to partition functions for dense systems is the use of distribution functions [26]. This is explained in the next section.

3.2 Statistical Mechanics for Dense Systems

The concept of distribution functions was developed to describe dense systems. Distribution functions give time-averaged spatial configurations of the molecules. The use of distribution functions in the gas phase is not useful because the molecules are very far apart and therefore no density inhomogeneities occur. Radial distribution functions are accessible in three ways, by X-ray and neutron scattering, from integral equation theories and from computer simu-lations.

In section 3.2.1 different distribution functions are introduced. Section 3.2.2 shows how inte-gral equation theories can be used to determine distribution functions. In section 3.2.3 it is presented how thermodynamic quantities can be calculated using distribution functions. In section 3.2.4 perturbation theories are introduced, which can be used for determining thermo-dynamic quantities in dense systems, too.


3.2.1 Distribution Functions

For isotropic fluids like the Lennard-Jones fluid the radial pair correlation function g(2)(ri,rj), depends only on the distance between the molecules. The distribution functions of non-isotropic fluids additionally depend on the Euler angles i of the molecular sites. In this section only isotropic fluids are considered. The radial pair correlation function g(2)(ri,rj) specifies the spatial distribution of pairs of molecules. It describes how the density of surrounding matter varies as a function of the distance from a distinguished point. The super-script (2) indicates that pairs of molecules are considered. This subscript is usually left out for brevity. In Figure 3-1 the radial pair distribution function for a Lennard-Jones fluid at constant reduced density * and different reduced temperatures T* is shown. This plot was generated from molecular dynamics simulations [28].

0

0.5

1

1.5

2

2.5

3

3.5

0 1 2 3 4 5r/

g(r)

T*=1.095 T*=0.936T*=0.591

*=0.88

Figure 3-1: Radial pair distribution function of a 12-6 Lennard-Jones potential fluid from molecular dynamics simulations [28]. The lines are drawn to guide the eye.

As can be seen in Figure 3-1 at a reduced distance /r of 1 the probability of finding another molecule is very high, due to collisions between the molecules. For large distances the probability approaches the value of 1. For the ideal gas g(r) is always 1. For the Lennard-Jones fluid at the given T* and * in Figure 3-1 the radial pair distribution g(r) approaches 1 for a reduced distance of r/ = 4. This reduced distance (here /r = 4) is an indicator for the minimum cutoff in a molecular simulation.

Another distribution function in statistical mechanics is the total correlation function h(r). It is closely connected to the radial pair distribution function g(r) and for an isotropic fluid defined as


( ) ( ) 11212 = rgrh (3-11) The third correlation function which needs to be introduced in this work is the direct corre-lation function c(r). It is related to the total correlation function h(r) through the Ornstein-Zernike equation3, which reads for an isotropic fluid as

( ) ( ) ( ) ( ) 332131212 rdrhrcrcrh += (3-12) The first term in equation (3-12) gives the direct part and the second term in (3-12) gives the indirect part of the total correlation function h(r12). As was already indicated, correlation functions can be determined by means of experimental methods, molecular simulations, and by integral equations. The concept of integral equations is explained in the next section.

3.2.2 Integral Equation Theories

Integral equations are a fast method for the determination of distribution functions needed to calculate thermodynamic quantities in dense systems.

Figure 3-2: Task of an integral equation [26].

Figure 3-2 shows the connection between the pair potential and pair correlation functions by integral equations.

Some integral equation theories are based on solving the Ornstein-Zernike equation (3-12), e.g. the mean spherical approximation and the Percus-Yevick equation. The resulting radial pair distribution function of such an integral equation theory is not exact, because an assump-tion is needed for the solution of the Ornstein-Zernike equation. This assumption is called closure and denominates the Ornstein-Zernike based integral equation theory. For the Percus-Yevick equation the closure reads as

3 Some textbooks perceive the Ornstein-Zernike equation as a definition of the direct

correlation function c(r), e.g. [26]. Others, (most prominently Hansen and McDonald [29]) argue that c(r) is independently defined as a derivative of the intrinsic free energy, and that consequently the Ornstein-Zernike equation only relates the direct correlation function c(r) and the total correlation function h(r).


( ) ( )

=

kTu

rgrc exp1 (3-13)

For a hard sphere fluid, see equation (3-4), where g(r) = 0 for r < and u(r) = 0 for r > , the Percus-Yevick closure simplifies to

( ) 1=rh for r (3-14)

Another integral equation approach for solving the Ornstein-Zernike equation is the mean spherical approximation (MSA) [30]. For a fluid with hard repulsion the pair potential can be written as

( ) =ru for r ( ) ( )ruru soft= for >r

(3-15)

For a fluid with hard repulsion the radial distribution function must be zero for distances smaller than the contact radius. The second approximation is developed from a cluster expan-sion of the direct correlation functions and is only exact for infinite intermolecular distance r. The MSA closure can be written as

( ) 0=rg for r

( ) ( )kT

rurc

soft for >r

(3-16)

For a hard sphere fluid the MSA coincides with the Percus-Yevick approximation. From equations (3-15) and (3-16) it becomes apparent that the MSA is restricted to fluids with hard repulsion.

Integral equation theories were used to determine pair correlation functions for many model fluids, see for example [31,32]. The results can be compared to the pair correlation functions from molecular simulations. In the next section it is shown how integral equations can be used for the determination of thermodynamic quantities. However, one has to keep in mind that distribution functions from Ornstein-Zernike based integral equations are not exact, and therefore the thermodynamic quantities which are determined from them can be inconsistent.

3.2.3 Microthermodynamics

Microthermodynamics is the interpretation of macroscopic thermodynamic properties in terms of molecular forces [26]. There are three major routes to obtain thermodynamic quantities from molecular distribution functions and pair potentials, which are the energy equation, the pressure equation, and the isothermal compressibility equation.


Figure 3-3: Idea of microthermodynamics.

The idea of microthermodynamics is illustrated in Figure 3-3. For the calculation of the energy the pair potentials between the particles according to their spatial distribution are summed up. For an isotropic homogeneous fluid the energy is determined from

( ) ( )

+=0

2422

3 drrgrurkTNkT

Upi

(3-17)

For the pressure equation the virial ( ur ) is summed up, and gives

( ) ( )

=

0

346

1 drrgdr

rdur

kTkTp

pi

(3-18)

The isothermal compressibility can be written as

( )( )

=

0

2 141 drrgrkTKT pi (3-19)

Equations (3-17), (3-18), and (3-19) show that knowledge of the radial pair correlation g(r) for a given pair potential u(r) is sufficient to determine macroscopic properties. However, the radial pair distribution function g(r) needs to be known for all state conditions. Consequently, the experimental determination of g(r) for one or two state points is not sufficient. Unfortunately, thermodynamic correlation functions determined from Ornstein-Zernike based integral equations are not exact. Therefore, thermodynamic quantities determined by equations (3-17), (3-18), and (3-19) considering g(r) from, say, the Percus-Yevick equation or the mean spherical approximation, are also not exact. This is illustrated by the following example. If the Ornstein-Zernike equation with the Percus-Yevick closure is solved for a hard sphere fluid, and the pressure is calculated on the one hand with equation (3-19) for the isothermal compressibility and on the other hand with the virial equation (3-18), the results are different. The two pressures pv from the virial equation and pc from the compressibility equation were compared to molecular simulation data [33]. The best agreement with the molecular simulation data was found by combining the two results, with weights of 1/3 and 2/3, respectively, so that


kTp

kTp

kTp cv

32

31

+= (3-20)

This is the well-known Carnahan-Starling equation [33] that provided an impetus for many dense fluid theories.

Besides these main routes to thermodynamic quantities (equations (3-17), (3-18), and (3-19)) statistical mechanics also provides equations for the chemical potential i and the Helmholtz energy A [26,27].

3.2.4 Perturbation Theories

The idea of perturbation theories is to describe the properties of a fluid as those of a simpler fluid, the so-called reference fluid plus a correction. In perturbation theories the intermole-cular pair potential is split in two parts, one is the reference part u0(r) and the other one is the perturbing part u1(r).

( ) ( ) ( )rururu 10 += (3-21) The reference part u0(r) can be the pair potential of a fluid with repulsive interactions only, e.g. a hard-sphere fluid. Alternatively, the reference fluid with the reference potential u0(r) can just as well be a fluid with repulsive and dispersive interactions, e.g. a Lennard-Jones fluid. The conditions for the perturbation theory to be effective is, that the Helmholtz energy A0(,T) and the radial pair distribution function g0(,T,r) for the reference fluid, are accurately known. The Helmholtz energy A of the fluid can be determined from a power series in terms of 1/kT, and can be written as

( ) ( ) ++= 33

2210

62 kTkTkTkTA

kTA

(3-22)

where 1, 2 and 3 are the perturbation terms of first, second and third order, respectively. The perturbation term of first order 1 equals [27]

( ) ( )= drrgruV 012

1 2

(3-23)

From (3-23) it is apparent that the radial pair distribution function of the reference fluid g0(r) needs to be known. For the determination of the higher order perturbation terms approxi-mations are used. This is explained in detail by McQuarrie [27]. Perturbation theories were used extensively for the development of equations of state. Impor-tant examples are the SAFT theory by Chapman et al. [34] which describes thermodynamic properties of hard-sphere chains using the thermodynamic perturbation theory of first order with the hard-sphere fluid as the reference. This work served as a basis for many equations of state for small molecules up to polymers, because it takes the effect of non-sphericity on the thermodynamic properties into account. One further development of the SAFT equation of


state is the PC-SAFT equation of state [10] that describes the dispersive interactions of chain molecules using the second-order thermodynamic perturbation theory of Barker and Hender-son [35]. The reference fluid for the dispersive interaction of the PC-SAFT is the hard-chain fluid.

3.3 Theories for the Description of Electrostatic Interactions

Electrolyte systems are not only dominated by short-range interactions but also by long-range attraction and repulsion, see Prausnitz et al. [5], p. 524 et seq. This is due to the electrostatic interactions, i.e. pair potentials decreasing slowly with increasing intermolecular distance, e.g. with r-1 for point charges (r denoting the radius). In addition, the electrostatic interactions at short range are very strong, what makes modeling and molecular simulation difficult [6].

For simplification of the modeling, most theories focus on describing the electrostatic inter-actions of ions that are surrounded by a dielectric continuum. This dielectric continuum repre-sents the solvent implicitly and screens the charges of the ions. Models that do not account for the molecular structure of the solvent are classified as implicit solvent or primitive models. However, such simplified models fail in describing characteristics of electrolyte systems, e.g. the ordering of solvated molecules around the ions. This is the reason why this work focuses on theories that explicitly account for the molecular structure of the solvent, so called explicit solvent or non-primitive models.

The structure of this chapter is as follows. In sections 3.3.1 an introduction into electrostatic pair potentials in the MKS and Gaussian system of units is given. Section 3.3.2 gives an over-view about two established theories, namely the Debye-Hckel and the Born theory. In section 3.3.3 perturbation theory approaches and in section 3.3.4 integral equation theories for modeling of electrostatic interactions are introduced.

3.3.1 Electrostatic Pair Potentials

For modeling and simulating electrolyte systems suitable pair potentials need to be chosen. These pair potentials need to represent the interactions of the real fluid properly, but also need to be simple enough to allow a theoretical description to be developed. Usually, the ion-ion electrostatic interaction is considered to be the key interaction in electrolyte solutions. The force between two ions is described by Coulombs law. If the ions are immersed in a solvent their charge is screened. The screening effect of the solvent can be described by its dielectric constant .

If the solvent is explicitly considered as a molecular species, the electrostatic interactions between the solvent molecules and between solvent molecules and ions need to be considered. All liquids that dissolve salts in a considerable amount are polar and polarizable. The dipolar interaction between the solvent molecules and the ion-dipole interaction between the ions and solvent molecules contribute significantly to the thermodynamic behavior of such mixtures. Non-primitive models account for these interactions in addition to the Coulomb interaction


between the ions. However, the restriction to these three interactions implies that effects due to higher multipole moments as well as polarizability effects are neglected.

If electrostatic potentials are studied, two different systems of units are used in literature. One is the SI or MKS system of units and the other one, used almost exclusively in electrolyte literature, is the Gaussian or CGS system of units. MKS stands for meter, kilogram and second, and CGS stands for centimeter, gram and second. The inconvenience of converting between the two systems becomes apparent if Coulombs law is considered; it reads

MKS CGS

r

r

ezzF ji 3

2

041

pi= r

r

ezzF ji 3

2

=

(3-24)

Equation (3-24) shows the force acting on two charges in vacuum. In the MKS system the charge is measured in Coulomb and the permittivity of free space 0 has the value of 8.8541910-12

C/Nm. In the CGS system there is no constant of proportionality. The charge is

given in units of 23

scmg

or esu. The conversion for the charge between the MKS and

CGS system is given by 1 C = 2.998109 esu. Further details about the two systems can be found in Purcell [36], p. 7 et seq. This work, in accord with the bulk of the literature on the subject, uses the CGS system exclusively. Table 3-1 summarizes the pair potentials for the electrostatic interactions in ion-dipole fluids, using the vector notation and the GGS system of units.

Table 3-1: Pair potentials for electrostatic interactions [6] in the CGS notation

Pair potential Equation

ion-ion ( )ij

jicc

r

qqru = with ezq ii =

ion-dipole ( ) 3,,,ij

jijiDiDij

cd

r

qrru

=

dipole-dipole ( ) ( )( )5 ,,3 ,,,, 3,,ij

ijjDijiD

ij

jDiDjDiDij

dd

r

rr

r

ru

=

The long-range nature of the electrostatic potential becomes apparent if the potentials are inserted into equation (3-17). Using this equation one can verify that for the ion-ion potential the potential energy assumes infinity.


In order to obtain a simple model for ions and dipolar molecules, the potentials in Table 3-1 can be combined with the hard-sphere potential (3-4) or with the Lennard-Jones potential (3-3). The next sections give an overview over statistical mechanical models for the electrostatic ion-ion, ion-dipole and dipole-dipole interactions. The main focus lies on models for the ion-ion and the ion-dipole interactions because the dipole-dipole interactions are not exclusively relevant in electrolyte systems. Furthermore, the dipole-dipole potential decreases faster with increasing intermolecular distance compared to the ion-ion and ion-dipole potential. Therefore, integral equation theories and perturbation theories for dipolar interactions can be found in literature, see for example [31] for an integral equation theory and [37] for a perturbation theory.

3.3.2 Historical Theories

The Debye-Hckel and the Born theory are starting points for many thermodynamic models for electrolyte systems. Therefore, they are introduced briefly in this section.

Debye-Hckel Theory

The Debye-Hckel theory [4] is a milestone in electrolyte theories [26] because it was the first model which describes properties of dilute electrolyte solutions. Debye and Hckel considered the ions to be point electrical charges in a dielectric continuum, representing the solvent, e.g. water. This assumption is valid if the amount of solvent is much greater than the amount of ions [18].

In the theory the electrical potential at a point in the solution is calculated in terms of the concentrations and charges of the ions and of the dielectric constant of the solvent. The starting point of the theory is the Poisson equation which links the symmetrical electrical potential with the electrical charge density e by

pi e42= (3-25)

Debye and Hckel assumed a Boltzmann distribution for the charges with

=

i

iiie kT

ezez exp (3-26)

where zie denotes the charge of the ions of species i, and i is the number density of ions of species i.

If the potential energy of the ion zie is very small compared to its thermal energy kT, which comes true for very dilute solutions, because the ions are far away from each other, equation (3-26) can be linearized. This yields


=

i

iie kT

ez 22 (3-27)

If the Poisson equation (3-25) and the linearized Boltzmann equation (3-27) are combined, a second order differential equation for the potential as a function of r is obtained

=

=

2

22

22

41

pi

iii zkT

e

drd

rdrd

r (3-28)

where is the Debye length in units of inverse length. By solving the linearized Poisson-Boltzmann equation the potential is determined. From the potential thermodynamic properties can be calculated, e.g. the internal energy and the chemical potential. Only the dielectric constant , the temperature T, the valence of the ions zi, and the number density of the ions of species i i = Ne,i/V need to be specified. The assumption for linearizing the Boltzmann equation (zie


solvent. Henderson considered perturbation terms up to fourth order. The perturbation term of first order is canceled due to the electroneutrality condition. The perturbation series converges very slowly. In addition, some integrals in the perturbation series are divergent because of the long range of the Coulomb potential. The divergent terms need to be summed skillfully. The perturbation series gives as the first term the Debye-Hckel term and thereafter the primitive MSA term. The main deficiency of the theory is the simplifying assumption that the solvent is regarded as a dielectric medium.

The theory was later extended to the non-primitive framework. Non-primitive framework means, that the solvent is explicitly considered as a molecular species. In the non-primitive perturbation theory hard-sphere ions with point charges are surrounded by hard-sphere solvent molecules with point dipole moments. Henderson et al. [7] expanded the Helmholtz energy up to the third order. The first order term cancels because the system is electroneutral and because of angle averaging. Like for the primitive perturbation theory, the difficulty is that most of the integrals diverge, because the electrostatic potentials are long-ranged. However, Henderson et al. sum up the divergent terms in order to cancel the appropriate divergent terms and receive the limiting cases, namely the perturbation theory results for the primitive model and for the pure dipolar hard-sphere system.

The ideas of Henderson et al. [7,40] were adopted by Jin et al. [41], and by Wu et al. [42] who modified some perturbation terms. In conclusion, it can be stated, that the main drawback for the application of perturbation theories is the slow convergence of the expansion series and the divergence of most perturbation terms. The divergence is caused by the long-range character of the electrostatic interactions. Furthermore, it needs to be pointed out that the primitive and non-primitive perturbation theories are based on a hard-sphere reference fluid.

3.3.4 Integral Equation Theories

Like for perturbation theories integral equation theories were developed for the primitive and the non-primitive framework. The most influential integral equation theory for electrolyte systems is the mean spherical approximation (MSA), see paragraph 3.2.2. The importance of the MSA for the hard-sphere ion and the hard-sphere ion-dipole system results from the fact that the solution is algebraic [8,9,43,44]. In contrast, the distribution functions for the hypernetted chain integral equation theory for hard-sphere ion-dipole mixture [45] need to be integrated numerically. However, the closure of the MSA rests on the condition that the molecules have hard repulsion, see equations (3-15) and (3-16). Primitive MSA

The primitive MSA was solved by Blum and Hye [43,44]. It gives thermodynamic properties and the distribution functions of a system of hard sphere ions embedded in a dielectricum. Hard sphere ions are described by their charge zie and their diameter i. The whole system is electroneutral.

The solution of Blum and Hye is demanding. Therefore only some key ideas are elucidated in this section. Blum and Hye solved the Ornstein-Zernike equation (3-12) using the method


of Baxter [46,47]. This leads to an algebraic equation, whose numerical solution gives a single screening parameter . The excess internal energy can be expressed using the energy equation (3-17) and the direct correlation functions c(r). The Helmholtz energy A can be determined from the energy U by integration, viz.

UA =

(3-30)

Because this integration is tedious Hye and Stell [48] developed a sophisticated method that combines results of the excess pressure by the energy (3-17), the virial (3-18) and the compressibility (3-19) route, and accounting for the inconsistencies of the MSA. The method Hye et al. [48] was used by Blum et al. [44] to determine the excess internal energy U, the excess Helmholtz energy A, the excess pressure p over a hard-sphere system, as well as the mean ionic activity coefficient * (m). In addition, the paper of Blum and Hye [44] presents radial pair distribution functions for the ions. The excess Helmholtz energy, which can be used to calculate all other thermodynamic quantities, is given [44] as

( )

+

+

+=

pi

pi

321

1 3222 kT

Pze

kTNkTA

Ni i

ii

total

MSA

(3-31)

In equation (3-31) , and PN are auxiliary quantities in the primitive MSA. The equations to calculate them are given by Blum et al. [44]. The iterative determination of the screening parameter requires that the pressure and the chemical potential according to equations (2-2) and (2-4) respectively, need to be determined by numerical partial derivation. If the ions have the same hard core diameter the solution of the so-called restricted primitive MSA is simplified because the screening parameter does not need to be determined iteratively.

Non-primitive MSA

The underlying model system for the non-primitive MSA (np MSA) is a mixture of charged hard-sphere molecules and dipolar hard-sphere molecules. The electroneutrality condition prevails. Blum et al. [8,9,49] used a similar solution method for the np MSA compared with the primitive MSA. However, because the ion-dipole and the dipole-dipole potentials are orientation dependent, see Table 3-1, an expansion is carried out in order to get rotation invariant parts of the total correlation functions h.

The solution of the np MSA yields three energy parameters that need to be determined iterative by solving a set of three nonlinear algebraic equations. All thermodynamic quantities can be expressed in terms of these energy parameters. Three different variants of the np MSA are known. For the most general case, the diameters of the ions are different from each other and different to the one for the dipolar solvent molecule. This is called non-restricted np MSA. If the ion diameters are equal but different from the solvent diameter, the np MSA is called semirestricted, and if the diameters of the ions are equal to the diameter of the solvent it is called restricted np MSA.


Like for the primitive MSA the excess thermodynamic quantities respective to a hard-sphere system are determined by the energy route, see equation (3-17). This leads to a surprisingly simple equation for the excess internal energy U over a hard-sphere reference system. The other thermodynamic quantities detailed in the papers, are the excess Helmholtz energy and the chemical potential. The excess internal energy and the excess Helmholtz energy are linked through the compressibility factor by the simple relationship

NkTA

NkTUZ = (3-32)

The papers of Blum et al. [8,9,49] are difficult to understand from a chemical engineers per-spective. The model equations for the np MSA are independently rederived only for the re-stricted case by Hye et al. [50,51]. Moreover, there are several misprints in the papers of Blum [8,9,49], that are only partly corrected. Nevertheless, the np MSA is a promising theory because the excess thermodynamic quantities of a hard-sphere ion-dipole system can be determined by algebraic equations. In order to apply the results of the np MSA to real electro-lyte systems the model equations and the equation for the excess Helmholtz energy for the semirestricted np MSA are presented in chapter 4. Strategies for verifying the model equations and the equations for the thermodynamic quantities are also presented in chapter 4.

The results of the statistical mechanical theories can be compared to molecular simulation data. For polar fluids and ion-dipolar mixtures this was done by Liu et al. [52]. The concepts and applications of molecular simulation in deriving and comparing statistical mechanical fluid theories are explained in the succeeding chapter.

3.4 Molecular Simulations

Molecular simulations can be referred to as computational statistical mechanics [53] because they deliver exact macroscopic properties (subject to only statistical uncertainties and possible finite size effects) for a statistical mechanical system. In this section, the idea and implemen-tation of molecular simulation methods are explained. More details about molecular simula-tions in general can be found in the textbooks of Allen and Tildesley [54]4, Frenkel and Smit [55], and Sadus [53].

The chapter is structured as follows. In section 3.4.1 computer simulations are introduced and the connection to fluid theories is explained. In section 3.4.2 some of the basic concepts of Metropolis Monte Carlo simulations are given. In section 3.4.3 Monte Carlo simulations in the grand canonical ensemble are presented. This ensemble is chosen in order to simulate phase equilibria. Subsequently the method that was used to analyze the grand canonical Mon-te Carlo simulation data is explained in sections 3.4.4 and 3.4.5. Paragraph 3.4.6 addresses the determination of the critical point from Monte Carlo simulation data. The last section 3.4.7

4 Several examples and pieces of FORTRAN code, explained in the book of Allen and

Tildesley, are provided on the CCP5 website, http://www.ccp5.ac.uk/librar.shtml.


introduces the Ewald summation that was used to compute the long-range electrostatic interactions.

3.4.1 Introduction

Molecular simulations are computer experiments carried out in order to get quasi-exact results for macroscopic properties (e.g. the potential energy or the pressure) of a defined model fluid. Defined model fluid means that the pair potentials between molecular sites are specified (often neglecting 3-body and higher-body potentials as well as quantum effects). The term molecular simulation stands for Molecular Dynamics and Monte Carlo methods. Molecular Dynamics simulations are based on the solution of Newtons equations of motion. This method investigates the evolution of a statistical mechanical system on a time scale. Therefore, with Molecular Dynamics it is possible to determine time-dependent quantities, e.g. diffusion coefficients. Monte Carlo simulations in the canonical ensemble are based on random displacement of molecules in a simulation box, which are accepted with a certain probability. The sampled states in a Monte Carlo simulation are not linked on a time scale. Therefore, Monte Carlo methods can be only u

Dissertation StefanieHerzog

Documents

Transcript of Dissertation StefanieHerzog