TN E5333

“J LIBFEARY a

% Pavia: E13593 Sm:

- I
V

This is to certify that the

dissertation entitled

Ion Transport Processes

presented by

John H. Leckey

has been accepted towards fulfillment
of the requirements for

Ph 0 D . degree in ChemiS try

 

M

Major professor

Date 7/494,

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

MSU

LlIlARllS

 

“-

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

 

 

ION TRANSPORT PROCESSES

By

John H. Leckey

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Chemistry

1981

ABSTRACT

ION TRANSPORT PROCESSES

By
John H. Leckey

Ions flow due to electrochemical potential gradients
so as to minimize the production of entropy. In a homo-
geneous substance the electrochemical potential can be
separated into a chemical potential gradient, the driving
force for chemical diffusion, and an electrical potential
gradient, the driving force for electrical migration.
The development of a macrosc0pic set of partial differential
equations from nonequilibrium thermodynamics, hydrodynamics,
and electrostatics, which describes the isothermal trans—
port of ions, is reported here. This system of equations
contains a set of transport coefficients which result from
the partial separation of the electrical and chemical ef-
fects involved in ion transport. The characterization of
these coefficients and the solution of this system of equa-
tions for liquid—liquid Junctions and double layer relaxa-
tions is also reported here. The solutions are in terms
of numerical finite difference computer simulations and

approximate analytical equations. All activity coefficients

John H. Leckey

have been retained in the numerical solutions and whenever
possible in the analytical solutions.

The time range of the dependence of the cell potential
on single-ion activity coefficients has been delineated
for several Ag/AgCl concentrationcells containing a
single chloride electrolyte. A simple formula is obtained
for the cell potential as a function of time, and in terms
of Onsager coefficients, ionic valences, and the permittivity.
The results have two major implications: (1) They can be
used as guidelines for choosing appropriate electrolytes,
time domains, and boundary conditions for future experi-
ments to determine single-ion activity coefficients; (2)
They suggest complicated interdependences between elec-
trical transport processes, on the one hand, and single—
ion activity coefficients, on the other.

An analytical solution in terms of a perturbation
series has been derived for the equilibrium double layer
in the presence of a charge on the electrode. Activity
coefficients appear explicitly in the solution. An ap-
proximate analytical solution has been obtained, using
Laplace transforms, for double layer relaxation following
rapid electrode charging. The solution results in a
simple formula for the potential difference across the
double layer as a function of time. The numerical solu-
tion of this problem has been used to delineate the range
of validity of the time dependent analytical solution.

Results are compared to previous studies.

To Pam

ii_

ACKNOWLEDGMENTS

I wish to thank the Department of Chemistry for the
financial support I have received from teaching assist-
antships. I also wish to thank the General Electric
Foundation and Amway Corporation for awarding me fellow-
ships through the Department of Chemistry for the summers
of 1979 and 1980, respectively.

I am most grateful to Professor Frederick H. Horne
whose teachings of nonequilibrium thermodynamics have
been the cornerstone of this work. I appreciate his en-
couragement and useful suggestions, which aided in the
completion of this dissertation.

I would like to thank Dr. Thomas H. Pierce for helpful
discussions about computer programming, and especially for
making available the matrix inversion package used in this
work. I wish to thank Dr. Tom Atkinson for providing his
curve plotting package MULPLT used to prepare the figures
in Chapter 2, and Folim Halaka and Rob Thompson for help-
ing me use this package on a PDP-ll computer.

I am indebted to Bruce Borey for reading and comment-
ing on several manuscripts, and to Dr. Stanley R. Crouch
for critically reviewing this dissertation on very short

notice.

iii

Most of all, I thank my wife, Pam, for her support,
encouragement, and unbridled optimism. Her understanding
throughout difficult circumstances has made the completion

of this work possible.

iv

TABLE OF CONTENTS

Chapter Page
LIST OF TABLES. . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . x
I. INTRODUCTION. . . . . . . . . . . . . . . . . . l
A. Phenomenology of Ion Transport. 1
B. Objectives of the Research. 3
C. Plan of the Dissertation. 5
2. ION TRANSPORT EQUATIONS . . . . . . . . . . 8
A. Introduction. . . . 8
B. Continuity Equations. 9
C. Nonequilibrium Thermodynamic
Equations . . . . . . . . . . . . . . . . . 9
D. Electrochemical Potentials and
Activity Coefficients . . . . . . . . . . . II
E. Electrostatic Equations . . . . . . . . . . 1A
F. Combined System of Equations. . . . . . . . 15
G. Transformation of Concentration
Variables . . . . . . . . . . . . . . . . . 17
H. Discussion of Transport Coefficients
for Transformed Equations . . . . . . . . . 20
3. ION TRANSPORT BOUNDARY CONDITIONS . . . . . . . Al
A. Introduction. . . . . . . . . . . . . . . . Al
B. Derivation of Kinetic Boundary
Conditions. . . . . . . . . . . . . . . . A3

Chapter Page
C. Applications of Kinetic Boundary
Conditions. . . . . . . . . . . A8

A. NUMERICAL FORMULATION OF THE
TRANSPORT EQUATIONS . . . . . . . . . . . . . . 50

A. General Finite Difference Approach
for Parabolic Partial Differential

Equations . . . . . . . . . . . . . . 50
B. Difference Equations for Ion
Transport . . . . . . . . . . . . . . . . . SLl
1. Space and Time Grids. . . . . . . . . . 5A
2. Application to Ion Transport
Derivatives . . . . . . . . . . . . . . 57

3. 'Boundary Conditions . . . . . . . . . . 58
A. System in Matrix Form . . . . . . . . . 62

C. Computer Simulation Discussion. . . . . . . 63
5. LIQUID-LIQUID JUNCTIONS AND SINGLE-
ION EFFECTS . . . . . . . . . . . . . . 68
A. Introduction. . . . . . . . . . . . . . . . 68
B. Initial and Boundary Conditions . . . . . . 73
C. Simplifying Approximations. . . . . . . . . 7”
D. Analytical Solution . . . . . . . . . . . . 75
E. Full Cell Equation Results. . . . . . . . . 79
F. Numerical Solution. . . . . . . . . . . . . 81
G. Conclusions . . . . . . . . .i. . . . . . . 95

6. DOUBLE LAYER RELAXATION AT BLOCKED
ELECTRODES. . . . . . . . . . . . . . . 97

A. Introduction. . . . . . . . . . . . . . . . 97
B. Equilibrium Solution. . . . . . . . . . . . 100

C. Time Dependent Solution . . . . . . . . . . 120

vi

Chapter
D. Computer Simulation Results and
Discussion. . . . . . . .
E. Comparison with Previous Results.
7. SUGGESTIONS FOR FUTURE WORK
A. Future Thermodynamic Studies.

B. Program Improvement Suggestions

APPENDIX A - TRANSPORT COEFFICIENTS
APPENDIX B - DIMENSIONLESS VARIABLES USED BY
PROGRAM "IONFLO" . . . . . . . .
APPENDIX C - TRANSPORT EQUATIONS IN FINITE
DIFFERENCE FORM. . . . . . .
APPENDIX D - NEWTON'S METHOD AS USED BY
PROGRAM "IONFLO" ' .
APPENDIX E - LISTING OF PROGRAM IONFLO.
APPENDIX F - SOLUTION OF SYSTEMS OF LINEAR
PARTIAL DIFFERENTIAL EQUATIONS
APPENDIX G - ANALYTICAL SOLUTION OF D FOR
TIME DEPENDENT DOUBLE LAYER
RELAXATION . . . . . .
REFERENCES.

vii

Page

133
1A6

15A
15A
156
158

161

163

167
170

197

200

205

Table

6-2

6-3

LIST OF TABLES

Comparison of Experimental Cell

Potentials (EEXP) for Cell (A) with

Cell Potentials (ECALC) From Equa-

tion (5-20)

Comparison of Equation (6-72) and

Computer Simulation for 0.1 M KCl.

A¢* Indicates Equation (6-72) Results.
-1

a:

a:

6 x 10

1 cm.

6

S

of = 10‘7 C/cm2.

Comparison of Equation (6-72) and

Computer Simulation for 0.1 M NaCl.

A¢*

O.

a

Indicates Equation (6-72) Results.
-1

6 x 10

1 cm.

6

S

Q° = 10'7 C/cm2.

Comparison of Equation (6-72) and

Computer Simulation for 0.1 M HCl.

A¢*

a:

a

Indicates Equation (6-72) Results.
-1

6x106

1 cm.

S

Q° = 10‘7 C/cm2.

viii

Page

82

139

. 140

lUl

Table

6-A

6-6

Page

Comparison of Equation (6-72) and

Computer Simulation for 10"3 M KCl.

A¢* Indicates Equation (6-72) Results.

o = 6 x 106 s‘l. Q° = 10'7 C/cm2.

a = 1 cm. . . . . ... . . . . . . . . . . . 143
Comparison of Equation (6-72) and

Computer Simulation for 10'3 M NaCl.

A¢* Indicates Equation (6-72) Results.

3 6 x lo6 s’l. Q° = 10‘7 C/cm2.
a = 1 cm. . . . . . . . . . . . . . . . . . 1AA
Comparison of Equation (6-72) and

Computer Simulation for 10'3 M HCl.

A¢* Indicates Equation (6-72) Results.

a 6 x 106 8'1. Q° = 10"7 C/cm2.

a 1 cm. . . . . . . . . . . . . . . . . . 1A5

ix

LIST OF FIGURES

Dependence of YSS
at B equals zero.
(0) NaCl, (d) KCl
Dependence of YDD
at B equals zero.
(c) NaCl, (d) KCl
Dependence of YSD
at B equals zero.
(0) NaCl, (d) KCl
Dependence of YDS
at B equals zero.
(c) NaCl, (d) KCl
Dependence of YSS
for NaCl. (a) B

o, (c) B = 0.5 M'

Dependence of YDS

on ionic strength

(a) HCl, (b) LiCl,

on ionic strength

(a) HCl, (b) LiCl,

on ionic strength

(a) HCl, (b) LiCl,

on ionic strength

(a) HCl, (b) LiCl,

on ionic strength

-005 M-l, (b) B =

on ionic strength
1

for NaCl. (a) B = -o.5 m‘ , (b) B =

0, (c) B = 0.5 M-

Dependence of YSE/S on ionic strength.

(a) HCl, (b) LiCl,

(0) NaCl, (d) KCl.

Page

25

27

29

31

3A

36

38

Figure

2.8

5.2

5.3

Page

Dependence of YDE/S on ionic strength.

(a) HCl, (b) L101, (0) NaCl, (d) KCl. . . . A0
Interfacial region at a phase boundary. . . A2
Uniform grid mesh for finite dif-

ference equations . . . . . . . . . . . . . 52
Nonuniform grid mesh for improved ac-

curacy near boundaries and for more

efficient approach to steady-state

equilibrium . . . . . . . . . . . . . . . . 56
Flow diagram for program IONFLO . . . . . . 65
Predicted charge density surface near

the liquid-liquid interface at short

time for a NaCl concentration cell

operating between concentrations of

0.11 M and 0.09 M. B is zero for

this simulation (y+ = y_). The time

row at t = 0 has been omitted since

it would obscure the lower portion

of the surface. . . . . . . . . . . . . . . 86
Constant time planes from Figure 5.1.

(a) 5 x 10.98, (b) 1 x 10-88

10-85 . . . . . . . . . . . . . . . . . . . 88

, (c) 2 x

Charge density X§° cell position for
different values of B at t = l x 10-85

for a NaCl concentration cell operating

xi

Figure Page

between concentrations of 0.11 M and
1

0.09 M. (a) B = —l.OO M- , (b) B = O,
(c) B = 1.00 M.1 . . . . . . . . . . . . . 90
5.A Cell potential vs time for cell.

(A) with M = Na, SL 0.11 M, and

1.00 M'l, (b)
l.

S = 0.09 M. (a) B

R

B = o, (c) B = -1.00 M’ 92

5.5 Coalescence time Kg YDE for 0.1 M

solutions, with 5% criterion. (a)

HCI, (b) BaCl2, (C) KCI, (d) NaCl,

(e) LiCl. . . . . . . . . . . . . . . . . . 9A
6.1 Contributions to s' from S2 and s3

for a 0.1 M NaCl solution with an

6

electrode charge of 2 x 10' C/cm2.

b = -0.1; (A) s', (B) s (C) s . . . . . . 113

2’ 3
6.2 Contributions to d' from d1, d2,

for a 0.1 M NaCl solution with
6

and (13

an electrode charge of 2 x 10- C/cm2.

b = -0.1, (A) dl, (B) d1. (C) d2.

d . . . . . . . . . . . . . . . . . . .
(D) 3 115

6.3 Dependence on a of D at x = -a/2 as
a function of time as predicted by

Equation (6-65) for a particular NaCl

2 -6

solution. SB = 10- M, Q0 = 10

C/cm2, (A) a = 3 X 107 8-1, (B) a =
l x 107 s_¥, (c) a = 6 x 106 3-1. . . . . . 126

xii

Figure Page

6.A Concentration dependence of the
relaxation of D at x = -a/2 for an
electrode charging time constant of

107 s'1 as predicted by Equation
2

(6-65). (A) SB = A x lo' M, (B)
SB a 2 x lo‘2 M, (c) SB = 1 x lo‘2 M. . . . 128
6.5 Effect of time on D versus z/a profiles

as predicted by Equation (6-67).
Q0 = -lo‘7 C/cm2, a = 6 x lo6 s‘l,

= 5 x 1075, (A) lo'7s, (B) lo'8s,

YDB
(c) lo‘9s . . . . . . . . . . .-. . . . . . 131
6.6 Computer simulation of charge density

difference versus distance from

6s for two 0.1 M

NaCl solutions. a = 6 x 106 8-1,

electrode at 10‘

Q0 = "10-7 C/cm2, 9 corresponds to

B = -1 M‘1

B = 1 M'1 . . . . . . . . . . . . . . . . . 135

and p* corresponds to

6.7 Effect of B on Na+ and C1- activity

coefficients. (1) yNa+ if B = 1 M“1

and yCl_ if B = -l M“1, (2) y: for

NaCl as measured by Brown and MacInnes

1

[1936], (3) YCl- if B = 1 M‘ and

yNa+ if B = -1 M‘l. . . . . . . . . . . . . 137

xiii

Figure

6.8

6.9

6.10

Page

Electrode charge versus time profile

assumed by Newman [1973] for analytical
solutions to the Nernst-Planck equa-

tions . . . . . . . . . . . . . . . . . . . 1A8
Comparison of the experiments of Anson

_£‘_l. [1969] (A), simulations using

IONFLO (B), and Newman's [1973] ana-

lytical predictions (C) for the relaxa-

tion potential of HCl at a blocked

electrode. a = 6 x 106 3-1, Q0 =

10'5 C/cm2. Time axis is logarithmic . . . 150
Comparison of the experiments of

Anson EE.§l- [1969] (A), simulations

using IONFLO (B), and Newman's [1973]

analytical predictions (0) for the

relaxation potential of KCl at a

blocked electrode. a = 6 x lo6 s‘l,

Q0 a lo‘5 C/cm2. Time axis is

logarithmic . . . . . . . . . . . . . . . . 152

CHAPTER 1

INTRODUCTION

A. Phenomenology of Ion Transport

Application of a current impulse to an electrolyte
solution causes the flow of ions. Positive ions move in
the same direction as the current and negative ions move
in the opposite direction. When the ions reach an elec-
trode one of two possible events will occur. One possibility
is the transfer of the ions' charge to and/or across the
electrode-solution interface. This is done by either
specific adsorption, an electron transfer reaction, or, as
in the case of some intercalated electrodes, the simple
migration of the ionsixux>the electrode (Armand 33 3;.
[1978]). The other possibility is that the ions will re-
tain charge and be somehow blocked from interacting with
the electrode (Chang and Jaffé [1952]). In this case the
current pulse is usually referred to as a charge impulse or
coulostatic pulse (van Leeuwen [1978]). If the current is
turned off and the interface remains blocked, the concen-
Itration of the ions at the electrode will not increase
indefinitely. This is so despite the fact that the charge

on the electrode will be attracting all of the ions of the

opposite charge to the electrode. The balancing effect,

of course, is that the ions will tend to diffuse from the
high concentration near the electrode to a lower concentra-
tion in the bulk of the solution. Eventually these effects
result in a fixed concentration profile of each ion and a
time invariant electric field.

Electrolyte solutions can also generate an electric
field without any type of electric pulse. A diffusing
electrolyte solution creates an electric field and,there-
fore,a potential difference between any two points in the
solution along the direction of the diffusion. This is
because the positive and negative ions do not diffuse at
exactly the same rate. Different masses, volumes, and
interactions with the solvent cause different rates of
diffusion. Of course, one type of ion does not get very
far ahead of or behind the other type. The electrical force
created by their separation keeps them from separating
very far in the solution. The potential difference caused
by this effect is nevertheless significant. This is the
origin of the liquid Junction potential present in certain
types of electrochemical cells.

Although these two events seem only remotely related,
they can be treated theoretically using the same non-
equilibrium thermodynamic formalism. In the language of
nonequilibrium thermodynamics, the ions will flow due to

electrochemical potential gradients so as to minimize the

production of entropy in the system (Onsager [1931]).
Besides the above effects, the same set of nonequilibrium
thermodynamic, hydrodynamic, and electrostatic equations
describes many other important processes. The flow of ions
through membranes, the movement of electrons and holes in
semiconductors, and the performance of completely solid

and polymer batteries can all be described macroscopically
using thermodynamics, hydrodynamics, and electrostatics.
The unifying concept is the competition between electrical
migration and chemical diffusion. In a general way this
dissertation deals with all of these effects, and in par-
ticular with the detailed theoretical study of double layer
relaxation at blocked interfaces, liquid Junction poten-
tials, and how single-ion activity coefficients are related

to both of these.

B. ObJectives of the Research

Although ion transport has been treated extensively in

theoretical research (Turq and Chenla [1977], Vallet and

Braunstein [1977], and Brady and Turner [1978] are some

of the relatively recent studies.), it has not previously

been approached using a full, rigorous nonequilibrium thermo-
dynamic analysis. Rarely have activity coefficients or cross-
diffusional effects been included in ion transport equa-
tions. Their omission is valid for many types of non-

electrolyte transport processes because in that case

their influence is often small. In electrolyte solutions,
however, they are significant for concentrations as low
as 0.01M. In a 0.01M solution of BaCl2 in water, for
instance, the cross-Onsager coefficient is 20% of the Ba
Onsager coefficient. The mean molar activity coefficient
of HCl at 0.01M is 0.906, £1 departure of nearly 10% from
ideality. Both the cross-Onsager coefficient and the mean
molar activity coefficient rise rapidly as the concentra-
tion increases above 0.01M. In an effort to determine the
effect of these terms on measurable quantities, such as
potential differences, they have been retained throughout
the numerical studies of the dissertation. They have also
been retained whenever possible in analytical derivations.
In a related obJective, the role of single-ion activity
coefficients in electrolyte transport processes has been in-
vestigated theoretically. Although they have not been
determined experimentally yet, information about them is_
essential for further interpretation of many phenomena.
As Conway [1978] has noted, detailed explanation of phase
transfer reactions, electrodialysis, and ion-solvent inter-
actions all require properties of individual ions. Hall
[1971], [1973], [1977], and [1980] has extensively studied
single-ion activities, and their relationship to the
equilibrium thermodynamics of charged interfaces. The
fact that single-ion activity coefficients deviate sig-

nificantly from unity at concentrations larger than 0.01M,

yet are undetermined, indicates that conventional steady-
state experiments do not depend on them in an obvious way.
With this in mind, most of the research has focused on time
dependent phenomena. One of the obJectives has been to
discover situations which might lead to experimental de-

termination of single-ion activity coefficients.

A fourth obJective has been to solve transport equa-
tions analytically whenever possible. The analytical re-
sults, although less accurate than the numerical ones,
have the advantage that the final equation specifically
reflects the effect that the input parameters (ionic
valences, activity coefficients, Onsager coefficients, egg.)
have on the final result.

It is hoped that the results presented here will lead
to a greater understanding of ion transport in general, in
whatever type of system it occurs. More specifically,
it-is hOped that these results will stimulate research in
experimental and theoretical studies of time dependent ion

transport phenomena.

C. Plan of the Dissertation

Chapter 2 begins with those fundamental equations of
hydrodynamics, nonequilibrium thermodynamics and electro-
statics necessary to formulate the ion transport equations.
The rest of Chapter 2 deals with formulation of suitable

ion transport equations and interpretation of the transport

coefficients. Because of the tremendous importance of
boundary conditions on ion transport phenomena, Chapter 3
is devoted entirely to their derivation and application.
Chapter A describes the combination of the ion transport
equations and boundary conditions into a general finite
difference computer program. This program is the first of
its type to retain all thermodynamic terms.

Chapters 5 and 6 describe the application of the ion
transport equations to two specific problems. Both chap-
ters present analytical as well as numerical results.
-Chapter 5 follows very closely Leckey and Horne [1981],
and deals particularly with single-ion activity coefficients
and their dependent liquid Junction potentials. The time
range of the dependence of the cell potential on single-
ion activity coefficients is delineated for several Ag/AgCl
concentration cells containing a single chloride elec-
trolyte. A simple formula is obtained for the Cell poten-
tial as a function of time, and in terms of Onsager co-
efficients, ionic valences, and the permittivity. Chapter
6 describes transient and equilibrium double layers created
by a rapid current pulse at a blocked electrode. The
equilibrium ion distribution equations are solved using a
perturbation technique. At infinite dilution the solution
approaches the Gouy-Chapman result. An approximate equa-
tion for double layer relaxation is solved using Laplace

transforms. From this a simple formula is obtained for the

potential difference across the double layer as a function
of time. This result and computer simulation results

are compared to previous theories and the experiments of
Anson 23 al. [1969].

The need for improved computer programs which account
for concentration dependence of the dielectric coefficient,
specific ion adsorption, and more sophisticated boundary
conditions is discussed in Chapter 7. The possibility of
studying liquid Junction potentials using thermal dif-

fusion is also discussed there.

CHAPTER 2

ION TRANSPORT EQUATIONS

A. Introduction

The macroscopic equations necessary for the description
of ion transport are assembled from the fields of hydro-
dynamics, electrostatics, and nonequilibrium thermodynamics.
The systems dealt with in this dissertation are isotropic,
nonreacting (except possibly at interfaces), isothermal,
binary electrolyte-solvent mixtures not subJect to magnetic
or gravitational fields. All concentration and electric
field gradients are assumed to be in one dimension. Electro-
neutrality is not assumed.

The fundamental equations employed are not new. There-
fore, more consideration will be given to simplifications,
transformations, and in latter chapters solutions, rather
than to derivations of these equations. Details of the
derivations and assumptions involved in obtaining the start-
ing equations are available from de Groot and Mazur [1962],
Fitts [1972], Hasse [1969], Duckworth [1960], and Lewis and
Randall [1961].

B. Continuitnyguations

The concentrations of any two of the three components
(solvent, cation, anion) are independent in a charged
system. For the positive and negative ions as the two
independent species, the relevant continuity of mass

equations are

(ac+/at) -(a/ax)(J+H + c+vo), (2-1)

(ac_/at) -(3/3x)(J_H + c_vo) . (2-2)

Here t is time, x is the space coordinate, c+ and c_ are
the positive and negative ion concentrations, J+H and J_H
are the Hittorf diffusion fluxes, and v0 is the velocity
of the solvent.

A Hittorf diffusion flux has as reference velocity the
velocity of one of the components. In many cases it is
useful to choose the velocity of the solvent for the

reference velocity since it is often small enough to be

neglected.

C. Nonequilibrium Thermodynamic Equations

 

The nonequilibrium thermodynamic equations which relate
the fluxes in Equations (2.1) and (2.2) to electrochemical

potential gradients are the Onsager flux equations for

10

coupled diffusion,

l
L,
I

+ - 2++(an+/ax) + 2+_(an_/3x), (2-3)

I
c_.
ll

£_+(3u+/3x) + £__(3u_/3x) . (2-A)

u+ and u_ are the electrochemical potentials of the cation

and anion respectively, and the t are the Onsager co-

13
efficients.-

These flux equations have been used in preference to
the simpler Nernst-Planck flux equations for two main
reasons. In a very dilute aqueous electrolyte solution,
the ions are so far apart that long range Coulombic forces
are the principal interactions, and specific effects are
negligible. In this circumstance the flow of an ion is
proportional to gradients only of its own properties and
is not influenced by the presence of other distant ions.
This reasoning leads to the development of the Nernst-
Planck flux equations. At higher concentrations activity
coefficients become specific, ionic conductances are no
longer additive, and the Nernst-Hartley equation for the
diffusion coefficient fails. Thus, intuitively, the flow
of an ion must be influenced by the other ions present and
the gradients of their properties as well. The macroscopic
theory that includes these effects is nonequilibrium thermo-

dynamics, which leads to the flux equations, (2-3) and (2-A).

11

A second reason for choosing these flux equations is the
RiJ themselves. They are more fundamental than the common
transport coefficients for special cases, such as the
equivalent conductances or diffusion coefficient. It has
been shown explicitly by Miller [1966] that any isothermal
electrical transport process, no matter how complex, can

be completely characterized once these A are known as

13

functions of concentration.

D. Electrochemical Potentials and Activity Coefficients

 

By strict thermodynamic definition u+ and u_ in Equa-
tions (2-3) and (2-A) should be termed chemical potentials
instead of electrochemical potentials. Since the only work
terms being dealt with are electrical, however, Guggenheim's
[1929] precedent will be followed and the term electrochemi-
cal potential will be used for the "entire" potential.

The electrochemical potentials for the ions are

m + RT tn(c+y+) + z+F¢ , (2-5)

1:
+
II

11+

u = u_m + RT £n(c_y_) + Z_F¢ (2-6)
where,

ui0° denotes the infinite dilution standard state, R is the
gas constant, T is the absolute temperature, Z1 is the

charge number, F is Faraday's constant, and ¢ is the electric

potential. The activity coefficient yi is on the molarity

12

scale, and c1 is the molarity of component i.
The chemical potential derivatives of Equations (2-3)
and (2-A) in the variables electric field, concentration,

and ionic strength are

(3u+/3x) = (RT/0+)(ac+/8x) + RT(v/2u+)(M+B)(BIc/3x)
+ z+F(3¢/3x) , (2-7)
(3u_/8x) = (RT/c_)(30_/3x) + RT(v/ZV_)(M-B)(3Ic/3x)

+ z_F(3¢/3x) , (2-8)

where IC is the ionic strength on the molarity concentra-

tion scale,

Ic = % (z+2c+ + z_2c_) , (2-9)
and M and B are defined by
M = (axing/316%,P , (2-10)
B = (BAnGi/aIC)T’P . (2-11)
yi has its usual definition
v v+ v-

Y+ = y+ 'y 3 (2‘12)

13

and 6+ is the mean molar activity coefficient deviation

first introduced by Frank [1963].

/y (2-13)

The exponents v and v_ are the positive and negative ion

4.
stoichiometric coefficients, and v = v+ + v_.

The motivation for this way of treating activity co-
efficients is to take advantage of the fact that yi has
been experimentally determined and tabulated (Lewis and
Randall [1961]), whereas 6i has not. From Equations (2-12)

and (2-13) it follows that

(yioi)v/2V+, (2-1A)

v/2v- (2_15)

<<:
ll

(ye/5:)

Thus knowledge of 5+ at a given ionic strength would yield

knowledge of both y+ and y_. Following Goldberg and Frank
[1972] it is assumed that 5+ depends on ionic strength

according to

in 6: = BIc . (2-16)

For purposes of approximate analytical solutions B is

assumed independent of ionic strength. For numerical

1A

solutions B is given functional dependenciestx>match various
activity coefficient theories.
In a NaCl-H20 solution at an ionic strength of 1M,
for example, a value of B of 0.5M-1 corresponds to
y+ = 1.13A and y_ = 0.688; at B = 0, y+ = 0,833 and y_ =

0.833; at B = -0.5M’1, y = 0.688 and y_ = 1.13u. Thus,

4.
at a given ionic strength, increasing B corresponds to

increasing the "stability" of the positive ion relative to

the negative ion. Clearly, knowledge of B would be very
useful in determining the structure of an ionic solution in

water on a molecular scale.

E. Electrostatic Equations

An additional relation between electric potential and

concentration is provided by Poisson's equation

(32¢/ax2> = -(Aw/e)D(X) , (2-17)

where g is the permittivity of the medium and p(x) is the

charge density,

p(x) = F(Z+c + Z_c_) . (2-18)

4.

Thus,

(a2o/ax2) = -(AnF/e)(Z+C+ + Z_C_) . (2-19)

15

For computer simulation purposes (Chapter A), it is
useful to replace the electric potential with the electric

field E,

E = -(ao/ax) , (2-20)

and to replace Poisson's equation with the displacement
current equation (Duckworth [1960]),
H

H)

(aE/at) = (An/e)(l - z J+ — z J

+ (2-21)

This conveniently introduces the current I as an inde-

pendently adjustable variable.

F. Combined System of Equations

 

Combining Equations (2.1) - (2.13) yields Equations for

the time derivatives of c+ and c_ in terms of more practical

space derivatives than the continuity equations provide:

(3/8x)[a++(3c+/3x) + a+_(3c_/3x)

(3c+/3t) =
+ 8+E(3¢/3x) - c+vo] , (2-22)
(ac_/3t) = (3/3x)[a_+(ac+/8x) + d__(ac_/3x)

+

B_E(3¢/3X) - C_vo] , (2-23)

16

where
B+E = F(z+2++ + zl£+_) , (2-2u)
B-E = F(z+i_+ + z_i__) , (2-25)
and
a++/RT = (2+,/c+) - %(z+/z_)(z+ - z_)[z+2+,(M+B) -z_2+_(M-B)].
a+_/RT = (£+_/c_) - %(zd/Z+Xz+-z_)[z+£++(M+B)-z_£+_(M-B)],
a_+/RT = (z_+/c+) - %<z+/z_>(z+-z_)[2+2_+(M+B)-z_2__(M-B)l,
a__/RT = (£__/c_) - %(z_/z+)(z+-z_)[z+2_+(M+B)—z_2__(M—B)J.

(2-26)

The stoichiometric coefficients and ionic valenchxsare re-

lated by
v_ = z+[v/(z+ - z_)], v, = -z_[v/(z+ - z_)l
(2-27)
which are easily obtained from
v+ + v_ = v, v+z+ + v_Z_ = 0 (2-28)

'The quantity [v/(z+ - z_)] is unity except for 2-2, 3—3,

etc., electrolytes because in those cases v_ + 2+ and

17

v. + 12-1.
Equations (2-19), (2-22), and (2-23) form a system of
nonlinear partial differential equations. To the extent
that vO can be neglected there are three dependent vari-
ables, c+, c_, and o, and two independent variables, x and
t. Of course, ¢ can only be determined to within an addi-
tive constant. Because of the complexity of these equations,
an exact solution for this system cannot be expected for any
physically meaningful set of boundary conditions. Part
two of each of Chapters 5 and 6 deals with simplifying ap-
proximations that can be made for this system, after a
suitable transformation of variables, in order to obtain

analytical solutions. Chapter A contains a numerical tech-

nique used in solving this system of equations.

G. Transformation of Concentration Variables

In the absence of chemical reactions, the isothermal
transport of electrolytes can be thought of as the result
of competition between an electric potential gradient and
an activity gradient. These, of course, may oppose or
reinforce each other. In order to separate the mass dif-
fusion effects from the electrical migration effects, the

following transformation is defined,

8 = %-[(c_/z+) - (c+/z_)] , (2-29)

18

% [(C_/z+) + (c+/z_)] . (2-30)

U
II

S can thus be thought of as a sum or neutral electrolyte
concentration, changing mainly due to diffusion. D can be
thought of as a charge concentration, changing mainly due

to electrical effects.

In terms of these new variables,

0+ = z_(D-S),
c_ = z+(D+S), (2-31)
0 = 2z+z_FD
Equation (2-17) becomes,
(32¢/3x2) = -KD , (2-32)
where
K = 8an+z_/e . (2-33)

Besides the fact that Poisson's equation takes on a
simpler form, this transformation has numerical advantages
‘when integrating the Poisson equation. Since the most
important experimental measurement for this system is the

puntential difference between two points, the calculation of

19

Ac from concentration profiles is very important. Because

small changes in (c+ - c_) produce large changes in the po-

tential difference, integrating Equation (2-19) can be risky.

The use of Equation (2-32) eliminates this undesirable

numerical effect.

The differential equations for S and D are, from Equa-

tions (2-22) and (2—23)

(BS/8t)

(ab/3t)

where

Yss

(3/3x)[yss(aS/ax) + YSD(aD/3x) + YSEK-l(3¢/3x)-SVO],
(2-3A)

(a/afiErDS(BS/3x) + yDD(3D/3x) + YDEK'1(3¢/ax)-Dvol,
(2-35)

1/2(a+++a__) - 1/2[(z+/z_)d+_ + [z_/z+)d_+],
-1/2(a++-a__) - 1/2[(z+z_)a+_ - (z_/z+)a_+],

-1/2(a++-a__) +1/2[(z+/z_)a+_ - (z_/z+)c_+],

l/2(d+++a__) +l/2[(Z+/Z_)a+_ + (z_/z+)a_+], (2-36)

(Ans/e)(z_e_E-z+e+E)=-(AnF2/e)(ziz++-zfz ),

(“NF/8)(z_B_E+Z+B+E)=(unF2/€)(Zi£+++2z+z_£+_ +

2

Z_£__) = (An/e)x = ‘1

TD ,

20

where A is the Specific conductance and TD is the dielectric
relaxation time.' In Equations (2-36) and throughout this
dissertation, 2+_ = 2_+. This is consistent with the use

of the flux equations, (2-3) and (2-A).

In terms of these same variables Equation (2-21), the

displacement current equation, becomes
(BE/3t) = (“w/E)I + KCYDS(BS/3X) + YDD(3D/3X) +
-1 R
YDEK (3¢/8X)]. (2_37)

Thus the flux term in the displacement current equation is
the same as the flux term in Equation (2-35).

Equations (2-32), (2-3A), (2-35), and (2-37) form the
set of general partial differential equations which are
solved to various extents in Chapters 5 and 6. To complete
the problem, a set of boundary conditions which can be used

to model many types of processes is presented in Chapter 3.

H. Discussion of Transport Coefficients for Transformed

Equations

 

The rather long expressions for YSS’ YSD’ YDS’ and
YDD in terms of the 213, M, and B appear in Appendix A.
For uni-univalent electrolytes in neutral solutions they,

along with the expressions for YSE and YDE’reduce to:

21

RT
YSS = §§E(£++ + 2£+_ + t__) + Ms (2++ + 22+_ + 2__)
+ BS (£++ —t__)] , (2-38)
' - BE[( ) + MS (2 A ) - BS (i -
YDS - 2S 2-- - 2++ —- - ++ ++
2£+_ + 2 )l,
— BI
YDD ‘ 2s (£++,+ 8--) .
(2-39)
_ BI
YSD ' 28 (2-- 2++)’
=(uF2/m -£)
YSE " 8 ++ —- ’
(2-A0)
YDE = (“UF2/€)(2++ " 2£+_ + £__)

The quantities YSS and YDD can be thought of as diffusion
coefficients for "salt particles" and "charge particles"
respectively. Thus, YSD and YDS are the cross diffusion co-
efficients. The corresponding fluxes for these "particles"

are

S - YSS(3S/3x) + YSD(3D/3x) + YSEK_1(3¢/ax), (2-u1)

I
c_,
I

D ’ YDs(3S/3X) + YDD(3D/3x) + YDEK'l(a¢/ax) . (2—u2)

I
c_,
I

22

These are not Onsager flux equations since YSD # YDS'
They are quite satisfactory for describing experimental
flows, however. The conventional Fick's law salt diffusion
coefficient (Miller [1966]) is proportional to (788 -

-1)
YSBYDsYDE '

From the positive definite nature of the entropy pro-

duction equation (Hasse [1969]) it can be shown that the

quantity (£++ - 2% + 2__) is always positive, as are

+-

1 and t__: From this it can be seen from Equations

++
(2-38) and (2-39) that 788 and YDD must then also be posi-
tive. Thus salt and charge always tend to diffuse in a
direction such as to eliminate the gradients of S and D,
respectively. The signs of YSD and YDS depend on the rela-
tive magnitudes of t++ and t__. Thus, for instance, a
positive D gradient causes a flux of charge in the negative
direction. This can be accomplished by negative ions dif-
fusing in the negative direction and positive ions diffusing
in the positive direction. If 2++ > t__, positive ions are
diffusing faster than negative ions, YSD is negative, and
there is a flux of salt in the direction of the D gradient.

If 2_ < 2++, then YSD is positive, and salt and charge

will diffuse in the same direction.

Figures 2-1 through 2-A are plots of respectively, YSS’
YDD’ YSD’ and YDS versus the square root of ionic strength.
The substances are HCl, LiCl, NaCl, and KCl. B has been

set to zero for these comparisons.

23

Figure 2-1 shows that the highly mobile H+ ion increases
the diffusion rate of HCl significantly and the highly
solvated Li+ ion decreases the diffusion rate of LiCl.

For KCl YSS is slightly larger than for NaCl, presumably
because K+ is larger than Na+. The upturn in all of the
YSS plots around 0.5M is due to the increase in their mean
molar activity coefficients near this concentration.

Figure 2-2 shows that with increasing ionic strength
YDD becomes less than YSS' There are two reasons for this.
For l-l electrolytes in neutral solutions there are no
cross Onsager coefficients in the YDD expression as there
are fer YSS' There are also no activity coefficient terms
for these YDD expressions. Thus the upturn around li/z =
0.5Ml/2 is absent. These effects combine to produce the
result that the ratio of charge to mass diffusion increases
with increasing ionic strength.

Figures 2-3 and 2-A show that YSD and YDS are generally
smaller than YSS and YDD‘ This is because they contain a
relatively large term which involves the cancellation of
diagonal Onsager coefficients. The cross coefficients,

YSD and YDS’ are negative only for HCl because it is the
only electrolyte shown for which £++ > 2__. There is no
upturn for YSD because it does not involve activity co-

efficient terms. For KCl YSD and YDS fall significantly

below the curves for the other electrolytes because of the

almost complete cancellation of t++ and t__

2A

Figure 2.1. Dependence of YSS on ionic strength at B
equals zero. (a) HCl, (b) LiCl, (c) NaCl,
(d) KCl.

 

 

25

. H . N mazwfim

Ectsmiyésm
ad 06 3

 

 

 

 

 

 

 

 

9N

 

 

SSA

17__()|, I‘SZwO /

26

Figure 2.2. Dependence of YDD on ionic strength at B
equals zero. (a) HCl, (b) LiCl, (c) NaCl,
(d) KCl.

 

27

m.m opswfim

 

 

 

 

Azctemlyécm
9— 06 md 16 90
i i I _ n n l
I! . n :13
//
illlttlllllr, o Illllllltlllllllllltttr-
o /h.. 3

 

 

 

‘-

db

‘-
-

q—
d)—

d-

 

00/.

9.0L |'_szulo /

28

Figure 2.3. Dependence of YSD on ionic strength at B
equals zero. (a) HCl, (b) LiCl, (0) NaCl,
(d) KCl.

 

 

m.m onswfie

Azctem2.o:tem
0.0 ad 3

 

29

c-(r-

 

 

 

 

 

 

nu)—
unh-
d)-
d)—
du—

31
0.»:
A
a:
nu
. u /
o N . 3
w
7v
8
_
o..- L
o.
n:
0.0

o.—

30

Figure 2.A. Dependence of YDS on ionic strength at B
equals zero. (a) HCl, (b) L101, (0) NaCl,
(d) KCl.

 

:.m madman

AzctemiyEcm
Q6 06 *6

 

-r
aL-g
1.

fi—

31

 

 

 

db

 

qr—
‘-
.4-

 

 

0;:

90

Q-

80A.

9-0L l'_Szl.l..|0 /

32

Figures 2-5 and 2-6 show the effect of B on YSS and
YDS respectively-for NaCl - H2O solutions. As B increases
788 decreases because 2__ > £++ for NaCl. When 2++ > t__
the Opposite is true. Qualitatively this means that if the
more mobile ion experiences a significantly higher activity
coefficient at lower concentrations, there will be less of
a tendency for the electrolyte to diffuse in its own con-
centration gradient. As B increases YDS decreases sharply
with increasing ionic strength, not only for NaCl, but for
all electrolytes. Because of the complex nature between
the coupling of S and D there appears to be no simple ex-
planation for this behavior of yDS.

The coefficients, YSE and YDE’ as expressed in Equation
(2-A0), are shown in Figures 2-7 and 2-8, respectively, for
HCl, LiCl, NaCl, and KCl. The YDE/S curves are higher than
the corresponding ySE/S curves because the YDE expression
involves a sum of diagonal Onsager coefficients, whereas
the YSE expression involves a difference. This reflects
the fact that a charged particle will move faster in an
electric field than a neutral particle.

The complete, formal Onsager thermodynamic treatment of
charged systems has not yet been accomplished. Presumably
the Onsager flux equations corresponding to Equations (2-A1)
and (2-A2) would contain coefficients similar to those
discussed in this section. The phenomenological approach

here is sufficient for the applications described in Chapters
5 and 6.

33

Figure 2.5. Dependence of YSS on ionic strength for NaCl.

(a) B = -0.5 M‘l, (b) B = o, (c) B = 0.5 M'l.

3A

m.m opswfim

Azctcmiyéem
ad Gd *6

od

 

——

db

u

 

cub

 

 

--

‘-

db

 

N;

H—

w;

9—

Q—

fi—

88"

9-0L L‘SZLUO /

35

Figure 2.6. Dependence of YDS on ionic strength for NaCl.

(a) B = -o.5 m‘l, (b) B = o, (c) B = 0.5 M‘l.

36

m.m osswfim

Azctcmioztem

90 *6

«d

 

 

 

 

 

 

37

Figure 2.7. Dependence of YSE/S on ionic strength. (a)

HCl, (b) LiCl, (C) NaCl, ((31) KC1.

38

e.m opswam

Azctsmiyésm
0.0 *0

 

 

 

 

 

II 0.7-

 

 

 

0.—

38A)

0L0L L-'°wL—52w° / ‘5’

39

Figure 2.8. Dependence of YDE/S on ionic strength. (a)

HCl, (b) LiCl, (c) NaCl, ((1) KCl.

 

 

A0

. w.m opswfim

133310.25
GO Q0 to

_ . _

 

p n
_ . _ . _ . _

I Q—

! 9N

 

 

 

Q5

30/6)

OLOL Lnow ,5ng / (S/

CHAPTER 3

ION TRANSPORT BOUNDARY CONDITIONS

A. Introduction

Any mathematical or physical description of a system
is, at best, limited without a description of the events
taking place at the system's boundaries. For the one-
dimensional systems discussed in Chapters 5 and 6 there
are two boundaries. At each boundary it is necessary to
make a statement about the concentration and/or gradient of
the concentration of each ion present in the system.. Since
these boundaries correspond to the intersection of two
phases, as shown in Figure 3-1, it is convenient to employ

kinetic boundary conditions of the form,

1 = kifcio ‘ kibcié ’ (3'1)

for the flux Ji of ion 1 through an interface of thickness
6 at each of the interfaces. The forward and reverse rate
constants are kif and k1b (the directions are indicated in
Figure 3-1), and Ci is the ion concentration. Boundary
conditions like these have been used before by Chang and

Jaffé [1952], Brumleve and Buck [1978], and have been

A1

A2

 

 

P7
H
H)
W
‘F
27
I-‘°
C‘

 

 

Phase L Phase R

Interfacial Region

Figure 3.1. Interfacial region at a phase boundary.

A3

discussed by Perram [1980]. The purpose of this chapter
is to determine the assumptions required for validity of
boundary conditions of this type. This discussion fol-

lows closely that of Horne, Leckey, and Perram [1981].

B. Derivation of Kinetic Boundary Conditions

At equilibrium, the fluxes vanish and the Nernst
Distribution Law (Prigogine and Defay [195A]) obtains for
the ratio of activities ai in the absence of an electric

field,

_ _ _ e-
aio/aio ' (Cicyio/cioyio) ‘ KiN ‘ eXp(‘AGi/RT) (3'2)

where yi is the molarity based activity coefficient, K1N

is the Nernst Distribution Constant, and AGG- the standard

1’
Gibbs free energy difference, is the difference in standard

chemical potential across the interface,
8’ 11.3- 8 . (3-3)

AGi g 6 - 1110

If an electrical potential difference, A¢ = ¢6 - c is

O,
maintained across the interface, then Equation (3-2) be-

comes

(aid/aio) = KiN exp(-ziF¢/RT) . (3-A)

AA

Equations (3-2) to (3-A) are readily derived from the

equilibrium requirement that ”id = U10: with the electro-
chemical potential given by

.9.
pi = “i + RT 2n a1 + ziFo . (3-5)

If Ji = 0 in Equation (3-1), then

(Cid/cio) = (kif/kib) ° (3'6)

ThuS, by Equation (3-2), (kif/kib) = Kf for ideal solutions
at equilibrium in the absence of an electric field and, by

Equation (3-A), (kif/kib) = K exp(-ziFA¢/RT) for ideal

iN
solutions in the presence of an electric field.

If cross terms are neglected, or if only one ion is
moving through the interface, the flux of an ion in one

dimension is given by

-Ji = tii(3ui/3X) , (3-7)
where iii is the Onsager coefficient and pi is the electro-

chemical potential. The electrochemical potential gradient

is, by Equation (3-5),

 

___= i i+ZFdl+——£ (3'8)

A5

with

d an yi
Pi = 1 + ——————— . (3-9)

d in C1

Although Equation (3-8) appears quite different from Equa-
tions (2-7) and (2-8), the only real change is the inclu-
sion of the general standard state term, (dug/dx) in Equa-
tion (3-8). This is necessary because here the ions are
passing from one phase to another, so the standard state
must change .

Substitution of Equation (3-8) into Equation (3-7)
yields

J = D dci + ziFDiCi Q9 + D ci 1 (3-10)
1 1 dx RT dx RT dx ’

 

where
D = _ ’ (3‘11)
and where activity coefficient terms have been neglected.

Equation (3-10), rearranged, is a differential equa-

tion for c1,

——+————=-——, (3-12)

A6
with
A1 = of + ztib . (3-13)

If A1 is linear in x and if J1 and D1 are constants across

the interface, then the solution of Equation (3-12) is

xAA xAA

_ i GRT 1
Ci - cioexp(- m) + J1 ———D1AAi[exp(- ——ORT) -1] , (3-1’4)
with
-e-
= + .. -
AAi Aci ziFdo . _ (3 15)

The interfacial boundary condition, Equation (3-1),
is obtained by putting x =(5in Equation (3-1A) and solving

for J . This results in

 

i
D AA AA AA
_ i i i i -l
kif - _6RT [exp(- WJEl-exph RT” 3
D A (3-16)
AA A
_ i i i -l
kin ‘ “(S-RT [1 ' exp(“ RT )1
The ratio of kif to kib is
kif Au? + ziFAq)
k1b - exp( RT ) , (3-17)

“7

as required by Equations (3-A) and (3-6) for ideal solu-
tions.

If Z1 and Ac are positive and if ziFAo > Aug; then kif <
kib' This is as expected since a positive ion would be
moving into a region of higher electric potential as it
moves in the positive x direction. Similarly, if Z1 is
negative, Ac is positive, and ziFA¢ > Aug; then kif > kib'
The other term in the exponent of Equation (3-17), Aug.
can be approximated by the not very accurate Born equation,

as Perram [1980] has done.

z2F2
'0 i 1 1
An = (——---- , (3-18)
i 231 SR eL

where ai is the radius of the ion and ER and EL are the di-
electric coefficients of the two phases. For Aug'> |ziF¢|,
the Born equation predicts a larger rate constant corres-
ponding to transfer to the higher dielectric phase. This
is also as expected, since the higher dielectric will
"stabilize" a charged particle more. Clearly the parameter
6, the thickness of the interface, is not readily accessible.
It may be possible to estimate it from experiments which

determine kif’ k and the distribution coefficient KiN

ib
for the ion between the two phases.
Thus, Equation (3-1) is strictly valid if: (i) the

neglect of cross terms in the flux equations is valid,

A8

(ii) A? and ¢ are linear functions of x, and (iii) J1

and Di are constants across the interface.

C. Applications of Kinetic Boundary Conditions

One of the main features of Equation (3-1) is that it
can be used for a wide variety of problems. It has been
used in program IONFLO to increase the program's versatility.
Without it, it would be necessary to rewrite large sections
of code for each different set of boundary conditions.

For liquid-liquid Junctions, discussed in Chapter 5
of this dissertation,the rate constants are set equal to
zero because there is no flux of ions through the ends of
the cell. For double layer relaxation studies, described
in Chapter 6, the rate constants are set equal to zero at
the blocking electrode since by definition no ions cross
this boundary. At the reversible electrode the rate con-
stants are set equal to infinity and the bathing concentra-
tions are set equal to the equilibrium system concentrations.

As mentioned in the introduction to this chapter, this
type of boundary condition formulatiOn.has been used pre-
viously by other workers. Chang and Jaffé [1952] used
similar equations to study the effect of increasing rate
constants at the electrode on the conductance. They did
this for alternating current studies. Brumlene and Buck
[1978] have used these boundary conditions to study num-

erically a variety of effects. Among other systems, they

“9

have studied charge steps at a blocked interface and mobile
site ion exchange membranes. For the blocked interface
simulation, they set the rate constants at the blocking
electrode to zero, and observed the effect of current
reversal on the ion distributions near the electrode. For
ion exchange membranes they considered three ionic species,
one trapped in the membrane by zero values for its transfer
rate constants, the other species initially on either side
of the membrane with finite rate constants. This simula-
tion was used to determine selectivity properties of ion
exchange membranes. Steele [1980] is using these boundary
conditions to simulate the steady-state performance of
various types of solid batteries. Perram [1980] has had
considerable success in describing biological membrane

transport with these boundary conditions.

CHAPTER A

NUMERICAL FORMULATION OF THE TRANSPORT EQUATIONS

A. General Finite Difference Approach for Parabolic Partial

 

Differential Equations

The transport equations and boundary conditions des-
cribed in Chapters 2 and 3, respectively, form a system of
nonlinear parabolic partial differential equations. Such
systems, in general, are not soluble by exact analytical
methods. To get an accurate solution it is necessary to
resort to numerical techniques, which, although approximate,
can be used to get as accurate a solution as desired. The
numerical technique of choice to solve the system of equa-
tions of interest here is the finite difference method.

The general presentation given in this section is based

on Smith [1978] and Rosenberg [1969]. A review of numeri-
cal simulations applied to ion transport is included in
Section C of this chapter.

To employ the finite difference method, continuous
variables are replaced by discrete ones, and partial deriva-
tives are replaced by finite differences. The resulting
equations are algebraic rather than differential. To

obtain discrete variables, the time-space solution domain

50

51

is replaced by a discrete grid as shown in Figure A-l.
The solution is calculated only at the grid points. Any
dependent variable U can be specified in time and space
with its appropriate indices Ui,n
A finite difference approximation to the first and
second space derivatives of U at constant time is obtained

and U in Taylor's series about the

by expanding Ui+l,n i-l,n
point Ui,n’ and taking the sum (for the second derivative)
and difference (for the first derivative). If terms con-
taining (Ax)3 and higher powers of Ax are omitted, the
results are

(3U/3x)i,n = (U - Ui_l,n)/Ax , (A-l)

i+1,n

(BU/8X)i,n = (Ui+l’n - 2Ui,n + Ui_1,n)/(AX)2 . (A-2)
A finite difference approximation for the time deriva-
tive of U can be made by combining Taylor's series expan-
sions about any of the points +, o , or o as indicated
in Figure A-l. The choice is mainly a matter of computa-
tional efficiency. The use of point + (forward difference
method) produces the simplest set of algebraic equations
to solve, but suffers from a severe drawback. ’For this
method to be stable the size of time steps, At, relative
to the size of space steps, Ax, is very restricted. The

use of point 0 (Crank-Nicholsen method) results in the most

52

1 IT+'I o o o o o -+ o o o
.1
t f] O O o I (o C) O o o
l1-I o o o o o o o o o o
H l i+l
)(-——->

Figure A.1. Uniform grid mesh for finite difference
equations.

53

accurate time derivative approximation, but also produces
the most complicated algebraic equations. A Taylor's
series about point 0 (backward difference method) produces
a time derivative analog that is unconditionally stable,
but not as accurate or as cumbersome to use as the Crank-
Nicholson method. After truncating terms containing

(At)2 and higher powers of At (actually the (At)2 term
vanishes for the Crank-Nicholson method), the time deriva-

tive analogs become,

Forward difference (aU/at)l,r1 = (Ul,n+l-Ui,n)/At ,'
(4-3)

Crank-Nicholson (BU/at)i,n+l/2=(Ui,n+l'Ui,n)/At’
(A-A)

Backward difference (BU/at)i,n+l=(Ui,n+l-U1,n)/At'
(14-5)

Since these time and space derivative analogs can be
applied at every interior grid point for each dependent
variable, there will be a total of Jx(R-2) unknowns for
each time row n, and Jx(R-2) angebraic equations. Here
J is the number of dependent variables in the original
partial differential equation system. Similar treatment
of the boundary conditions produces a total of JxR equa-

tions and JxR unknowns for each time row. If the original

5A

partial differential equations are nonlinear, the system
of algebraic equations will be so also. It then requires
some sort of iteration procedure (Fletcher-Powell, Newton-
Raphson, egg.) to solve them. The right hand sides of the
equations in any case are formed from the solution to the
system of equations at the previous time row. Starting
from the initial conditions, it is then possible to step
through time, calculating the solution to the partial dif-
ferential equation system at each grid point, one row at

a time.

B. Difference Equations for Ion Transport

The finite difference method described in Section A of
this chapter can be applied to the system of electrolyte
transport equations and boundary conditions detailed in

Chapters 2 and 3 with several modifications.

1. Space and Time Grids

It is well known that near charged interfaces quanti-
ties such as concentration and electric field change rapidly
over small space dimensions (<10"5 cm). To reduce trunca-
tion error in the Taylor's series expansions, it is neces-
sary to increase the number of space points in these
regions. Since it is still necessary to transverse the

entire solution domain With grid points, it is necessary

55

to use a space grid like that depicted in Figure A-2.

This means that in regions where Axi and Axi-l are not
equal, the finite difference space analogs are more com-
plicated than previously stated. Equations (A-1) and (A-2)

then'become,

2

_ 2 2
2
- (Axi) Ui=l,n}/[AxiAX1-1(AX1+AXl-l)J , (A-6)
(82U/3X2) = [2Ax U -2(Ax +Ax )U
i,n i-l i+l,n i 1-1 i,n
+ 2Aini_l,nJ/[AXiAxi_l(Axi+Axi_l)1 . (A-7)

\
‘o

Although McAfee [1976] has theoretically develOped fairly
general procedures for optimizing nonuniform space grids,
a less cumbersome empirical approach is used here.

Since the dependent variables change slower as steady-
state is approached, it is useful to employ a nonuniform
time grid to reduce the number of time steps necessary to
reach a steady-state or equilibrium. A typical time grid
of this sort is also shown in Figure A-2.

Choosing the appropriate time derivative analog re-
quires consideration of the relative computation times,

complexity of the computations, and ease of finding and

56

n+loooooo o. o o 0000
t [1000... o o o 00.0
I
n—[oooooo o o o 0000
I l-I l i+l

Figure A.2. Nonuniform grid mesh for improved accuracy
near boundaries and for more efficient ap-
proach to steady-state equilibrium.

57

eliminating errors. It has been found empirically, that
the backward difference analog performs best for this

system of equations.

2. Application to Ion Transport Derivatives

Equation (2-37) can be converted to finite difference
form by a direct application of Equations (A-5) and (A-6).
Equations (2-3A) and (2-35) contain derivatives of the
form (a/ax)[P(U,V)W] and (a/ax)[P(U,V)aU/ax], where W is
the electric field, U and V are the sum and difference
concentrations, and P is a transport coefficient, 788’
YSD’ YDS’ YDD’ YSE’ or YDE‘ The discrete forms of these
derivatives are best derived by repeated applications of
Equation (A-6). They become,

{(a/BX)[R(U,V)W]}i = P<U1’Vi){(Axi-l)2wi+1 +

[(Axi)2-(Axi_l)2]wi - (Axi)2Wi_l}/[AxiAxi_l(Axi+Ax1_l)1

+ Wi{(AXi_l)2P(Ui+1,Vi+l) + [(Axi)2 - (Axi_l)2]P(Ui,Vi)

" (Axi)2P(Ui_lavi_l)}/[AX1AX1_1(AX1+AX1_1)J , (ii-8)

58

{(a/ax)[P(U,V)3U/3X]}1 = {2AxiP(Ui+l/2’Vi+l/2)Ui+l

‘ [2AX1P(Ui+l/2’Vi+l/2) + 2Axi-lP(Ui-l/2’Vi-l/2)]Ui

+ 2AX1-1P(Ui—l/2’V1-1/2)Ui-l}/[AxiAxi-l(Axi+Ax1-l)J’

(4-9)

where the first application of Equation (A-6) to (3/ax)

[P(U,V) U/ax] has been over the interval x
)2

i-l/2 to xi+l/2

and terms of order (Ax - Axi-l have been neglected in

i
Equation (A-9). In practice this introduces little error
since grid points near interfaces are actually spaced at
equal increments while the dependent variables are changing
rapidly. From 20-30 grid points at each interface are
typically spaced at equal increments for an 80 grid point

network. 3

3. Boundary Conditions

Treatment of the first and last space points requires
special consideration. If the first space point is
labeled 1 and the last space point labeled R, direct ap-
plication of Equations (A-7) and (A-8) requires knowledge
of the dependent variables at points 0 and R+l and evalua-
tion of some of the transport coefficients at points 1/2

and R+l/2. Elimination of these so-called "image points",

59

which are outside the system boundaries, requires use of
the boundary conditions. Since application of Equation
(A-6) to the boundary conditions, Equation (3-1), also
generates dependent variables outside the system boundaries,
it is possible to eliminate these variables between the two
sets of equations. This can be illustrated by considering
Equation (2-3A) at the first grid point. In finite dif-
ference form, after neglecting the v0 term, Equation

(2-3A) becomes

T + T

1

2 + T3 + TA = Sl /At (A-IO)

,n

where

T1 = -[YSS(l/2)/(Axl)2]SO,n+l + [1/At + 783(1/2)/(Axl)2

Q

+ YSS(3/2)/(Axl)2131,n+l - [YSS(3/2)/(Axl)2]s2’n+1,
(A-ll)

T2 = -[YSD(l/2)/(Ax1)2]Do,n+l + [YSD(l/2)/(Axl)2

+ YSD(3/2)/(Axl)2]Dl,n+l - [YSD(3/2)/(Axl)2]D2,n+ls
(n-12)

T3 = -[YSE(l)/(2KAX1)1E2,n+l+[ySE(1)/(2KAxl)]E0,n+l,
(A-13)

60
Tu = -[ySE(2)/(2KAxl)-YSE(O)/(2KAxl)1E1,n+l, (A-lA)

and the fact that Axl = Ax2 has been used. PUV(UJ’n+1,
VJ,n+l) has been abbreviated as PUV(J), where J represents

the grid point at which P is to be evaluated.
The flux boundary conditions in finite difference form

at the first grid point are
' [YSS(l)/(2Axl)]so,n+l + [YSS(l)/(2Axl)]s2,n+l
- Emlyn/(2AM)1D,,n+1 + [YSD(1)/(2Ax1)]D2,n+1
+ [YSE(1)/K]El,n+l = (kaSl,n+l'kaSL) +

L
+ (kaDl,n+l ‘ kDfD ) ’ (“'15)

-[yDS(l)/(2AX1)]SO,n+l + [YDS(l)/(2AX1)]S2,n+l

_ [YDD(l)/(2Axl)]Do,n+l + [YDD(1)/(2Ax1)]D2,n+l

_ L
+ [YDE(l)/K]El’n+l - (kasl’n,l - kDfS )

L
+ (kaD1,n+l - kaD ) , (A-l6)

iMhere’ka, ka’ ka, and kDf are appropriate combinations

61

of the rate constants described in Chapter 3, and SL and
DL are the concentrations of S and D beyond the system's
boundary at the first grid point.

Since YSS’ YSD’ YDS’ and YDD are only weakly concen-
tration dependent it is a good approximation to evaluate
them at grid point 1 instead of grid point 1/2 in Equations
(A-11) and (A-12). This permits elimination of SO and D0
in terms T1 and T2 by Equation (A-15).

In finite different form, at grid point 1, Poisson's

equation is

(1/2Axl)E2,n+l = (1/2Ax1)EO’n+l = (22+Z_F/e)D1,n+1 _
(A-l7)

Thus E0 in Equation (A-l3) can be eliminated using
Equation (A-17).

Elimination of (0) in Equation (A-lA) is somewhat

YSE
more complex. This requires solving the finite difference
boundary condition equations, (A-15) and (A-l6), simul-
taneously for S0 and DO, and then substituting these values
into the expression for YSE' The result, along with the
entire set of finite difference equations, is listed in

Appendix C.

62

A. System in Matrix Form

The system of nonlinear algebraic equation is put in
suitable matrix form by placing all dependent variables
from time row n+1 on the left hand side of the equations
and all dependent variables from time row n on the right

hand side. The system can then be written in the form
AY = B , (4-18)

where the solution vector i contains all three variables,
S, D, and E, evaluated at every space point of time row
n+1. The vector B is generated from the solution to time
row n.' For the first time row B is generated from the
initial conditions. The matrix A is a function of Ax,
At, and the concentration dependent transport coefficients
evaluated atstime row n+1. These coefficients have been
fitted to functions of S and D from tables of Onsager co-
efficients and activity coefficient data compiled by Miller
[1966], Jalota 33 El. [1972], and Agnew g3 a1- [1978]. Of
course YSS and YDS still contain the adJustable parameter B.
Newton's iteration procedure (Appendix D) is employed
to solve Equation (A-l8). This procedure is used at each
time row, with the solution used to form B for the next
time row. This routine is repeated until the dependent

variables no longer change with time. At any time row the

electric field components of X may be integrated to calculate

63

the potential difference between any two points in the
system.

As can be seen from Equations (A-15) and (A-l6), the
rate constants for transport across the phase boundaries
enter the system through the boundary conditions. They
will appear in Equation (A-18) as part of B. Each may be
varied separately to simulate any type of interface trans—

fer consistent with Equation (3-1).

C. Computer Simulation Discussion

The finite difference procedure discussed in Section B
of this chapter has been coded as a FORTRAN program for
use on a CDC Cyber 750 computer. A general flow diagram
is shown in Figure A-3. A listing of the program, which
includes numerous comments about its use, is presented in
Appendix E. ‘For each iteration using Newton's method, the
A matrix of Equation (A-18) is inverted using a program
package contained within LINPACK [1979]. Since a computer
solution approach requires dimensionless variables, a
complete listing of those chosen for this system is as-
sembled.in Appendix B.

It is not surprising that computer programs similar
to this one have been developed before for electrolyte
transport. Schryer [1977] has noted that virtually any
classical problem can be solved numerically. The theoretical

details of finite difference methods have been available

6A

Fi ure A. I
g 3. Flow diagram for program IONFLO

65

IONFLO

 

 

 

 

SPACE
GRID

"HTWAL
STATE

 

 

 

 

 

 

 

TRANSPORT
COE.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

L BEGIN
TIME
LOOP
[ NEW TIME
[STEP SIZE
$ ASSEMBLE .
a "
f - l
ASSEMBLE
A
A V _
SOLVE
Ah 3 WRITE
OUTPUT
" CULALI: ATE
VOLTAGE

 

v I

Figure A.3

66

for some time, but only in the last 15 years have computers
become efficient enough to handle the solution of coupled
nonlinear partial differential equations routinely. Cohen
and Cooley [1965] were the first to simulate electrolyte
transport on a large scale in electrochemical systems. Al-
though they used simplified flux equations and omitted
activity coefficients, they recognized the advantage of
using the displacement current equation for simulations.
Gummel [196A] and DeMari [1965] numerically solved steady-
state and transient equations for charge distributions in
semiconductor Junction devices. DeMari used flux equations
similar to Cohen and Cooley's, but applied them to holes

and electrons. Kahn and Maycock [1966] also used simplified
flux equations to simulate numerically the electron diffusion
in polar crystals. Feldberg [1970] has simulated double
layer relaxation using the Nernst-Planck flux equations and
Poisson's equation instead of the displacement current
equation.

Recently the trend has been toward more general programs.
Newman [1973] has developed general programming techniques
for the solution of ordinary differential equations as-
sociated with the steady state of a number of electro-
chemical systems. Sandifer and Buck [1975] have developed
an algorithm for simulation of transient and alternating
current electrical properties of conducting membranes,

Junctions, and finite Galvanic cells. By far the most

67

general program has been written by Brumleve and Buck [1978].
They used simplified flux equations, but were the first to
employ numerically the generalized boundary conditions con-
ceived by Chang and Jaffe [1952] and used here. They used
a nonuniform space mesh and space derivative analogs of the
form defined in Equations (A-1) and (A-2), instead of the
more appropriate space derivatives defined in Equations
(A-6) and (A-7) for a nonuniform space mesh. Although this
seems incongruous, it does not cause serious errors if
Ax1 = Ax1+1 for regions near interfaces, where dependent
variables are changing rapidly. However, for systems with
extended double layer effects, for instance, their results
are subJect to error. Although their results have not
been compared to experiments quantitatively, they have
simulated a number of solid state, membrane transport, and
electrochemical systems, and have obtained qualitative
agreement with experiment.

The program described in this section contains all of
the essential features of the above programs, plus more ac-
curate flux equations, space derivations, and options for

the use of activity coefficients.

CHAPTER 5

LIQUID-LIQUID JUNCTIONS AND SINGLE-ION EFFECTS

A. Introduction

 

The significance of single-ion activity coefficients
has not diminished despite their persistent experimental
unattainability. Conway [1978] has listed several areas
where information about single-ion activity coefficients
and other individual solvated ion properties is required
for further interpretation of observed phenomena. Among
the more important areas are ion binding at proteins,
specificity of ion-solvent interactions, and kinetics of
electrode processes. Buck [1976] has emphasized the de-
pendence of interfacial potential on single-ion activities.
Frank [1963] and Goldberg and Frank [1972] have reviewed
the importance of single-ion activities and the apparent
futility, of their experimental determination. Following
Goldberg and Frank [1972] concentration cells with liquid
Junction have been investigated with the aim of determin-
ing the feasibility of measurement of single-ion activity
coefficients. In this chapter both analytical and numerical
results are presented which delineate the time intervals

during a concentration cell experiment in which single-ion

68

69

activity effects might be measured.

A concentration cell can be made by bringing into
contact two solutions of different electrolyte concentra-
tions and containing identical electrodes. The cell emf

(E = ¢R - ¢L) that results will depend on the electrode

cell
potentials and the charge distribution near the region of
contact between the two solutions. At t = 0, when the solu-
tions make contact, there is no separation of charge in the
liquid-liquid interface region. After the solutions make
contact the ions will begin to diffuse at varying rates.
This creates a charge distribution, which produces the
liquid Junction potential. If the Junction is formed

from two solutions containing the same binary electrolyte,
the Junction potential, and hence the cell potential, ob-
tain a constant value very rapidly. Eventually the cell
potential will decay to zero as the entire cell reaches
equilibrium. The plateau region begins when the liquid
Junction potential has risen to its maximum value and ends
when diffusion begins to decrease significantly the elec-
trolyte concentrations at the electrodes. During this time
period, the cell potential is constant.

The relationship between plateau region cell potentials
and activity coefficients, for concentration cells with
liquid Junctions, has been described theoretically since
Guggenheim [1929]. Bearman [195A] derived this relation-

ship from rigorous nonequilibrium thermodynamics. For the

70

plateau region, most electrochemistry texts such as Mac-
Innes [1939] and Bockris and Reddy [1977], state that the
cell potential for concentration cells containing a single
electrolyte is independent of the activity coefficients of
the individual ions and depends only on some mean activity
coefficient. Although the time regions before and after the
plateau region are not often discussed, MacGillivray [1970]
has used singular perturbation techniques to obtain analytic
solutions to the time-dependent Nernst-Planck equations for
ion flux across a membrane under voltage-clamp conditions.
The results presented in this chapter rely on the equations
developed in Chapter 2 of this dissertation and not on the
approximate Nernst-Planck equations. Horne [1981] has shown
that the Nernst-Planck equations are correct only at
infinite dilution.

Hickman [1970] used singular perturbation techniques
for investigating the Junction potential created by diffus-
ing electrolytes. His approach is not directly applicable
to electrochemical cells because he used unrealistic boun-
dary conditions. Goldberg and Frank [1972] and Chen and
Frank [1973] also attacked the question of concentration
cell potential dependence on single-ion activities. Their
approaches were entirely numerical, however, and they were
mainly concerned with the plateau time region. Goldberg

8

and Frank report early time behavior (<10' 5) for a par-

ticular KCl concentration cell which exhibits significant

71

cell potential dependence on single-ion activities up to

about 10’9

s, when significant effects vanish.

The principal goal here is to relate the liquid Junc-
tion potential as well as the full cell potential, to single-
ion activity coefficients in a time dependent analytical
equation. This has the advantage that the final equation
specifically reflects the effect that the input parameters
have on the final result. In this case the input parameters
are the transport coefficients defined by Equations (2-36).
These in turn are functions of the more fundamental parameters,
activity coefficients, Onsager coefficients, and ionic val-
ences.

After applying suitable approximations, Equations
(2-32), (2-3A), and (2-35) are solved for the cell po-
tential as a function of time and the single adJustable
parameter B. In keeping with the goal stated in the last
paragraph, it is desirable to determine whether there 18‘

a time range in which the cell potential depends strongly
enough on B to determine y+ and y_ experimentally.

The analytical expression for the liquid Junction

potential derived here has the form

EJ = tPS exp(-St) , (5-1)

for short times. Here EJ is the liquid Junction potential,

t is time, PS is the pre-exponential factor, and S is the

72

short-time time constant. PS depends on single-ion activity
coefficients, but S does not. During this time range the
full cell potential depends on single-ion activity co-
efficients. At about 10-85 the expression for the liquid

Junction potential becomes
EJ = PM a (5'2)

where PM depends on single-ion activity coefficients, but
not on time. In this time interval the single-ion activity
coefficients do not, however, appear in the expression-for
the full cell potential. At much longer times, the liquid

Junction potential is
EJ = PL exp(-Lt) 3 (5-3)

where PL depends on single-ion activity coefficients, but
L does not. Single-ion activity coefficients do not appear
in the expression for the full cell potential.

The results presented in this chapter confirm analyti-
cally that cell potential dependence on single-ion activity
8

coefficients vanishes very rapidly (<10- 3) for concentra-
tion cells operating with a single binary electrolyte.
Previous results are extended by quantitating this effect
for a series of chloride electrolytes. These results sug-

gest interesting and complicated interdependences between

73

electrical transport processes, on the one hand, and

single-ion activity coefficients, on the other.

B. Initial and Boundary Conditions

The initial conditions for Equations (2-3A) and (2-35)

for the liquid Junction problem are

S(x,o) = SE, D(x,o) = 0, ¢(x,o) = 0, -a/2<x<0,
(S-A)
S(x,o) = 8%, D(x,o) = 0, ¢(x,o) = 0, 0<x<a/2
where SE and SE are the electrolyte concentrations in the

two cell compartments at the beginning of the experiment
except for 2-2, 3-3, etc. electrolytes. The length of the
cell is "a" and the electrodes are at a/2 and -a/2.

The boundary conditions are

(aS/ax)

0
II

(BS/8x)

-a/2,t = a/2,t ’

O
I

(ab/ax) (ED/3x)é/2 t , (5-5)

-a/2,t =

(3¢/3x)_a/2,t = 0 = (3¢/3X)a/2,t

The gradients of S and D vanish at the electrodes because the

fluxes of S and D vanish there; the electric potential

gradient vanishes because the electrodes are not charged.

7A

C. Simplifying Approximations

For the systems of interest here (0.01 to 1.5 molar),
the valuescfi‘the diagonal Onsager coefficients, 2 and

++
11 2-J"1'cm"l-s"l and are

£__, are approximately 10' mol
approximately proportional to the ionic strength. The

off diagonal coefficient, 2+_ = 2_+, ranges in value from

0 to 25% of the diagonal coefficients and is approximately
proportional to the ionic strength to the 3/2 power. [MI
and |B| range up to about 1000 cm3-mol'l.

The coefficients YSS and YDD are dominated by the rela-
tively constant term (RT/2I(2++/c+) + (2__/c_)] and are
approximately 10.5 cmg-s-l in value. The coefficients
YSD and YDS are smaller, but are also dominated by the
similar constant term (RT/2)[&++/c+) - (t__/c_)]. The
coefficients YSE and'YDE are both proportional to concen-
tration and are approximately 109 s-1 in value.

Since the coefficients depend on concentration,
Equations (2-3A) and (2-35) are non-linear, as stated
in Chapter 2. Horne and Anderson [1970] and Ingle and
Horne [1973] solved similar systems of non-linear partial
differential equations with the help of perturbation ex-
pansions. The perturbation technique provides useful prac-
tical results only when the expansions converge rapidly,
which they do for relatively constant transport coefficients.

In any case, the zeroth order term of the perturbation pro-

vides the general form of the solution. The zeroth order

75

solution corresponds to solving Equations (2-3A) and (2-35)
with constant coefficients. It is the zeroth order solu-
tion presented in this chapter. As [32 - SE| decreases,
the results approach the exact results. Since the fpgg
of the analytic solution is of prime importance, the zeroth
order solution will suffice. In order to obtain explicit,
accurate numerical results, the computer program described
in Chapter A is used to determine ion concentrations and
cell potentials as a function of time. I

The solvent velocity terms in Equations (2-3A) and
(2-35) are neglected at this point since V0 is very small
when there is neither appreciable bulk flow nor a sig-
nificant density gradient.

These approximations, along with direct substitution
of Equation (2-32) into Equations (2-3A) and (2-35) result

in the pair of linear partial differential equations

(aS/at) Yss(328/3x2) + YSD(32D/8x2) - YSED ,
(5-6)

_ 2 2 2 2
(aD/at) - st(a 3/3x ) + YDD(3 D/ax ) - YDED

D. Analytical Solution

The solutions to Equations (5-6), subJect to conditions

(5-A) and (5-5), are (see Appendix F for details)

76

_ o o o o -1 m -l
s(x,t) - l/2(sR + SL) + (sR - sL)n égi [(2k-1)le

[P_kexp(-R_kt) + P exP(-R+kt)]Sin[(2k-1)UX/a]:

+k
D(x,t) = -2UYDS(S; - Sin.2 é:£[(2k-l)/Qk][exé(—R+kt)
-exp(-R_kt)]Sin[(2k-1)nx/a] , (5-7)
Where
R+k = l/2{yDE + (Yss+YDD)](2k-l)t/a12 - Qk} ,

R_k = l/2{YDE + (YSS+YDD)C(2k-1)n/a]2 + Qk} .

- 2
P_k - Qk - (YSS-YDD)[(2k-1)n/a] + YDE ’

"U
I

2
+k .. Qk + (YSS—YDD)](2k-1)TT/a] " YDE ’

= {[( - )2+A ][(2k-l)n/a]u+2 ( )
k Yss YDD YSDYDs YDB YDD‘Yss

+ 2stySE][(2k-l)w/a]2 +YDE2}1/2 . (5-8)

Explicit equations for c and c_ are easily found by com-

4..
bining Equations (5-7) according to Equations (2-31).

With the last of Equations (2-31), the charge density

77

is

p(x,t) = -unz+z_Fst(S§-s§)a'2 Agi [(2k-l)/Qk][exp(-R+kt)
- exp(-R_kt)]sin[(2k-l)nx/a] . (5-9)

By integrating Equation (2-32), the liquid Junction potential

is calculated,

_ o o -2 - i
EJ(t) = ¢R - ¢L — 2hKyDS(SR - sL)a ‘1' dx 3!;2 D(x,t)dx,

where K is given by Equation (2-33). The result is
- o o m k -l
EJ(t) - (A/n)KyDS(sR-SL) £2& (-1) [(2k-1)Qk] x
[exp(-R+kt) - exp(-R_kt)] . (5-10)

This expression may be summed, with the help of a computer,
to obtain EJ for any t.

For three particular time intervals Equation (5-10)
assumes simpler forms: (1) for t<(lIyDE)'l (242;: for t 5

10

10- s), to better than 1% accuracy,

0 O
EJ(t) = -tKYDS(SR - sL) exp(-l/2YDEt) . (5-11)

(ii) for t>(6a2/YSS) (i.e., for t3105 s for a 1 cm cell),

78
to better than 1% accuracy,
EJ(t) = -(u/w)K(vDS/YDE)(8§ - s§)exp<-R,lt> , (5-12)
with
R+1 E (w/a)2(YSS + You) = (n/a)2(a++ + c-_) . (5-13)

(iii) For (5/YDE)<t<(105a2/ SS) (i.e., for 10'9 sstgl s),
to better than 1% accuracy, the liquid Junction potential is

independent of time,

0

E = -K(YDS/YDE)(s§ - sL) (5-1A)

J

The expression for the lower limit of Equation (5-1A)
is particularly important. (5/YDE) is a measure of the on-
set of a constant liquid Junction potential. From Equation
(2-33), YDE is proportional to the conductance A, so (5/YDE)
= (SE/Ark) = 5tD. In the next section it will be shown
that YDE is also directly related to the time interval in
which the cell potential depends on single-ion activity co-
efficients.

From Appendix A, Equation (5-lA) becomes, for 1-1

electrolytes, for (5/YDE)<t<(105a2/YSS),

EJ = -(RT/icF)(sg-si)[(t+-t_)(1+Mlc) + BIC] . (5-15)

79

E. Full Cell Equation Results

The liquid Junction potential is, of course, not measur-
able by itself. Any measurement will also include poten-

tial differences at the two electrode-solution interfaces.

For definiteness and simplicity the full cell potential for
the following single electrolyte concentration cell is cal-

culated:

Ag(s)|AgCl(S)|MC1(SL,aq)IMC1(SR,aq)lAgCl(S)|Ag(S)
(A)
’The cell potential is

where
RT ’
EN = TTE£n(sR/SL) + An(y_R/y_L)l , (5-17)

and EJ is given by Equation (5-10). SR and SL are the
electrolyte concentrations at the right and left hand
electrodes. Strictly speaking the chloride ion concentra-
tions should be used, but SR/SL is equivalent to [Cl-JR/
[Cl-JL. y_R and y_L are the chloride ion activity co-
efficients at the right and left hand electrodes. It is
straightforward to put Equation (5-17) in terms of y1 and

B. By Equations (2-10) and (2-11),

80
E - BE[2n(S /S ) + 2n( / )1 ( 18
N ‘ ‘3' R L yiR yiL) ' B(SR ‘ SL ’ 5' )

where the fact that IC = S for cell A (1-1 electrolyte)
has been used. The results are not changed significantly
for other electrolytes.

Combining Equations (5-10), (5-16), and (5-18), and
simplifying the pre-exponential factors using the expres-
sions in Appendix A, results in the following full cell

potential equation,

ECELL = -%;[£n(SR/SL) + 2n(y R/y L) _ B(sR - SL)]

+ (A/n)KyDS(s§-sg) £1 (-1)k[(2k-1)Qk]-l[exp(-R+kt)

- exp(-R+kt)] . (5-19)

SR’ SL, yiR’ and yiL will depend on time only after dif-
fusion has begun to change the electrolyte concentrations
at the electrodes.
The implications of Equation (5-19) are clearest during
the time periods when Equations (5-11), (5-12), (5-13),
and (5-1A) are applicable. It can be shown that:
-l

(i) For t < (AYDE) , ECELL depends on B

(ii) For t > (5/YDE), E is independent of B

CELL

81

5 2
(iii) For (5/YDE) < t < (10 a /YSS), the cell po-
tential can be expressed, after considerable

algebra, in the particularly simple form

ECELL=(-2RTt+/F)Itn(sg/s§) + tuna/mu . <5-20>

where t+ is the transference number of the positive ion
(see Appendix A). This is the classical electrochemical
cell potential derived by Guggenheim [1929] assuming con-
stant transference numbers across the cell.

Table 5-1 shows a comparison between Equation (5-20)
and some experimental results for concentration cells
operating with liquid Junctions and moderate electrolyte
concentration differences. The differences are quite small
as long as the cell electrolyte concentration differences
do not become too large. Although these cell potential
comparisons are all in the constant cell potential time
interval, they give some idea of the error induced by the
linearization for longer and shorter times.

The next section describes the use of program IONFLO
to check the validity of the analytical solution for times

less than (5/YDE).

F. Numerical Solution

In an effort to validate the analytical results at

short times, program IONFLO was used to solve the system

82

Table 5-1. Comparison of Experimental Cell Potentials

 

 

 

(EEXP) for Cell (A) with Cell Potentials (ECALC)

From Equation (5-20).
Cell sg/mol-dm‘3 sL/mol-dm'3 ECALC/10'3-v BEXP/lo'3-v
KCl 0.10 0.040 21.16 21.17a
KCl 0.10 3.00 -80.37 —80.A2a
NaCl 0.10 0.080 _ A.13 - A.05b
NaCl 0.10 0.040 -l7.09 —16.8ub
HCl 0.10 0.078 9.9A 9.95a
HCl 0.10 0.0u0 35.39 36.21?
H01 0.10 0.005 103.u 118.8a

 

 

aShedlovsky and MacInnes [1937].

bBrown and MacInnes [1936].

83

of Equations (2.3A), (2.35), and (2.37) numerically. The
V0 term was omitted from Equations (2.3A) and (2.35) for
the reason discussed in Section C. For convenience the
electric potential was replaced by the electric field in
Equations (2.3A) and (2.35) according to Equation (2-17).
The current I was set to zero since it vanishes for this
experiment. Because the dependent variables S, D and E
change rapidly near the liquid-liquid interface (X = 0)
at short times, a non-uniform space mesh with a larger
number of grid points near the interface was used. The

7

grid points were 10- a units apart near the interface and

gradually expanded to lO'la units apart at the electrodes.
Ninety grid points were used. Time steps of lO’lOs were
used initially. The step size was increased as the cell
approached equilibrium.

Numerical stability was checked by comparison of results
obtained using different step sizes. The numbers of spa-
tial and temporal grid divisions were varied by factors of
3 and 10 respectively without change in result. For cell
(A) operating at concentrations of SR = 0.11 mol-dm-3
and SL = 0.09 mol dm'3, non-linear computer simulation
results and analytical solution results agreed to within
5% for all time ranges and all electrolytes. Since the
only approximation in the computer simulation is neglect

of the solvent velocity, 5% is a good estimate of the ac-

curacy of the analytical solution for cell (A) operating

8A

as indicated.

Figure 5-1 shows the numerically predicted charge den-
sity surface near the liquid-liquid interface for a NaCl
concentration cell with an assumed value for B of zero
(y+ = y_). Figure 5-2 shows the same cell at three dif-
ferent times. As time increases the maximum value of the
charge density decreases and moves away from the interface.

Figure 5-3 is the same cell at 10-8s for three values
of B. As can be seen, changing B changes the magnitude of
the curve, but does not change its shape. This corresponds
in the analytical solution to B appearing in the pre-ex-
ponential factor but not in the time constants.

'Figure 5-A is a graph of cell potential versus time
for cell (A) with M = Na, SR = 0.11 moi.dm‘3 and sL =
0.09 mol-dm-3. As a benchmark, the criterion that the
cell no longer depends on B when the ECELL versus time

3

curves for B = 1 dm -mol-1 and B = -1 dm3-mol'l differ

from the plateau region cell potential by less than 5%

was chosen. With this somewhat arbitrary criterion, Figure
5-A indicates that from a time of 1.6 x 10793 on, the cell
no longer depends on B. Figure 5-5 shows a graph of the

coalescence time (tc) versus YDE for several electrolytes.

Figure 5.1.

85

Predicted charge density surface near the
liquid-liquid interface at short time for a
NaCl concentration cell operating between
concentrations of 0.11 M and 0.09 M. B is
zero for this simulation (y+ = y_). The time
row at t = 0 has been omitted since it would
obscure the lower portion of the surface.

 

 

 

86

 

H.m ossmflm

 

 

 

 

 

 

 

 

 

87

Figure 5.2. Constant time planes from Figure 5.1. (a)
5 x 10'9s, (b) 1 x 10'8s, (c) 2 x 10‘8s.

88

m.m ensues

. sob.)
m o s N o m- T o- n-

 

 

 

 

 

Figure 5.3.

89

Charge density YE cell position for different
values of B at t = 1 x 10'8s for a NaCl con-
centration cell operating between concentrations
of 0.11 M and 0.09 M. (a) B = -1.00 M'1
(b) B = 0, (c) B = 1.00 m‘l.

3

90

m.m ossmfla

sob.)
o N-

c.

mwl

mwt

 

 

_ _

_

 

 

91

Figure 5.A. Cell potential ys time for cell (A) with M =

Na, SL = 0.11 M, and SR = 0.09 M. (a) B = 1.00
M'l, (b) B = 0, (c) B = -1.00 M‘l.

92

ON

0..

:.m oszmwm

mmIO_\~

0..

0.0

 

 

 

_

 

 

n.mI

m.m..

93

Figure 5.5. Coalescence time ys YDE for 0.1 M solutions,
with 5% criterion. (a) HCl, (b) BaCl2, (c)
KCl, (d) NaCl, (e) LiCl.

9A

 

l6- '

tc / I0"°s

 

l l l

 

 

I .
4 8 I 2 I 6
the / loss"|

Figure 5.5

24

95

G. Conclusions

The cell potential dependence on single-ion activity co-
efficients vanishes very rapidly(t < 10785) for the systems
investigated here. From Figure 5-5 it is clear that YDE
is a good indicator of this effect for an arbitrary elec-

trolyte. For the electrolytes studied, the relationship

tC = (1'9/YDB) = (1.9 /AwA) = 1.9TD (5-21)
approximates the coalescence time to within 5%.

These results can be interpreted as follows. Initially
in the interfacial region the electrical potential gradient
is zero and the activity gradient is infinite. Consequently
the positive and negative ions begin to diffuse from the
region of high activity to the region of low activity.
Since the aiJ are not equal, the diffusion begins at dif-
ferent rates for the two types of ions. This soon creates
a separation of charge and an electrical force that pulls
the positive and negative ions back together. Thus a
steady-state balance of diffusive and electrical forces
is established. Before this balance is established the
cell potential depends on the positive and negative ions
separately. After this the forces and hence the ions them-
selves are coupled in such a way as to remove their indi-
vidual identities. Buck [1976] has similarly described

the behavior of the Junction potential in terms of

96

differential ion motion.
This model is consistent with the expression for

yDE (the last of Equations (2_36)),

2 2 _
YDE = (HUF2/E)(Z+Z+++2Z+Z_£+_+z_2__) = (“fl/€)A = TD].

(5-22)

Large values of t++ or £__ (HCl for instance) indicate rap—
idly diffusing ions, large values of YDE and A, small TD’
and consequently a rapid vanishing of single-ion effects.
The same is true for large values of 2+ and z_. For BaC12,
the 2: factor offsets the relatively small values of 2++
and £__, which leads to relatively large values of YDE

and A, small TD, and a small coalescence time. In terms
of this model, the valence number of +2 for Ba++ creates

a larger electrical force sooner. This apparently over-

+ to diffuse slowly, and thus

whelms the tendency of Ba+
the coalescence time is relatively fast.

For large values of e, YDE is small, TD is large, and
the single-ion effect time is prolonged. This is consistent
with the idea that a large value of 5 indicates that the

positive and negative ions are relatively independent of

one another.

CHAPTER 6

DOUBLE LAYER RELAXATION AT BLOCKED ELECTRODES

A. Introduction

 

Double layer relaxation studies, both theoretical and
experimental, have interested electrochemists to various
degrees over the last 35 Years. Physicists and electrical
engineers have also been interested in double layer relax-
tion (also called space-charge decay), especially in semi-
conductors and solid electrolytes. Laser and Bard [1976]
note that for an intrinsic semiconductor the relaxation of
free charge carriers resembles that of the cations and
anions in the diffuse double layer at a metal/electrode
interface following charge inJection to the metal. This
is not surprising since the same thermodynamic and electro-
static equations govern both processes.

In an effort to study the kinetics of fast electrode
processes, electrochemists have often been concerned with
eliminating the double layer effects on electrode processes.
Barker [1959] suggested the use of short duration current
pulses as a means of eliminating the double layer charging
current. Despite some experimental success by Reinmuth

and Wilson [1962], Delahay and Aramata [1962], KooiJman

97

98

and Sluyters [1967], Weir and Enke [1967], Daum and Enke
[1967], and Yarnitzky and Anson [1970], van Leeuwen [1978]
has concluded that short duration pulse methods have not
achieved the success expected of them because they in
principle provide no more information than periodic tech-
niques. Anson, Martin, and Yarnitzky [1969] have conducted
one of the few experimental studies whose main goal was to
study nonequilibrium double layers in the absence of
electrode reactions. They studied the potential-time -
behavior of several electrolyte solutions following coulosta-
tic charge inJection into a hanging mercury drOp electrode.
Their main conclusion was that the thickness of the double
layer was not constant during rapid charge inJection, and
that transient potentials calculated on the basis of a
constant double layer thickness were in error.

The research reported in this chapter is a theoretical
thermodynamic study of systems similar to the ones investi-
gated by Anson e£_al. Although the kinetics of fast elec-
trode processes is not dealt with here, it is hOped that
a better understanding of double layers in the absence of
electrode reactions will lead to new ways to eliminate and/
or correct for them during electrode kinetic studies. More-
over, the equations and their solutions are intrinsically
of interest to both the macrOSCOpic and the molecular
theorist.

In an attempt to analyze the response to a coulombic

99

impulse applied to an electrode at equilibrium, Reinmuth
[1962], starting from absolute rate theory and Fick's law,
derived potential-time relaxation response expressions.
These expressions were intended to cover the range from
charge-transfer dominated processes to completely blocked
electrodes. For Reinmuth's results to be valid, however,
the solution must contain a large excess of supporting
electrolyte. As Francescetti and MacDonald [1970] and
Anson, Martin, and Yarnitzky [1969] have pointed out,
there are many electrochemical and semiconductor systems
where the use of supporting electrolytes is not feasible.
In this chapter relaxation of the double layer at com-
pletely blocked electrodes is considered in the context of
the theory presented in Chapters 2 and 3. In an approach
somewhat similar to that used in Chapter 5, both analytical
and numerical results are presented. Section B deals with
the equilibrium double layer in the presence of a charge
on the electrode. The solution for the equilibrium case
is in terms of a perturbation series with all activity
coefficient effects included. Section C details an approxi-
mate time dependent analytical solution using Laplace trans-
forms. Section D deals with numerical simulation results

using program IONFLO.

100

B. Equilibrium Solution

The equilibrium double layer has been studied in great
detail by chemists since the historic treatment of the
subJect by Gouy [1910] and Chapman [1913]. Grahame [19A7],
Delahay [1965], and Payne [1972] have presented extensive
summaries of the topic. The only thermodynamic treatment
which ineludes single-ion activity coefficients has been
attempted by Hall [1977]. However, he presents no
calculations. This section contains a systematic perturba-
tion solution to the double layer problem. All single-ion
effects are retained.

At equilibrium the fluxes of all ions vanish, and
therefore JS and JD, defined by Equations (2-A1) and (2-A2),-

do also. This leads to the system of equations,
3E
J = 0, J = 0, SE = KD , (6-1)

for the ion and field distributions in the system in the
presence of a charge on the electrode. Due to the presence
of activity coefficients, Equations (6-1) cannot be solved
exactly. Instead, a perturbation approach is used.

In terms of S, D, and E, Equations (6-1) become

YSS(dS/dx) + YSD(dD/dx) - YSEK‘lE = 0 -a/2:X:a/2a
(6-2)

101

YDS(dS/dx) + yDD(dD/dx) - YDEK-lE = 0 -a/2ix5a/2,
(5-3)

dE/dx = KD -a/2:x:a/2.
(6-A)

Immediate simplification results from simultaneously
eliminating the dD/dx term from Equation (6-2) and the

dS/dx term from Equation (6-3). The results are

dS/dx + (NS/H)E = 0 , (6—5)
dD/dx + (ND/H)E = 0 , (6-6)
where
Ns = YSEYDS ' YDEYSS’ N0 = YDEYSD ‘ YSBYDD ’
(6-7)

and H = K<YSSYDD ‘ 180103)

By Appendix A, the coefficients of E in Equations
(6-5) and (6-6) simplify considerably for 1-1 electrolytes.
For other electrolytes the essential results are the same.

The simplified equations for 1-1 electrolytes become

dS/dx + (F/RT)(LSDE) = 0 , (6-8)

102

 

dD/dx + (F/RT)(LDSE) = o , (6-9)
where
2 2
LS = ITMEIBB and LD = S+M<§ 'D ) . (6—10)
, S+MS -BDS

This simplification has been fruitful for at least two
reasons: (1) there are no transport coefficients, only
activity coefficients remaining in the equilibrium equa-
tions for S and D; (2) for ideal solutions, where M = 0
and B = 0, Equations (6-A), (6-8), and (6-9) yield, among
others, the Poisson-Boltzman equation (Moore [1972]).

The boundary conditions for the double layer relaxation

are

E = (An/E)Q x -a/2 , (6-11)

x = a/2 , (6-12)

where SB is the bulk concentration of S (=10), and Q is the
charge on the electrode placed at x = -a/2.

Before introducing the perturbation technique the
system is next transformed to dimensionless form by the

following substitutions,

s = 8/88,

e = (aF/RT)E,
k2 = (8HFaga2/ERT) ,
m = SBM ,

103

II
U)

D/S

B,

x/a ,

(-AnaF/ERT)Q,

(6-13)

These transformations are slightly different from those

used by IONFLO and are used here primarily for simplifica-

tion. In terms of these reduced variables, the system of

equations becomes

ds/dg + Asde = 0

d -
dd/dg +,A se - 0

2
de/dg + k d = 0
e = -q
s = l, d = 0

where
s = 1
A 1+ms-bd and

A

-1/25§:l/2 ,

‘1/25531/2 ,

-l/2:§:l/2 ,

d _

-1/2

1/2

s+m(s
s+ms

3

9

2-d2)
2--bds

(6-1A)

(6-15)

(6—16)

(6-17)

(6-18)

(6-19)

10A

The perturbation ordering is defined by

g = 1 + 681 + 6232 + 5353 + 0(5") ,
- . 2 3 A
d 6dl + 5 d2 + 6 d3 + 0(5 ) , (6-20)
A _ 2 3 A
e - 6el + 6 e2 + 6 e3 + 0(6 ) ,

where 6 is a bookkeeping parameter which characterizes

the departure of s from one and d from zero. When 5 = 1,
then B = s, d = d, and 8 = e. The parameter 5 can be
thought of as turning on the perturbation (the charge on the
electrode at 6 = -l/2) to the system. At 6 = 0 the system
is not perturbed and the system has the trivial solution

3 = l, d a 0, e = 0, for all values of E. At 6 = l the

full effect of the charge on the electrode is felt by the
system. In practice 6 is used to order subsequent manipula-
tions and is always set equal to one after the equations are
actually solved.

In order to effectively apply the perturbation order-
ing to the coefficients, A5 and Ad, they are first expanded
in Taylor's series about s = l and d = 0.

The reduced coefficients in Taylor's series form

are

s _ s s s l 2 s ‘s
A (s,d) - A0 + [dAd + (s-1)Asl + [5d Add + d(s-1)Asd

+ %(s-1)2A:S] + ... , (6-21)

105

d
A (s,d) = is

d d 1 2 d d
+ [dAd + (S-1)AS] + [2d Add + d(s-1)ASd

+ Asset-gs] + , <6-22>

where the derivatives in the Taylor's series are

AO = 1%5- 1% = l

A: = -m(A:)2 A: = 0

A38 = 2m2(A:)3 128 = 0 _ (6-23)
82d = -2mb(>.:)3 Aid = -mb(A:)2

A: = b(l:)2 A3 = bl:

Agd = 2b2(A:)3 )gd = [-2m(m+1)+2b2](A:)2

Substituting the expansions of s and d, Equations (6-20),
into the coefficient expansions, Equations (6-21) and (6-22),

and grouping according to powers of 6 yields

106

= S S
A A0 + 6{Assl + A: d1} + 6 2{AS s32 + Add2 +
1 s 2 1A8
2‘sssl I AsdS 1dl I 2 dddl } I 5 ”{A sS3I Add3 I
(6-2A)
s s l s 3
+ +
Asssls2 I Isd(sld2 d2 5 l) I Idddld2 6Issssl I
1 s 3 l s 2 M 2 2 A
6Iddddl I 2Issdsldl I 2Isddsldl} I 0(5 )’
d _ d d d
A - AO + 5(ASsl + Addl } + 6 ”(A Ss2 + Add2 +
1 d 2 d 1 d 2 3 d
2*sssl I Asdsldi I 2 dddl} I 5 {Ad s33 I Add3 I
(6-25)
d d d 1 d 3 '
Asssls2 + Asd(sld2 + d2sl) + Addd1d2 + EASSSSI +

ExgdddI I 2Igsd81dl I Bxgddsl I u

The perturbation approach outlined here is similar to
that used by Horne and Anderson [1970] for thermal dif-
fusion with one notable exception. The convergence of their
solution relied partially on the smallness of higher order
derivatives of their transport coefficients. Here no such
smallness is expected. Instead the convergence properties
of the solution rely on the magnitude of the charge on the
electrode, q. If q becomes large enough it is probable

that the solution, Equation (6-20),will diverge.

107

Application of Equations (6-20), (6-2“), and (6—25)
to the system of Equations (6-14) - (6-18) and subsequent
ordering of the terms according to the power of 5 yields

the following system:

dsn/dg + fn = o - 1/2 5 g g 1/2 , (6—26)
d2dn/dE2-k2dn+dhn/d€ = 0 - 1/2 1 g g 1/2 , (6-27)
den/d5 + k2dn = O - 1/2 3 g E 1/2 , (6-28)
sn = 0, an = o g = 1/2 , . (6—29)
en = -qn, ddn/dg = gn - hn g = -1/2 , (6-30)

n: 1,2,3’°°' 9

where a second derivative of Equation (6-15) has been taken
to permit substitution of Equation (6-16) into Equation

(6-15). The fn, gn, and hn are defined by

fl = 0

r2 = Agdlel (6-31)
_ s s s 2

f3 ‘ Ao(d191 I d2el) I Assldlei I Addlel ’

108

81 ' q

gn = O n = 2,3,... , (6-32)

h1 = O

h2 = (A: + Ag)slel + Agdlel

h3 - (Ag+kg)82el + (AS+A:)sle2+kg(dle2+d2el) + (6-33)
<x:.+xa>s.d.s + <A2+%x2.>sie1+%xfiddiel

Higher orders of h and f are omitted for brevity.

This system of equations can now be solved if equations
of order n are solved before equations of order n+1, and
if for a given n Equation (6-27) is solved before Equation
(6-28). As in any perturbation method (Bellman [l96u]),
there is no guarantee that Equations (6-20) will converge
for 6 = 1. In this case convergence is not expected to be
rapid for large charges on the electrode. Since this same
system is solved numerically in the next section, the per-
turbation solution is presented here only through third
order.

Through second order s, and first order d and e, the

solutions are

 

 

s 2
IO q 1 h2[k(€ 1/2)J (6 3”)
S‘I'S= SI'l " a -
l 2 2K2 cosh2(k)
_ 2g
d1 - k cosh(k) sinh [k(€-l/2)] , (6-35)

109

el = 33§%TE7 cosh[k(£-l/2)] . (6-36)

At 298K for a 0.1 M bulk concentration and cell length
of 1 cm, k is about 1x107. This means that, for instance,
-exp(k) is an extremely good approximation to 2sinh (k).
With this type of approximation for the hyperbolic func-

tions, the results through third order are

s' = l + (qu2/2k2)exp[-2k(£+l/2)J +
. (6-37)
(qu3/3k3){exp(-2k(€+l/2)]-2 expE-3k(€+l/2)]},

d' = -(q/k)exp[-k(g+1/2)] + (qu2/3k2) x

{2exp[—2k(£+l/2)J-eXPE-k(€+l/2)J}-

(q3/k3){(3/16)[A:+2(A§)2+Agd]exp[-3k(g+1/2)J (6-38)

_ %(Ag)2exp[-2k(a+1/2)J-(1/16)EA§-%#<A§)2 +

AngexpE-k(£+l/2)]} ,

110

e' —lep[-k(€+1/2)J + (qu2/3k>{exp[-2k(a+1/2)J

-exp[-k(€+l/2)]} - (q3/k2){<1/16)£A:+2(A§)2 +

(6-39)
Ad ]exp[-3k(g+l/2)] — 3(Ad)2ex [-2k(g+l/2)] -
dd 9 d p

(1/16)[A: - §§<A§>2 + AgdlexpE-k(a+1/2>J} ,

where the prime indicates the solutions are valid for k >>

The identity

A (6-40)

Om
CLO.

has been used to simplify the 8' equation.
The convergence properties of the s' and d' series are
best evaluated at g = -l/2. At a = -l/2, the series for

s' and d' become

3' 1 + (qu2/2k2) - (ASqu3/3k3) , (6-A1)

d!

-(q/k) + (qu2/3k2) - (q3/8k3)(l - f%Ag) , (6-u2)

where the identity

d 2

111

I has been used in the d' equation.

Since an additional power of (q/k) appears in each
progressive term of the series for s' and d', it is a good
measure of the rate of convergence of the series. In

terms of dimensioned variables (q/k) is
(q/k) = (2'ITQ‘2/6RTSB)I/2 . (6-uu)

Thus, small charges on the electrode and high concentra-
tions favor rapid convergence. Small charges on the elec-
trode indicate the system is not perturbed far from its
original state, so the solution is well represented by the
first two or three terms in the series. High concentra-
tions indicate that for a given charge on the electrode the
relative change in concentration at a given point in solu-

tion is less than for low concentrations. If, for example,

Q is 10-6 C/cm2 and SB is 0.1 M, then (q/k) is 0.29, and
s and d are accurately represented by s' and d'. If, on
the other hand, Q is 10"5 C/cm2 and S is 10-3 M, then

B
(q/k) is 291, and the series are apparently diverging.

Figures (6-1) and (6-2) are plots of (s'-l) and d',
as defined by Equations (6-37) and (6-38), versus distance
from the electrode, showing the relative contributions of

the terms in the series. The graphs are for a 0.1 M

6

NaCl solution with an electrode charge of 2x10- C/cm2,

and b assumed to be —O.l. This corresponds to an activity ,

112

Figure 6.1. Contributions to s' from 32 and s3 for a
0.1 M NaCl solution with an electrode charge
of 2 x 10.6 C/cm2. b = -O.l; (A) s', (B)
s2, (C) 53.

113

H.@ mpswfim

..o_>~\.-3
0N m. 0. n 0

d A d

 

 

KWO.OI
00.0

000
9.0

nNd
mmd

9.0

 

no.0

 

 

 

-mmd

13)

ssamun/ (I

11A

Figure 6.2. Contributions to d' from d1, d2, and d3 for
a 0.1 M Nagl solution with an electrode charge
of 2 x 10‘ C/cm2. b = -o.1; (A) d', (B)
d1, (C) d2, (D) d3.

m.w mpzuflh

scion».

ov mm 0» mm 8. m. o. n o.

. a q a _ u _ .

 

115

 

 

 

 

 

 

 

 

ssamun/p

116

coefficient for Na+ of 0.71 and for C1' of 0.86. Two fea—
tures of the graphs are readily apparent: (1) Even though
(q/k) is greater than one (1.08) for these conditions,
both series appear to be converging due to the numerical

constants in the denominators of the perturbation terms.

(2) The contributions from the terms containing Ag and
Igd’ which depend on single-ion activity coefficients, are

very small.

Since the values of s' and d' decrease so rapidly as
distance from the electrode increases, a more practical
experimental variable is the potential difference. The
potential difference in reduced form due to the equilibrium
double layer is calculated from the electric field using

Poissons equation in integral form,

1/
Aw = fe'dfi, (6415)
-l/2
where
Aw = (F/RT)A¢ , (6-46)

and A¢ is the potential difference due to the double layer.

The usual negative sign in Equation (6-A5) has been omitted

to conform with standard electrochemical convention.
Through third order the reduced potential difference

is

117

Aw = awl + Aw2 + aw3 , (6-u7)
where
Awl = -q/k , (6-u8)
Aw2 = -qu2/6k2 , (6-49)
aw3 = q3EA: + Agd + 4(Ag)21/2Ak3 . (6-50)

For the electrolyte solution whose d' profile is shown
in Figure (6-2), the corresponding potential difference
calculated using Equations (6-H5) - (6-A9) is -27.06 mV.

I If b is assumed to be +0.1 (This corresponds to switching
the activity coefficients from yNaI = .71 and yCl‘ = .86
sto yNaI a .86 and yCl’ = .71.) instead of -0.1, the cal-
culated potential difference is -27.1u mV. Thus, accord-
ing to these results the potential difference is the same
to 3 significant figures for this solution. The potential
difference seems insensitive to whether yNaI > yCl‘ or
yCl‘ > yNaI'

Although Hall [1978] has proposed that single-ion ac-
tivity coefficients can be estimated from equilibrium
measurements on electrolyte solutions at charged electrodes,
when the electrode charge is small and specific adsorption

is absent, it is difficult to confirm this from these

118

results. Exceedingly large or small single-ion activity
coefficients must be postulated for the terms containing
them to be significant.

With Equation (6-A7) it is now possible to compare the
results presented here with those of Gouy [1910] and Chap-
man [1913]. They solved the Poisson-Boltzmann equation
(Philip and Woodling [1970], Chen [1971], Olivares and
McQuarrie [1975]) for the boundary conditions used in this
chapter. Although their results only partially correspond
with experrmental observation, they provide a logical
starting point for comparison. In terms of the reduced
variables defined by Equation (6-13), they predict the

following charge-potential relationship

q = —2 k sinh(A¢/2) . (6-51)

Assuming for the present comparison that A: = l, Agd =

0, and A3 = 0 (i.e., m = b = 0, which is valid at infinite

dilution), Equation (6—A7) becomes
Aw = -(q/k) + (q3/2uk3) . (6—52)

In order to compare this equation with the Gouy-Chap-
man result, Equation (6-51) is expanded in a Taylor's series
about Aw = 0. (This series always converges.) The result

is

119

l 5
1920(AW ) + ...] 0 (6.53)

 

q = -kEAw + §fi(aw)3 +
Inversion yields
Aw: -(q/k) +' (q3/2uk3) - (9q5/1920k5) + 0(q7) . (6-5A)

Thus Equation (6-52) is identical to Equation (6-5H)
through the second term in the series. There is little
doubt that neglecting activity coefficients in the perturba-
tion expansion produces a result completely equivalent to
the Gouy-Chapman result. Another way of saying this is that
the thermodynamic solution derived in this chapter is
equivalent to the Gouy-Chapman result at infinite dilution.
It is not surprising that the Gouy-Chapman result does not
predict experimental results well at high concentrations
and large electrode charges. At high concentrations A:
begins to deviate significantly from one, and for large
electrode charges higher order terms in the perturbation
series, which almost certainly contain stronger activity
coefficient dependence, become important. Although the
failure of the Gouy-Chapman theory is generally attributed
to its neglect of ions adsorbed on the electrode, it seems
likely that part of its failure is due to the neglect of

activity coefficients.

120

C. Time Dependent Solution

In an effort to obtain the form of the analytical solu-
tionfku‘double layer relaxation, the partial differential
equation for D, Equation (2-35) has been Solved after assum-

ing that: (l) v is very small; (2) YDS << YDD3 (3) YDD

O

and YDE are constant. Neglect of the v0 term is valid
because there is no appreciable bulk flow during the relaxa-
tion. The error introduced by neglecting the YDS term
compared to the YDD term varies from system to system. As
can be seen from Figures 2.2 and 2.A, YDS is about 1% of

YDD for KCl and about 60% of YDD for HCl. Although for a
given electrolyte YDD varies only about 10% up to 1 M,

YDE is much more concentration dependent. Figure 2-6 shows
it to be nearly linear in ionic strength up to l H. Al-
though assuming YDE to be constant is a severe approximation,
the results derived here will be compared to more exact

computer simulation results in the next section. With

these approximations, Equation (2-35) becomes

(ED/3t) = YDD(82D/3x2) - YDED - a/2 : x g a/2
(6-55)
Since there is no flux at the blocking electrode JD
is zero there. Applying the above discussed approximations
to Equation (2-42) at zero flux yields the boundary condi-

tion at x = -a/2,

121
(ED/3X) = -(YDE/2FYDD)Q(’G) X = '3/2 , (6‘56)

where Poisson's equation has been used to put the boundary
condition in terms of the charge on the electrode, Q.

At x = a/2, the reversible electrode location, the
solution has its bulk prOperties, so the boundary condi-

tion there is
D = 0 x = a/2 . (6-57)

Initially the solution is not charged, so the initial

condition is
D(x,t) = 0 -a/2 i x i a/2 . (6-58)

Using a Laplace transform technique, as described in Ap-

pendix G, the solution of this system of equations is

D( t) (/ 1(2) f ”BRUNEI
x, = k aFw n 2 2
n=0 <-1> dj‘ Q(t-U /YDE){exp[-U
-<pfi/U2)J-expE-U2-<r§/U2)J dU}. (6-59)

where

pn = §<An+1-2x/a), In = §<un+3-2x/a) , (6-60)

122

and k is defined in Equation (6-13). Here k is actually

(YDE/YDDAl/B’ but this is equivalent to the k defined in

Equation (6-13) because YDD and YDE are assumed constant.
Since charges are typically injected onto electrodes by

discharging a capacitor, a practical functional form of

Q(t) is
Q = Qotl-exp(-at)] , (6-61)

where a is the time constant for the charge injection and
00 is the maximum charge on the electrode during the charging
period. Using this functional form for Q, Equation (6-59)

integrates to
D(x,t) = (-kQ°/MaF) z (-l)n{-exp(r )erfc(r /'r+'r) 4-
n=0 n n

exp(-rn(erfc(rn/T-T)+exp(pn)ePfC(Dn/T+T) -

exp(~pn)erfc(pn/T-T)-(l/w)exp(-at)[-exp(wrn)erfc(rn/T+wt)

+exp(-mrn)erfc(rn/r-wT)+exp(wpn)erfc(pn/T+w1) - (6-62)
exp(-mpn)erfc(pn/T-wr)} ,
where
T = (YDEt)I/2, w = <39§l9>1/2, <6-63)

123

and the complementary error function, erfc, is defined by

1/2 w -v2
erfC(Z) = (2/1r ) f e dv . (6-614)
z

For Equation (6-8) to be valid a must be less than YDE'
In practice, this does not cause a problem since a is rarely
larger than 107 s'1 and YDE is greater than 107 for most
systems down to concentrations as low as 5 x 10'“ M.
If a is greater than YDE the solution to Equation (6-59)
takes on a different form, but nothing dramatic happens
physically. It merely means that the charge has been applied
to the electrode so rapidly that the ions have not had time
to move far, until the electrode has been completely charged.
In the other extreme, if the capacitor slowly charges the
electrode, then a is small, w approaches one, and the solu-
tion moves through a series of equilibrium double layers.

Although this solution appears a bit cumbersome at
first, it has several advantages over the Fourier series
solutions derived in Chapter 5. Principally, it converges
very rapidly for small times. It also takes on simpler forms
for long times and for regions near the blocking electrode.
As t goes to infinity, the exp(-at) portion of the equation
goes to zero, and the remaining portion of the equation
can be summed. The result is Equation (6-35), the first
order solution for d. At x = -a/2 and all values of t, the

solution for k >> 1 (k is much greater than one for all

12“

systems except possibly membranes.) reduces to

D(-a/2,t) = (kQO/2aF)[erf(T) - w-lexp(-at)erf(wr)] , (6-65)
where the error function, erf, is defined by
z 2
erf<z> = (2/nl/2) f e‘V dv . (6-66)
0
Figure 6-3 shows the dependence on a of D,at x = -a/2

as a function of time predicted by Equation (6-65) for a
particular NaCl solution. As the electrode charging time
constant increases, D at x = -a/2 reaches its equilibrium
value much faster, as long as a is less than YDE' The
difference in shape between the error function and expon-
ential function is also evident from Figure 6-3.

Figure 6-U shows the concentration dependence of the
relaxation of D at -a/2 for an electrode charging time
constant of 107 8‘1. Lower concentrations produce longer
relaxation times and smaller charge densities at the elec-
trode. Physically, the ions must travel further to reach

the electrode for dilute solutions, so the relaxation

time is longer.
For regions near the electrode (x near -a/2) and k >> 1

it is possible to approximate Equation (6-62) accurately by

the following expression

Figure 6.3.

125

Dependence on a of D at x = -a/2 as a func-
tion of time as predicted by Equation (6-65)
2

for a particular NaCl solution. SB = 10- M,

Q0 = 10'6 C/cm2, (A) a = 3 x 107 s'I, (B)

a = 1 x 107 s‘l, (C) d = 6 x 106 s'l.

 

‘1')

m . m cLZHrTm

s.0..0\«
0.0 0.? 0.? 00 .00 0M 0.N 0.. 0.. 0.0 0.0

0N 0

 

 

 

 

 

Figure 6.4.

127

Concentration dependence of the relaxation of

D at x = -a/2 for an electrode charging time

constant of 107 s"1 as predicted by Equation
(6-65). (A) sB = u x 10‘2 M, (B) sB = 2 x
10‘2 M, (c) s = 1 x 10'2 M.

B

\

 

z.m mpzwfih

 

 

 

sbro:
0.0 0...? 0...? 0..m 0.0 0..N 0..N 0“. 0... 0.20 000.0
10.0
.. 0..
mm
W. a»
8 0 M
m 4 ..ON WI
u...
.mN 0.
g
.00

 

 

L .
:0
r0

129
D(z,t) = -(kQO/AaF){exp(kz/a)erfc[(kz/2ar) +
TJ-exp(-kz/a)erfc[(kz/2aT)-T] -
(6-67)
w-Iexp(-at)[exp(sz/a)erfc[(kZ/2at) +

wt]-exp(-kmz/a)erfc[(kz/2aT)-wr]]} ,

where z is the distance from the electrode and is related

to x by

z x + a/2 . (6-68)

Figure (6-5) is a graph of D versus z/a as predicted by
Equation (6-67) for three different times. For this com-
parison the charge on the electrode at equilibrium is--10'7
C/cm2, a is 6 x 106 s-I, and YDE is 6 x 107 8. At 10.7 s
the system is essentially at equilibrium.

The potential difference as a function of time due to
the double layer is calculated from Equation (6-67) by
integrating Equation (2-32), Poisson's equation. The

result is

O
A¢(t) = (aQ /ek){l-exp(-YDEt) + -;BE:;—

(6-69)

[eXp(-dt)-exp(-YDEt)]} ,

130

Figure 6.5. Effect of time on D versus z/a profiles as
predicted by Equation (6-67). QO = -10’7
C/cm2, a = 6 x 106 3'1, = 6 x 1078, (A)

E
10-7s, (B) 10-85, (C) lO-Bs.

 

131

m.© opzwfim

obSoE

 

 

 

0. 0.
V m
#0: 33.1113 - IOU-V0

0.
I0

Q
(0

 

 

132

where again a must be less than YDE' As t approaches

infinity Equation (6-69) becomes

A¢<w> = (aQO/kE) . (6-70)

This is equivalent to Equation (6-48), the first order
equilibrium potential difference.

If a << YDE’ then Equation (6-69) becomes

A¢(t) = (aQO/ek)£l-exp(-dt)] , (6—71)

1
E'
ionic distribution. The change in potential difference with

for t >> yB This situation indicates a virtual equilibrium
time is due entirely to the change in the charge on the
electrode.

A more useful situation is when ak >> YDE’ This is
true for most systems,except for membranes, extremely di-
lute solutions, or slowly charged electrodes. For ak

>> YDE’ Equation (6-69) becomes

A¢(t) = (aQO/e){k-I[l-8Xp(-YDEt)1 +
(6-72)

a(YDE-a)-I[exp(-at)-eXp(-YDEt)1}

In the next section this equation is compared to some

computer and experimental results.

133

D. Computer Simulation Results and Discussion

Program IONFLQ has been used to simulate the same type
of systems as discussed in Sections A and B of this chapter.
The goals in the simulations have been to: (1) delineate
the range of applicability of the analytical results of
this chapter, (2) determine the effects of single-ion
activity coefficients on double layer relaxations, (3)
compare the theory presented here with experimental results.

Figure 6-6 is a charge density difference plot at 10"6 3

generated from two 0.1 M NaCl double layer relaxation

simulations. The electrode charging time constant was

6 x lO6 3.1 and the charge on the electrode was -10'7

C/cm2. The only difference in the simulations was the

1

value of B. For B = -l M- (YCI‘ > YNaI) the charge

density is slightly greater near the electrode than it is

1 simulation. Beyond about 5 x 10'8

1

for the B = 1 M-
cm the charge density for the B = l M- simulation becomes
larger. Physically this can be understood in reference to
the activity coefficient curves for NaCl shown in Figure
6—7. Since the charge on the electrode is negative, the
Na+ ions will tend to increase their concentration at the
electrode relative to the bulk. If B is negative, then,
from Figure 6-7, the relative change in the activity co-

efficient for Na+ between the region at the electrode and

the bulk is greater than for positive B. (i.e., The slope

 

of curve 2 is steeper than curve 1 at 0.1 M.) Thus, there

Figure 6.6.

13M

Computer simulation of charge density dif-
ference versus distance from electrode at
10‘65 for two 0.1 M NaCl solutions. a = 6
x 106 s‘1
1

B = -1 M- and p* corresponds to B = 1 M"1

, Q0 = -10-7 C/cm2 0 corresponds to

 

135

00

w.m ohswfim

 

 

 

 

to. . Eo\x .

o... 0% Qw. o“. 80...

AS 0.0
- o._

. on

- on

9.

 

 

136

Figure 6.7. Effect of B on Na+ and C1- activity coef-

_ -1
ficients. (1) YNaI if B - l M and YCl’

1

if B = -1 M- , (2) Y, for NaCl as measured

by Brown and MacInnes [1936], (3) y01_ if

1 -1

B = 1 M’ and YNaI if B = -1 M

 

 

137

0N0 0&0 0...0 1.0 N..10

m.m omzmfim
2qu
_.o 00.0 moo v0.0 «0.0

00.0 .

 

 

 

 

 

 

l

00

0.0

CD
0
swmamaoo Kumov

 

 

0..

138

is less of a tendency for Na+ to diffuse to the bulk, where
it must increaSe its activity coefficient more than it would
for negative B. In simpler terms the more positive B is,
then the more Na+ enjoys the environment at the negative
electrode. The opposite is true for Cl‘.

Despite the fact that the charge density profiles are
different for these two simulations, the potential difference
across the double layer is the same for both of them to 5
significant figures. This result appears to be general for
all electrolytes at all times during the relaxation process.
Although different values of B produce different charge
density distributions, the potential difference varies by
‘a virtually insignificant amount.

Several simulations of relaxations due to extremely
rapid electrode charging were made for 0.1 M solutions of
NaCl, KCl, and HCl solutions. The time constant for the
electrode charging process was set at 1011 5'1, 00 was 10"7
C/cm2, and B was varied from 1 M21 to -1 M-1. None of these
showed any potential difference dependence on B, even for
very short times (<10-9 s). This is quite different from
the liquid junction potential behavior at short times dis-
cussed in Chapter 5. The differential ion motion is ap-
parently absent here.

Several computer simulations have been made in an
effort to assess the limitations of Equation (6-72).

Tables 6.1—6.3 contain potential differences from selected

139

 

 

 

Table 6.1. Comparison of Equation (6—72) and Computer
Simulation for 0.1 M KCl. A¢* Igdicates Equa-
tion (6-72) Results. a = 6 x 10 s‘1 Q° =
10-7 C/cm2. a = 1 cm.

t/S'10-7 A¢/V A¢*/V A¢/A¢*

0.125 ”3.9 “3.3 1.02
1.25 22.9 22.0 1.0“
2.50 10.8 10.“ 1.0“
5.00 2.42 2.32 1.0“
7.50 0.54 0.52 1.0“
10.0 0.122 0.117 1.0”
-2I -2 +

1A.l 1.32x10 1.13x10 1.16
-3I -3 +

18.3 2.AAx10 2.18x10 1.12
25.1 1.u2x10’3 1.u0x1o'3 1.01
36.2 1.39x10"3 1.u0x10'3 0.99

 

 

+Time steps too large for the computer simulations at these

time rows.

140

 

 

 

Table 6.2. Comparison of Equation (6-72) and Computer
Simulation for 0.1 M NaCl. 00* Indicates
Equation (6-72) Results. 0 = 6 x 1063-1.
0° = 10"7 C/cm2. a = 1 cm.

t/s-10’ A0 /v A¢*/v A¢/A¢*
0.125 51.8 51.6 1.00
0.625 40.2 38.2 1.05
2.50 13.0 12.“ 1.05
5.00 2.91 2.76 1.05
7.50 0.647 0.617 1.05

10.0 0.146 .139 1.05

_2+ -2 +
12.7 3.2lx10 2.72X10 1.18

-3I -3 +
16.0 6.18x10 5.15x10 1.20
21.3 1.63::10‘3 1.54.:10’3 1.06
30.0 1.40;:10"3 1.391(10‘3 1.01

 

 

+Time steps too large for the computer simulations at
these time rows.

1A1

Table 6.3. Comparison of Equation (6-72) and Computer
Simulation for 0.1 M HCl. A¢* Indicates Equa-

 

 

 

tion (6-72) results. a = 6x106 s-I. Q0 =
10'7 C/cm2. a = 1 cm. I
t/s-10'7 A¢/v A¢*/v A¢/A¢*
0.1 ' 1A.5 16.A 0.88
0.5 11.6 11.u _ 1.02
1.0 8.60 8.97 1.02
2.0 u.72 4.65 1.02
u.0 1.u2 1.u0 1.01
7.5 0.176 0.172 1.02
10.0 4.021(10"2 3.83x10‘2 1.05
14.8 3.921(10‘3Jr 3.5ux10‘3 1.10I
19.0 1.61::10'3 1.56::10‘3 1.03
26.0 1.39x10”3 1.391(10"3 , 1.00

 

 

1.
Time steps too large for the computer simulations at
these time rows.

142

time_rows from these simulations for 0.1 M solutions of KCl,
NaCl, and HCl, along with the corresponding results from
Equation (6-72). There are relatively large deviations
between the analytical and computer results for KCl and

NaCl after 10-6

s and before equilibrium because the time
step sizes were increased too rapidly during that part of
the simulation. It is interesting to note that the devia-
tions are about the same for all 3 electrolytes. This is
somewhat surprising since one of the approximations which
lead to Equation (6-72) was that YDS << YDD' At 0.1 M

the ratios of YDS/YDD for HCl, NaCl, and KCl are -0.59,
'0.21, and 0.018, respectively. Closer examination of the
simulations, however, reveals that d2s/dx2, the term which
multiplies YDS in Equation (2-35) is much smaller than
d2D/dx2 for all three electrolytes. Presumably at higher
electrode charges this would not be so, and the effect of
omitting the st term would be more pronounced.

These same simulations have been repeated at 10'3 M
for the same three electrolytes. The results for selected
time rows are listed in Tables 6.4 - 6.6, along with the
results computed from Equation (6-72). Again the computer
results are somewhat too high after 10-6 s and before equi-
librium because the time steps have been set too large
there. The relative deviations between the analytical and

computer derived potential differences is slightly larger

here than at 0.1 M. This is because the relative change

143

 

 

 

Table 6.4. Comparison of Equation (6-72) and Computer
Simulation for 10"3 M KCl. A¢* Indicates Equa-
tion (6.72) Results. 0 = 6x106 s'l. Q0 =
10'7 C/cm2. a = 1 cm.

t/s-lO"7 A¢/v A¢*/v A¢/A¢*

0.250 1u19 1555 0.91
1.25 2268 2288 0.99
2.50 1293 1241 1.00
5.00 299 283 1.06
7.50 66.9 63.1 1.06
10.0 15.0 1u.1 1.06
1u.1 1.111+ 1.20 1.17+
18.3 0.11Mr 0.107 1.35I
25.1 1.671110"2 1.63x10'2 1.02
36.2 . 1.1401(10'2 1.39'x10"2 1.01

 

 

fTime steps too large for the computer simulations for
these time rows.

144

 

 

 

Table 6.5. Comparison of Equation (6-72) and Computer
Simulation for 10"3 M NaCl. A¢* Indicates
Equation (6.72) Results. 0 = 6x106 3‘1.
00 = 10"7 C/cm2. a = 1 cm.
t/s-10“7 A¢/v A¢*/v Ad/A¢*
0.625 2502 2617 0.96
1.25 2628 2643 0.99
2.50 1594 1524 1.05
5.00 384 357 1.08
7.50 85.9 79.7 1.08
10.0 19.2 17.8 1.08
1u.1 1.82I 1.5a 1.18I
18.3 0.18144r 0.136 1.35+
25.1 1.77x10‘2 1.59110‘2 1.10
3.62 1.40x10-2 1.39x10‘2 1.01

 

 

+Time steps too large for the‘computer simulations at
these time rows.

145

 

 

 

Table 6.6. Comparison of Equation (6-72) and Computer
Simulation for 10"3 M HCl. A¢* Indicates
Equation (6-72) Results. 0 = 6x106 3'1.
00 = 10;7 C/cm2. a = 1 cm.
t/s-10"7 Ac/v A¢*/v A0/A¢*
0.500 1052 1096' 0.96
1.00 875 864 1.01
1.50 657 6A3 1.02
3.00 268 261 1.03
5.00 80.7 78.7 1.03
7.50 18.0 17.6 1.03
10.0 4.03 3.93 1.03
1A.8 0.276I 0.23u 1.18I
19.0 3.70x10‘ 3.16x10’2 1.17I
30.9 1.39x10' 1.39x10'2 1.00

 

 

1ITime steps too large for the computer simulations at
these time rows.

146

in YDE throughout the 10'3 M solutions is greater than it
is for the 0.1 M solutions. (Assumption 3 of Section C
of this Chapter is that YDE is constant throughout the
solution.) The potential differences for the 10"3 M

simulations are larger than those for the 10'1

M simulations
because the relaxation is slower and because the equilibrium

charge distribution is more diffuse.

E. Comparison with Previous Results

 

Anson 32 El: [1969] have experimentally measured poten-
tial differences as a function of time for double layer re-
laxations at dropping mercury electrodes. In a related study
Newman [1973] has analytically solved the Nernst-Planck equa-
tions for the boundary conditions discussed in this chapter.
Besides the previously discussed thermodynamic assumptions
inherent in the Nernst-Planck equations, Newman has made
several other approximations in the process of solving
them. He assumes that the charge on the electrode versus
time is always of the form shown in Figure 6-8. Moreover,
he assumes that while the electrode is being charged the
ions move only due to migration in the electric field, and
that after charging is complete, the ions move only due to
diffusion.

Figures 6-9 and 6-10 are comparisons among the ex-

periment of Anson £3 a;., Newman's analytical predictions,

147

Figure 6.8. Electrode charge versus time profile assumed
by Newman [1973] for analytical solutions to
the Nernst-Planck equations.

 

148

w.w mpswflm
.

0
“

 

 

 

 

 

Figure 6—9.

149

Comparison of the experiments of Anson e: §l°
[1969] (A), simulations using IONFLO (B),
and Newman's [1973] analytical predictions
(C) for the relaxation potential of HCl at a
blocked electrode. a = 6 x 106 8.1, 00 =
10-5 C/cm2. Time axis is logarithmic.

 

150

m.m opzwflm

 

 

 

 

 

 

3:.
0.0.. 0.0.. 0.0... 0..... ON...
. _ 4 _ . 4 N0

0
I 10.0
m

1V0

lllllllmlllll/llll
I 10.0

_ _ .

 

A/(IOH) 0V

Figure 6.10.

151

Comparison of the experiments of Anson 23 MM.
[1969] (A), simulations using IONFLO (B),

and Newman's [1973] analytical predictions
(C) for the relaxation potential of KCl at

a blocked electrode. a = 6 x 106 8.1, Q0

= 10'5 C/cm2. Time axis is logarithmic.

 

152

oa.c cpzwwm

 

 

 

 

 

 

C. c.
00... OK... 0.0 .. 0.0: 0.0. ..
. _ _ _ N0
I /
.. 0.0
1. ..v.0
/
1 0 .0
L _ . .

 

A/(IOXMV

153

and the simulations using IONFLO, for HCl and KCl relaxa-
tions respectively. Time has been plotted logarithmically
so that all results could be displayed on the same graph.
Anson's transients are much longer in duration than either
of the theoretical predictions. However, he reports varia-
tions of a factor of 10 in the lifetime of the transients
due to the shape of the capillary used in preparing the
mercury drops. Because of the large charge on the elec-
trode (10-5'C/cm2), the computer simulation results are
subject to truncation error for the region near the elec-
trode. The only definite agreement among the results is
that the potential difference due to KCl decays slower than

that due to HCl.

CHAPTER 7

SUGGESTIONS FOR FUTURE WORK
This chapter contains two sections: (A) suggestions
for future thermodynamic studies of ion transport processes,

and (B) suggestions for upgrading program IONFLO for future

use and for validating it in its present form.

A. Future Thermodynamic Studies

 

(1) It does not seem possible at present to make ob-
servations in the first 10.8 seconds of a liquid junction
potential experiment. If the balancing of forces hypothesis
is correct, it may be possible to slow the balancing time
by initially imposing a temperature gradient on the system.
This can be done by using constant temperature reservoirs at
the two electrodes. A measurable potential difference will
result (Agar and Breck [1956]). At steady state the concen-
tration distributions will be approximately linear in
space (Horne and Anderson [1970]). If the reservoirs are
switched to a third, mean temperature reservoir, the con-
centration gradients will begin to decay to zero, as will
the junction potential and the entire cell potential. The

decay time will be on the order of the relaxation time of the

154

155

thermal gradient in the cell (>1s). Besides the equations
discussed in Chapter 2, the additional equation needed to
analyze this experiment is the conservation of energy equa-
tion (Hasse [1969]). This system of equations is quite com-
plex, and a numerical solution may be required to analyze
the data.

(2) For many e1ectrode/electrolyte/electrode systems
the double layer relaxation boundary conditions are more
complicated than those used in Chapter 6. For these systems
it is necessary to account for specific adsorption of ions
at the electrode. To analyze these systems thermodynamically
requires adding another region to the solution domain with
its corresponding boundary conditions. This new region
would be too narrow to allow for conventional ion transport
within it, so the main equation to solve there would be
Poisson's equation. The solution to the equations in this
region and in the previously discussed region would be
joined by kinetic boundary conditions as described in Chapter
3. This is also a complex problem and would probably re-
quire a numerical solution.

(3) Turq 22.21: [1979], using radioactive tracers,
have reported the appearance of transient concentration
gradients from initially homogeneous species in multi-
component electrolyte solutions. This experiment could be
analyzed analytically by a direct application of the matrix

method described in Appendix F. Other more complicated

156

electrochemical and biological multicomponent systems could
be analyzed numerically. Miller's [1967] ternary system
data and n-component projections would be useful in these
extensions.

(4) By necessity, Onsager coefficients have been cal-
culated from measurements on essentially electrically neutral
solutions. Thus, it is not known how they depend on 0+
and C_ separately. Pikal [1971] and others have postulated
various models from molecular considerations. It should be
possible to test these models for macroscopic manifesta-
tions using program IONFLO.

(5) It is generally accepted that the dielectric co-
efficient is not a constant, particularly near charged inter—
faces where the water dipoles become oriented in preferential
directions. Although the electrostatic implications of a
non-constant dielectric coefficient are more fundamental
(Perram and Barber [1974]) than the subject of this dis-
sertation, the equations necessary to test various dielectric
models are contained in this dissertation. Again, a numeri-

cal approach may be necessary.

B. Program Improvement Suggestions

(1) IONFLO has only been used for the boundary condi-
tions discussed in Chapters 5 and 6, and for a narrow range
of electrode charges and concentrations. It should be

validated for other applications such as membrane transport

157

and solid electrolyte transport.

(2) At present IONFLO reads in expansion coefficients
of Onsager coefficients and activity coefficient deriva-
tives. These are used in function subroutines to form the_
Y transport coefficients. These functions are called col-
lectively about 5000 times per iteration per time row by
subroutine COEMAT and RHSVN. It would be more efficient to
form expansions of the Y transport coefficients as functions
of S and D directly from the input data at the beginning of
the simulation. These direct expansions of the transport
coefficients could then be transferred to streamlined func-
tion subroutines for use by COEMAT and RHSVN.

(3) IONFLO should be redimensioned to permit simula-
tions of multicomponent systems. Temperature could be
treated as a component with some minor adjustments for en-
thalpies of mixing.

(4) The dielectric coefficient should be programmed to
have the capability to depend on the dependent variables.

(5) The space grid and boundary conditions should be

expanded to permit use of multi-region domains.

APPENDIX A

TRANSPORT COEFFICIENTS

The y in Equations (2-3A) and (2-35) are

13

788 (RT/2z+z_)(ae - bg)

YSD = (RT/22+Z_)(-af - bh)
(A-l)
YDS = (RT/2z+z_)(-ce + dg)
YDD = (RT/2z+z_)(cf + dh) ,
with

m
ll
0‘
u

(Z_/c+)(z+£++-Z_£_+), (z+/c_)(z+2+_-Z_£__).

O
I
Q.
I

- (z_/c+)(z+2+++z_£_+), - (z+/C_)(z+2+_+z_£__),
e=1-%(z+/z_Xz+-z_)2c+(M+B), f=1-%(z+/z_)(zi-zf)c+(M+B),

g=1-%(Z_/z+)(z+-z_)2c_<M-B), h=l+%(z_/z+)(zf—zf)c_(M-B).

For l-l electrolytes,.
158

159

_ -1 _ '1
a1 — -c+ (£+++2_+) , bl - C_ (£+_+£__) 9
-1 -1
Cl = _c+ (2++-g_+) , d1 = c_ (2+_-2__) ,
(A-3)
el = 1 + c+(M+B) , f1 = l,
81 = 1 + c_(M-B), hl = 1

and in that case

YSS = %RT{[(£+++2_+)/C+1+[(£+_+£__)/c_] + (2+,+21+_ +
+ 2__) M + (2,1 - 2__)B} .
YSD = - %RT{[(£++ + £_+)/c+] _ [(g+_ + g-_)/c_]} ,
(A-4)
YDS = - §RT{[(£,+ - £_+)/c+1 + [(2,_ - 2__)/c_] +
+ (2++ - 2__)M + (2++ - 2z+_ + 2__)B} ,

YDD = %RT{[(£++ - 2_+)/C+] - [(£+_ - £_-)/c_]}

In terms of transference numbers ti and specific con-

ductance A, c and d in Equations (A-2) become

0 = (z_/z+c+)(vt+/F2) d = (Z+/a_c_)(At_/F2), (A-5)

160

and for 1-1 electrolytes,

%(ART/F2){[(t+/C+) - (t_/c_)J + (t+ - t-)M + B}
(A-6)

YDs

YDD = %(1RT/F2[(t+/c+) + (t_/c_)]

For electrically-neutral regions in a 1-1 electrolyte,

c = c = S = Ic’ and

1 2
-§(ART/Ic)[(t+ - t_)(l + IC) + 1c] ,
(A-7)

_ 1 2
‘ §(ART/Ic) :
and

(YDS/YDE) = (eRT/BNI§)[(t+ - t_)(l + Ic) + BIC]
(A-8)

APPENDIX B

DIMENSIONLESS VARIABLES USED BY PROGRAM "IONFLO"

Dimensionless variables are signified by a line above

M is the mean value of S in the system at the

the symbol. S
beginning of the simulation. Other variables are defined

in the body of the dissertation.

EN = (SM)
YgE I YDE(SM’OA
s = (s - SM)/sM
B = D/sM
E = EME/(2F2+z_aSM)
i=x/a
T = 2Fz+z_YDgEt/eM

ZI'= EMI/(4aF2zizEYgE)
_ M

70 - Q/(2F2+z_aS )

ET'= EMA¢/(2a2FZ+Z_SM)

161

m
rd
U)

9.153%:

3

M
e/e

M
E Yss

M

=€YDD

M
8 YSD

M
E YDS

M
e kij

162

2 M
/(2Fz+z_a YDE)
2 M
/(2Fz+z_a YDE)
2 M
/(2FZ+Z_a YDE)
2 M

/(2FZ+Z_a YDE)

M
/(2Fz+z_ayDE)

TRANSPORT EQUATIONS IN FINITE DIFFERENCE FORM

In order to render the transport equations into a

APPENDIX C

compact finite difference form, it is useful to make the

following definitions.

the system.

U and V

1/(AXi + axi_1)
2/[AX1(AX1 + AX1_1)]
2/[AX1_1(AX1 + AXi_1)]

YUV(si,n+l"Di,n+l)

D

YUV(Si+l/2,n+l, i+l/2,n+l)

YUV<si-l/2,n+l’ Di-l/2,n+l)

(S )

YUE i+l,n+l’ Di+l,n+l

D

YUE(Si=l,n+l’ i-l,n+l)

(S D

i,n+l)

YUE i,n+l’

can be either S or D.

163_

R is the number of space points in

(C-l)

S-Equation

l6u

2 5 1 5 R-l

- - +
(MYSS)Si-1,n+l + (1/At + MYSS + PYSS)Si,n+l

+
-(PYSS)S

+
PYSD)Di,n+l ‘

+

D-Equation

1+1,

2

+ -
n+1 ‘ (MYSD)D1-1,n+1 + (MYSD +

(c-2)
+ ..
(PYSD)D1+1,n+i + (AYSE)Ei-l,n+l

i,n

5 i 5 R-l

- - +
-(MYDS)Si-l,n+l + (MYDS + PYDs)31,n+1 ‘

+ _.
(PYDs)Si+1,n+1 5 (MYDD)Di-l,n+l

_ + + -
MYDD + PYDD)D1,n+1 ' (PYDD)D1+1,n+1 + (AYDE

+ (l/At + (C-3)

)Ei-l,n+l '

+
’(AYDE)E1+1,n+1 = (l/At)Di,n

E-Equation

 

2

R-l

IA
p.
IA

(AYDS/€)Si-l,n+l ' (AYDS/€)Si+l,n+l + (AYDD/€)Di-l,N+l

'(AYDD/€)Di+l,n+l + (l/At ‘ YDE/€)Ei,n+l = (C'”)

(l/At)Ei,r1 + (In+l/e)

ll
H

S-Equation i

 

(l/At + MAkSO +

(uAkDO + PYSD

(”AYSE)E1,n+l =

II
H

D-Equation i

(”AkDo + PYDS Sl,n+l '

PYss)31,n+1 ‘

YSE/€)Dl,n+l '

(l/At)51,n + AA(kSOSO + k

165

(PYss)32,n+1

+

(PYSD)D2,n+1 ‘

DODO)

(P DS)SZ,n+1 + (l/At + MAR

(C-S)

SO

+ PYDD ‘ YDE/€)D1,n+1 ' (PYDD)D2,n+1 ‘ (“AYDE)E1,n+l

= (l/At)Dl,n + AA(k

II
P

E-Equation i

 

—(kDO/€)Sl,n+l

= (l/At)El,n + (In+1

S-Equation 1

II
III

 

‘(PYSS)SR-1,n+1

”(PYSD)DR-1,n+1

+ (”AYSE)ER,n+1

+

+

DO 8o +

+ (kSO/€)D

/e) +

ksoD o)
1,n+l + (l/At)E
(kSDDO/e) - (k

(l/At + AAkSa + PYSS)S

(UAkDa

(1/At)sR

l,n+l

DO SO/e)

R,n+l

+ PYSD ' YSE/€)DR,U+1

+AA(k

88.83.

+ kDaDa)

(C-6)

(C-7)

(c—8)

166

D-Equation i = R

‘(PYDS)SR-1,n+1 + (”Akpa + PYDS)SR,n+l ‘ (PYDD)DR-1,n+1

+ (l/At + “Aksa + PYDD - YDE/€)DR,n+l + (”AYDE)ER,n+l

 

= (l/At)DR,n + “A(kDaSa + kSaDa) (C-9)
E-Equation i = R
(kSa/€)SR,n+l ‘ (kSa/€)DR,n+l + (l/At)ER,n+l ( )
C-lO

= (l/At)ER,n + (In+l/e) + (kDaSa/e) - (kSaDa/e)

APPENDIX D

NEWTON'S METHOD AS USED BY PROGRAM "IONFLO"

The matrix equation,

e<i> X = a <D-1>
can be represented by the equations
LJ(X) = BJ j = 1,2,...RR , (D-2)

where each value of J represents a row of the matrix.

Expansion of L (2) in a Taylor's series about the RR

J
variables yields,

1 RR RR 32L
- ( ) A AY AY + (D-3)
2
RR 3 L
1 J 2
- ————) (AY ) +
2 i=1 3Y12 1 ’

167

168

where

L = LJ(?O) and AYi = Y O (D-A)

o
J 1'Y1
Truncation of the Taylor's series after the AY1 term and

substitution of Equation (D-3) into Equation (D-2) produces

RR 31;
1-1 (§§i)§_§OAYi = R3 - L0 J = 1,2,...RR . (9—5)

This can be written as the matrix equation,

£0 , (D-6)

IN)
I

QLY=

where the matrix elements of g are defined by,

L}?

) . (D-7)

G = ( A A
13 1 X‘XO

m
K:

This suggests the following iteration procedure:

1. Choose 20.

2. Calculate g and 20 using Equations (D-4) and (D-7).
3. Solve Equation (D-6) for 4:,

A. Use the last of Equations (D-A) to calculate 2.

5. Compare £2 and 2.

6. If IAYi/YII i = 1,2,...RR is less than some specified

tolerance, the iteration is terminated and the

program proceeds to the next time row. If not

169

+ 20 and the procedure is repeated starting with Step

Ir<>

Of course, there is no guarantee that this procedure
will lead to convergence for all initial choices of 2°.
Program IONFLO has two options for choosing go. Option A
chooses 2° to be the solution to the previous time row.
This is simple to use and the time step can always be
decreased until the procedure converges. Option B projects
the solution at time row n to the solution at time row n+1
using space derivatives and transport coefficients evaluated
completely at time row n. This projected solution is then
used as go. These values of 20 are almost always more ac-
curate than the ones from option A, so larger time steps
may be used. Option B is usually more efficient for rapid

approach to the steady-state of the system.

APPENDIX E

LISTING OF PROGRAM IONFLO

The following input deck executes program IONFLO:

ATTACH,MAIN,IONFLOBINARY.
ATTACH,X,JHLLINPACKBINARY.
ATTACH,Y,JHLBLASBINARY.
ATTACH,T,AZCOlBINARY.
ATTACH,U,INTEGRATIONS.
ATTACH,W,EQNlBINARY.
ATTACH,V,GRIDPRBINARY.
ATTACH,R,JHLPLOTBINARY.
LOAD(MAIN)

LOAD(X)

LOAD(Y)

LOAD(T)

LOAD(U)

LOAD(W)

LOAD(V)

LOAD(R)

EXECUTE.

The following is a listing of program IONFLO:

170

OOCOOOOOOOOCOOOOOOOO 000000000000OOOOOOOOOOOOOOOOOOOOOOOO000000000

 

171

PRCXSRAM IONFID (INPUT,GJT‘PUT,TAPE GIIOIJ'I'PUT‘HIAPE 50=INPUT,TAPE 7)

PRERAM IONFLO SOLVES THE SYSTEM OF EQJATIOI‘B GENERALIZED NERNST-
PLAMIK, POISSCN DISPLACDIENT‘ CURRENT FOR THE TRANSPORT OF A
BINARY EIECTROLTT'E IN A FINITE SPACE IXNAIN. DERIVATIVES ARE
APPROXIMATED USING FINITE DIFFEREMZE.

2T'YPESOFBOUNDARYVALUE PRCBLEMSAREPESIBLE. KEY=1 ISE‘CRT'HE
DIFFUSION POTENTIAL KEY=2 IS FOR KINETIC ELECTRODE MODELIM;

BOI‘H EMPLOY A NON-UNIFCIM SPACE MEH WHICH DEPENE ON THE

MEAN CQCENIRAT‘ION OF THE EIECT'ROLYTE AND THE INPUT PARAMETER

INPUT

ALL INPUT IS READ THRGBH SUBROUT'INE AIO. THIS SIBROUT'INE SHOULD
3E INSPECTEDT‘ODET‘ERMINE PROPERGRDERAI‘DFGRMATOF INPUT‘

THEFOLLONING ISALISTOF meVARIABLESAmmEIRMEANm.
TITLE DESCRIPTION OF SIMULATION

KEY 1 FOR DIFFIBIQI POTENTIAL SIMULATION
2 5m KINETIC ELECTRODE MODELING

R NLNBEROFSPACE POINT‘SUPTOIEB
Q NLHBEOFTIMERGIS

ISEI‘ NIMBER 05‘ W SPACE POINTS AT EACH ELECTRODE
on SIDE OF ACE

IJK FRACT‘ICN OF TIME ROVS TO BE WRITTEN

IST'OP NINBER OF TIME ROVS T‘IME STEP SIZE ARE‘CCIBT'ANI‘

M14 SPECIFIES A TIME ROI HR WHICH MANY INTERNAL PARAMETERS
BE WRITTEN IN UNFORMAT'I‘ED OUT‘HJT. THIS GENERATES A EDT
WPPUT‘ AND SHOULD ONLY BE USED FOR DELGGM. WHEN NOT
USE SEI' W GREATER THAN Q.

NPLT‘ FRACTIQIOFTIMERONSTTRT‘AREWRIT'IENT'HATARET‘OBE
PLOT'I'ED

WILL
OF
IN

IOFF 1 Fa? C(MPLETE ITERAT‘ION mm
a PG? no ITERAT‘ICN

IW GWRITESOUPPUTCNLYFCRIASTITERATIWOFEPCHTIMEROI

SPECIFIED BY IJK
lWRITESCXJ'I‘PUT‘FCREVERYITERATIwOFEACHT‘IMERm
SPECIFIED BY IJK

ISTUPIAST‘TIMEROIITT‘RT'CURRENT‘ISTURNEDINFE

RAEE DETERTINES RATE AT WHICH GRID SPACING EXPANIB FROd SMALL
STEP SIZE TO LAEE (NE.

171‘ INITIAL TIME STEP SIZE IN SECQ‘JE
F DETERMINE RATE OF TIME STEP SIZE INCREASE.

TOL MAXIMIMAILOIEBLE FRACT‘IQIAL CHAEEFROJICNEITERAT‘IW T‘O
T'HENEDCI‘FORPRERAMTOPREEDETOT‘HEWTDIERGV

ZP CHAEE NIHBER Em POSITIVE ION
ZM CHAEE NIMBER Fm NEGATIVE ION

0 0000 0000 0000 0000000000000000 00000 000 0000 00000 000 0000 0 00000000 000 000

172

20 FINAL m (N ELECTRCIE IN C/OI"*2

ALF TIME Cm Fm CHAEE INJECTIW

A CHARACTERISTIC LEMSTH IN (N

GDE MEAN 2m VALUE IN MOL*S*C/G*CM'*3

SIZE LEICT'H w SYSTEM IN 04

EPM MEAN DIELECTRIC CIEFFICIENT IN C*C*S*S/G*QT**3
84 MEAN ELECTROLYTE CGCENT‘RATION IN MOLE/04”“

[030,169 KSLNDA,KSAI<DRKSR RATECQTSTANTSFORTRANSPORTACROSS
ELECTRmE-E LECTROLTT‘E'DTPERFACES.

SR,SL,m,DL, BATHIM; CQCEM‘RATIODS.
C(1) THROIEH CIEISI'g)AND EXPANSION COEFFICIENTS Fm ONSASER
COEFFIC ANDACTIVITY COEFFICIENT DERIVATIVES.
C(1)LAND c 2) FOR L-H-, C(3) AND C(4) FOR L—, C(Sg AND C(6)
AREF7%*AND {-19 SAMOL‘vMof/Cégw CIGRfiGER
MANSIONSARE
L-H- : C1*CP + C2*CP**1. 5
L— = c3*c4+ C4 «14*an
L-I: . cs*s**1. 5 + C6*S 2
E=w+c
NAME(1,2,3) TITLES FOR PLOTS
JNAME(1,2,3) TITLES Fm PLOTS
PomT(1.2.3) SYMBOLS Fm PLOTS
OUTPUT
ALLOFTREINPUT PARAMETERS THEDEBYELWI‘HOFTHESYSTEN
AND THE TRANSPORT PCOEFFICIEN'I'S IN DIMENSIONLESS FORM ARE WRITTEN
BEFORE THE FIRST TIME ROI IS EXECUTED
S SLMCONCENTRATION. WRITmNMEVERYTIME RONAS INDICATED
BY IJK AND EVERY SPACE POINT AS INDICATED BY 1w. LNITS
ARE MOLES/owa
SM IS SUETRACTED FRO! 5 FOR THE ourou'r
D CRAROEDENSITV IABELEDRROONourRJT. WRITTENKREVERY
TIMEROAIAS INDICATEDBY IJKANDEVERYSPACE POINTAS
INDICATED BY Iw. LNITS ARE C/CM**3
SAME comm AS FCR D. LNITS ARE VOLTS/OT

PSI SAMECQMEN'TASFCRD. LNITSAREVOLTSRELATIVETOTHE
POTENTIAL AT X39.

DIMENSIWIES POSITIVE ION CQCENTRATION
DIMEIBIONLESS NEGATIVE ION COICENTRATION
SPACE POSITICN IN INIT‘S OF A

TIME IN [HITS OF SECNE

CELL POTENTIAL DIFFEREMZE IN VOLTS

I131

8*“98

000 0000000 0000 0000000

000 000 000 000 000

173

21 CURRENT mm IN AMPS/04“?
CLE TOTAL CW8 IN SYSTEM IN COULGVIBS

RCCND SEE LINPACK DWLMENI‘ATION
IF SUBROUTINE SGBFA IS USED BY IONFID RCGID IS MEANINGLE‘SS

INTEGERgR

REALKD RbKSR
DIMENSION 7:59:32? Prsi'i'émmeé , E30 6;
DIMENSION x 39m ,DELCS BO) ),XER 300 ,YY(3DO)
DIMENSIONAB
DIMENS mfimpggflm (1O£,NAMEém)
D,KSO,KD ,KSL, ,KSA,KDR,KSR,SR,SL,m,DL
CO‘MON /ELEC/
gfl$§§WEC2 zRQDTFTOLRANGE ISET
I I I
CCMVIQJ /UNITS/ SM, GDE ,Em,A,2P,2M,CLN
CCMMCN / BOwa/ NPLT, IJK
COMMON /PLTPTS/POINT(1O)
CMMON /CHRG
SUBROUTINE AIo READS THE INPUT DATAAND ANDCONVERTS VARIABLES To
DIMENSIONLESS FORM. MOST OF THE INPUT DATA IS ALSO WRITTEN
CALL AIO (IOFF,IS'IOP,SIZE, DEBL,NP!‘, NAME,JNAME,MM)
SUBROUTINE GRIDxx, WHERE xx Is PR FOR mam BLOCKING ELECTRODE,
RR FOR Two BLOCKING FORL'IgID JUNCTION
S%MTIONS{<EYGENERATE 3 THE NONfiNIFORM SP PACE H ACCOTSDUIJLSTGM TO R,

GA AND DEL ISTHE mm AND XERH
SPACE PCBITIONS OF THE GRID PTS.

CALL GRIDPR (DEL,XER,SIZE,DEBL)
Slim ORIGIN FCR KEY=1 PR‘BLEM
KEYg. 1.38. 2) GO TO 13
C(NTI)=XEfi(I)-g. .5
CQITINUE
RR=3*R
TIMff-OJ
WW1.0
ZI = DJ
M =
KOUNT IS THE NIMBER OF ITERATIOIS PERFORMED AT EACH TIME Rm.
KOW
SUBROUTINE IC ASSEMBLES THE INITIAL CCNDITIOIS FCR S,D, AND E.
CALL IC (KEY,Y,R)
OUTPUTI‘IM; THE INITIAL CWDITIONS
CALL IANW (Y,TD4E,KOUNT,KQD,XER,NPI‘,MME,JNAME,M,KFLAG)
SAVIMSFASFFFGRUSEAPI'ERCHAMIM; ISCMPLETE.
FF = F
BEXSINNIM; OF TIME LOOP

DO
17MB ALF*2&EXP(—ALF*TmE)

0000 00000 000 000

...a
co

000000E00 00

19

000000 000000 0000

N
E

000000

17”

IN GENERAL TIME(N,F,UT)8DT*((1+F) **N-1)/I"

TIME-TIME-i-UT

IF Ream LESS 'H'IAN LEE-25 THEN MATRIX WILL BE NEAR SMUIAR

IF (m .LE. 1.03-25) Mo“
SUBROUTINEREBVOCOTHHESTHEPCRTIWOFTHERPB-WGENERATED
BYfliESOLUTIClBFROflTEEPREVIOlBTIMEROV. RHSVOISCAILEDONLY
GEE FCR EACH'TIME RON. THERHS-VEC‘TCR ISRETLRNEDASB.

CALL REBVO (RAEL UT ,Y,B,M,Mfl)

THEWMPORTICNOFBISSAVEDASBBKJRIBEWTTHEACH
SUBSEQJEN'TITERA

ID 18 IBlfRR

5%331

SUBROUTINE C(IMAT ASSEMBLES THE CWFFICIENT MATRIX, ABD, IN BAND
STORAGE ch

CALI. COEMAT (R,WL,UT,Y,ABD,M,W)
SIBROUTINERTBVNCCNPUTESTHEPORTIWOFTHERWTORUJETO

THE SOLUTION! TO THE PREVIOUS ITERATICN. THIS m IS RETURNED AS
B. 'IHEHW B+BB IS THECGTPLETE RIB-VECTOR PCB ANY GIVEN

CALL. QRHSVNRARpDELmTflpB)
B(J)1=B(J)J-BB(J)

TEEN

IDAISTHELEADDBDIMDBIONOFABD,ML IS'H'IENLMBEROFDIAmNAIS
ABOVE‘IHEMAIIN,MLIS'IHENUIBER EROF DIAGNAISBWTHEMAIN

LEA-16
MIPS
MUIS

SSBCOFWETHEMA'IRIX ABD. MWISRETURNEDASl/COND.

° MERE C(ND IS THE C(NDITION OF THE MATRIX.

IF mam IS DDT m, SLBRWTINE SCBFA MAY BE [BED TO DECREASE
EXECUTICN TIME

Juana

CALL SGBFA (ABD,LDA,RR,ML,m,IPvr,INFO)

SGBSL SOLVES THE MATRIX EQUATION ABD*Y=B. WHERE ABD Is TEE BAND

MATRIX FACTCRIZED BY SGBCO B IS FORMED FRm RRSVO Am RHSVN,

ANDYISTTTESOLUTIONVECTLSRRETURNEDASB. TREENIERING
VECTOR, B, IS DESTROYED.

SQSL ABD, LONELMLJIJ, IPVT,B,JOB)

KFIAG IS
NEXT TIME RCN. OTHERWISE KFLAGBI AND AMI-ER ITERATION IS
PERFORMED

0000 0

0000

N
\l

000 0

000 00

AH

12

0 0000 0000 00 0000 0

I7
16
10

175

CALL (NPR (Y,YY,TOL,R,KFLAG)

THE NEXT CARD REDlfll'B EXECUTION TIME FOR S‘IALL PERTURBATIOIB
BY ELIMINATIM; THE ITERATICN MEDJRE

KFLPG=IOFF*KFIAG

KOLNI‘BKOUNT-i-l

IFKOlNTREPCHESllITERATIWISPRCBABLYDIVEFCINGSOPRCBRAM
ISTEFMINAT‘ED.

IF noémerR .1530. 11) GO TO 16
Y(J)=YY(J$RR

cmrmus

YERR—Zg =- 0.0

IFMISNOTE T‘OlORAMULTIPIEOFIJK,TT-IECXJTPUT
ATTHATTIME WIILNOTBEWRIT'I‘EN

L-mommx)
IFEL. ..Né 9 ..AND M.NE.1)GOT014
IF IW..EQG)GO'I'011

mummcummsmcmmemmsmmsysm nil-2mm
ASAFLNCTIONOF EVERYIJKmeERmAmpLo'Is
s,D,Am3ASAnJN6rION or x.

CALL IANDVI (Y, TIME,KOLNT,KOID,XER,NPT,MME,JNAME,M,KFLAG)
C(NI‘INUE

. GOTO
IF L.NE. EQ .AND.M M.NE. 1) GO TO 12

IF IW .NE. g0 GOTO 12
IANDN (Y, GO,KOINI',FCCND,XER,NPT,NAME,JNAME,M,KFLPG)
CQJ‘TINUE

UPISINCRENENTEDATEACHTIMERCWACCORDIMSTOTREVALUEG‘F
THAT IS SPECIFIED IN THE INPUT DATA

KEEPS TIME STEP SIZES CW WHILE ELECTRON
IS BEIK; CHAMED

IF (M .LE. ISTOP) F I 0.0

AFTERELKTRQJECHAmE ISCOdPLETEDTIMESTEPSIZESARE
GRAHJALLY INCREASES IF F IS SET POSITIVE

IF (M .GT. ISTOP) F = ABS(FF)

IF F IS SET BETWEENI AND -1 TIME STEP SIZES WILL BE
DECREASEDB YABS FF) ’Fm TIME R00 ISTOP+1 AND THEN
GRADUALLY EDTHEREAFTER AS US

IFEM. .ISTOP.AND.FF.LT.O.D DT=ABS *DT
11258 F+1?8)*DI‘ ) (FF)
C(NTINUE

STOP

CCNI‘INUEl
”Rm (611%
WT (1X, 2 HSOLUTION VECTCR IS DIVEFCII‘B)

000000

Ian- HHHO—‘H

N NNNN N
‘0 QOUl NON” bWNI—‘Q

WU
NH

33

176

SUBROUTINE AIO (IOFF, ISTW,SIZE,EBL,WI‘,NAME,JNAME,M)

SUBROUTINE AIO R“ IN ALL THE DATA NEEDED BY IONFID. ALL
INPUTDATAISAISOWRIT'I‘ENBYAIO. ALLINPUT'FGIMATIS EITHER
I5 (R E10.3 EXCEPT NAME AND JNAMEWHICH AREAS AND TITLEl

AND TITLEZ WHICH ARE A10.

DIMENSIONkS NP'I' II) .JNAME(1I) .NAME(1I)
INTEGERR
%wm/%Q3;ik?bL,kL,m,mflm,m,52,311,111 DL

comm /ELEC/
comm /GREEK/EPS 335 mags?

CWQI/UN EPMA,ZPZM cm
comm /p11%é6% Dr' E,TDL.RAmE.ISET
cmMCN /PL
IéNPLT,i
comm /CRRG/ ALE
EI .37 T2 LE1 TITLE2

READ
READ 50,10 KEY,R, ,ISET, IJK ,RIPLT,IOFF,IW,W,ISTOP,RAM3E,DP,F
READ 50,13 (EéEm'zoéiALFSIEPS'A'mL

 

 

READ 5I.13

READ SG,1§1§DWKS0'ELSKS,LHEAIGA,MIGR
READ 5I,1RaéL.DDE

READ Sl,13 SCiI).I--1.1Ii:::

READ SI.11 pomigmgl

READ 50,12

READ SI 12 1 4.181

WRITE (61,3 ) TITLE ITLE

WRITE 61,28 QKEY

WRITE 61.35 DT..ETDLRANGE,ISET
WRITE 61,36 SR 'sL.DR,DL

WRWRITEITE 61 'EZI JKTISE IDEE. ISTOP, NPLT
SC=SE§(EPEEPI+ JKEE 13.31312?”
DEBL=1.631*58RT(EPM/

WRITE 61.4I DE EPM .SM.zo,SIZE
WRITE 61,22

WRITE 61.27 C1.C

WRITE 61,26 c3 .C4

WRITE 61,25 c5 .C6

WRITE 61.14

WRITE 61.15 C7 ,C 8%,

WRITE 61.16 C9 .C1)

WRITE 61.34

WRITE 61,38 KDI.KSI.KDL.KSL
WRITE 61.39 KDA .KDR.

WRITE 61 29 DEB

FORMAT III5.3E1I.3)

EGRMAT 1E

EGRMAT 8 )

A1
8E}0. .3)
FORMAT // ,2X,198MANDBED(PAEBIOIB%’Q
31mg =.E1I.3 .SM/S 131.810.3118 s**1. 5./)

ROME fizE“? 31E EQSAEERH +0EFFI CIENT EEEANSIDNSJ)

 

 

 

 

FORMAT 3X,SHLPM =-;.E1I.3. I * S)**l. 5 + III. 33 73* S)** 2 /;
mflmfl ’“sammp: 01.15133 3 33* 2 +81. 3‘ :E18°3 .3 1331 E 3’ 81131? 73
FORMAT /x1x.14 spAcE1ns..4x.1410H'mEEavS. °

 

12x 63 .12 .//
11111111; é/ésixopfinm DEBYE LENGTR=1I.3.//
ECREAT' 2X4 E13

3;,31/ .ES

EGMAT 1H1. 6Ix .21111I./Emmi)?!”NFL '//)

36

37
38

39
40

0000

S0

000

000

00000

177

T 3x 1411RATE GGNS'IRNTS
T géé/Igufirggyh 3. '3335'33 FIIIQLG. 3. 311.511 T0L-.E1I.3. 311.711 RANG

I

1E= E1I.3x

mT £11,111 SR=.EII. 3. 3x.411 SLa. EM. 3. 2x. 411 1383.810. 3. 4x. 411 DL-.
0)

35:51! 0103308. E10.3, 2X 51-] KSO-.Ell.3 1X SH KDLB,El0.3,
EBEDAF, E103: 21,51! KSAI, E10. 3 131.511 KDR- ,E10.3.

13x 511 R-. EII. 63$} )3
WT 1I.3. 211.511 EpMn. EM. 3. 111.511 314- -.E1I.3.
13x. 511 202x -IE1I.3.5K 5x.711 SIzE -.E1I.3I///)

C(NVERSIQJ OFUI‘ 2AM MAN!) THE C-EXPAINBION CCEFFICIENI‘S TO
DIMENSIONLESSQIANI' IESECR (BEINTHEMAINPRQERAM

cm 8 Ema.“ “0*96485 0*ZP*ZM)
1271' I GDE*
=- QCWIZ%A*91*EPM)

58 CW

CSI) I)=1. I23§E-09*C(I)

13x 511T
Edw.

 

C1 I5C1 “(3&4

C 2 8 C 2 0.513“ /(GDE*A*A)
C3 31C3

C4 8C4 *0.5* GDE*A*A
C2 IICS 05 *GDE*A*A
C 8C6

C7 8C7 SM

C8 =C8*SQRT(91)

C9)='C9

CONVERSION OF BATHIM; CQCEN'IRATIOBS 'I‘O DIMENSIONLESS FORM
SR: SR- 91/

SL: S

[1513919d

DL- aIDL/91

DIMENSIONLESSIZATIDN 0E RATE CCNSTAN'I‘S
DIM . $0W6485.I*A*ZP*W)
... D

SUBROUTINE IANDN (Y, TIME ..KOWP,M(ND,XER,NPT WE,JNAME,M,KFLAG)

WCAICUIA'IES'I‘HEPOTENPIALASAFIMIONOF TIIETGTAL
CHAISEANDMASSDITTIESYSTEI.ANDWRITESTTIE THEOUTEUIEGI
TIME RCW SPECIFIED BY IJK.

31118111111811.1111 “31%81’133111133’

BMWIE 1.1111115 33321111113311“)

0000 000 0000

000

0000

N
U1

0000000

178

DIMENSION PD 150 ,PUI‘(15@)

COMMON #3ch 21

com ONP/ KEQéD R ,,Q,DI‘ F, 'IOL, RAMEJSET
comm /UNITS/SM EPM, A' ,zp,5M, cm

com 1 /Box/1w 1w,NPLT EIJK

ccmou R/PLTP'I'S/ ’Pom'fua)

INTEGER

DATA IIO/61RRMgX/75/,ISCALE/l/,NF/3/
3335*}: :XARA /225*0. a/

STRIPPIM; OUT THE MASS CHARGE, AND ELECTRIC FIELD VENRS FCR
INTEGRATION AND FORMING THE + AND- ION COCEN'I‘RATIODS Fm PLOTS.

Q

 

D) 16 K=1,R
ZARAY K,1 = 814* YE3*K-2+ - Y1 3*K-1 ;*1.0E+fl6

ZARAY SM* *K-2 *1 . 0E+¢6
ZARAY — Y(3 )g 1-.$E 11*Aj3945

‘11 K; = Y 3*K-1;

YJ K 3 ¥ gglf-Z

E53565

WERANDUPPERBCXJNIEFORAVINT

XLO=XER1)
XUP=XERR)

CALCULATDK; THE TOI'AL CHAME AND EXCESS MOLES 0F IODS IN THE
SYSTEM

'ICS= AVINTXER,YI, RWXLOXUPIND

'IMS= ==,AVIN'I‘XERYJ,R::,Xl’..OXUl=’ IND

CAI£ULATIM3 THE POTENTIAL AS A FUNCTION OF X

m 26 J=(1‘R

PSI (J) = -AVINT(XER,E, R,XLO, XUP, IND) *l flE-G7*A*S'4/CDN

gDEWLS'K; = PSI (R) + 1.0E—07*A*SM*XER(R)*E(R—l)/CON

couvsa'rmcm'rovoms TIMETDSDcoNDs, ms'rochoMBs, ms'ro
MOLES, AND 21 To AMPS PER (MHz
RIME = TWCON/GDE
mm + M/IJK; =
RI = A*SM*GDE EPM*ZI/CQJ**2
ms = TCS*2.I*ZP*ZM*96485.0*SM
= ms*sm*2. a
WRITE W31 ,24 ngéxowr, .1115
1525/ 1m<'1ox,3asu~1,9x 31mm 9x, 2112+ lax, 211c-, ll-IE
1,612.57x w132x7 K6Nr=157x,6 =,Em 013% '
(86X,4H EJ=,E10.3,1X,5H ms=, E19. 5, 1x, 41;iMS 21=,Em. 5)
11133114 comm-non LESS THE BULK CMEN'IRATION
DENSITY IN couwMBs cc

VOL'I‘S
IS THE DIMENSIONLESS WITIVE ION CQCENTRATION
IS THE DIMENSIONLESS NEGATIVE ION CQCEN'I'RATION

1 *SM*ZP*ZM*2. 0*96485. a
+2) *1 . IE-G7*A*Sfi/CON

000000

00 0000000 0000 0000 000000000 000000

17
15

27
28

179

at; : $5:- 1 - 2214-1;

WRITE 23) SXER((I+2)/3) ,AY,BY,CP, 01,CY,PSI((I+2)/3)
MT 7152

CONTINUE1

PLOT-rm; ELECTRIC FIELD POS. ION can. AND NEG. ION caVC. vs.

LE ELECTRODE. SCALE IS CORRECT ONLY FOR
THEFIRST ISTOPNIMBEROFP'IS. CPANDCMARE INM ILLEMOLES ITER
RELATIVE TO THE BULK. ELECTRIC FIELD Is IN VOLTS*1E**(4)/

.NE. 0 GOTOl7
IPL’A‘SI JKT*NPL )

JPLT 8 MG) M, IPLTAND

IF JPLT. M .NE. 1) GO TO 17
NPT l =MAX

NPT 2 =MAX

NPT 3 =MAX
XUP = P.74 LOTW‘XEREZ 5
Y,ZXARAY, ISCALE,JNAME,MME,XLO,XUP,NPT, NF,MAX)
IF‘MP .EQ. (0) GO TO 15

31+

WRITETAGl 22:11:5ng“ ,PD(I),I=1,M)
FORMAT //,1X, 46H) POTENTIAL DIFFEREICE CORRECTED Fm OHIO WP,”
END

SLBROJTINE COEMAT (R,ML,UI',Y,ABD,M,MN)

SUEROUTINE COEMAT ASSFMELES THE CORE MATRIX FOR USE m THE

NFm'ON—RAPHsm ITERATION PROCEDURE. THE MATRIX IS ASSFMELED

DIRECTLY INTO BAND STORAGE FORM FOR USE IN THE LINPAC SUEROUTINE

SGECO. R IS THE NLMEER OF SPACE GRID POINTS. RR IS THE TOTAL

NLMBER OF DEPENDENT VARIABLES. DEL IS THE GRID SPACING VECTOR.

DT IS THE TIME STEP. ABD(LDA RR) IS THE MATRIX 1N BAND STORAGE
ML—NLMEER OF DIAGONAIS

m. LOAF-=2*ML+MU+1 ABOVE THE MAIN.
MJ—NIMBER OF IAGQIAISDBELW'IHEMADI. YIS'IHESOLU'I‘IONW
FRO! THE PREVIOUS TIME R01 (R PREVIOUS ITERATION. CEMAT CALLS

'H'IEBADDA'MAGE'MATRAX' IAPASSEDTOTHEMAINPROGRAMW
THE comm BLOCK COE. THE RATE CQISTANI‘S ARE PASSED To C(EMAT
momHTHEccmmBmCK

RCNS'IS.
IFTHEREWERE WERElSSPACEGRIDPTS. 'fliEAANDABDMA‘IRICESWQJLDLCDK
LIKE'IHEFOLWDB B.LCX3KS

E

QQEQQ m>~ooo Q05
“8 “SE :9me 0301
EQEQQ‘OQQO‘U‘ID
aawromméw
SQ ’KDCDQO‘U‘thN
EGOQQO‘U‘MNH
"DMO‘UIIBW SQ
\OQQGUDUNQQS
mummwat-‘aas
\lmthUQQGQQG
O‘UIMNQQGGQG
thwNO-‘GEGGQQ
waQQQQQGQQ
UNGQQQQQGGS
NHSGQGGQEQG
ass SQQBSGGQG
OMQSQQQGGGQ
mm bGGQQGQGQ
awkwaaaasaa
mthNEQSQQS
GthUNl-‘QQQSG
O‘U‘bNQQQQSQG
mmbwmaaaaaa
OSU'IDWNI-‘QGGQQ
O‘U‘lwaQQQGGQ
GUIWNQQQGGG
O’SUlchUNI-‘QSQQG
GMMQQQSQQQ
GMhUNQQQQQG
O‘thWNt-‘Q&QQQ

000000

000

0000

MM

180

secs
aaaa

as
aaaa
scam
see»
$330
013%
omqm
mqmm

5
9

GUI-DU

DION

0 90

WDKDM

REAL KDO,KSO,KDL,KSL,KDR,KSR

Q>DW

seem

WFWQQ

“>0“

QQ‘OM

wymmq
a>mmq

EEOC)“

mywcoq

$3,me

QQ‘OGQ
QEQWQ

QEEQQ

INTEGER R RR RSTOP
80M wl%%}a%éy BSGAEIlaRSSGRSAgASA KDR KSR SR SL m DL
C /ELEC/ EPS'Z£ I I I I I I I I I I

$333“!
IBflP=IMI-5

DENIUZDBTHEMmflmmHXTOZfim

DO 130 91:1m .16

A3%I?J 865.5

MEDBUNSflfllflflflXEUflflflSAflEDflImfllfimimmmMN

ammmmusxrww

1352'=D1§S£11{*DEL(1)/2. 0

2%}

F1 = 1.91/01
=KSH/DzD

F3 = zSS(s D)/D22

T1 = F1+ 5‘2 + F3

T2 = zSSP S,D)*Y(1)/D22

%; 3

T2. /2 2mm)
T5 = i/Aégg °

T6=

ABD(11 1)=

F1 : 2DS(S 2D)/D22
F1+ 52

zmpés, Wig/022

-ZDEPPSé:DLY

-HEP *D
-5354; :g; M” 2.
-2DDPS
5“€

KISS/D
SDQSD

Z
22

-233; S 2,} 5:3}

655352333333
I

MIIINIIIIIIII

HH

D

afidfifig
S§:$

  

.‘DZ
10,2/D g/(z 0 )

4ZDDS:I%¢R2

IIIIIIIAIIIIIII'IIIIIIIIIII

:SSEEHESR

2+T3+T4+T5+T6

2+T3+T4+T5+T6

T3+T4+T5+T6+T7+T8+T9

00000

181

T =ZDSDS,ZD*Yl/D2
T = musi' 133* vf23/D2 22
T7 = -2DSD S,D D*Y 4 /D22
T8 = -ZDDD S,D *y 5 /D22
T9=-ZDEDSD*Y2/EPS
Tm -ZDE (S, )* ( )/(2.a* D2)
ABDéll, =T1+T2+ T3+T4+T5+T6+T7+T8+T9+T1z
2 2 = KsyEPS
1 = -is S D?) (2.9MD2?D
T2 = -ZSE SP P)/§2. 2* 2)
,3?) = Ti + T
T1 = —z E S D))/(2.D*D2)J
T2 = -2DE SP P)/+2. 3* 2)
D 91,3) = T1 + 2
ABD 11,3) =
ABD8, 4 = —zSSS,D/D22- a. 5*ZSEPSP,DP/D2 *Y3
$8 9644):” -2DS S, D /D22 - a.5*2DEP SP, DP /D2 *2 3
ABD 7,5) = -zSD S ,D)/D22 —. 0.5*ZSED SP,DP /D2 *2 3)
ABD 8,5 = -ZDD S ,D)/D22 — 0.5*ZDED SP,DP /D2 *Y 3
ABD 9,5 =
ABD6 6 = 0.9)
ABD (3,6) = 3.9
ABD(8,6) = 2.0
LING THE INTER

ASSEMB IOR MATRIX Em
IT IS CQNENIENT T0 LOOP IN MULTIPLES 0F 3 BEAUSE THERE ARE 3
EQJATIOI‘S AT EVERY GRID POINT

m 1Q1K=4 ,RSTOP,3
1013“; 1)/3

fK-33
K-3
K+3

3:13: 2&3 +1{)/2.a

K
DP= = WEY yKfi + Y(K+l))/2.l
=DEL KM + DELéKK)

+ Y(K))/2.a

 

 

SPP:Y

T2 = -ZSSSM*:1§K-3/DEM
T3 = *YK-Z
T4 = ZSSP ”A“? éfim
a ZSEPSMQ RA {Dézmm/AD
T4 + T5 + T6
:3" {53“‘3mm
ZDDP 'm *YiKMD
2138 9451.113“) RAMK+
81313913 =I=T1+ +334. +)44+T5+T6

ABD 15,K-

g? 16,;333 = zm(S, TD)*RAT/(EPS* AD)
T2 = -ZSS 134W K-3 /DEN\1

T3 = -2

T4 = ZSSD

T5 = 25145ng K+l)/DEM

 

 

182

a 2S? Sm mo*RAT*Y K+2
T1D(13£11§D -§4='T1/&§42 4344 + T5 4- T6
T2=- ZDSDSM,DM*YK-3/Dem
>0 9‘6?“ ”Maw
DDD 94:34 *Y K+1)/DEm

a T1 + T 4 + T5 + T6
51135, D) *RAT/(EPS AD)

919

0.0

 

 

 

T2 = zss SP,DP /DENK

T3 = 255 sum *le

T4 = -zs S,D nap“

T5 =( -zsap D*S K4m +2)/mm

ABD11K)=- ST1+T2+T3+T4+Ts
zfiSSM m/DENK

 

*YK/+l /1-:ps
T4 = -ZDEP=§ §*Y K+2 *SUB/MUL
ABD(12 K T4
T1 a 266 ES,D 1*Y K-Bi *RAT EPS*AD
T2 = zDDpi S'D *Y K-2i BPS*AD
T3 = -zDSpé /PS*AM
T4 = -2DDPS*YK+1S*S/(RPS* )
T5 = -ZDEP D*Y K+2 /EPS
T6 = -ZDSP S' D *Y K+3 / RA'IV‘EPS*AD
T7 = -ZDDP S' D *Y K+4 / RA'I'*EPS*AD
T8 = -ZIB(= .151 SUB/(E AMUL
ABD(13K)= +T2+T3+T4+TS+T6+T7+T8
T1 = zSD SP, E)/D15:KK
T2 = 281835;?
T5 = l-ZSED?S?DE*SLSB*Y(K+25/1§MUL
M(g+ ='=1+T2+T +T4+T5
T2 = ZDD SP, DP m
T ..T T $6, TT)
T4 = -z Sé 'AD 395m:
T5 = -
T6 = -2DED{S' D2.*S+KT *Y(K+ 2W
A13D(1%ls‘5 +T +T4+T5+T6
T1 2 SM K-3 i*RAT/ EPS*AD
T2 = D*Y K-Zi *RAT EPS*AD
T2 = -ZDSS 1“P1101803 PS*AM m)“.
T = (EPS*AM )
T5 = -ZDED 5' D3 2/EPS
T6 = -2DSD S, D *Y K+3 / RAT*EPS*AD
T7 = -2DDDS D*Y 54(4/ /RAT*EPS*AD
T8 = -zw( ,b) AMUL)
ABD(12K+1'=T1+T2+T3+T4+T5+T6+T7+T8
$2 a 2% $95” DPP ”5ARAQBAD)
«r3 = -zsaism sub Timur.
ABD(9 K+2 2D 1 + T2 + T3

= ins: Dm *RAT/AD

T2 : :Tégafgseéagéwm
36 K - T1 + T2 4- T3

533,5 2’4"” '
§DZ112§§P 1 + T2
T2 =- zssp§SR,Dp;‘*Y(K)/DENK

O 000

101

183

Tit-5335:; 592?" %§& mflR‘E‘RT

T6 = -ZSEPS 113/DEM“ Y(K+2)/(AD*RAT

ABD(8,K +3)=PT% +T4+T +T6

T1 = -2Ds SR,

T2 = 205? SRIDR *Y K) IEmix/Diana?”
T3=ZwPSPDP*YK+1yD

T4 = -2DS R mm

T5 . -ZDDP SR, DR *Y K+4 /DENK

T6 = -ZDBP SRR DRR Y(K+2)/(AD*RA

ABD9K+3T1+2+T3+T4+T+T6

ABD1RéK+ £21): -ZDSMD)/(E:*RAT*AD)

T 1
T2 = ZSSD SR, DR *YK
$31 a ZSDD SR' {:DP *Y ”135%
T5 = -ZSDD SR, DR *YK
T6 = -ZSED SRR DRR/)1.2 (K+2)/(AD*RA
ARD(D7,K+4)=T1+ +T3+T4+T +T6
T1 = -ZDD SR, DR /D
T2 =- ZDSDS RIDR *YK )K/DE NK
T3 = ZDDD SR, DR *Y K+1)/DENK
T4 -2DSD SR, DR *YK K+33/DENK
g =- -Z-I uSR'D

a -2DBRS SRRD Rzz (K+2)/(AD* RATE

ABD8,K+4 TiD+ 2+T3 3T6+T4+T+
ABD 9,K+4

521.30 (S, D) / ( EPS*RAT*AD

5
3‘
2‘.
5 a;
IIIIIIIO)

0..

ASSEMBLIMSTHEMATRIXEWARISINGFRCMIHEW
CCNDITIOINB AT X=A

B§2==D§EL RRR11{§RBL(R-1)/2. a
S a Y 3(RR-
D . Y -1 ‘

-zss S, D/D 2+ *ZSEP $14,114 /D2 RR
-Zm28: DI/flRz 2+ 5 ZDEP SM ,mI/ngwi RRI
-23DS,/DD 2+ *zsman,m/D2*YRR

is: 62% 22+ 2 i an?

.5
.*5
.5
.*ZIED5 Slhm/DZWRR

a

U
R E
lllllllllllllllll

 

F
éSE:(S§:fi D D:*Y1RR)Z£2. 6*D2)
§D(11,RR-2$ + T3 + T4 + T5 + T6
52 : 21% ”(ESéD) /D22
T2=zDRRSD*YRR/ 2.a*D2
T3 - -2w (1%, R) 4 (RR-iszz I

1
102

000000000000 0000

1'3

1814

_—_- * _
$4 = g8§PES,D ,D YéRR 2)/D22
$3 = IRl PP(R2* Y

= +
T6D = -2DSR S ,R) ;R(RR;1R)/D22

WI g= +B+M+fl+fi

D13I RR-z =
T = 2
$3 3 2:3 SéD/ 22
T4 = -ZSSR R, /*Y RR-Sg/DZZ
T5 = SDDS *Y RR-4
$9 3 ZSSERS,6E*YRRR- RR 22
T8 = -zs (é, )* (RR- )/ERS
T9 = ZSED(S D)*YiRR)/ 2. 2 *D2)
fiW-Qf 15 +n+w+m+m+m+m+m
T2 =
$2 3 29D SzD/ /D22
T5 = -2DSR(S,/EPS )*Y(RR-5)/D22
$9 3 EESDDDD, S bDlyYéRRR RR-2 ‘4’ 232
T8 = 22 MR *Y RR-l /D22
T9 = (R, ;* ERR- )/ERs

T10= 2DED s, D *Y lRR)r/(2.E; D2)
gmz’m ‘3: 2+T +M+fl+m+w+m+m+na
ABD 9 éRR = 25 S E/ 2. E*D2) + 23 SM, DM / 2. 0*D2
1,RR) RDR(R,R) }(2. 2*D2) + ZDR(SM, DR)}(2.E*D )

IF M .NE. MMM [GO To 123

D01

WRITE2 (61,1 iRBD(IM,JM)
%Tbé2x,

CONTINUE

RETURN

END

SUBROUTINE RIBVN (R,IEL,U1‘, Y ,B)

WIEN IFLAG=0 TIE FIRST TIME RIBV IS CALLED AT A GIVEN TIMER
SUBmmE CQIPUTES 'IIERHS VECTOR FORTHE H%TION) AT
'I'HEN-I-lTHTmERONFROI'I'I-IESOLUTIONTOTHEMA'IRIX

TIE N TH TIME RON. WHENIFIAG=1, RHSV COIPU'I'ES'IIE 'IIECURRENI‘
VALUEOF'ITEAUXIIARYVEC'IOR. INTIEMAINPRCXSRAM'I‘I‘IISVEC‘IOR
WILLBEAIJDEDTOTIERHSVECTORTOFQQM'HECURRENTRHSVEXHOR
FOR A GIVEN ITERATION. R IS THE NLMBER 0F SPACE PTS.

TH IME RON IF IFLAG=G AND
PREVIOUS ITERATION STEP IF IF LAG=1. B RR) IS THE RIGHT HAND SIDE
VECTOR IF IAG=0AND ANDTIE AUXILARY VEC IFIF LAG=1.
THE (NLY DIFFEREI‘CmE‘S’o BE'IWEEN THE RIB VECTORAND THE TEAUX VECTOR
AREmTI-IE'I'EMSDNO HEVINGDI‘ANDINTIEBGJNDARY MS.

REALIGDOKSEJQDLIISSLWKSAKDAKDRKSR

INTEGER
DTMENSIDN' DE 32 Y 62 B 322

/Rc::TR/ KRéiRSé,I<15L(I(SL),KDA,KSA, KDR,,KSR,SR Siam, DL
CDMMDN/ /ELEC/ EPS,ZI

000

00000

185

RR=3*R
RS'IOP3RR-5

vm ELEMENTS ARISIMS FRCM THE Bum CWDITIOBB AT X30

D2 =- DEL(1 /2.0
D228 ay 811??“ )*DEL(1)/2. 0

DaYé
353—{5/02
ms-Kéa 302
P3a-z

1‘1: 21322442323?

$3 45 41363.8: ASZ %:NZ%SE(§§%§5 912859432 2)*) (31m)

 

 

T5 = ZSD S,D /D22 *Y 54
B 1) 8 T1 + T ‘1' T5
‘1' =- (-KD02 - +zm(s, D)4/D22)*Y(1)

2‘2 Ksa/Dz
S2 " 2%.”?232965’ 22
'12: 214/22 3+F4)*Y2

 

 

$2 = guns /2. 6*02) + 2&(39. DP)/(2. amznwm
. . ... s 4.22;. .

$32; =1032/ saw? 5 T4 T5

T . -§:ng/ 5M5 )

3 3): 21 + +1-

VEC'KR Em FRO! THE INTERIOR SPACE POINTS
THE m-LmP IS IN MULTIPLES OF 3 SINCE ‘IHERE
ARE 3 EGJATIOBB AT EACH SPACE POINT

91 K=4 ,RSTOP,3

 

 

 

 

 

 

no 1

22 " 2:2 53

34 = Y K-3 + Yam/2.2
em = Y K-3

5 =- Y n+3

39': 5:) x+3 YK) /2.a
w= YK—Z YK+))/20
B": 2 2.2?

ngan 1;: + Y(K+1))/2.a
Ian-02:01:24 +03 ADéKK)
222 =22 $22212 3‘
mu. 2 DE *DEL

an a DEL 81:“) - 0:1}.(131)

F1 1. 0
F2 3 -2 SP, DP /DENK

m-4$m
T4 - (n + Film + 3)*Y(K)

000

000

000 '

000

186

F122: Iéé'sig'fiwfifi
2F131>+S'Fé/EPS+ F3 2(*Y( +1

= ZSD sp Pinpgwi Rig/Dam

VECTOR ELEMENTS mm s summon

BiK)=T1+T2+T3+T4+T5+T6+T7
T =- we SM,m)*Y(K-3 BN4

)
$2. : 22. 22.22222 3.“...
T4 T5; Exzaoép, W13? +4i/6Dm

T63-ZDDSP,DP*YK+

ZDE(S D5
(F1 + ’F ”Yaw-1%
zusuésp pr *Y [(+33 )
-z 13.335022. Ma
21238 pp ,Dpp)/(RAT*AD)
208 s, ,Défs UBAMULé
T10 = (F 2 + 3)*Y(K+2)
m swarms mm D EQUATION
8991) =T1+T2+T3*+ mgv+TS+ +T6+T7+T8+T9+T10
a -2Ds s, D: K-3 EPS*AD;
= -zmé§lg Diwlgx-zz/ 3332;422:929}
T3 = zns *MULLL

=21!) ,D* *S
3%; -13' 126312 m/éps)m .1)

T68 m S, D *YéIzG-Bg/ S‘EPS PS*RA'NAD
Wflwa,D*YKM/EPS*RAT*AD

WWMBEGIATIW
BK+2 =T1+T2+T3+T4+TS+T6+T7

351339
“I

{38333333
I II II II II II II

101 C

vacrm Em ARISIMS E'Rm 'IHE BGJNDARY CQIDITIOIS AT XSA
D 8 =DEL éR-l 0
éR-1{£DEL(R-1)/2. 0

S 3 RR.

§L YEAR 1 ’

3132235131ng (RR -5;/022

T2= 28D 3933' *Y -4/D22

F1 = -1.

F2 = -

F3 ... -z 88(2DS

T3 ... (F i’éz + m:3)*Y(RR—2)

F1 = - 2

F2 . -ZSD s D)

F3 = 233( .15

$15: a (1.25st ,D) 2. 3)2'132)“ + 4zsr: SM 130/(2. 0*D2))*Y(RR)
322-... 2, .352 2.43" ..2
T a 2m ”(ERR-5 T/D22T4

T2 = no s'D *Yf ran-4 /D22

5» ‘ ’2 1?ng

F1:

15+ 8&2 *Y(RR-2)

JUN

0000

100

00000

187

-Zw 82D)2

3222138 (F1(+ D5fi3+g4w
§3£WI(ZDB3S,D)432.B*D2;+izw3é‘13M/(2oWDZHWmR)

gmgrfRR-l)

13% ‘22 ' "

SIBRUJTINE RIBVO (R,IEL,UI‘,Y,B,M,M4)

SUBROUI‘INEREBVOGENERA’IES'IHEPORTION OFTHEBVEII'IURTEETIDES
MGMWITHEACHITERATION.

REAL KDflfiKSO, Lm'II'DKSL' ,KDA ,KSA,KDR,IGR

mmm
81113125 ION EgméapsgfiéixgggkL ,9ch ,1<DRI®A,1<SA, ,xsmsmsmmmr.

‘2’ 36
T35- L*3SL + 1.fl)/EPS
T = KSL D

Bf=n+m+n+u+fi
WI mm

EPf“. mommun,nmn,1

T'( .314. 6, 15)

mm(1
IF {Tum .éo. a) B(I)-B(I)+ZI/EPS
3mg 0.3

quama

T4

m4

Mm=n+m

DID

SIBRQITINB GRIDPR (DEL,XER,SIZE,‘L)

GRIDPR GENERA‘IES A NONINIFGW SPACE MESH WITH IICREASING
DISME BETWEEN GRID POINTS AS X II‘CREASES. IT IS
DESIM FOR SIMULATIONS USING A POLARIZABLE ELECIRGE
AT XIO AND A REVERSIBLE ELEC'IRGE AT X81.

000

10

20

30

4E

56

60

0000000

10

\IH

188

THEVALUECFFAMSESHOULDBESEIEC'IEDWITHCARETO
AVOIDLNWANDOVERFWFM‘IHEWIALFINCTION

DIMENSION DELUIG) ,XER(166)
INI'EGBR R

H/UNITS/ EPMA ZPL 214 cm
0%: /p SKI-3‘9, n.6,ur IF I'mf. ,RANGEJSET

 

no
DEL L)= DEL(L)*SIZE/(A*sm)
c mun
mag?” =- a. 6R

axm'zm-n + DEL(M-l)

c

m6 8
flNgN =- A*kER(N)
REI'URN

mm AVINI' (X,Y,N,XLO,XUP, IND)

AVINT IS USED TO mm THE Nommmm SPACE MESH
ma cmer'rmN OF THE porENrIAL DIFFEREMZE m1. CHARGE
IN THE sysma AND 'I'OI‘AL moms IN THE sysmfi
IFINDISRBI'ORNEDAs-l xup IS LESS 111mm

AND m cum'A'rION WILL NOT moot-ma

DIME'NSIQQ X Y
I (N): (N)
IF N iLT. N3) Rm

IEWI‘MZ .IE.X X(I-1)) RETURN
SWO. 0

IF (XLO .LE. XUP) GO TO 5
SYIFXUP

IX) l 181,}!
ggiXU) .GE. XLO) GO TO 7

IB=MAXG 2 ,IB)
IBIMINO IB,N-1)
II) 2 I3 ,N

0 00000

0000

14

15

170
160

189

§F_1(XUPJ .GE.1X(J)) GO TO 8

m3

X1=X EJM-1)J
X2=X JM)
X38 JM+

TERM}- Y JM+1 3-Xl) * (X3-X2) )

gal—x X2+X3 W14X$+X3EW- ‘XI+X2)*TER‘|3

IF (JM3 .GI‘. IB)GO

T‘s-gig 3%} 35.32 019%,?“
MM

,,,;gfiii

15

M. 5* RICA

C830. 5* B-l-CB

CW. 5* C+CC

SWSWOCM (X2**3-SYL**3) / 3 . 0+CB*0 . 5* (X2**2-SYL**2)+CC* (XZ-SYL)
CBIB
CC=C
SYLFXZ
A WiJSJW3-SYLW3V3. ..+CB* 5* (XUP**2-SYL**2)-I-CC*(XUP-SYL)

(IND .EQ. RETURN

SYL=XUP
XUPBXLO
XLO=SYL
AVIN'b-AVINT
END

SlBROU'I‘INE IC (KEY, Y,R)

SUBROUTINE IC ASSWIES 'IHE INITIAL CODI'I‘IODB m: S AND D
IFKEYISIASTEPFLNCTION NFCRSIS INITIALIZED
ANDDALDEAREINITIAIJEZEDASZERO EVERM‘IERE
HKEYISZS,D,AWE1mEsmm2m0Wm

DIMEINBIQI Y iS300)
INTEGERRR

RR-B*R
II) 160 I=1,RR,3

Y I+l 8% 0 \
Y 1+2 3%. 0

I (I .13 . W2) GOT0170
EAIEIJ

@493
INUE

PORANODDVALUEOFNI THIS INITIALIZES‘IHEVAIDECFSATTEE
mwmcmmbm

 

IF
IF
m 40 I=:-1,RR,§
Y(I

”Elem? .Né? so”) HER/$353
)-

14B
150

00000000000000

000

100

000 0 000

101

000

00000 0 00

102

190

ETA; . A a
CMINUE
CWTINUE
RETURN

END

SIBRQJTINE O‘lPR (Y,YY,TOL,R,KFLAG)

SUBROUTINE OHPR COTPARES THE CURRENT ITEEATION VECTOR TO THE
CURRENT SOLUTION VECTOR mm. (n+2 IS RENAMED YY DESTROYINC TEE
OLD SOLUTION VECTOR THE CURRENT ITERATION VECTOR IS TESTED TO
DETERMINE IFITISSMALLENOLKEHFGR‘I'HEPRmRAMTOPROCEDETOTI-E
NEXT TIME STEP. SINCE SOTE OE THE SOLUTION VECTOR cmPONENTS ARE
MUCHSMALLERANDHENCE LESSACCURATETHENOTHERS THEYARENOT
INCLUDED IN THE TESTING. EACH OE THE 3 COMPONENTS MAKING UP THE
SOLUTION VECTOR S D AND E ARE CHECKED SEPARATELY. MEANS
THE PRESENT ITERATION VECTOR ENOUGH PCB THE
PROCEDE TO THE NEXT TIME STEP. KELAGsI MEANS TEE PROGRAM WILL
ITERATE AGAIN WITH THE DERIVATIVES EVALUATED AT YY+Y. TOL IS THE
RELATIVE CHAbGE TOIERENCE ALLOWED.

DIMENSION! 3m YY 3m
mmmg )' ( )

RR=3*R

mm; THE NEW SOLUTION VECTOR

no me I=1 RR
YY§I gunman)

M-I S'IRIPSOUI‘THBSCW, M-2THEDCOTPONENTS, M-3THEE

CGIPONENTS.

mAx=AES(yV(M))

EINDINO THE MAXIMLM VALUE OE 5 (IE Mal) ,D (IE M-2) , m E _(IF M-3)

DO 101 I=M, m3

YMAX=AMAX1(ABS(YY(I)) ,YMAX)

C(NTINUE

FINDING THE MINIMLM VALUE THAT WILL BE CHECKED

YMIN= DEB-09
II) 132 JIM, ,3

TESTIMTODETEMINE IFAVALUE SHCXJIDBECHECKED
IF (ABS(YY(J)) .LT. YMIN) GO TO 1'2

CEECKIMSTOIE’IERKINE IF'IHECHABGEINTHESOLUI‘IONVECMIS
WMATAGIVENJ

THE FOLLOWED 2 SEW PREVENT DIVISION BY ZERO
AY=ABS(YY(J))

RAT10=Y.I(‘T7YY.a§-99) GO TO 102

IF (TEST(.G1‘. L) GO TO 104

C(NI‘INUE

191

103 CQITINUE

RETURN

104 KFLAGII
RETURN
END

 

.CCC . ..

mw: a aflflhflfia
mpwmmmmnlllxmm

C12 8362

 

mom (m.
mmdsm

***

o);
GT135P l

 

mm

 

192

3; - C(8)/SQRT(S + 1.3)
A
+
D +
’61

E?

gauge:
mmmmmm

mm

 

 

m
C O O
11 1 .
+1.9 + )
o ) a
DUI MW 5. L
. ++ ... 1 +
sss a... .m.«. a... + a” n
()1! 11“.» 1 D 1 —
mmmw ++a. + + + S
. C DDZ D S D (
* * ) (
.m“+) .++5 + w/ . MW
1((\la SSII. S .\I S ...m
CC3 . ((9 . l\ M ( Z
+ ... 1 .I\ / M /
air, .... m + m +
L90 .8 ... 1. +m + ...m”
. o oll‘ \I\IS o — + )2
)311+ 24(1 J“ J(
$4313 mm + m 1. M 1+
DIN/DD“) 55mgmw + + + +
S Cu . + mw 1191. P Dm WP DZ
P 1..SS ++(8C( mm + Z“ . _
S * .. t . )1. C 8+ _ S
w L.CCCIxGCClLGWWIJ 2NW3JW4
S:==CC===:.Z UIZU. U1.U
(W ma=mmw Masgaasasssaa
l 12 3 34 4
cmzmzmmz 2W UWTUWflUWTUWT

193

2881”3 T1‘+ T2 - T3 - T4

MHWW
am

 

 

 

 

 

”3 3H “3
*u u* *u
902 m 002 m 002
o o o o 0* o o 0*
11 1 11* 1 11*
++J + ++J + ++J
DDl $7 DD1 5 D01
.++ * .++ * .++
SSSJJW SSSJJW SSS JJ
((illm (Nill ((5 11
ammm++a a 6++0 0mm“ ++
l o 1 o 1
. SCDDZ) . SCDDZ) . SC DD
PP+_++W ) WPV+.++W ) PV+ .+
024 m 024 m 024 )
1E)SS\I~L lEgS\lz 1E)“ $8
+ * o+ + *T o+ + *1
a++* 1 . 0++P 1 + E++V+
021g o‘l‘l o)\-l
1aa.w$+p m loamm$+ m looms ..
.LL1”4(J 2 Z 2 .LL124(W 2 Z 2 .LL1(.24
Qh++wwm W a us a flh++hmmmul u N u Qw++wm mm
um 835% m m m n. (2 Bug L m mm mm new. Bu) .3
S m.+ 1131 W1 S m.+ 113+ I“ S J.+ mnll
P 188 ++(+m+n+ +W lss ++( +n+ +Wfi P 138 ++
_ C 2 . . CD2 m .
w ***\I\I*D D D+D ***\l\‘ D D-D + ***- )\.I
2 oil- 5 - + +1 2 Giza—.3. .... a? + z 0;) 13
Tang c. L. . 1+ m+2m.. wage. c.1.s .1+m+2m ($8,511))9c. c:
.CCC =.$MSUSJ$W. .CCC (WSUSdfiw. m .CCCCC..

ml\ 3 ..P
mmmmmmmmm

SW

_/
mmwg==h gssaszsspm mmw:===aah7ausaaxaam
mmmmmmmmmmmnmwmmm m mmmmmmmmmmmmnmnmnz m

191-l

. . * * O
ZPMP 1.§E?é?l?§Q§T éofiiig 332 0 C(6) (S + 1 0)

W805

ABM.G

Ul=-ZPP+ZH‘1

V1= +ABP+(AM+AB)/(S-D+1.$)
T1801 1

V2 = -ZPPP + ZPMP + (ZPP - ZPM)/(S - D + 1.0)
T2=UZ*V2

U3=ZPM-Zm
V3=MP-ABP+(AM-AB)/(S+D+l.0)
T2=U3 3

v4=ZPMP-zm1>+ zm-ZPM)/(S+D+1.0)
T4=U4*V4

ZIEPSTl-i-TZ-T3-T4

END

9; **3
6 **3
2

 

G
\

+ 2PM
0/ S - D + 1.0)**2

8.53%

NO

END
§%§S+D+1.m+m-AB

NCl-‘NCD-‘NC HNOO II II

SEESSSSSSSSSEZ
“llllllllllllllllllllllllll

fl??? D + 1 awrz
u4*v" 4 '
ZSSD=T1+T2+T3+T4
RETURN

33'"?
:r'6
.1.

C NN
“55
II ll
O
War-(.3 1' l- t
N V
I
HH
0 o b-b-s
GU“ A
V mm
.1328
M

555355

195

$303 9%1'45’35/"2
RENRN

 

COMQJ/
IF (S .13. -1.a) S a -1.0 + 1.08-la
ZPP=C1*S-D+1.0g+C$2} s-D+1.a;**3
zm=c3*s+o+1.a +C4* S+D+1.0**3
ZPM=C5* ((S+1.0;**3)+C6) s+1.a)**2
ZPPD=-C1)-1.5*C§2;* Tés' + .0)
zmn=-C()+1.5*C$ (+D+1.a)
Vls-PP-ZPM
'r = U1+v1)
%- S+D;m.a)*m
w2= S+D+ .m**2
T2: UZ+V2
zsw-T1+T
RETURN
END
WIONZIEDSD
coum/ 6(4)
ZPP=C1*S-D+1.0 +c * S-D+1.0**g
zm=c3*s+o+1.I +C *S ‘r s+n+1.a**
ZPM=CS* ((S+1.¢£**3)+C6) S+1.0)**2
ZPms—C§1)-1.S*C523‘* és- + .0)
m=C()+1.§*C) (+D+1.I) ~
AM=C7 -C(8) (3+1.a)
fiac9
name“
3%”?ngng lfl)**2

= o - + 0
33:0?“
¥2=$.%SS+D+1.0)**2
znsna 1+T2+T3+T4
END

mm 2339 s D)
cm Nun's éq GDE,EPM,A,ZP,ZM,CCN
ccmou {G (3&4)

IF (S c To -10”) S a -100 + 1.03-10
zppp . cu) + 1.s*C(2)*SQRT(s - n + 1.91)

196

E'é‘z‘fi : §$3§4E$駧$filifl§$§ 1’ ‘z’pfap 178211

END

mum 21189 S D) x
comm mum's én GDE ,EpM,A,zp,zn,cm
ccmw c; C(é4

IF [gs .mcs = -1.a+1.aE-1a
299 ='cTi +1. S*C2*SQRTS-D+
map: + 1. 5*c 4 *SQRTS
szp 1. 5

21189 : 3. 894é—éGSgfi*A*A*(ZPP + 21149 pitgzfitfiflm

END

SED S D
Cwmm% ITS}C S36”: ,EPM,A,ZP,m,C(N

COMO! /GREEK/

ZPPD= :ngZéJ.)+-11(:(S *Cé £21§$LEDS£ D+1+1.0)
mamas-03914111111 lZDPPD
RETURN

END

nuc'rxou 2). -s D)
comm Nun's: éa éGm, EPM ,A,zp,zm,cm
comm /GREEK/

ZPPD=-C1:S§*&i£2*ST +115“)
man-+1.5C DD++)
ZDEDS 3x8 4E-06*A1A* ZPPD

RETURN

END

APPENDIX F

SOLUTION OF SYSTEMS OF LINEAR PARTIAL
DIFFERENTIAL EQUATIONS

In order to solve Equations 5-6 for S and D, they are

expressed in matrix form as
pg = o , (F-l)

where E is the linear operator matrix and Q is the solu-

0

tion vector.

If w is defined such that
DET(g)w = o (F-2)
then it can be shown that,
Ui = COFJi(P)w , (F-3)
where J can be any i (the cofactor can be expanded along

any row).

For our system,

197

198

(3/3t) - YSS (a2/ax2) 'YSD(32/3x2) + YSE

(32/3x2) (a/at)-YDD(aZ/ax2)+yDE ,
(F-u)

‘YDs

C}
II

The equation for w is

[A(3u/8xu)+B(33/3t X2)+C(32/3x2)K32/3t2)+G(3/3t)1W = 0,

(F-S)
~where
A = YssYDD ‘ YSDYDS’ B = “(Yss + YDD) ’
(F-6)
C G:

= YSEYDS ' YDEYss ’ YDE '

Since W, in this case, must be an odd function about x = O,

the solution is expanded in a Fourier sine series,

W(x,t) = f, bk(t) sin (kflx/a) . (F-7)
k=1

Substitution of Equation (F-7) into Equation (F-S) produces

the ordinary differential equation

199
bg<t> + [G - <kn/a)2BJ bi(t) + [<kw/a>”A -
- (kn/a)203bk = o (F-8)
The solutions of Equation (F-S) are
bk(t) = L+k exp(-R+kt) + L_k exp(-R_kt) , (F-9)

where R+k and R_k are given by Equation 5—8, and L+k and
L-k are determined by using Equation (F-3) and the initial
conditions. The result is Equation (5-7).

APPENDIX G

ANALYTICAL SOLUTION OF D FOR TIME DEPENDENT

DOUBLE LAYER RELAXATION

The time dependent system of equations in Section C of
Chapter 6 is solved using Laplace transforms. Because this
is a time dependent problem the following set of reduced

variables is used in this appendix for simplicity.

d = D/SB E = x/a
” = YDDt/32 k = a2YDE/YDD (G‘l)

a<n) = Q(t)/(2aFSB)

In terms of these variables the system becomes

3% = iégf-- k2d -1/2 S g 5 1/2 , (G-2)
%% = -k2§(n) a = - 1/2 , (G-3)
d = 0 g = 1/2 , (G-u)
d(€,o) = o - 1/2 5 a g 1/2 . (G-S)

200

201

Duhamel's Theorem (Carslaw and Jaeger [1959]) states

for the above system that
D(€,n) = §%{4E U(€,A,n-X)dx} (G-6)

where U is the solution of:

%% = §‘%T - k2U - 1/2 5 g g 1/2 , (G—7)
35

3% = -k2a<1) = — 1/2 , (c_a)

U = o g = 1/2 , (G-9)

U(€,o) = o - 1/2 5 g g 1/2 . (G-lO)

The Laplace transform L of d is defined as

(D

G(a,s) = L{d} = Jr d(€,n)e'snon . (G-ll)
o
The inverse transorm, L-l, of G is
-1 1 F+iw ns
d(€,n) = L {G} = EFT e G(£,s)ds , (G-l2)
T-im

where P is chosen to be so large that all of the singulari-
ties of G lie to the left of the line (P-iw, P+im). In

practice many Laplace transforms and their inverses have

202

been tabulated. The ones used in this appendix appear in
either Carslaw and Jaeger [1959] or Roberts and Kaufman
[1966].

Transformation of the reduced system of equations

yields
2G
1.5 _ (s + k2)G = 0 - 1/2 5 g 5 1/2 , (G-13)
36
3G - k2 _ 1
SE — - E? q(1) g — - 1/2 , (G-l )
G = o g = 1/2 . (G-15)

G can now be easily solved as a function of s and

g. The solution is

-§(1)k2sinh[(s+k2)1/2(g-%)J
s(s+k2)l/2cosh[(s+k2)l/2]

 

G(€,s) = (G-l6)
If this function is transformed back to a function of

n using the inverse Laplace transform or using Laplace

transform tables, the result is identical to that which

would have been obtained if the problem were solved using

Fourier series. Thus there would be a practical difficulty

in computing the value of the series for small values of

n . Instead the Binomial Theorem (Beyer [1976]) is ap-

plied to Equation (G-l6) before G is inverse-transformed.

After some algebra and before inversion the expression is

203

 

2-
_ -k q(1) 2 1/2
G(€,S) - {exp[(s+k ) ( -3/2)] -
S(s‘gk2)172 5
(G-17)
expE-(s+k2)1/2(g+l/2)]}ZO(-l)nexp[-(s+k2)n]
n:

This is then inverse-transformed using several sets of
identities provided by Roberts and Kaufman [l966] and
evaluation of two integrals which appear in Chapter 7 of
Abramowitz and Stegun [1970). The expression for D in

terms of x and t becomes

an

t
D(x,t) = g% {(k/uaF) JFQ(1) O(-1)n{-exp(2rn) x
O

n:

1/2(t-1)‘1/2 1/2(t-1)1/2

erchrnyDE +YDE 1+8Xp(-2Pn) X

erfcfrnyfié/2(t-A)-l/2-y%é2(t-l)1/2]+exp(2pn) x (G-18)

-l/2 -1/2 1/2

J-exp(-2pn) x
erfCEPnyaé/2(t-A)-l/2-y%é2(t-A)l/2]}d1} :
where

rn = k(un+3-2x/a)/A, pr1 = k(un+l-2x/a)/u , (G-l9)

and the complementary error function, erfc, is defined by

204

m 2
erfc(z) = (2/nl/2) .f e.V dv . (G-20)
2

For this system, Leibnitz's rule (Greenberg [1978])

can be expressed by the equation,

t

t .
595% r(1,t)cu] = 43' (sf/atm + f(t,t) . (3-21)

Application to Equation (G-18) yields

m t
D(x,t) = (kYDE/Zale/Z) 20 H)“ f cum
n= 0

-[YDE(t‘A)J-l/ZGXPE‘riYB%(t‘A)-1'YDE(t-A)1 (G'22)

'l/2expt-PiYBE(t'A)-1-YDE(t—A)]}dx

+[YDE(t-A)]
This equation can now be put in a compact form by transform-
ing the dummy integration variable A to U with the rela-

tionship

= YDE(t — 1) . (G-23)

The result is the solution for D given in Section C of

Chapter 6.

REFERENCES

REFERENCES

Abramowitz, M., and I. A. Stegun, Handbook of Mathematical
Functions (Dover, New York, 1970).

Agar, J. N., and W. G. Breck, Trans. Faraday 800., ii,
167 (1957). ‘

Agnew, A., and R. Paterson, J. Chem. Soc. Faraday I, 13,
2896 (1978). .

Anson, F. C., R. F. Martin, and C. Yarnitzky, J. Phys. Chem.,
11;. 1835 (1969).

Armand, M. J. J. M. Chabagno, and M. Duclot, 2nd Int.
Meet. Solid Electrolytes, Scotland (1978).

Barker, G. C., Transactions of the Symposium of Electrode
Processes, E. Yaeger, Ed., (Wiley, New York, 1961).

Bearman, R. J., J. Chem. Phys., gg, 585 (195A).

Bellman, R., Perturbation Techniques in Mathematics,
Physics, and Engineering (Holt, Rinehart, and Winston,
New York, 196“).

Beyer, w. H., Ed., Standard Mathematical Tables, 2Ath
Edition (The Chemical Rubber Co., Ohio, 1976).

Bockris, J. O'M., and A. K. N. Reddy, Modern Electro-
chemistry (Plenum, New York, 1977).

 

Brady, J. F., and J. C. R. Turner, J. Chem. Soc. Faraday
I. 13, 2839 (1978).

Brown, A. 8., and D. A.MacInnes, J. Phys. Chem., El,
1356 (1936).

Brumleve, T. R., and R. P. Buck, J. Electroanal. Chem.,

Buck, R. P., Critical Reviews in Analytical Chemistry, i,
323 (1976)?

Carnahan, R., H. A. Luther, and J. O. Wilkes, Applied
Numerical Methods (Wiley, New York, 1969).

 

205

206

Carslaw, H. S., and J. C. Jaeger, Conduction of Heat in
Solids (Oxford University Press, London, 1959).

Chang, H.-C., and G. Jaffé, J. Chem. Phys., 20, 1071 (1952).
Chapman, D. L., Phil. Mag. 25, A75 (1913).

Chen, H.-C., and H. S. Frank, J. Phys. Chem., 11, 15A0
1540 (1973)-

Chen, J.-S., Ph.D. Dissertation, Michigan State University,
1971.

Cohen, H., and J. W. Cooley, BiOphys. J., 5, 1U5 (1965).
Conway, B. E., J. Sol. Chem., 7, 721 (1978).
Daum, P. H. and C. G. Enke, Analyt. Chem., 31, 653 (1969).

deGroot, S. R., and P. Mazur, Non-Equilibrium Thermodynamics,
(North-Holland, Amsterdam,II969).‘

Delahay, P., Double Layer and Electrode Kinetics (Inter-
science, New York, 1965).

Delahay, P., and A. Aramata, J. Phys. Chem., 66, 2208
(1962).

DeMari, A., Solid-St. Electron., 11, 33 (1968).
DeMari, A., Solid-St. Electron., 11, 1021 (1968).

Dongarra, J. J., J. R. Bunch, C. B. Moler, and G. W. Stewart,
Linpack (SIAM, New York, 1979).

Duckworth, L., Electricity and Magnetism (Holt, Rinehart,
and Winston, New York, 1960).

Feldberg, S. W., J. Phys. Chem., 13, 87 (1970).

Fitts, D. D., Nonequilibrium Thermodynamics (McGraw-Hill,
New York, 1962). .

Franceschetti, D. R. and J. R. MacDonald, J. Applied Phys.
i9. 291 (1979).

Frank, H. S., J. Phys. Chem., 61, 155A (1963).

Goldberg, R. N., and H. S. Frank, J. Phys. Chem., 16,
175 (1972)-

Gouy, G., J. Phys. Radium, 2, A57 (1910).

207

Grahame, D. C., Chem. Rev., 51, “H1 (19u7).

Greenberg, M. D., Foundations of Applied Mathematics
(Prentice-Hall, Englewood Cliffs, 1978).

 

Guggenheim, E. A., J. Phys. Chem., 33, 8H2 (1929).
Gummel, H. K., IEEE Trans. Electron Devices, 11, “55 (196M).
Hafemann, D. R., J. Phys. Chem., 69, U226 (1965).

Hall, D. G., Trans. Faraday Soc. 11, 2516 (1971).

Hall, D. G., J. Chem. Soc. Faraday II, 62, 1391 (1973).
Hall, D. G., J. Chem. Soc. Faraday II, 13, 101 (1977).
Hall, D. G., Faraday II, 16, 125“ (1980).

Harned, H. S., and B. B. Owen, Physical Chemistry of
Electrolyte Solutions (Reinhold, New York, 1950).

Hasse, R., Thermodynamics of Irreversible Processes (Ad-
dison-Wesley, Reading, MA, 1969).

 

Hickman, H. J., Chem. Eng. 801., 26, 38l (1970).

Horne, F. H., and T. G. Anderson, J. Chem. Phys., 63,
2371 (1970).

Horne, F. H.,‘in Physics of Superionic Conductors and
Electrode Materials, J. W. Perram, S. deLeeuw, and
B. C. H. Steele, Eds., (Plenum, New York, 1981).

 

 

Horne, F. H., J. H. Leckey, and J. W. Perram, in Physics
of Superionic Conductors and Electrode Materials,
J. W. Perram, S. deLeeuw, and B. C. H. Steele, Eds.,
(Plenum, New York, 1981).

Ingle, S. E., and F. H. Horne, J. Chem. Phys., 62, 5882
(1973).

Jalota, S. K., and S. Paterson, J. Chem. Soc. Faraday I,
£2. 1510 (1973).

Kahn, D., and J. N. Maycock, J. Chem. Phys., 56, ““3“
(1966).

Kaufman, H., and G. E. Roberts, Table of Laplace Transforms
(Saunders, Philadelphia, 1966).

 

KooiJman, D. J., and J. H. Sluyters, Electrochim Acta,
13, 1579 (1967).

208

Laser, D., and A. J. Bard, J. Electro. Chem. Soc., 123,
(828 (1976). ———

Leckey, J. H., and F. H. Horne, J. Phys. Chem., 66, 2504
(1981).

Lewis, G. N., and M. Randall, Thermodynamics (McGraw-Hill,
New York, 1961).

 

MacDonald, R., Phys. Rev., 91, 912 (1953).
MacGillivray, A. D., J. Chem. Phys., 63, 3126 (1970).

MacInnes, D. A., The Principles of Electrochemistry (Rein-
hold, New York, 1939):

 

Miller, D. G., J. Phys. Chem., 19, 2639 (1966).

Miller, D. G., J. Phys. Chem., 71, 616 (1967).

Miller, D. G., J. Phys. Chem., _l, 3588 (1967).

Moore, W. J., Physical Chemistry (Prentice-Hall, Englewood
Cliffs, 1972). ‘

 

Newman, J. S., J. Phys. Chem., 73, 1843 (1969).

Newman, J. S., Electrochemical Systems (Prentice-Hall,
Englewood Cliffs, 1973).

 

Olivares, W., and D. A. McQuarrie, BiOphys. J., 16, 143
(1975).

Onsager, L., Phys. Rev., 66, 2265 (1931).
Onsager, L., Phys. Rev., 61, 405 (1931).
Payne, R., J. Electroanal. Chem., 111, 277 (1973).

Perram, J. W., and M. N. Barber, Molecular Physics, g6,
131 (1974).

Perram, J. w., J. Math. Chem., Q, 37 (1980).

Philip, J. R., and R. A. Wooding, J. Chem. Phys. 6g, 953
(1970).

Pikal, M. J., J. Phys. Chem., 16, 3124 (1971).

Prigogine, I., and R. Defay, Chemical Thermodynamics, D.
E. Everett, Trans., (Longmans, London, 1954)}

 

209

Reinmuth, W. H., and C. E. Wilson, Analyt. Chem., 66,
1159 (1962).

Reinmuth, W. H., J. Analyt. Chem., 63, 1446 (1962).

Roberts, G. E., and H. Kaufman, Table of Laplace Transforms
(Saunders, Philadelphia, 1966). '

Rosenberg, D. U., Methods for the Numerical Solution of
Partial Differential Equations (American Elsevier,
New York, 1969).

 

Sandifer, J. R., and R. P. Buck, J. Phys. Chem., 12, 384
‘(1976).
Sch6nert, H., J. Phys. Chem., 73, 62 (1969).

Schryer, N. L., J. Phys. Chem., 81, 2335 (1977).

Shedlovsky, T., and D. A. MacInnes, J. Phys. Chem., _6,
503 (1937).

Smith, G. D., Numerical Solution of Partial Differential ,
Equations: Finite Difference Method (Oxford University
Press, Oxford, 1978).

Steele, B. C. H., Current Research (1981).

Turq, P., and M. Chemla, J. Phys. Chem., 66, 485 (1977).

Turq, P., J. M'Halla, and M. Chemla, Ber. Bunsenges. Phys.
Chem., 66, 982 (1979).

Vallet, C. E., and J. Braunstein, J. Phys. Chem., 66,
2438 (1977).

vanLeeuwen, H. P., Electrochim. Acta, 66, 207 (1978).
Weir, W. D., and C. G. Enke, J. Phys. Chem., 1;, 280 (1967).

Yarnitzky, C., and F. C. Anson, J. Phys. Chem., 16, 3123
(1970).