THERMODYNAMICS AND STATISTICAL MECHANICS OF
SYSTEMS OF "REACTIVE" COMPONENTS WITH
APPLICATIONS TO STRONG ELECTROLYTES
BY
RANDAL LEWIS PERRY
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1980
... to my parents, for their continual love, encouragement and
support throughout the years.
ACKNOWLEDGEMENTS
The author would like to thank Dr. John P. O'Connell for his aid
and guidance in this research, and John C. Telotte for useful discus-
sions and his help in preparation of the manuscript. Thanks also goes
to Dr. Hans C. Anderson for his questions and comments which opened
doors to the final results, and to the National Science Foundation for
financial support.
TABLE OF CONTENTS
PAGE
ACKNOWLEDGEMENTS...................................................iii
KEY TO SYMBOLS......................................................vi
ABSTRACT........................................................... xii
CHAPTER
ONE INTRODUCTION.................................................1
TWO SOLUTION THERMODYNAMICS FOR "REACTIVE" COMPONENTS...........6
Introduction ...............................................6
Mathematical Basis ......................................... 7
Construction of Projectors--Material Balance.............. 10
Properties of Solutions................................... 11
Application to Analytic Solution of Groups: UNIFAC and
UNIQUAC.............................................. 26
THREE STATISTICAL THERMODYNAMICS OF "REACTIVE" COMPONENTS IN
THE GRAND ENSEMBLE .......................................... 33
Introduction .............................................. 33
Components and Species in Complex Systems.................34
Grand Ensemble for Systems of "Reactive" Components.......39
Composition Fluctuations................................... 45
Complete Reactions......................................... 54
FOUR SOLUTIONS OF STRONG ELECTROLYTES ............................61
Introductions............................................. 61
Material Balance Relationships among Salts, Solvents
and Ions................................... .......... 62
Formulation in Terms of Direct Correlation Functions......68
Thermodynamic Properties.................................. 72
Single-Salt, Single-Solvent System........................ 76
Double-Salt, Single-Solvent Case.......................... 81
Single-Solvent, Common Ion Case; Harned's Rule............ 86
FIVE ELEMENTS OF A MODEL FOR ELECTROLYTIC SOLUTIONS ..............91
Introduction.............................................. 91
Primitive Ion Model ........................................92
Calculation of Solution Properties.......................94
Electrostatic Contributions to Direct Correlation
Function Integrals...................................97
The Form of the Complete Model .......................... 100
Preliminary Results...................................... 101
SIX AN EXPERIMENT FOR VOLUMETRIC PROPERTIES OF SALT SOLUTIONS
AT HIGH PRESSURES.........................................104
Introduction............................................104
Thermodynamic Basis .....................................105
Experimental Design......................................107
Operating Procedures...................................... 112
Conclusions.............................................117
SEVEN CONCLUSIONS...............................................118
APPENDICES
A DERIVATIONS OF PRINCIPAL RELATIONS ........................122
B EXAMPLES OF THE PROJECTION PROCEDURE......................125
C EXPRESSIONS FOR THE UNIQUAC AND UNIFAC MODELS.............139
D CONSTRUCTION OF MATRICES FOR GENERAL MATERIAL BALANCE.....143
E DERIVATIVES WITH RESPECT TO COMPONENT CHEMICAL
POTENTIALS................................................146
F PROJECTION OF INVERSE MOLE FRACTIONS OF LIMITING
REACTANTS IN THE LIMIT OF COMPLETE REACTIONS..............151
G ALTERNATIVE FORMULATION OF FLUCTUATION PROPERTIES FROM
INTEGRALS OF THE RADIAL DISTRIBUTION FUNCTION............. 155
H INDEPENDENCE OF SALTS..................................... 160
I ELECTROSTATIC EXPRESSIONS FOR MULTISALT SYSTEMS........... 164
J HARD-SPHERE FORMULAE......................................167
REFERENCES........................................................ 170
BIOGRAPHICAL SKETCH.................................................. 176
KEY TO SYMBOLS
A species-composition fluctuation matrix
A matrix of composition derivatives of activities for components
a.,a activity of species i, component a
B. name of reactant j
Bik kth coefficient of activity coefficient expansion for ion i
C matrix of direct correlation function integrals
C matrix of integrals of nondivergent portion of direct corre-
lation functions
Co salt-salt "direct correlation function" integral
ctB
c.. "centers" direct correlation function
ij
c.. nondivergent portion of direct correlation function
ij
D.ik kth coefficient in direct correlation function integral ex-
pansion for ions i and j
ENq energy of system in state defined by [N,q]
e magnitude of charge on an electron
e unit vector with nonzero ath element
-ca
f(1) one-particle density function
f(2) two-particle density function
G,AGMIX Gibbs free energy, Gibbs free energy change on mixing
G transformation of A to Lewis-Randall system variables
G' transformation of A to density variables
gij pair correlation function
H,MI X enthalpy, enthalpy change on mixing
H matrix of integrals of total correlation functions
h.. total correlation function
IJ
I ionic strength
I,I identity matrix, identity matrix in component space
i vector of ones
Kk equilibrium constant for reaction k
K part of projection matrix for direct correlation functions
k Boltzmann's constant
L nonsquare matrix related to identity matrix
M,AMMIX general property, general property change on mixing
m molality
N total number of moles or particles
N total number of systems in grand ensemble
Ni,N number of moles of species i, component a
i oct
n,n ,n number of species, components, solvents
nNq number of systems in state defined by [N,q]
0 lowest order of remaining terms
P pressure
P. name of product species j
P projection matrix of component subspace in species space
PNq probability of system being in state defined by [N,q]
Q heat of reaction
Qi surface parameter for group i
q set of variables describing quantum state of system
q surface parameter for molecule a (UNIQUAC eqns.)
aq measure of charge strength of salt a (Salt eqns.)
R gas constant
R. volume parameter for group i
vii
vii
R coefficient in test for independence of salts
a
r number of reactions
r volume parameter for molecule a
r position vector
MIX
S,ASI entropy, entropy change on mixing
S ,SvSk coefficients for Debye-Hueckel limits for mean ionic activ-
ity coefficient, partial molar volume, apparent compres-
sibility for single salt system
S1 S1 1
S v,Sk same as S v,S except for multisalt system
T temperature
U internal energy
U total interaction energy of groups (UNIFAC eqns.)
U orthogonal left inverse of W
U.. interaction energy of groups i and j
u molecular total interaction energy
u orthogonal left inverse of v
u t interaction energy of molecules a and B
V,AVMIX volume, volume change on mixing
v,v molar volume, partial molar volume
W material balance matrix
X diagonal matrix of mole fractions
i.,x mole fraction of species i, component a
1 O
x set of mole fractions at reference state of component a
-(a)
Y constant left inverse of W
Yoa density-based activity coefficient for component a
constant left inverse of V
Z matrix of stoichiometric constraints in "reactive" system
Z constant portion of Z
-o
viii
Z nonconstant portion of Z
-1 _
z matrix of charge and stoichiometric constraints on system of
strong electrolytes
a coefficient from Mean Spherical Model for electrolytes
ay Harned coefficient for salt B in solution of salt y
B = 1/kT Lagrange multiplier for energy
B.. physical interaction energy parameter
1J
y mole fraction based activity coefficient
Ym molality based activity coefficient
6(r) Dirac delta function
6.. Kroniker delta
E,E1 dielectric constant of solution, solvent
C quantity for hard sphere contributions
1k weighted average of ionic size parameters
0. surface fraction of group i
O.. local surface fraction for groups i and j
0. surface fraction of group i in molecule a
i
6 surface fraction of molecule a
Ok Lagrange multiplier for constraint k
9a local surface fraction for molecules a and B
K inverse Debye length
K ,K1 isothermal compressibility of solution, solvent
X absolute activity of component a
1iUoa chemical potential of species i, component a
V matrix of stoichiometric coefficients
5 intensive extent of reaction
P,po total density of species, components
pi.,p density of species i, component a
(1) ensemble-averaged one-particle density
i
(2)
(2) ensemble-averaged two-particle density
ij
C. size parameter for species i (hard sphere diameter)
T r interaction parameter for molecules a and 3
Dk apparent compressibility for multisalt systems
D. volume fraction for group i
(k apparent compressibility for single-salt systems
(a volume fraction for molecule a
.ij interaction parameter for groups i and j
Q volume of orientation space
w vector of Euler angles
Subscripts
c chemical contribution to AGMIX
C set of components
D set of dependent species
i,j,.. species or group quantity
k, reaction number
L set of limiting reactants
n normalized quantity
o properties in component space
P set of products
p physical contribution to AGMIX
S set of solvents
a,B,... component or molecular quantity
E
elec
f
hs
MIX
T
_ (as in A)
(as in v.)
< > (as in )
o (as in o )
(as in v.)
(as in m )
in (as in ZnX)
Special Symbols
vector or matrix quantity
partial molar property
ensemble average
reference or standard state
infinite dilution value
mean ionic property
vector of logarithms of a quantity
Superscripts
excess property
electrostatic contribution
final state for integration of thermodynamic properties
hard-sphere property
reaction set in sequential reactions
mixture property
transpose of vector or matrix
Abstract of Dissertation Presented to the Graduate Council of the
University of Florida in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy
THERMODYNAMICS AND STATISTICAL MECHANICS
OF SYSTEMS OF "REACTIVE" COMPONENTS WITH APPLICATIONS
TO STRONG ELECTROLYTES
By
Randal Lewis Perry
December, 1980
Chairman: John P. O'Connell
Major Department: Chemical Engineering
The task of quantitatively describing physical properties of
solutions, whether for use in engineering design and optimization, or
for understanding natural processes, grows increasingly difficult as
more accuracy is demanded and more complex systems are studied. To
meet these demands, engineering thermodynamics now looks more thoroughly
into the intermolecular forces in a solution, and rely more on the the-
oretical, rather than empirical, bases for their quantitative descrip-
tions. This molecular basis is particularly important when "reactive"
systems--those systems where component molecules interact so as to form
different chemical species--are considered. This work develops a sta-
tistical mechanical solution theory for such systems, examining the
constraints caused by their "reactive" nature, anda model which may be
applied to correlate and predict their thermodynamic properties.
The thermodynamic consistency requirements of models for "reactive"
solutions are expressed. These are used to formulate a procedure, using
matrix projections, which builds these requirements into the solution
model. The statistical mechanics of such a solution in the grand canon-
ical ensemble is explored in light of the "reactive" constraints and the
consistency requirements. By using the projection procedure mentioned
above, the derivative, or fluctuation, properties of the components in
a general "reactive" system are related to equilibrium relations and
molecular correlation functions.
The case of total dissociations as the only "reactions" in the
solution is explored, with specific applications to strong electrolytes.
First, the singularities of the long range interactions are shown to be
removed by the projection procedure. Then a model is formulated based
on a composite of hard-sphere and charge effects. This model, though
presently not complete, shows considerable promise in correlating the
properties of salt solutions. An experimental apparatus for volumetric
behavior of liquids was completed and methods are described for using
data from it to evaluate the parameters in the solution model, as well
as develop a better understanding of density effects in electrolyte
solutions.
xiii
CHAPTER 1
INTRODUCTION
Design and optimization of industrial processes and equipment,
as well as quantitative understanding of natural or biological pro-
cesses, require quantitative descriptions of the physical properties
of the components of those systems. A very important need is the cor-
relation and prediction of P-V-T-x relationships and phase equilibria
in multicomponent liquid systems. Meeting this need can be quite dif-
ficult, especially if there are complicated intermolecular interactions
or tight constraints on accuracy. The situation is complicated further
in "reactive" systems, where the relative amounts of the species in a
given system are constrained by equilibrium conditions or stoichiometry.
These "reactive" systems include, in general, solvation, dissociation
and association--including processes such as micellization or ioniza-
tion of salts. A large number of real systems exhibit this "reactive"
nature, although it may proceed to such a small degree that it is
negligeable. In other situations, however, the additional species and
constraint relations contribute appreciably to the solution properties.
Although specific types of systems have been treated previously (e.g.
Renon and Prausnitz, 1967), there has not been a general approach to
these systems which could serve as a guideline for the consistent use
of solution models.
Solutions involving electrolytes are a class of "reactive" systems
which are of particular interest. Electrolytes appear in geological
liquid systems (oil reservoirs, geothermal energy sources), in systems
of biological interest (e.g. in active transport in membranes and in
using salts to separate proteins), and have important applications to
industrial systems (extraction, distillation and reverse osmosis).
These systems have an additional complication, though, because of the
presence of long-range coulombic interactions in addition to equili-
brium and stoichiometric constraints.
Statistical mechanical approaches have been successful for many
non-"reactive" systems, as well as some specific "reactive" systems
(Friedman and Ramanathan, 1970). A recent approach using fluctuation
solution theory (Mathias and O'Connell, 1980a,b) has shown application
for gases dissolved in various solvents at high pressures. This pro-
cedure is based on molecular correlation functions and works particu-
larly well in systems which contain components of very different na-
tures. Hence, it should be useful for salt solutions.
This work aims, firstly, to treat "reactive" systems in a manner
that is applicable to any system, and so to present a general and
unifying thermodynamic framework in which correlations can be consis-
tently developed. This framework is then applied to the above statis-
tical mechanical method to develop a formalism for expressing the
properties of "reactive" solutions in terms of these molecular correla-
tion functions. Finally, this formalism is applied to solutions of
strong electrolytes and a solution model is examined.
Chapter two develops the classical thermodynamics of "reactive"
systems. There, the "reactive" nature is defined mathematically in
terms of the material balance (stoichiometric) relations and intensive
extents of reactions. The total amounts of chemical species in a sys-
tem are related to the amounts of components (where the components each
represent a degree of freedom in composition) by a matrix relation which
is, in effect, a projection of the species space into the component sub-
space. The projection operators are completely defined by stoichiometry
and equilibria, and are used to relate solution properties in terms of
species to the required solution properties in terms of components.
These relations also show restrictions'which must be placed on solution
models for such systems, and two specific examples are detailed.
Chapter three begins with the definition of a grand ensemble with
overall constraints on the mole numbers of the species. It then derives
a solution theory based on total and direct correlation functions,
between species that yields component properties. This functionality
is important because, although the ensemble-averaged species composition
vector must lie in the component subspace, this is not necessary for
the species composition vectors for any of the systems in the ensemble.
The chapter then gives the component density derivatives of the com-
ponent activities in terms of the reaction extents and the integrals
of the species direct correlation function integrals. Lastly, a de-
tailed analysis of the matrix equations for a system which undergoes
complete dissociations is given, and the matrices are shown to col-
lapse to a very simple form.
Chapter four considers a specific case of complete dissociation--
a system of strong electrolytes dissolved in a non- "reactive" solvent.
First, it reformulates the results of chapter three using dimensionally
smaller matrices, explicitly accounting for the absence of the undis-
sociated salts. In this formulation, it is shown that the long-range
portions of the direct correlation function integrals (which diverge
in integration) are exactly removed by the projections. Isothermal
compressibilities, partial molar volumes and composition derivatives
of the activity coefficient are given in terms of special integrals
of the ion-ion, ion-solvent, and solvent-solvent direct correlation
functions, and the Debye-Hueckel limits of these properties are given.
Two specific cases are examined, the simplest case of one salt and
one solvent and the case of two salts with a common anion. The latter
example gives the Harned coefficients in terms of the correlation func-
tion integrals.
Chapter five deals with the development of a model for the elec-
trolyte solutions. Results of previous theory and experimental results
are used to determine the type of model to be used. The direct cor-
relation function between molecules i and j is given:
hs hs
Cij(r) = C.. (r) + {Cij(r) C.. (r)}
where the terms superscripted hs are the direct correlation functions
for a system of hard spheres. The bracketed term, which corresponds
primarily to electrostatic effects, is approximated by a composite
5
expression, based on a low-density expansion such as the mean spherical
model, with correction to the proper infinite dilution properties in-
cluding hydration of ions. The model is formulated in terms of several
universal (not salt dependent) constants, and functions of ionic size,
charge and density. Some rough, preliminary calculations are given.
Chapter six examines an experimental method by which volumetric
behavior of salt solutions can be determined over a sizeable range of
of pressures and temperatures. The range of operation allows the ef-
fects of total density and composition on free energies to be examined
individually by allowing available activity data at low pressures to
be extended to higher pressures (at constant density). It also pro-
vides a basis for determining the model parameters for the electrolyte
direct correlation functions. The design and operation of this ex-
periment are described.
Finally, further applications of the results of this work, as
well as suggestions for extensions, are given.
CHAPTER 2
SOLUTION THERMODYNAMICS FOR "REACTIVE" COMPONENTS
Introduction
Fundamental to the development of thermodynamic properties is the
mathematical rigor that yields basic relationships such as Maxwell rela-
tions (Van Ness, 1964), identification of partial derivatives (e.g.,
heat capacities), and Legendre Transformations (Beegle et al., 1974).
Provided the requirements of the specification of thermodynamic state
are met, natural phenomena appear to follow all the mathematical conse-
quences of the First and Second Laws and subsequent definitions, even
for systems of many components and modes of interaction with the sur-
roundings (Redlich, 1968).
The most complex systems to deal with are those involving chemical
reactions because the compositions of the species are connected at equi-
librium. In regard to this fact, the purpose of this chapter is to
display a mathematical formalism for expressing the solution thermo-
dynamic properties of mixtures made of components which can react to
form additional species. Their concentrations can be related through
equilibrium constants or there can be complete dissociation. The
particular situations of current interest include electrolytes, both
strong and weak, the "chemical theory" of solutions (Prausnitz, 1969)
and the "solution of groups" method for activity coefficients (Derr and
Deal, 1968, 1969, 1973; Fredenslund et al., 1977a). Our motivation is to
6
provide a degree of unification not heretofore attempted so that deri-
vations can be streamlined and any potential inconsistencies can be
avoided, particularly in the development of molecular and statistical
thermodynamic theories and correlations. The results obtained for the
first two cases are the same as found previously. They have been ex-
tended to a more complete set of thermodynamic properties. In the
solution of groups method we have found precisely how analytical solu-
tion of groups methods such as UNIFAC (Fredenslund et al., 1977b) yield
different results from the molecular equation upon which they are based
(e.g., UNIFAC from UNIQUAC, Abrams and Prausnitz, 1975).
The results of this analysis form the framework for the more
detailed analysis of the statistical mechanics of solutions of "reac-
tive" components, and in particular for the modelling of properties of
solutions of strong electrolytes, which are the subjects of subsequent
chapters.
Mathematical Basis
The mathematical situation is one in which the thermodynamic
functions of interest (enthalpy, entropy, activities, etc.) are most
conveniently expressed in terms of a set of variables related to the
number of moles of species, N, or their mole fractions, x, but actu-
ally must be functions of the set of independent values related to
the number of moles of the original components, N or their mole
0
-o
fractions, x The fact that the number of independent variables in
-O
the former case is different from the number in the latter case re-
quires that the functions be projected from one space to another.
Figure 1 is a graphic representation of the projection process.
Here, we treat the case of three species, with two independent mole
fractions, in a two component system with only one independent mole
fraction. The surface represents the values of the thermodynamic
property, M, for all values of the species fractions satisfying
Zx =1
ii
However, the "allowable" compositions of species are controlled by
the stoichiometry of the "reaction" system. These values are con-
strained to the projection of the species fractions, which here is a
line in the base plane. The attainable values of the system thermo-
dynamic property, Mo, are the projection of M, the curve formed by
the intersection of its surface with the vertical plane containing the
line of "allowed" species fractions. In multicomponent-multispecies
systems the system thermodynamic property values will be found in the
intersection of the property surface in species space and a hypersur-
face for the allowable species fractions. Our purpose is to establish
a general procedure to obtain this intersection. The mathematical op-
erators for this procedure in the present case can be straightforwardly
defined and are easily used as standard matrix manipulations.
This method has been implicitly used in a wide variety of sit-
uations such as electrolyte solutes, the "chemical theory" of solu-
tions and "analytical solution of groups" correlations. While the
1
In certain statistical thermodynamic developments the method is
more complex but represents the only technique we know to remove cer-
tain apparent divergences (Perry et al., 1980a)
O
LL
C)
LL
0
0
0
Z.
O2
0
CxJ
x
!i.- .' ;--.. .. -..... ,.
.. ::.,
.-.-: :.: K2'.3 ,.
..::..... ..t/ ,,
., :
.- .. :.. ... . :, ,r -,.;. .
.. ::. .. '.. . 1 ..', ..
Figure 1. Projection of a thermodynamic property surface.
I
I
I
I
I
xO
0
x
I
/
/
I
I
I
/
/
I
10
number of "new" results is not large, our objective is to unify appar-
ently diverse approaches and present a single procedure for use in all
problems where the properties of a solution are to be obtained from
model expressions for the properties of the species.
Construction of Projectors -- Material Balance
We will begin by examining a closed system at fixed T and P made
from the mixing of n components, with the mole numbers of the com-
o
ponents given by the vector N This system may now undergo successive
-o
associations, dissociations or reactions leading to an equilibrium sys-
tem of n species, with mole numbers given by the vector N. In order to
relate the species vector N to the component vector N we will use a
-o
process which is chosen for mathematical convenience and is not intended
to reflect the chemical evolution of the system. The general develop-
ment is given in Appendix A; two complete examples are given in Appendix
B.
The general result is
N = WN (2-1)
o
where W (defined in Appendix A) is a matrix involving the stoichiomet-
ric coefficients and extents of reaction of the independent reactions
which change components into species.
The importance of this relationship lies in the fact that N has
-O
n degrees of freedom, and since the extents of reaction are functions
of N only (at constant T and P), N also has only no degrees of free-
dom. However, it is in a space with n degrees of freedom. Thus W is
a projector from a space of dimension n into a space of dimension n.
o
11
As shown below, use of this and related matrices allows us to project,
without any loss of information, expressions for thermodynamic proper-
ties into the desired no- space that are naturally written in n-space.
It must be emphasized that the apparently simple relationship of equa-
tion (2-1) contains a considerable number of subleties. The next few
sections show the ramifications that this relation has on the thermody-
namics of solutions. An example involving 2 components, 4 species, and
2 reactions is given in Appendix B.
Properties of Solutions
Partial Molar Properties
We will now see how solution properties, especially partial molar
properties, are projected using the W matrix defined above. First,
assuming our mixture to have reached equilibrium and for T and P to be
constant, we know
dG) = 0 where G = NT (2-2)
T,P
BG
and i =
i T,P,Njfi
Since G is really only a function of the N with T and P constant,
-o
we can use Euler's theorem to also write
T 'G I
G = N with 1 (2-3)
-o--o o o T,P,N
12
Reaction equilibrium can be expressed as
T
v = 0
(2-4)
where v is an n x r matrix with (v) being the stoichiometric coef-
ik
ficient of species i in reaction k and r being the total number of
independent reactions (see Appendix H). Then, after some manipulation
(see Appendix A), we obtain
T
S -o
(2-5)
A consequence of this is the well-known result that for any component
a which is also found as a species in the solution,
(2-6)
i = 1
(2-7)
c
jc
ccT,P,N jfa
Gc
Po a
T,P,N a
o $#cx
Then, since all the partial molar quantities of interest are linear
combinations of P (or P) and its derivatives, all of the other partial
-o
molar properties of components are also unchanged by the reactions.
That is, for a general partial molar quantity (M = volume, energy,
entropy, log of fugacity, etc.) we have:
a = 1,...,n
oM M
(2-8)
13
Equation (2-7) is surprisingly simple since the derivative in n -space
is complex. However, as long as M is a linear combination of G and its
derivatives, equation (2-8) holds because of the linearity of (2-2).
This is valuable in the chemical theory of solutions.
Further, these same linear relations, when applied to equation
(2-4) also yield
v=M = 0 (2-9)
Mixing Properties
For mixing properties, the differences in reference states of the
components and species appear. Looking, for example, at the change in
Gibbs free energy on mixing,
T-mix T To T
AGMIX = NTAGIX = NT NT0 RT NT an (2-10)
while
MIX T -MIX T T o T
AG = N AG N N_ RT N Zna (2-11)
O -O -- -- 0-0 -0 ---
where the relationship between N p and N p depends upon the choice
of component reference state composition. From equation (2-5), equa-
tions (2-10) and (2-11) and the Gibbs-Duhem equation
-MIX = w MIX (o ) (2-12a)
AGO = ) (2-12a)
-0 -
In terms of the initial components the two forms for the chemical poten-
tial are
(p ) p = (po) + R T {(nx ) + (kny) } (2-13a)
-oa oa --o --o oa
= (wT o) + R T {(W nx) + (W Tny) } (2-13b)
Equation (2-13a) is the normal form for iou. Equation (2-13b) is the
general expression for this quantity in terms of the species variables:
it is applicable to all cases regardless of the reference states chosen.
The same relationships exist for all partial molar properties
since they can all be found by linear operations (including differen-
tiation and multiplicative functions) on the chemical potentials. Thus
for
M = f(u.)
i 1
then
-MIX MIX T o) (2-12b)
All =wAMIX f(Pi0- W P
-o --
Figure 2 shows a representation of the projection method for
AMmIX when the reference state composition for component 1 is pure,
x = 1, and "infinite dilution" (or a hypothetical state of pure com-
01
ol
ponent whose fugacity is Henry's constant) for the others. Such a
choice is made with electrolytes and supercritical substances. Here,
it has been assumed that species 2 and 3 are contained only in component
z
0
-J
LU
--J
D
(9
C)
LU
Li
(r
LU
Projection with "infinite dilution" reference states.
C\
0
X
Figure .2.
2, so an ideal solution occurs when o = x = x = 0. The value of
o2 2 3
IIo2 in the reference state is oM which corresponds to no particular
o2 o2
pure species value. It is merely the extrapolation to (hypothetical)
pure component 2 (x02 =1) of the tangent to the line of the projection
at x The line of extrapolation is the ideal solution line which
ol
becomes horizontal in the plot of AMMIX vs. x Note that we have
o o2
terminated the plane at an arbitrary point where the physical state
might change due to a solid phase or a critical point being encoun-
tered. The last term in equation (2-12b) would be the difference be-
tween the line shown and the line of intersection of the plane connec-
o
ting all of the M? = lim M and the (vertical) stoichiometric plane.
1 x.l
i
These two lines need not be parallel.
On the other hand, Figure 3 indicates the results for a case where
the reference states for the components are pure component. Thus
o o= 1 and
ol 1 and x2= 1. This graph of AMIX vs. xo2 is more likely to
show a sign change even though the surface in species space might not
indicate such behavior when examined by itself. Slight undulations of
the surface can yield sigmoidal projections.
Among the many ways of dealing with equations (2-13), several
common ones have been chosen. The method of Prigogine and Defay (1954)
has the same reference state for any initial component which may also
be present as a species in the solution. That is
0 0
o o
Pocx = cx
(2-14a)
CO
LU
LJ
LJ
U-
-r
E
CL
Figure 3. Projection with pure component reference states.
cO
0
>
z
0
k-
--9
0
0r
n-
o To
U = L (2-14b)
-O
This is most useful for solutions where the components are considered
pure, but react when mixed with other components (e.g. solvation).
It is this case, and only this one, for which the most common relation
holds
nnx oy = Jnx ~ (2-14c)
T
A second choice is the case where WT is independent of composi-
tion. It is normally used when complete reactions are assumed. Ex-
amples include strong electrolytes and solutions of groups, which are
discussed in detail later. Here, it is assumed that
Po = T o (2-15a)
o
with
T T
Jnx + ny = W Tnx + W Jny (2-15b)
O
In such a case W is a constant and there can be no successive reactions
so that the upper n xn elements of W are either unity (unreacting
O O
solvents) or either zero or ratios of stoichiometric coefficients (com-
pletely reacting components) and the lower (n-no)xn elements form the
matrix V. It is thus more convenient to reformulate the expressions in
terms of a matrix which has only nonzero rows as discussed below.
Finally, there are cases in which the species standard states bear
no relation to the component standard states. An example of this is an
associating substance which even when "pure" is made up of a number of
species (Figure 3 is an example of such a system). Then the full equa-
tions (2-13) must be used. The technique here is to subtract equation
(2-13b) for the component reference state from the expression at the
solution conditions. For equation (2-13b) this is written as
0 T o -MIX T -MIX
PO p W{ + AG } lim W oP + AG
0O Ox
T_ MIX T-MIX T T o
W AG lim W G +[W lim W ] j (2-16a)
O -
X~ x(
x(a) -(a)
where x0 is the vector of species mole fractions at the reference
(a)
state of component a. In the case of systems where the components are
present as species and the model for AG includes all the contribu-
MIX
tons to AG ,this becomes
S- =LTAGMIX LT lim AGMIX (2-16b)
0 x -
(a)
If the standard state for species is chosen such that (x0). = 1 and
the phase is the same as that of the solution, the last term in equa-
tion (2-14) removes any composition independent terms arising from
the differences in phase or structure in the expression chosen for
AG An example of this procedure is given in Appendix B.
20
Often, the reactions are described in terms of equilibrium con-
stants which are defined as
ZnK T -VT y /RT (2-17a)
which, with equations (2-4) and (2-10), yields
knK = TV na (2-17b)
Equations (2-1) and (2-17a) allow us to rearrange (2-12) to a
relation explicit in the extents of reaction and equilibrium constants.
The general result is given in Appendix A. Equation (A-ll) for one set
of reactions (no sequential reactions) reduces to
T -MIX -MIX T
W AG M = AGI + RTT nK (2-18a)
This relation, or equation (A-ll), can be used in place of equation
(2-12a) in any of the above derivations, including an alternative to
equation (2-16a).
The other mixing properties have similar relationships. For one
set of reactions, we can in general write
M. = f (.)
1 1
where f is a linear operator (including derivatives and multiplicative
functions). Then
T-MIX -MIX T
W AM = AM + RT f(TZnK) (2-18b)
-o
For example, using standard operations on equation (2-15) for
definitions of p. not functions of pressure
1
dinK
T--MIX -7MIX T
WAH = AH + RC (2-18c)
o dl/T
wT MIX = A-M IX (2-18d)
0
T-MIX -MIX T d(TXnK)
SAS = RS __ (2-18e)
o dT
and so on. Thus equation (2-18c) gives the familiar relation for the
heat effect associated with mixing the components expressed in terms
of the activities of the species and the equilibrium constant
MIX T-IX T T 1/T T dTnK
Q = AH N H = N W RN
o -o--o--o 1/TpN -o dl/T
T nyv dinK
= RN T -R NT T (2-19)
D/T.pN -o dl/T
If the temperature variation of the activity coefficients is ignored,
the result is that usually found for reactions in the standard state.
This relation is also appropriate for the chemical theory of solutions.
Excess Properties
Excess properties differ from mixing properties because an addition-
al composition term appears in some excess properties (GE,AE,SE). Using
the usual notation we have for component values:
GE/RT = NT na
0 0
- N Tnx N Tny
0 o 0 0
aTn3
E T T
S /R = NT -o +N Tnx
o -o 3T P,N -o--
-o
For the species values, however, we have
GE/RT = NT na N TnxE NT ny
E T
- S /R = N
Tn
P,N
-o
Use of equation (2-5) on the above relations gives
GE = WTGE RT(knx WT nx) (,o WT o) (2-20a)
-0 -- -0 O
T -MIX o To
= WTG RTinx (p wp ) (2-20b)
For the case of equations (2-14), the activity coefficients become
nnyo = iny + kn(x /xo )
oa a a oa
etc.
+ NTnx
+ N Znx
(2-21a)
For the case of equations (2-15) the relation is
kny = (W Tny) + (W Tnx) enx (2-
oa -c a - a oa
For the case of equation (2-16) the relation is
kny = (L GMIX) (L lim AG~IX) }/RT knx (2-
oa a o a oa
In terms of the equilibrium constants of a single set of reactions,
T-E -E T T
WG =G G + RT InK + RT(knx W Tnx) (2-
o o
or in terms of activity coefficients
W Tny = ny + CT nK + onx W Tnx (2-
-Fin fo -e
Finally, for the entropy,
aTinK
T-E -E T TnK
WS = S + R T
-- -o 3T P,
-o
T
+ R(Znx W Tnx)
o
Of course, for successive reactions, CT nK must be replaced with the
general expression.
Unfortunately, due to the presence of the logarithm, the excess
function expressions do not have the kind of cancellation there is for
the partial molar and mixing properties.
21b)
22)
23a)
23b)
(2-23c)
An apparent simplification of equations (2-20) is normally made
for electrolyte solutions. A set of components is defined for which
T
In x 1 W In x (2-24)
where WT is WT with its elements suitably normalized.
n
Then,
T
In y = W In y
-- -n -
(2-25)
However, this results in a set
which could be used to make up
x # x =
a oaC
This difference requires extra
partially dissociating systems
of components which are not those
the system. That is,
n
N / Nr
oa =1 og
care in expressing the properties of
(weak electrolytes).
Complete Dissociation
An important special case of the above formulae, applying
especially to strong electrolytes and "solutions of groups" methods,
is the case of complete dissociation. Here, we simplify W by elimi-
nating rows of zeros, making it into a matrix V. That is, the reactions
are of the form
B -- > Xv). P.
j j-' :ij
25
The V matrix contains the stoichiometric coefficients of only the
products, and subscripts for reactions have the same index as those
for the corresponding components. Solutes or non-dissociating compo-
nents are considered to undergo the reaction
B > B.
j J
so that they appear in the final expressions.
Our equilibrium equations now become
= v-T (2-26a)
0
T_
M = v^M (2-26b)
0
and the material balance is written
N = v N (2-27)
-O
Now, however, the species vectors do not include any dissociating
components, so are of dimension n n + n where n is the number
o s S
of components that do not undergo dissociation (solvents = n < n ).
In the case of ions and groups, the actual standard states are
left to be defined, so we choose 10 as for equations (2-15)
o To
p = (2-28)
0
with no loss of generality. This leads to a projectable vector of
activities,
vTnia = kna (2-29a)
T T
or v ny = kny + nx v Tnx (2-29b)
o o
As noted above, the "mean" component set can be chosen equation (2-24)
so that for completely dissociated salts equation (2-25) becomes
iny = v iny (2-30)
- n
with mean ionic activity coefficients y+ and mean ionic concentra-
tions for x+ (Harned and Owen, 1958; Robinson and Stokes, 1965).
Application to Analytic Solution of Groups: UNIFAC and UNIQUAC
There are two general situations where group contribution methods
are used in correlation of thermodynamic properties. One is when the
parameters of a molecular correlation are obtained by summing contribu-
tions from the atomic groupings in the molecules. For example, energy
and volume parameters have been calculated in various models from sums
of group energy and volume contributions (Reid et al., 1977). The method
of projection operators as we discuss above does not add anything to
these developments because the groups are not considered to have any
thermodynamic property values. The projection of parameters does not
require the same care as the projection of properties.
The other situation where group contributions are used is in such
cases as the "analytical solution of groups" methods (ASOG) (Derr and
Deal, 1968, 1969, 1973). Here the groups are considered to have activity
coefficients which are summed to obtain molecular activity coefficients.
The present analysis leads to significant results because it allows a
detailed exploration into the relationships between such theories for
groups (species) as UNIFAC and their complementary molecular (component)
theories such as UNIQUAC. Projection operators show how it is that the
two equations do not yield the same composition dependence for system
partial molar properties. As shown below, this difference can yield
complete miscibility using UNIFAC but only partial miscibility for
UNIQUAC when the two binary parameters for the latter are obtained only
from the infinite dilution activity coefficients predicted by UNIFAC.
The UNIQUAC model (Abrams and Prausnitz, 1975; Maurer and Prausnitz,
1978) for mixtures gives an expression for the excess Gibbs free energy
of a system of molecules. In the UNIQUAC Functional-group Activity
Coefficients model (Fredenslund et al., 1975, 1977a,b; Skjold-Jorgenson,
et al., 1979) (UNIFAC), the residual portion of the excess Gibbs free
energy is written in terms of the group contributions using the identi-
cal expressions for the composition relations, but involving projections
of the group parameters.
This expression was chosen by physical arguments after UNIQUAC was
developed and the question can be raised about whether equation (C-3b)
represents a proper projection of (C-4), which would then lead to iden-
tical results for both equations. If not, the activity coefficients
can be matched at only one composition other than pure component.
Since it is preferable to project a complete thermodynamic property
instead of a portion, we have examined the energies, instead of the free
energies, as given in the derivations of the models. This eliminates
the combinatorial term which yields the same results for both equations
and focuses only on the residual terms where the differences arise.
Figure 4 shows how projections can be applied to solutions of
groups. A three dimensional representation is given of the excess
enthalpies, HE, of a solution of hydroxyl, methyl and methylene groups
HE UE U X x. lim U
i x.-i (2-31)
where equation (C-6) has been used with parameters AUki from Skjold-
Jorgensen et al. (1979). Vertical planes representing mixtures of
n-pentane with methanol, ethanol or n-pentanol are shown with the
projections being the heavy lines in the curved surface. The ideal
solution line for these binary mixtures is the straight line in the
plane which connects the values of HE at the two pure component posi-
tions. A portion of the methanol-n-pentane line is shown. The desired
E
value for the system, H is the vertical distance from the ideal solu-
0
E E
tion line to the H surface. As can be seen, the shapes of the H
o
curves will vary because curvature of the HE surface varies with chain
length.
While the plot represents an exact projection, the ASOG method
also relates the group and molecular properties. In fact, because the
same analytic expressions are used, expressions for the partial molar
energies, U and U., look very similar.
on
EXCESS ENTHALPIES FROM UNIFAC FOR n-ALCOHOL-n-PENTANE SYSTEMS AT 298K
n-PENTANOL
ETHANOL
I TE THAriO
n-PENTANE
HE
JOULES
SMOLE
400
200
XCH
Excess entlinlpies from UITPAC for n-alcohol-n-pentane
systems at 298K.
XCH
Figure 4.
However, the test can be made whether the true projection
U = v U (2-32)
exists as required by equation (2-26b). The actual expressions are
given in Appendix C. They show that due-to the complicated interrela-
tions of the indices in what are otherwise very similar summations, it
is clear the molecular energy parameters cannot be related to the group
energy parameters in order to make the two expressions become equivalent
This is because the weighted group surface fractions, Oij, cannot be
projected to give a set of weighted molecular surface fractions, 6,
consistent with the projections of the group surface areas, Q, to the
molecular surface areas q, and N to N. It is possible to achieve
O -
numerical equivalence of the energies at any particular solution compo-
sition besides unity, but the parameter relationships will vary with
concentration. This means that the composition dependence of the com-
ponent activity coefficients calculated from UNIFAC will be different
from those calculated from UNIQUAC even if the values match at any one
point, say infinite dilution. Thus, UNIFAC should be considered as
different from UNIQUAC as, for example, the Wilson equation is.
If the numerical differences are not large, then the procedure
suggested by Fredenslund et al. (1977b) can be used to minimize computa-
tion time by finding UNIQUAC component binary parameters from UNIFAC
infinite dilution activity coefficients and calculating vapor-liquid
equilibrium using UNIQUAC equations.
On the other hand, a particularly significant problem can occur
if the difference is large. Since the slope of the binary UNIQUAC
activity coefficient is always less than that predicted by UNIFAC at
low concentrations, when the infinite dilution values are matched,
situations can arise where UNIQUAC predicts liquid-liquid immiscibility
even though UNIFAC does not. For instance, this happens in the
isopropanol-water system at 1.013 bar (see Figure 5). This problem
can be avoided by fitting the UNIFAC infinite dilution values to obtain
Wilson equation parameters. In fact, this may be the best procedure
in any case for the miscible systems.
1Multiple solutions for UNIQUAC parameters can be found by this
procedure in some systems (Fredenslund, 1980), but we have not found
this in the present case.
0 0 0
Sx i
0Hr
Free energies of mixing for isopropanol(l)-water(2) at 1 atm.
Figure 5.
CHAPTER 3
STATISTICAL THERMODYNAMICS OF "REACTIVE"
COMPONENTS IN THE GRAND ENSEMBLE
Introduction
There are two general methods for relating thermodynamic properties
of solutions to statistical mechanical correlation functions. One arises
from the canonical ensemble, usually assuming pairwise additivity of the
intermolecular potentials, which expresses excess energies or free ener-
gies in terms of differences of integrals involving the species pair
potentials and radial distribution functions (Hill, 1956; Reed and
Gubbins, 1973). The other method, from the grand ensemble, relates con-
centration derivatives of the chemical potential to integrals over the
total correlation function (Kirkwood and Buff, 1951; Hall, 1971;
O'Connell, 1971) and the direct correlation function (O'Connell, 1971).
While modelling studies have been done for both procedures, the latter
has recently been shown to be a powerful technique for solutions con-
taining supercritical compounds dissolved in liquids at elevated pres-
sures (Mathias and O'Connell, 1979, 1980a,b). The correlation of elec-
trolyte solution properties (DeGance and O'Connell, 1975) is another
application of interest and will be treated in subsequent chapters.
Furthermore, the direct correlation functions from the grand ensemble
can be used in the formalism of Chapter 2. The purpose of this chapter
is to apply the projection formalism in such a way as to establish
fundamental statistical mechanical relations in systems with arbitrary
reaction schemes.
Components and Species in Complex Systems
We consider the closed system of Chapter 1, constructed from n
neutral components that may undergo reactions, dissociations and asso-
ciations to arrive at an equilibrium system of n chemically identifiable
or hypothetical species,, neutral or ionic. As before, the material
balance which relates the vector of the particle numbers of the species,
N, to that of the components, N ,is
N = W N (3-1)
-o
where W is as defined previously.
In this more detailed analysis, we must use the structure of W
and related matrices to obtain the desired results. We look here at
the case of no successive reactions, and leave the general case to
Appendix D. In this simple case, we have, as in Chapter 1,
W = L + v (3-2)
1 In this chapter, the vector N will have two meanings. When used
in an obviously macroscopic sense, N is the vector of the particle num-
bers of the species in the system of interest. However, in the statis-
tical mechanical derivations beginning in the next section, N will de-
note the composition vector of a specific system in the ensemble, and
will denote the average of N over all the systems of the ensemble,
and so corresponds to the N used previously. Later, we will return to
using N instead of , which is cumbersome.
The L matrix is of dimension n x n where the first n rows have the
SO
form of an n x n identity matrix and the remaining n-n rows contain
only zeroes. L can then be written as
L = (3-3)
The v matrix is of dimension n x r and can be decomposed in the
following manner. We call the first n rows v (n x r) which contain
0 M-C 0
the stoichiometric coefficients of the components in the reaction
scheme. We next note that the reaction equilibria are defined by
vT = 0 (3-4)
where V again is the vector of species chemical potentials. Since the
n component chemical potentials are independent and equation (3-4) puts
r contraints on p, a given system (fixed p ) has n +r total contraints
-o o
on P. This leaves n-n -r "degrees of freedom" in p.1 That is, n-n -r
-- O
chemical potentials can be varied independently without changing the
values of P and without violating equation (3-4). Of course, the
--o
other r chemical potentials will be functions of the chemical potentials
which are chosen to be "independent".
Following these identifications, we now define matrix v (r x r)
which contains the stoichiometric coefficients of the "dependent"
species, and the matrix v which contains the stoichiometric
-I
1These are not degrees of freedom in a normal sense in that a
fixed T,V,po fixes all of p, but if yo is fixed there may be multiple
solutions to vP = 0.
2The identification of these species is not unique. There will be
some species which cannot be independent in this sense, but in general
there are many ways to divide the independent and dependent species.
coefficients of the "independent" species. Because of the manner of
-1
construction, vD will be nonsingular, so vD will always exist.
Using these identifications, we can rewrite v and W as
-C
v= N (3-5)
--D
-1
and
I + v6
W= D (3-6)
The matrix W has dimensions n x n with rank n provided all the
O O
original components appear in the final system. This means that W maps
N onto an n -dimensional "component subspace" in the n-dimensional
-o O
"species space," and therefore there exists a space of dimension n-n
in the species space which is orthogonal to the component subspace, and
thus to the columns of W. We define Z as an n x(n-n ) matrix made up
of a set of vectors which span this space orthogonal to W. Then
ZTW = 0 and so, ZTN = 0 (3-7)
Of course, any set of vectors which span this space can be used to
construct Z, and in general each element of Z is a function of the
lAnother situation arises in the case of completely dissociated
electrolytes, for example. Then the rank of W can be less than no.
The formulation must be slightly modified so that the above analysis is
still applicable. The details are given below.
37
extent of reaction. For convenience, we will construct a specific Z
consisting of two parts: Z which is n x(n-n -r) and has elements
-0 0
which are independent of the k 's, and Z which is n x r and has
elements which are functions of the 's. Z can be written explicitly
k -
as
Z =
-0
0
-(T -1
-)
I
T
Y
(3-8)
with 0 being n x(n-n -r) and I being (n-n -r) x (n-n -r). It is clear
Sthat
that
zT = T (I + v5) \v v + Iv -v + v = 0
-o- -C -I-DD -D --I I -I -
as desired.
The Z matrix is defined by the following relations:
-1
Z W = 0
-1 _
T
SZ = 0
-1 -o -
T
T
(3-9a)
(3-9b)
Z Z = I (3-9c)
-1 -1
which completely define Z The Z matrix can now be written as
-l 1
z= (z )
--0 -1
(3-10)
which does span the (n-no)-dimensional subspace orthogonal to W as
desired.
We also must construct two n x n matrices, U and Y, which project
in the direction opposite to W. First, U is important since it spans
the same space as W, whereas Y is important because its elements are
independent of the extents of reaction. The defining relations for U are
UW = I (3-11a)
UTz = 0 (3-1lb)
and it has the property that
UN = N (3-12)
-o
We can write Y explicitly as
I
T -1
Y = -( ~1) C (3-13)
0
with I being n x n and 0 being (n-n -r) x n Also
O O O O
YTN = N (3-14a)
but-o
but- -
but
YTZ # 0
(3-14b)
Grand Ensemble for Systems of "Reactive" Components
To establish the statistical thermodynamic relations among the
components and species, we begin by establishing an ensemble with NT
systems, total energy N , and N particles of each species i.
The number of ways, t, of having a given distribution of nNq 's, where
n is the number of systems in quantum state q with particle number
vector N is
NT!
t = (3-15)
H H n N
t=lNq'
q N
In Stirling's approximation,
knt = N TnN N (nNqnn n ) (3-16)
Nq
Using the assumption of equal a priori probabilities, knt is maximized
to find the most probable distribution, subject to the constraints of
constant total size and energy and the overall material balance
equations. These constraints, with their Lagrange multipliers, are
SNq = NT a (3-17a)
Nq -
Sn N E = N -6 (3-17b)
N q
40
n U N = NU T = N N
SNq - T T-o
N q
I I q ZTN = N ZT = 0
N Nq - T-
Nq
where knX and Wne are vectors of dimension no
whose elements are nX and nO k. The equation
d(knt) + a I Idn IENq dnNq +
N q N q
+ (n8))T I IzTNdNq = 0
N q
and n-n respectively and
0
for the maximum in knt is
(enX) T I UTNdq
N q
(3-19)
Equation (3-16) gives, upon differentiation
d(knt) = nnNq dnNq
Nq Nq Nq
N q
which with equation (17) leads to
-ZnnNq + a BENq + (nX) UN + (n) TZ TN = 0
Nq Nq -- - -___ -
(3-20)
since in the method of Lagrange multipliers the nNq are treated as
independent.
Solving equation (3-20) for n we have, for all values of N and q
nNq = exp {a SENq + (n) TU TN + (nO) TZ N}
Nq Nq-- -- -
(3-21)
PRn X
Zne
(3-18a)
(3-18b)
Defining
5E(,6,8,V) = exp {-E + (nm)T UN + (nO) Z N} (3-22)
Nq Nq
the probability of finding a system with particle number vector N in
state q is written as
Nq exp {-EENq + (CnA)TUT + (nen)TZT N}
PNq = -- (3-23)
T ( (A,B3,v)
In the usual case of nonreacting components, the partition function
would be written
(,B,V) = exp {-SE + TN}) (3-24)
NNq
where N is the same vector as in equation (3-22), B is 1/k T and p is
the vector of species chemical potentials. By comparing equations (3-22)
and (3-24) we should expect that nA and nO can be identified by the
relation
By = U knX + Z tno (3-25)
This identification, if correct, has several interesting results.
Using equation (2-5) and equations (3-11)
3p = wT = nA (3-26)
-o
Furthermore, equations (2-2, (2-3) and the Gibbs-Duhem equation
lead to, after some manipulation and use of equations (3-11)
{(nX)TdUT + (nne)TdZT} = 0 (3-27)
provided the differentials are taken along an equilibrium path.
-l
Equation (3-26) identifies -lZnX with the chemical potentials pf
the components. Equation (3-27) will be useful in later manipulations.
The identification in equation (3-25) is only a proposition, and
it must be shown to be a solution. We will assume it is a solution and
then test the validity of this assumption. We begin by examining the
average entropy for a system in the grand ensemble,
~~ = -k B C PNq kPNq~~
Nq -
or, substituting equation (3-23) for pNq in this relation
~~/ = B (nAX)TN + knE (3-28)~~
B -o
For the Helmholtz free energy, A, we have
BA = a ~~/ = (AnX)TN + knE (3-29)~~
kB -o
The chemical potential of component a in component space is
aA ( n_)r
oa 9No N + nX n- (3-30)
Toa aN o =aN ,,-o a ,N
T,V,N o a T,V,N oa T,V,Noy
oy~a i'oy~a ^"oy~a
43
where it is important to note that it is the amounts of the other com-
ponents (Noy) as initially charged to the system before reaction occurs,
and not the amounts of the species, that are held constant. We assume
that a change in a component concentration is done with all reactions at
equilibrium. Thus the species concentrations are known from the com-
ponent values, N and the material balance based on equilibrium extents
-o
of reaction (W). With these restrictions all the derivatives in equa-
tion (3-30) are taken along an equilibrium path, making equation (3-27)
applicable.
Since knX and nO are functions of N equation (3-30) can be
-o
written
a(Rnk) (RnX) aU
a = N + ZnX u ( nk)
oa N -o a 8N -N <-
T,V,N T,V,N TVN
y^ a 00y
(Zne))T azT
Z T (nA)T (3-31)
oa oa
o T,V,N oa T,V,N
The first and third terms on the right hand side of equation (3-31)
cancel by equation (3-12), the fourth and sixth are zero by equation
(3-27) and the fifth term is zero by equation (3-7). This leaves us
with the desired result
gpJ = S.nX
oa a
or
(3-32)
p = nX
-o -
as proposed in equation (3-25) and (3-26). Thus we find that these
equations are indeed a solution for the Lagrange multipliers ZnX and
knO.
Equation (3-32) not only defines the Lagrange multiplier ZnX
but also shows us what the independent variables for the grand canon-
ical ensemble are. We already know that a macroscopic (canonical) sys-
tem can be determined by fixing T,V and N a total of n +2 variables.
-o 0
Similarly, in the grand ensemble, we fix the same number of variables
T, V and Vo (or ZnX) to completely determine the thermodynamic state.
This fact is very important in taking derivatives of ensemble averages,
since this means there is a relationship between the elements of tne
and the elements of nnX. This dependence and important thermodynamic
derivatives of the ensemble are given in Appendix E.
Composition Fluctuations
Macroscopic thermodynamic properties of the ensemble can be re-
lated to microscopic properties through the density fluctuations in the
systems of the ensemble. We define the fluctuation matrix, A, by
A- ( ) (3-33)
1T
with = N exp {-1E + NT(U Zn + Z enO)} (3-34a)
N q
= B (3-34b)
1
= NN exp {-BE + N (U nX + Z Zne)} (3-34c)
NNq
= N qNq
Equation (3-34a) can be differentiated with respect to a component
chemical potential along the equilibrium path to obtain the relationship
between the fluctuations and macroscopic properties (details are given
in Appendix E). This differentiation leads to
1 D = N1 T I +
(Y) + Z
8By a -oai
o T,V, oa T,V,p
oyfa oyja
1 T
(Y) (3-35a)
a
= A (Y) + Z
a -o ap
T,V,y
oyfa
where the notation (Y)
--
We can use equation
equation into component
(A) =(A) <
1 T T
I- (Y)
-B -
denotes the a-th column of Y.
(3-14a) to project the left hand side of this
space. Defining the matrix A
ocN
40
1 T a
(Y)
-B B3 o I
T,V, N oa T,V,p
oy#a oyfa
1'(Y) +z
I a polT,V,p 1Ja
S(Y) (Y) (3-36)
Using equation (3-7) we can rewrite A as a projection of A,
() = (Y) A
=(Y) +
(Y) + Y)
ST
-- Z A (Y)
4u -o a
The next section shows how the A matrix is related to the
microscopic properties of the species, thus forming a bridge between
them and the macroscopic component properties.
(3-35b)
(3-37)
Correlation Function Integrals
The method of Kirkwood and Buff (1951) can be used to relate the
fluctuation properties of integrals of the radial distribution function
and direct correlation function. We define specially and orientation-
ally dependent one and two particle densities by
N.
(1)r = ( ) 6(w p ) i, = 1,...,n (3-38a)
i --1 -1 p -1p
p=l
N. N.
1 1
f(2) !2 = I ) p q( p-)
fij2)(r r ; 6(r -r ) 6(r -2) 6(
ij -1,-2 2q=1 1 1 q 2 -p 1
p=l q=1
(p#q if i=j)
x 6(w -q2) (3-38b)
-q -2
where 6(r-r ) is the three dimensional Dirac delta function, r is the
position vector and w is the vector of Euler angles of molecular orien-
tation.
By definition, we have
f(l)(r ;w ) dr d, = Ni (3-39a)
UJV0 i 1-1 1 11
J. VV32) rlr2 2) r dr dw dw = N N N 6. (3-39b)
SV2 ij1- 1 2 1 -2 i j i ij
We now define ensemble average number densities by
(1) (1)
P -(rl;) (3-40a)
S -1 1 i -1 1
(2) (2)
Pij (r r ; 2) ij
Equations (3-39) are integrated to obtain
S ( r ;o) dr dw =
I J VJ Pi (r i;2 1,2) d d d d = N -
VV P (r, ; '2) dr dr2d d = ..
ij 1 'E2; 15 -1-2 1 2 i 1 13i
(3-40b)
(3-41a)
(3-41b)
For fluid systems, we define the mean density and the pair correlation
function by
Pi p (r ;) = /V
2 (2)
PiPjij rl2;l' 2) Pij (r2;l1' 2
(3-42)
(3-43)
where I=E dw.
The pair correlation function can be integrated over its arguments
and related to the concentration fluctuations using equation (3-39b), as
Pi
6 = -1a V g(r ;g ,rW ) dr dw dw (3-44a)
1 3 1 1 R2 fV 1ij 121 -12 -2
6..
1 J 1 13
= Pi pjV 2dr1
= 11, 2 -12
The angle-averaged (or "centers") total correlation function, h..(r ),
is defined by
is defined by
hJ ij ) = l, 2-1
(3-45)
Then using matrix notation
1 3
(3-46)
kJ (H).
2 ij
Sij 1 J<
S 6..
1J
(3-47)
letting x. = /, equations (3-33) and (3-47) yield
1 1
(A).. = x.6.. + x.x.(H).
-1i] ij 1 ij
or, defining a matrix (X).. = x.6..
-13 1 1J
A = X(I + H X)
(3-44)
we find
1'
=f h, )dr N
()ij V V (E12 -12=
(3-48)
A "centers" direct correlation function can be defined by the
equation
c ( 12) E hij 2 -
Integrating over r12
(3-49)
n
S1 cI j (r13)P hi ( r23)dr3
and multiplying by /V we find
fhi .
V jd -12
f ,dr
V < IN + -V2
fV Vc i(r 13)Ph
V )p h (r23)
(3-50)
x dr dr1
-3 -12
For systems where the matrices A, H and C are all bounded and non-
singular Fourier transform theory (Pearson and Rushbrooke, 1957) can
be used to obtain
H = C + H X C = (I + H X) C
(3-51)
where
= I )d
(C) V c.. (r )dr
ij V 1i -12 -12
Some manipulation of equations (3-48) and (3-51) leads to the relation-
ship between C and A,
iHere it has tacitly been assumed that the integral of cij(r;p) is
not divergent. This is not the case for certain systems, such as those
with coulombic interactions. In such a case, the derivation of
Chapter 4 must be used.
A(X-1 C) = I (3-52)
Since we have only the components as degrees of freedom, we want to
project this into component space. Following equation (3-37) we pre-
multiply by
T
T aI T
(Y) + Z
-a ap -o
ST,V, 1-
and post-multiply by W, giving
T __T T -1 T
(Y) + Z A(X C)W = (Y) W
a 2o -o
oa
T
T,V,poy
o oya
or
T
ST,V,p oy#
"oyfO
where e is a vector such that (e ) =6 a Differentiation of
equation (3-1) leads to
T 0 VoW + N (3-54)
oa aB oa -o YBaoa
TVdc p TV,0(- TVOa y #
TVAoy#a oya oyA
We define the matrix K such that
-1p
no a(W)
(K )i N iy OY
y=l oa V, oy
oy~a
(3-55)
Then we can write the general expression in component space by com-
bining equations (3-37), (3-53b), (3-54) and (3-55)
AT(X-1-C)W + K (X-IC)W = I
V- -1 _
(3-56)
-1
where I is the n xn identity matrix. Multiplying through by A we
have the desired result
A-1 = (W + K)T(X-1-C)W
where
o a(W).
(K)ia = Noy No
y=l Nct
T,V,N
oyfa
n n
o o
= N Y
y=l 8=l
oB (W).
oa w)p 's
o T,V,N o T ,
oyoYc
n
0 ~-l l(Klli
=l
S(-1K T
- -1 ai
Equation (3-57), therefore, relates the derivatives of the component
chemical potentials in component space to the species direct correlation
(3-58)
(3-57)
53
functions and the component composition derivatives of the extents of
reactions (through W).
Because of its complexity, it is not apparent that a single optimal
procedure exists for obtaining changes in chemical potentials (and also
pressure) by integrating equation (3-57). However, in the case of to-
tally dissociating systems this is straightforward.
Complete Reactions
One kind of system of particular interest is that which, for each
reaction, there is a limiting reactant which appears only in that re-
action and whose concentration goes to zero in the solution. In such
a case, the species composition vector and the extents of reaction can be
uniquely determined from v and N using only algebraic relations (no
-- -o
equilibrium calculations). Examples of this kind of system include
solutions of strong electrolytes and "solutions of groups."
Here we make a different partitioning of the matrices of interest.
The components are divided into two groups: the limiting reactants which
are subscripted L, and the other components which are subscripted S
(solvents). The species which are not found as initial components form
only one group, which is subscripted P (products). The W matrix is then
I 0 v
-S
W= o0 + v (-s -) (3-59)
0 0 YVp
-P
Since a system of only complete reactions with the above restrictions
must have the same number of reactions as it does limiting reactants, we
will number each reaction by its corresponding limiting reactant. If
n is the number of solvents, the partitions of the W matrix have the
S
following dimensions:
1 In the case of where two components of a reaction would go to
zero concentration, one is chosen as the limiting reactant and the other
treated as an ordinary component. This situation could only occur if
these two components had initial compositions with the proper stoichio-
metric ratio. In the derivation given, only non-sequential reaction
schemes are considered. This is general because any sequence of com-
plete reactions can be written as a non-sequential scheme.
55
v is n x (n -n )
vp is (n-n ) x (n -n )
2 is (n-n ) x (n -n )
S is (no-ns) x ns
L is (no-ns) x (no-ns)
To see the simplifications that result in the case of complete
reactions, equation (3-57) can be written
A-1 = (W+KT (XI-C)(W+K) (3-60)
This equation may be obtained from equation (3-52) by following a pro-
cess identical to equations (3-53) through (3-58) plus postmultiplying
equation (3-52) by K to obtain
T -1 T aT
A(W+KT (X-1C)K = Y K+ ZTK = 0 (3-61)
yoa TV -o-
T,V
The second equality comes from
YK = 0 and ZK = 0
-.- -0- -
and is the Jacobian matrix of the derivatives
T,V T,V,po
Thus,
(W+K) (X-1-C)K = 0 (3-62)
Addition of equations (3-58) and (3-62) gives equation (3-60) as desired.
The projection matrix W + K is a very sparse matrix compared to either
W or K.
For a complete reaction, written
n
SB.v. = 0
B 0t
i=l1
with component a being the limiting reactant, we know that
N = 0 = N + )v E N (3-63a)
a oa aa a or (3-63a)
a
where r is the reference component for reaction a. In the limit of
complete reaction, the extent of reaction can be written
N
1 oa
S= 1 N (3-63b)
aa or
The general term in the W matrix is
n
o N
() =6 + (-- 1 oN) (3-64)
ia v N Br
a=n +1 aa or a
s a
Similarly, for the K matrix
n n 6 N N 6
o o v. a or oa Br
(K).- = N ( la) a a
n.oy L v yr 2
y=l a=n +1 aa a N
s or
a
n n
o v. o v. N
la la oa 6
6 + v N- (3-65)
a=n +1 aa a=n +1 aa or
S S a
Addition of equations (3-64) and (3-65) yields the general element of
the projection matrix
(W+K) = 6
iBs 1i
n
0 v).
v aB
a=n +1 aa
s
If index B is between 1 and n (solvent), we have simply
(W+K) i 6 =
1 : :n ;1 a n
1~B n; i'5
which defines the first n columns of (W+K). If 8 is between n +1
S - S
and n (limiting reactant) we have three possibilities: For i being a
solvent (1 i 5 n ),
(W+K)
iB
1 i < n ; n +l f B < n.
s s
(3-68a)
For i being a limiting reactant (ns+1 s i S no)
S 0
(W+K) = 0
is
(W+K)
(W+K) = 1 0
--iB v B
n +1 < i,B < n
s o
n +1 < i,B n
s o
(3-68b)
(3-68c)
Finally,for i being a product (n +1 5 i n)
0
(W + K)
iB v g
n +1 < i < n; n +1 < 8 < n
o s 0
Since vL (made up solely of the V 8's) is diagonal, equations (3-65) and
(3-69) can be written in matrix form.
-I
I -v v
-S-L
,W+K 0 0 (3-70)
-1
0 -v 1V
PL
(3-66)
(3-67)
(3-69)
Vi
Partitioning (X-1-C) in a similar manner yields
-1
X -C =
-1 T T
X -C -C -C
-S -SS -LS -PS
T
-C -C -CpL
LS -LL PL
^PS 4PL X _PP
-PS -PL -PP
0 0
0
0 -L
0 0
(3-71)
The first matrix on the right hand side is well-behaved in the limit of
complete reactions. However, the second matrix diverges and must be
operated on by (W+K) before the limit of complete reactions is taken.
This procedure is carried out in Appendix F and leads to the equation
0
limit
complete (W+K) 0
reactions
0
0
0 (W+K) = 0
0
Equations (3-60), (3-70), (3-71) and (3-72) lead to the following
equation for A-1 in the complete reaction limit.
-1
A
-1 -1 T -1
X -C {-(X -C )v +C v }v
-S -SS S SS -S PS-P -L
-1 T T -1 T
(V ) {v (X -CssC )v -S pC
-1 T T -1 T T T T -1 -1
(v ) {-v(X -C ) + v c} -v C + v (Xp -C )v } _
-L -S -S -SS -PPP p _PS-P -P -P -PP -P -L
(3-73)
Although equation (3-73) appears somewhat complicated, it contains
some important results and has two important limiting cases. Most im-
portant of these is that there are no contributions to A-1 from direct
(3-72)
correlation functions involving limiting reactants. Also, it can be
-l
seen that A- is symmetric as required.
If we look at a system with only solvation or association of com-
ponents, then all reactions are of the form
SV aLB + v B L B
L VBa VLLBL P
a=l
where the B 's are solvents, B is the limiting reactant and Bp is the
a L P
solvated product. Thev matrix is the identity matrix and
-P
-1
X -C
-S -SS
-1 T T -1
( ) (xs -C) +c PS
L -S S -SS PS
-1 T -1
{-(X -C )v + C }
S SS -T S -PS -L
-1 T T -1
L _S XS SS S -PS)
T T -1 -1
-v c + ( -c )}1
S-S -PP L
Equation (3-74) simplifies still further for systems with only either
solvation or association. For the former, v is the identity and
-L
-1
A =
-1
X -C
S -SS
T -1
-v (x -C ) + C
iS -6 SS PS
For the latter case, v = 0 and
-S
-1
X -C
S -SS
-1
A
-1 T
(v_ ) C
L -PS
S S SS S PS-S
TT -1
-_V sCps + (Xp Cpp)
-S PS _P -PP
T -1
C v
PS-L
-1T (X- -1
L _P PP -L
-1I
A=
(3-74)
(3-75)
(3-76)
The other case of interest, and the case to be pursued more fully
in later chapters, is that of complete dissociations, that is, reactions
of the form
n
0
B + v B
B L kLBk
k=n +1
o
where BL is the limiting reactant (dissociating component) and the
Bk's are product species. The solvents do not take part in the re-
actions (v =O). In this case, v is the identity matrix and
-S -L
-1 T
X -C C vp
S -SS -PS-P
-1
A =-- (3-77)
T T -1 -
-P-PS P P -PP P
These results are quite simple to use, in fact evaluation of
activities requires only simply quadratures, not the solution of partial
differential equations as in the case of partial reactions.
CHAPTER 4
SOLUTIONS OF STRONG ELECTROLYTES
Introduction
O'Connell and DeGance (1975) give a formalism for expressing the
fluctuation properties of electrolyte solutions (concentration deriv-
atives of the chemical potential, partial molar volumes and compressi-
bility) in terms of integrals of the radial distribution functions
(Kirkwood-Buff solution theory, Kirkwood and Buff, 1951; O'Connell,1971),
and of nondivergent integrals of the direct correlation function. The
formalism was based on a constraint on these radial distribution func-
tions (Friedman and Ramanathan, 1970) which is not necessary to arrive
at the desired results in the case of the direct correlation function
integrals, but significantly affects results in the case of the radial
distribution function integrals. The previous formalism was both rigor-
ous and practical for properties of multicomponent solutions expressed
in terms of radial distribution functions (assuming the constraint to be
valid), but for the multisalt case yielded intractable equations when
expressed in terms of direct correlation functions. This is primarily
because the multicomponent charge neutrality conditions were mishandled.
The method used here parallels that of the previous chapter and
yields equivalent results to those of O'Connell and DeGance for the sin-
gle salt case, and quite simple results for multisalt systems. In addi-
tion, this method does not require that the constraints on the radial
distribution function integrals be valid, although this situation is
addressed further in Appendix G.
Material Balance Relationships among Salts, Solvents and Ions
In the case of strong (completely dissociating) electrolytes a more
compact, but equivalent, set of matrices are used in the material balan-
ces and other relations than were used in the general formalism. This
set is used because of the sparseness of the W matrix and the fact that,
with the salts as reference components (which is so in this chapter),
the K matrix is a zero matrix. The new matrices are defined below.
We look at a system of strong electrolytes dissolved in n solvents,
where the total number of salts plus solvents is n0 and the total number
of ions plus solvents is n (n > n if independent salts are used see
Appendix H). The n components will be denoted by greek subscripts (a,
s,...) and the n species by i,j,k,....
We construct a matrix V defined:
for solvents
i = 1,...,s
(. = 6. 4-la)
ia la a = 1,... ,s
and for electrolytes
the number of ions i = s+l,...,n (4-1b)
vi = of type i in a molecule
S of salt a. c = s+l,...,n
with all other elements being zeroes. In the notation of Chapter 3, the
v matrix can be written
I 0
S= 0 v (4-2)
-I
Such a system has n-n constraints on its composition in terms of ions
and solvents, and they may be represented by an n by n-n matrix z
O
63
having the property
T
v z =0 (4-3a)
which gives
0
-1 T T
z= -(D )-I (4-3b)
I
The elements of z include the ionic charges as well as any additional
stoichiometric constraints. As these elements are not functions of com-
position, the z matrix corresponds to the Z matrix defined previously.
In order to perform the desired projections, two more nonsquare
matrices, u and y, are defined such that
T T
Tu = u z 0 (4-4)
T T
v = y z 0 (4-5)
where I is an n xn identity matrix; (I).. = 6.. i,j = 1,...,n
n
The u matrix spans the same no-dimensional subspace of R as V and is
orthogonal to z also. The form of u will be derived below. The y
matrix is a very simple matrix, but is not orthogonal to z. It has
the form
I 0
y = 0 -( T (4-6)
DD
0 0
and from equations (4-5) and (4-6) we know that
u = + zb
(4-7)
where b is an (n-n )xno coefficient matrix, the details of which are
not important.
Finally, we define the nxn matrix
T T T
P = vu uv = PI
(4-8)
where P is the orthogonal
ting perpendicular to the
property of P is that for
the vectors of z, i.e. Mz
MP = M
Equations (4-8) and (4-9)
Pu = u
PN) = V
Pz = 0
and with equation (4-7)
projector for the component subspace, projec-
kernel and parallel to the image of v. A
any matrix M with a kernel that contains all
= 0, then
(4-9)
also imply that
(4-10a)
(4-10b)
(4-10c)
(4-10d)
PY = y (4-11)
All of these properties of P will be used below.
The necessary matrices are easily constructed from the known v
matrix. Since u and V span the same subspace of the n-space, each
column of u must be a linear combination of the columns of v. In all,
2
the construction of u from v requires n coefficients which make up
the n xn matrix, c, such that
o O
(4-12)
u = Vc
Using this identification in equation (4-4) gives
T T
V u = V Vc = I
or (VT V)-1
orc = (V T)
This identification leads to
u = v( v)-1
I
u= 0
0
U
T T -1
2D (--DD + IVI)
T T -1
I D4D II)
S= v(vT )-1 T
T T -IT
o0 v(v +v )
--D -D---D --I-I --D
T T -1 T
0 I -D-D -I-I -D
0
T T -1 DT
v (VV V ) V
_I -D-D -I--I -D
T T -1 T
--I --D II -I
The material balance relationship between the species and compo-
nents is
N = vN
-- ---O
(4-16a)
where N is the vector of moles of solvents and salts (components).
--o
I (v v)-I = (v)- iff v is square (not possible here).
(4-13a)
(4-13b)
(4-14a)
and
(4-14b)
(4-15a)
(4-15b)
The
corresponding relations from species to components are
N = uN = yN (4-16b)
-o -
As before, the matrix expression relating species composition
fluctuations and radial distribution function integrals is
A = X + XHX (4-17)
where the matrices A,X and H have been defined previously. The species
composition fluctuation matrix restricted to n degrees of freedom, A,
is given by
T
(A) = ( + T A (4-18a)
alloC T,V, p
where the elements of A are given by
8N
kT oa kT DN 0
-) aB- N N (A)
T,V,i ]oyft T,V, poyf
a,8 = 1,...,n (4-18b)
Equations (4-18) can be derived in a manner analogous to equation (3-37)
using the methods of that chapter. Here, however, the matrices are
slightly different (as shown above) and we explicitly define the salt
chemical potentials by
o = v 1 (4-19)
where p is the vector of chemical potentials of solvents and ions. It
67
should be noted that although this is somewhat different from equations
(2-4) and (2-5), the procedure followed in Chapter 3 still applies.
As in O'Connell and DeGance (1975), A is nonsingular and can be
inverted to give
-1 N ooa
(A 1) kT
a3B kT aN
T,V,N
08#~
N oa3
kT aN
oCa
T,V,N
(-1)
The radial distribution function integrals are related to A by
T
z
V
(X + XHX)y
(4-21)
This equation is much less formidable than it seems, since
T + T
S T,V
T
D-11I
DI- 0
T,V
(X + XHX)z = 0
(E-11)
= -YIT(X + XHX)z{zT(X + XHX)z}-1
This allows us to rearrange equation (4-21) tol
A = yT{I (X + XHX)z[zT + XHX)z]-Iz(X + XHX)y
A = T (I (X + XHX)[Tz (X + XHX)z] zlT + XHX)X
(4-22)
-1
Of course, A is just the inverse of the right hand side of (4-23).
See Appendix G for certain simplifications.
(4-23)
(4-20)
The set of equations that equation (4-23) represents can be solved
by simple quadrature, as opposed to the more complex method required to
solve the system of partial differential equations given in equation
(3-57). However, even these quadratures are involved due to the mul-
tiple matrix operations (especially inversions) which need performed at
every quadrature point. This difficulty can be alleviated if the dimen-
sionality of the matrices are small enough to allow doing the matrix
manipulations analytically. An example of the analytical method is
given below. Thus, as before, the direct correlation approach should
be of greater practical value. It has been shown (Gubbins and O'Connell,
1974; Brelvi and O'Connell,1973,1975) that direct correlation function
integrals are insensitive to angle-dependent forces and yield excellent
correlation and prediction of the properties of gases dissolved in
various liquids, including water (Mathias and O'Connell, 1980a,b).
Formulation in Terms of Direct Correlation Functions
O'Connell and DeGance (1975) showed how expressions for single
salts which were simple functions of nondivergent direct correlation
function integrals could be attained. However, due to incomplete
expressions for charge neutrality and composition constraints, this
simplicity disappeared in the multisalt systems. An ancillary effect
of this oversight was that some terms were stated to be included, but,
in fact, are not. The following results are considerably simpler.
On pages 767 and 777 of O'Connell and DeGance (1975), it is noted
that all terms of order (M-1) and less in the divergent integrals are
preserved. In fact, only the first order terms remain.
Recognizing that the direct correlation function integrals of ions
are divergent, we return to the definitions of the direct correlation
functions themselves. For species of type i and j
n
cij(r2) hi(r12) 1 hik(r23)Pkckj(r3 )dr3 (4-24)
In order to perform the transformations necessary to attain the desired
matrix relations, both c..(r ) and h..(r ) must be integrable. How-
13 12 13 12
ever, c..(r 2) has a long range limit (O'Connell and DeGance, 1975) of
2
z.z.e
lim c..(r 2) = (4-25)
13 12 EkTr12
r-o 12
which has a divergent integral over volume. Here, e is the charge on
an electron and E,the bulk dielectric constant of the solvent mixture.
The presence of the z.z. factor in equation (4-25) causes these portions
of the direct correlation functions to lie in the kernel of the projec-
tion operator, P. Therefore, these divergences can be removed in the
process of projecting the species composition fluctuation matrix. We
now write c. (r) as the sum of a nondivergent part, c. (r), and the
divergent, long-range limit, giving
2
z.z.e
c.(r) = c. (r) + (4-26)
13 1J EkTr
Substituting equation (4-26) into (4-24) and operating from the right
with P gives
2
n n z.e n
Sco.(r 2)(P)jk + kr zj(P) h (r )(P)k
j1 j= kTr jk I j
j= j= c j=1
n ze
E kTr
j=1
n f
- ?
=1 '
z.(P)ij
dr
--3
(4-27)
Since the vector of species charges (the z.'s) is in the space spanned
J
by z, we have by equation (4-10d)
n
x z (P)k = 0
j=1
Use of this relation in equation (4-27) simplifies it to
c (r )P)
ijr2 jk
n
j=l
(4-28)
n n r
- (r 23)
j=1 =1
o
x co (r )dr (P)
cj 13 3 jk
or, after Fourier transforming,
COP = HP HXCP
where
S NV 1f c
(4-29a)
(4-29b)
(4-29c)
Pre-multiplying equation (4-29b) by X and rearranging gives, with
equation (4-17)
XHP = ACP
Adding P to both sides and using AX-IP = P + XHP, we find
Adding P to both sides and using AX P = P + XHP, we find
A(X-1 C)P = P
(4-30)
which is nearly identical to equation (3-52).
c (r3) (P)jk
Pj 13 -jk
hi (r23) P
ij (r12) (P)jk
71
We project to component variables by first pre-multiplying equation
(4-30) by
T
T + E
T,V
T,V
equations (4-10d) and (4-11).
(xT + z T A(X2
-o
T,V
Equation (E-11) makes it clear that
T
Tz the right hand side simplifying via
*
- C )P = uT
T
Z
+
T,V
zT A
Tj
(4-31)
has the
required characteristics of matrix M in equation (4-9). Thus we know
(4-32)
T,V A
T,V
Substituting this into equation (4-31) and using (4-8) gives
zT AuvT -1
AVz (X
,V
- C)P = uT
(4-33)
Post-multiplying by v and successively using equations (4-10a),(4-11)
and (4-18), equation (4-33) becomes
AvT(X C ) = I
A-1 = T(-1 C0)
(4-34a)
(4-34b)
T
z
V
Equation (4-34b) provides a very straightforward relationship
-1
between the desired macroscopic salt properties (elements of AI ) and
the molecular direct correlation function integrals for the ions and
-1
solvents. Note that both v and X- are both simple and well-known
matrices and that all solvent-solvent relationships are left totally
unaffected by these transformations. This result is identical to
equation (3-77), as required.
Thermodynamic Properties
The general matrix relationships between measurable thermodynamic
properties and the correlation function integrals are summarized here.
We choose to work with the mean ionic salt properties, such as the mean
ioni.c activity coefficient on the mole fraction scale,
I n n
n =y 1 v. n Y kT = I V ~~i kT in Ni/N)
V la 1 v kT = 1- \i-N i ( -5
a =s+l a i=s+
P 1j a n 0
V n y, n v N N (4-35)
-a t kT
i=s+l y=l
with
n
v = i ia a=s+l,...,no
i=l
Here we have tacitly assumed that p = vp which are the usual standard
-o -
yields (Friedman and Ramanathan, 1970)
Sinl By 8" n v. v.
nny N oa o~l 1
N -oI + v v (4-36)
a NN kT TN +
08 T,P,N ,0 T,P,N i=s+1 i
oy#B oy#B
where,
x. = N./N
1 1
If equation (4-36) is restricted to salts of two different ions (+ and
-) and to systems without ions common to different salts, it becomes
any
+a
Nv
a 8N
3 T,P,N
oy#f
N Poc
kT N
oB
T,P,N
The partial molar volumes are, for solvents
- v
voi 38N
oe#i
o
o0
1oi
DP
T,N
-o
-o anYoi
=v + RT
T,N_
-0
SN
6 + v v
N c6 a
oa
RT Zn x iY .)
DP
i=1,...,s
B=1,...,n
O
(4-37)
T,N
(4-38a)
(4-38a)
and for salts
v-
oa T,P,N
0 13oac
pct
P
T,N
'-o
o
(IJ + v RT kn x y )
oc a a+a'
P
T,N
-o
-o
v + v RT
oa ac
Uny.
T,N
-0
C=s-1,. ,no
B=l,. ..,n
(4-38b)
The partial molar volumes are related to the A matrix by the
relation
v
oc -1 T
SRT = (A u Xi)
K RT
T
(4-39)
PP/RT
Pp
oa
T,PoBio
74
where po N o/V, K is the isothermal compressibility and i is a
column vector of n ones. Also,
aP/RT
'-P T,N
1 T -1 T
K- = iX uA u Xi
pKTRT
(4-40)
where p is the
Further,
molar density for the ions plus solvents, p N/V.
G = A-1 1X + VT1 (4-41)
G=A -V-v~vvV(-1
where
Zny~a
(G) = Nv -aN
S T,P,N
oy/8
pv v
(y) ooa o
aGB KTRT
= (G) a a, =l,...,n
a,8=l,...,n
o
and 1 is an nxn matrix of ones.
All of the desired properties can then be obtained from correla-
tion function integrals using equations (4-34) and (4-39) to (4-41).
The expressions in terms of the (H).. are complex (O'Connell and De-
Gance,1975; Friedman and Ramanathan, 1970) whereas they are simple when
the (Co).. are used.
lj
1n n
= 1 i TXC Xi = [ ) x.x. [-(Co)
PK RT X-1 C ij
TRT i=1 j=1 i
v n n
- R V (I-C X)i = v xi[-(C)
T i=1 i j=1
(4-44)
=l,. . ,n
o
(4-45)
1 n n n n 0
PKTRT (G)aB L= VX x k ikiB Xk-( [1-(C)ij
T i=1 j=1 k=l R=i1
(4-42)
(4-43)
(4-46)
These results are in the same form as for nonelectrolytes (Mathias and
O'Connell, 1980a,b). Equations (4-45) and (4-46) are much less formi-
dable than they appear due to the sparseness of v and cancellations.
For salts of two different ions, the double summation over i and j leads
to, at most, four terms.
Since they can be modeled well in nonelectrolyte solutions, we also
believe it is useful to work with fluctuation matrices in terms of den-
sity rather than moles. We know that for both salts and solvents,
kT p T
N oat
kT 9N
T,V,N
0-yfQ
(A-1 -
-ciI -ci
Defining an activity
density,
coefficient, y, based on concentration or molar
0
u 0 + RT kn p y
oCoa oa oca oa
oct cT
= u + v RT n p y.
*oa + '+a-.
for solvents
(4-48)
for salts
We can write
-1 T -1
G' =A vX v
ninyot
(C') = p n
nT,p
and
= (') (VTC ) c
Bct _
Po = No/V
(4-47)
where
(4-49)
(4-50a)
(4-50b)
[1 (,ikl -(, jz 1
Note that for a= salt, we write V ony instead of ny Previous
success (Mathias and O'Connell, 1980a,b) in modeling (C ).. for non-
electrolytes may indicate that the same may be done for electrolyte
systems. In this case, the expressions should include short-range,
infinite dilution and Debye-Hueckel contributions. However, these
are ion-specific, not salt-specific.
Single-Salt, Single-Solvent System
The simplest example of salt solutions yields results obtained
by less elaborate methods. However, for illustrative purposes, we will
use the above method to derive these expressions. In this system, the
solvent is indexed 1 and the salt, with ions + and -, is indexed 2.
Since the choice is arbitrary, we chose + as the dependent and as the
independent ions. This choice gives, for equations (4-2) and (4-6)
1 0 1 0
S= 0 V y= 0 1/v+
0 v 0 0
The z matrix is, from equation (4-3b),
0
z= -V_/V
1
If we multiply z by z it becomes
0 0
zz = -v z /v = z
z z
(4-51a)
(4-51b)
which shows that z is indeed parallel to the vector of ionic charges
as required in equation (4-27).
First, we give matrix A in terms of (H)... Equation (4-23) gives
Sx+x {(H) -(H) }
L(H Xo0 + J 1+ 1- -
(A)ii = + x2 (H) xJ()- (4-52a)
-11 01 01 11 VXo2 + x+x ((H) ++-2(H) +(H) }
A)12= ()21 = xolo2 x(H)1+ + x_(H)1 + x+x_(H) {(H) __ (H)_}
+ x+x (H) I{(H)H) + H)+}- [xo2 + x x {(H) -2(H) +(H)
(4-52b)
(A2 = x o1 + x+(H) + x(H) + x (H) (H) (H) -
S+ xx (H) -2(H) +(H) } (4-52c)
where
xol = Nol/N = N1/N = xI (4-53a)
Xo2 = N/N = N =N x+/v+ = x /v (4-53b)
and = + + v_ (4-54)
Interchanging the + and subscripts does not affect the results. This
must be true since the choice of dependent and independent ion was
-I
arbitrary. The A elements can be determined by inversion.
For the direct correlation function integrals, the elements of
-1
A can be determined directly from equation (4-34b). From equation
(4-20), we now find
8N
ol
T,V,No2
o2
= 1 (C)
ol
2
x2 1 + x (H)+ + x_(H)
denom.
T,V,Nol
ol
N __02
k T DN ,
T,V,N
o2
= -{+(c)+ + (c)
1- 12
= olxo2{x+(H)++ x_(H)_+ x+x_[(H) -(H) ] (H) 1
+[(H) +-(H) ] (H)1 _}} (denom.)
N _o2
kT MN o
T,V,Nol
o1
1
Xo2
2 o 0 2o
- {)(C) ) -2v v (C) + V(C) 1
+ X + o2
S x (H)
Xol o201 -
2 2 1
- +Xox [(H) I-(H) i] (denom.)
-x+ X- 0 1 1(~+ 1-
2
denom. = {x + x o-(H)
ol ol0 -
2
}{1 + x (H) + x (H) + xx [(H) (H) -(H) 2]
11 -+ + + +-
2 2 2
- XolX2{o2 x (H) + x (H)1
0 1 2 o -2 + 1-+ 1- I-
2 2
+ x x [(H) (H) +(H) (H)
+ 1+ 1- ++
- 2(H) 1+(H) (H) ] }
(4-58)
1 2 0o o+2 x
pKTRT = Xo1(1-C1) 2VXolXo2( 12
+ v o2(1-C2)
02 22
(4-55)
+ x x [(H) -(H)2 ]
+ ++ +-
(4-56)
S2C
22
11}{ X22+
where
(4-57)
Then
(4-59)
v o/KRT = v{x 1(-C2) + xo (1-C2)}
o2 T 01 12 o2 22
Nv ny
pT RT 3No2
T,P,N
UZny+
pv o2 T
02 ,
TPol
= (vx )2(1-CG )(1-CG9-o-o ) 2 (4-61)
01 11 22 12
(4-62)
2 o
2V C
22
The direct correlation function integrals can be related to infin-
ite dilution values and the Debye-Hueckel limiting law expressions.
The correct expressions are
vkny = S (vxo2) + O(x2)
(4-63)
where
S = 2 2 N i3/2
S = I N2 e { .z2}]
y A clRT 1 ii)
(4-64)
with E1 the solvent dielectric constant, NA, Avagadro's number and e,
the unit charge. Then (Redlich and Meyer, 1964)
-oo 3 ny+
o vo2 = VRT P TN N
2 v 01 ol' o2
= ~ SV(PXol2) + O(Xo2)
Snel K1
S = S RT -
v P 3
(4-65a)
(4-65b)
(4-66)
where
(4-60)
80
and K1 is the solvent isothermal compressibility. The solution isother-
mal compressibility defines the apparent molar compressibility, k,(Guc-
ker, 1933).
KT = XolK P /Pol + Xo2ckP
(4-67)
where
xol + o2 = 1
and
k Sk(PXo2) + O(xo2)
k ko =3ko22
SK1
+ 2 a
- 6K P
+9( 2nE
+ 9
'8 nP '
(General expressions for the multisalt case are found in Appendix I.)
These equations yield
o 1 kXo2
1-Cl = 1 (1+x +
11 plK RT2 2 o2
P1K1 (K.RT
S
v
K RT
S
+ I
2
2 Sk
+3K
3 2
KIRT
3/2 2
P x3 + 0(x )
1 o2 02
-oo
S v2
1-Co = 2I
12 K RT
1
1-0
22
S
2v2
S
2 VK1
1
PnE
I ,n ,
3- P ) + 0(x2)
T
-1 + constant + 0(x2)
x 02
S021
with
S RT
S =
k 4
(4-68a)
(4-68b)
(4-69)
(4-70)
(4-71)
Double-Salt, Single-Solvent Case
Assuming the two salts contain different ions, we write the V, x
and z matrices, which are simple extensions of the one-salt case. Again
we choose the positive ions as dependent, and the negative ions as in-
dependent.
O 0
0- 0
V3+
V3-
0
l/v2+
0
0
0
0
0
I/v3+
0
0
0
-v2- 2+
0
1
L 0
0
0
-v 3/v3+
V3- 3+
0
1
(4-72)
Multiplying the first column of z by z2_, and the second by z3_, we
again get the matrix of ionic charges for salts 2 and 3. The A matrix,
in terms of total correlation function integrals, becomes quite com-
plex. In terms of component matrices, it is
Xol+ xol l11
xolxo2(H)12+
Xolxo3(H)13+
xolxo2(H)12+
xolxo3()13+
2
xo2++ xo2(-H)2+2+
xo2xo3(H)2+3+
2
xo3 /V3+ x 3 H)
o3 3+ o3 3+3+
xol2_{(H)12- (H)12+}
x2- /2++ o2x2- {()2+2-- ()2+2+}
x2-xo3{()3+2--() 2+3+}
xolx3-{(H)13--(H) 13+}
o2x3- (H)2+3-- (H2+3+
X3-
- + o3X3- (H3+3- 3+34
3+
Z =
82
x3 /v3+ + x 2_{ () 3+3+(H) -2 (H) }
3- 3- 3+ 3- H3+3+ 3-3- 3+3-
-x 3{ (H) +(H) _3- (H) (H) 3
2-3-- 2+3+ 2-3- 2+3- 3+2-
x -x x {(H) +++(H) (H) (H) }
I2-3- 2+3+ 2-3- 2+3- 2-3+
S_2- 2
- x + x {((H) +2(H) -2(H)2 }
V2+ 2- 2- 2+2+ 2-2- 2+2-
xol2- ) 12-- H 12+}
2-
X+ xX {(H) -(H)} -(
S2+ 2-- 2+2- -2+2+
x 2 { () 3+2--(H) 2+3+
03 033+2- 2+3+
XolX3- (H) 13-- -(H) 13+}
x o2x3- (H) 2+3-- (H) 2+3+
x3-
V3+ 3+3- -3+3+
2 2
2- -/V2+ 2-[(H) 2+2+(H)2-2- -2(H)2+2- 3-3- 3+ X3- (H)3+3+
+(H) -2(H)3+3]}- xx () +() -( ) -() 12
3-3-- 3+3- 2- 3 2+3+ 2-3- (H)2+3-- 3+2-
(4-73)
-1
The elements of A in terms of the total correlation function
integrals must be determined by carrying out the matrix operations in
equation (4-73) and then inverting the result. It is clear that this
would be a very complex expression.
However, equations (4-34b) and (4-72) allow us to write the ele-
-1
ments of A easily, in terms of the direct correlation function integ-
rals. These elements are:
N D ol
kT N No2'o3
T,V,N ,N
o2' o3
_ 1 o
Io C
x 11
ol
(4-74)
N ol
kT 8N o
T,V,NoNo3
N 02
k T 9N
o 02
kT N VNo'No2
N D1 o2
S 2C 12
-2 -2 C
x2 22
o2
- -v VC
2 3 23
where the C (a,B=1,2) above are defined in equations (4-56)
and
and (4-57)
o +(C o o +o
o 2+3+ ()2+3+ '2+3-(co)2+3- V2-3+( 2-3++ V2- 3-( 2-3-
23 V2 3
(4-78)
In the above equations, the other derivatives of the chemical potential
may be found by simply interchanging subscripts 2 and 3 and defining
Co and Co in a manner analogous to C and C Expressions for the
13 33 12 22
thermodynamic properties are the same as those of Friedman and Ramana-
than only when the constraint
(X + XHX)z = 0
(4-79)
is assumed (for the results in this case, see Appendix G). The expres-
sions in terms of the total correlations quickly grow cumbersome, but
the expressions in terms of direct correlations remain unchanged. Since
the matrix in terms of the (H).. must be inverted, and because of the
--I]
(4-75)
(4-76)
(4-77)
complexity of the relationship between A and the (H).., we prefer to
focus on the expressions using direct correlation function integrals,
which avoid both complexities. Using equations (4-44) to (4-46), we
find
1 2 o o2 o 22
p[KTRT olCl + 22XolXo2C12 3XolXo3C13 2o222
o 2 2 o
+ 2\) N) xx C 0 + \2x 2 C 0
+ 2 23o2Xo3C23 3Xo3C33
Vo2
So 1 T 2[XolC + 2Xo2C22 + 3Xo3C23]
K RT 12 2 22
(4-80)
(4-81)
pKTRT N2 T,P,
oy2
2 o 2 o o
= 2{[1-C22 [x 2(-Co ) + 2V3xol3(1-C3)
2 2 (l o 2
+ 3x o3(1-C 3)-[x ol(l-C )+ o3 (1-C23)]
3 o3 33 01 12 3 o3 23
(4-82)
2 NV2 '
pK RT
T)
9 ny2
8N o
03
T,P,N
oy#3
2 o
= V V{[-C 3[x (1-C) + xolXo (1-C
2 3 23 01 11 2 01x2 12
oo
+ 3xolxo3(1-C3) + \3Xo2 x3(l-C23)]
- [x(ol-C 2) + 2x 2(1-C22)][x o(l-C3)
01 12 2 o2 22 01 13
(4-83)
8,ny2
p2 o2
o2
2 o
2 22
(4-84)
T,p#2
oyf2
+ 2Xo3(-C33)]}
8anyi2 o
P2 p2 = 3C 23 (4-85)
03
2 o3 T,p
oy#3
Expressions involving salt 3 can be obtained from those in equations
(4-81) to (4-83) by interchanging subscripts 2 and 3.
The preceding examples show that it is advantageous to write the
thermodynamic properties in terms of salt-salt and salt-solvent func-
tions, taking linear combinations of the ionic direct correlation func-
tion integrals. If we note that
n
o
x. = x V (4-86)
i o t ia
a=l
equations (4-44) through (4-46) can be transformed into these salt
functions in the completely general case. We have
n n
o o
-1- = 0 v V x o(1-Co ) (4-87)
PKTRT a=l B=l
n
v o
oa = v x (1-C0) (4-88)
K RT a 6 ,o 1
T 3=1
n n
o o
1 0o
PKRT (G) = VB v V oyX2 [(-C 0)(l-C ) -
T y=l 6=1
(1-Co )(1-C ] (4-89)
ay
where
n1 n
CO X V [l-(C). ]/V V (4-90)
i=] j=l
and V = 1 if a is a solvent. Similarly, equation (4-50a) becomes
oa
M&ny
Pv -lt = V v C0 (4-91)
V a p ,p a B a3
Toy#p
Single-Solvent, Common Ion Case; Harned's Rule
The equations change somewhat when a two-salt system has a common
ion, e.g., the anion, now denoted -. Then
= [2 (Co) + 2v (Co + 2 (Co) i2 (4-92)
22 2+ 2+2+ 22+2- 2+- 2- 2(
C v (Co (C (Co) +v
23 = 2+ 3+() 2+3+ V2+ 3-()2+- 2- 3+(C) 3+-+ '2-_V3(C-O)
S2 V23 (4-93)
C [v (C) + 2v (C ) 2 (Co) _]/ (4-94)
33 3+- (C3+3++ 3+ 3-( 3+- 3- 3
and the number of different direct correlation functions has been re-
duced from fifteen to ten. We choose the two positive ions as the
dependent species.
The v, v and z matrices for this system are
1 0 0 1 0 0
0 2+ 0 0 1/v2+ 0
S= 0 0 =02 (4-95)
0 0 v3+ 0 0 1/3+
0 V2- 3- 0 0 0
and
z = 0 -v /v2+ /V3+ 1 IT (4-96)
--2- 2+ 3- 3+
where again, multiplication by z gives the vector of ionic charges.
Since the C o's use the v matrix directly, whereas the (H)ij's
must be transformed by z, inverted, further transformed by combinations
of X, H, z, and y, then inverted again, it is obvious that thermody-
namic properties in terms of the direct correlation functions are
much simpler mathematically.
Harned's rule for the mean ionic activity coefficients states
that the logarithm of the mean ionic
of one salt varies proportionally to
when the total ionic strength, I, is
ny m2
Dm3 T,P,l
= c23
molal activity coefficient, y ,
the molality of the other salt
kept constant. That is
(4-97)
where ym2 and I are defined by
mZ
U2/kT = V 22n [m + 02
I= 2 2 2 2
( 2+2+ + v2 z2)m2 + (3+ z3+ V z2)m3
2+z2+ 2- 2 3+z3+ 3- 3,
m is the molality of salt a and
salt 2.
m2 is the mean ionic molality of
V 2+ m2- I 2
m2 = 2+ -
(4-99)
From the definitions of density- and molality-based activity
coefficients in equations (4-48) and (4-98a), equation (4-97) can be
rewritten
(4-98a)
(4-98b)
~~
~~ |