Appendices

APPENDIX A: CONTACT
ADSORPTION OF IONS
Derivation of the triple layer model
If contact adsorption occurs, the double layer becomes a triple layer. The third
line can be drawn through the centres of the contact adsorbed ions. It is called
the IHP or Inner Helmholtz Plane. If it is assumed that there is no diffuse
double layer, the total potential difference across a triple layer interface (from
carbon (M) to bulk liquid (S)) can be written as follows:

M S = ( M IHP ) + ( IHP OHP )

A.1

Both potential drops can be expressed in terms of the integral capacities KM-IHP
= qM/(M IHP) and KIHP-OHP = qOHP/(IHP - OHP) that are independent of
potential and charge. Differentiating to qM gives:
d ( M S )
dqOHP
1
1
=
+
dq M
K M IHP K IHP OHP dq M

A.2

The electroneutrality constraint dictates that qM = qCA + qOHP. Differentiating

this constraint gives:
dqOHP
dq
= 1 CA
dq M
dq M

A.3

Combining Eq. A.2, Eq. A.3 and Eq. 2.2 yields:

1
1
1
=
+
C K M IHP K IHP OHP

dqCA
1

dq

1
1
=
+
K M IHP K IHP OHP

dqCA
1

K IHP OHP dq M

195

A.4

Appendices

This formula shows that the capacity of a system where contact adsorption of
ions take place does not depend only on the integral capacities of the Inner
and Outer Helmholtz regions, but also on the differential of contact
adsorption with electrode charge. The dependence of the capacity on qM
(which is directly related to the applied potential) can hence account for
asymmetric capacity-potential curves and especially for the non constant
values at higher potentials.
The contact adsorption of ions is examined using the following equilibrium:
Empty site + ion at the OHP Contact Adsorbed ion
The transition of an ion from the OHP to the IHP can be expressed using a
change in Gibbs free energy. Applying the law of mass action to this
equilibrium results in:
nCA
n empty sites

= n ion e Gr / RT

A.5

with n the amount per surface area and Gr the free energy change associated
with the reaction. The nion is related to the concentration of ions in the
solution and nempty sites to the fraction of the surface covered with ions. Because
we work under isothermal conditions the change in Gibbs free energy equals
the amount of work done. The total work when transporting an ion is the sum
of three contributions:
1) Chemical work arising from forces between electrode and ion.
2) Electrical work arising from interactions of the ion with the electrical field.
3) Electrical work arising from interactions of the ion with its surrounding
contact adsorbed ions.
The chemical work is accounted for by defining a change of free chemical
energy Gr,chem. The second contribution is related to the distance (x2 x1) our
ion is moved with respect to the field. The field inside the capacitor equals the
potential times the distance and the work done while transporting a charge e0
over this distance is:
196

Appendices

w OHP IHP =

q M (x 2 x 1 ) e 0
0r

A.6

Calculating the interactions of an ion with its contact adsorbed neighbours is

complicated. Bockris and Reddy (1970b) described the interactions of a
central ion with the surrounding ions by assuming the latter were pinned
down in a hexagonal grid (see Figure A.1) in order to simplify mathematics.
By smoothing out the charges of ions at the same distance into a ring of
charge o, they obtained a central ion surrounded by charged rings of increasing
radii. In the first ring 6 ions can reside, in the second 12 and in the oth ring 6o.

2r

r
central
ion
Area
(r/2)2
Site

Figure A.1: The adsorption sites for ions in the IHP according to
Bockris and Reddy (1970b).

The charge density in the oth ring is given by:

= 6oe 0

A.7

If r is the distance between the centres of two ions, the site area is (r/2)2 so
there are 4/r 2 sites per unit area. Because not all sites are filled with ions we
have:

197

Appendices

r=

4
nCA

A.8

The potential at the central ion due to the oth ring is equal to the total charge
on the ring divided by the distance to the ring and the value for the dielectric.
For the (repulsive) interaction energy we obtain:
e 0
2 e 0
=
4 0 r or 2 0 r or

A.9

This result needs to be modified due to the closeness of the carbon solid
phase. The adsorbed charges induce charges on the electrode surface. These
induced charges are mathematically represented by image charges lying a
distance 2rion inside the electrode (Bockris and Reddy, 1970b). The distance to
the reference ion is [(or)2 + (2rion)2]1/2. Thus the total interaction work taking
into account all rings both image and real becomes:
w LI =

e 0
2 0 r or
n =1

1
1

1 + (2rion / or )2

1
2

A.10

Using the following factorisation:

(1 + x ) = 1 21 x + 83 x 2 + ...
1
2

A.11

and incorporating the expressions for r (Eq. A.7) and (Eq. A.8) leads to the
final result:
1

w LI

2
3rion2 e 02 2 nCA
=
8 0 r

2
1 9rion4 e 02 2 nCA

2
32 0 r
n =1 o

o
n =1

A.12

where the two summations are equal to 2/6 and 4/90 respectively. With the
expressions derived for the chemical, electrical and interaction contributions
the total work is known. Substituting the results in Eq. A.6 gives the sought
after expression for the population density of contact adsorbed ions:

198

Appendices

Gr ,chem q M (x 2 x 1 )e 0
nCA
= n ion exp
+
n empty
0 r kT
RT

16 0 r kT 320 0 r kT
5
3
2
rion2 e 02 2 nCA

11 5
2
rion4 e 02 2 nCA

A.13

If the following additional relations are used (Bockris and Reddy, 1970b):
nCA =

N A 2ri
CA q max N A n empty
,
= 1 CA and n ion =
F
n total
1000n total

A.14

with qmax the charge density occurring at maximal adsorption, NA Avogadros

number and CA the fraction of the surface covered with ions. Incorporating
these relations in Eq. A.13, neglecting the (smaller) second term of the
interaction work and taking the logarithmic form, one obtains:
ln

Gr ,chem (x 2 x 1 )e 0
2rion F
CA
+ ln c ion
= ln
+
qM

q
RT
kT
1 CA
1000

max
0 r
2 e 2r 2 q N 2 3
0 ion
max
A
2

CA

16 0 r kT F
5

A.15

or simply:
ln

3
CA
2
= kCA + ln c ion + Aq M BCA
1 CA

A.16

with kCA, A and B the appropriate constants. Eq. A.16 can be used to calculate
the differential dqCA/dqM in Eq. A.4. Because CA and qCA are linearly
dependent the differential dCA/dqM can be used directly. However Eq. A.16
cannot be differentiated directly due to the two terms containing CA. Taking
1/(dqM /dCA) results in an equation containing only CA, yet an equation with
qCA is needed. The differential dCA/dqM must therefore be calculated
numerically.

199

Appendices

As can be seen from Figure A.2, the occurrence of contact adsorption can
lead to asymmetric potential capacity curves. It was not possible to reproduce
exactly the capacity curves calculated by Bockris. The inflection point in the
differential qCA/qM was found to occur at more negative potentials. This
results in the relatively high capacities (the left hump in Figure A.2) for
negative potentials contrary to the low values reported in literature.

C [F/cm2]

20

15

10
-0.2

-0.1

qM

0.1

[C/m2]

Figure A.2: Calculated differential capacity from the Triple layer

model (Eq. A.16) in series with the GC model (Eq. 2.X).
Chloride ions are contact adsorbing from a bulk containing 0.5 N
KCl.

200

0.2

Appendices

APPENDIX B: LISTING
GPROMS INPUT FILES

OF

Program listing 1. A single current collector is used that can be put at

the begin or the end of the packed bed. The electrical capacity is
constant.
#######################################################
##
##
##
gPROMS model source
##
##
model: Langmuir-Freundlich, ext resist, 1D
##
##
single current collector
##
##
December 2000
##
##
Vincent Fischer
##
##
##
#######################################################
DECLARE
TYPE
VOLTAGE
=
FLUX
=
CONDUCTIVITY =
CONCENTRATION =
LOADING
=
ADSCONST
=
SPEED
=
DISPERSION
=
FLOW
=
TRANSFER
=
SPECIFICSURFACE=

END

0.5 :
0
:
1.0 :
1
:
1
:
10 :
1e-6:
1e-9:
1e-6:
1e-4:
1E6 :

-1E15
-1E15
-1E15
-1E30
-1E30
-1E15
-1e15
-1e15
-1e15
-1e15
-1e15

:
:
:
:
:
:
:
:
:
:
:

1E15
1E15
1E15
1E15
1E15
1E15
1e15
1e15
1e15
1e15
1e15

UNIT
UNIT
UNIT
UNIT
UNIT
UNIT
unit
unit
unit
unit
unit

=
=
=
=
=
=
=
=
=
=
=

"V"
"mol/m2.s"
"1/Ohm.m"
"mol/m3"
"mol/kg"
"m3/mol"
"m/s"
"m2/s"
"m3/s"
"m/s"
"1/m"

# DECLARE SECTION

#####################################
##
##
## MODELS START HERE
##
##
##
#####################################
MODEL Tube
PARAMETER
Length AS
dbed
AS
dp
AS
rho
AS
Ebed
AS
qmax
AS

REAL
REAL
REAL
REAL
REAL
REAL

#
#
#
#
#
#

Reactor Length
Reactor diameter
Particle diameter
density particle
External porosity
Maximale belading

201

[m
[m
[m
[kg/m3
[m3/m3
[mol/m3

]
]
]
]
]
]

Appendices
kappaL
CAP
A
B
C
diff
SBET
lambda

AS
AS
AS
AS
AS
AS
AS
AS

REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL

#
#
#
#
#
#
#
#

Specific Conductivity
Differential Capacitance
interaction coeff. 1
interaction coeff. 2
interaction coeff. 3
diffusivity coefficient
specific surface area
Ohmic loss factor

DISTRIBUTION_DOMAIN
Axial AS (0:Length)
Axial2 AS (0:Length)

VARIABLE
q
AS
J
AS
cB
AS
V
AS
K
AS
kappa AS
cBin AS
Vzero AS
u
AS
D
AS
kL
AS
Abed AS
F
AS

[1/Ohm.m ]
[F/m2
]
[m3/mol.V2]
[m3/mol.V ]
[m3/mol
]
[m2/s
]
[m2/kg
]
[]

# coordinate for the concentration

# coordinate for the voltage

DISTRIBUTION(Axial) OF loading
# in solid phase
DISTRIBUTION(Axial) OF FLUX
# adsorption flux
DISTRIBUTION(Axial) OF CONCENTRATION # in reactor
DISTRIBUTION(Axial2) OF VOLTAGE
# Elect. potential
DISTRIBUTION(Axial2) OF ADSCONST
# Adsorption equil.
CONDUCTIVITY
# Bulk conductivity
CONCENTRATION
# Concentration in inlet
VOLTAGE
# External applied potential
speed
# interstitial velocity
dispersion
# Dispersion coefficient
transfer
# mass transfer coefficient
specificsurface
# surface area installed
flow
# volume flow

BOUNDARY
# @ x = 0
# At the inlet
cB(0) = cBin;
# closed boundary
# cB(0) = cBin + D/u * PARTIAL(cB(0),Axial) ; # open boundary
# V(0) = Vzero - lambda*Length* PARTIAL(V(0),Axial2) # for co-current
operation
PARTIAL(V(0),Axial2) = 0 ;
# for counter-current operation
# @ x = Length
# At the outlet
PARTIAL(cB(Length),Axial) = 0 ; # no flux condition
V(Length) = Vzero - lambda*Length*PARTIAL(V(Length),Axial2);
counter-current operation
# PARTIAL(V(Length),Axial2) = 0 ; # for co-current operation

EQUATION
# Specific interfacial area
Abed = 3*(1 - Ebed) / (dp/2);
# Interstitial liquid speed through column
u = F/(Ebed*3.1415*0.25*dbed^2);
# Dispersion coefficient

202

# for

Appendices
D = (2*Ebed/(3-Ebed)) * diff + (u*dp/2) ;
# Mass transfer correlation of Wilson and Geankoplis, 1966
kL = 1.09*diff/(Ebed*dp)*(dp * u /(Ebed*diff))^0.3333 ;
# Component B balance fluid phase
FOR x := 0|+ TO Length|- DO
$cB(x) = (D * PARTIAL(PARTIAL(cB(x),Axial),Axial)) (u * PARTIAL(cB(x),Axial)) - (J(x) * Abed) ;
# Diffusion, convection and adsorption parts
END # for loop
# Component B balance solid phase
FOR x := 0 TO length DO
$q(x) = J(x) * Ebed * Abed /(rho*(1 - Ebed)) ;
END # for loop
# Transport Equation to solid phase with LF isotherm
FOR x := 0 TO Length DO
J(x) = kL * (cB(x) - (q(x)/((qmax -q(x) )*K(x) ))^2 ) ;
END # for loop
# Langmuir-Freundlich Adsorption Constant from Electric Field
FOR x := 0 to Length DO
K(x) = (A * V(x)^2) + (B * V(x)) + (C) ;
END # for loop
# Effective conductivity from specific conductivity
kappa = ((2 * Ebed) / (3-Ebed)) * kappaL ;
# Electric potential differential equation
FOR x := 0|+ TO Length|- DO
$V(x) = (Ebed*kappa/(rho*SBET*(1-Ebed)*CAP)) *
PARTIAL(PARTIAL(V(x),axial2),axial2) ;
END # for loop
END # Model

#####################################
##
##
## SIMULATION STARTS HERE
##
##
##
#####################################
PROCESS ElecBuis
UNIT Elecbuis AS Tube
MONITOR
# The variables
#Elecbuis.J(50)
Elecbuis.cB(25)
Elecbuis.cB(50)
#Elecbuis.q(5)

that will be reported

;
;
;
;

203

Appendices
Elecbuis.q(25) ;
Elecbuis.q(50) ;
#Elecbuis.K(25) ;
#Elecbuis.K(50) ;
Elecbuis.V(13) ;
Elecbuis.V(25) ;
Elecbuis.V(37) ;
Elecbuis.V(50) ;
SET # Constants
WITHIN ElecBuis DO
Axial := [BFDM, 2,
Axial2 := [CFDM, 2,
qmax
:= 10
;
Length := 0.02
;
dbed
:= 0.7
;
dp
:= 5E-4
;
kappaL := 5.88
;
rho
:= 1146
;
Ebed
:= 0.45
;
CAP
:= 0.15
;
A
:= -0.1
;
B
:= 0
;
C
:= 0.2
;
diff
:= 9E-10
;
SBET
:= 1.1E6
;
lambda := 0.05
;
END # within

50]; # Discretization for convective flow

50]; # Discretization for dispersive flow
# [mol/kg
] monolayer adsorption
# [m
] of the bed
# [m
] diameter of the bed
# [m
] particle diameter
# [1/Ohm.m ] electric conductivity liquid
# [kg/m3
] particle density
# [m3/m3
] porosity of the bed
# [F/m2
] differential capacity
# [m3/mol.V2] 2th order influence pot
# [m3/mol.V ] 1th order influence pot
# [m3/mol
] 0th order influence pot
# [m2/s
] diffusion coefficient
# [m2/kg
] BET surface area
# [] external loss parameter pot

ASSIGN # more variables

WITHIN Elecbuis DO
F := 33.33E-6 ;
cBin := 5 ;
Vzero := 0.0 ;
END # within
INITIAL # Differential Variables
WITHIN ElecBuis DO
FOR x := 0|+ TO Length|- DO
cB(x) = 0;
END # for
FOR x := 0 TO Length DO
q(x) = 0;
END # for
FOR x := 0|+ TO Length|- DO
V(x) = 0;
END # for
END # WITHIN
SOLUTIONPARAMETERS
GPLOT := ON ;
REPORTINGINTERVAL := 2000 ;
SCHEDULE

204

Appendices
SEQUENCE
CONTINUE FOR 80000
RESET # Desorption starts here
Elecbuis.F := 33.33e-6 ;
Elecbuis.Vzero := 1.0 ;
Elecbuis.cBin := 0.0 ;
END # reset
CONTINUE FOR 320000
END # sequence
END # Process

Program listing 2. The packed bed is divided into n sub regions in

order to simulate the effect of n-1 current collectors in contact with the
carbon. The electrical capacity is considered constant.
#######################################################
##
##
##
gPROMS model source
##
##
model: Langmuir-Freundlich, ext resist, 1D
##
##
multiple current collectors
##
##
December 2000
##
##
Vincent Fischer
##
##
##
#######################################################
DECLARE
TYPE
VOLTAGE
FLUX
CONDUCTIVITY
CONCENTRATION
loading
ADSCONST
speed
dispersion
flow
transfer
specificsurface
END

=
=
=
=
=
=
=
=
=
=
=

0.5 :
0
:
1.0 :
1
:
1
:
10 :
1e-6:
1e-9:
1e-6:
1e-4:
1E6 :

-1E15
-1E15
-1E15
-1E30
-1E30
-1E15
-1e15
-1e15
-1e15
-1e15
-1e15

:
:
:
:
:
:
:
:
:
:
:

1E15
1E15
1E15
1E15
1E15
1E15
1e15
1e15
1e15
1e15
1e15

# DECLARE SECTION

#####################################
##
##
## SUBMODEL STARTS HERE
##
##
##
#####################################

205

UNIT
UNIT
UNIT
UNIT
UNIT
UNIT
unit
unit
unit
unit
unit

=
=
=
=
=
=
=
=
=
=
=

"V"
"mol/m2.s"
"1/Ohm.m"
"mol/m3"
"mol/kg"
"m3/mol"
"m/s"
"m2/s"
"m3/s"
"m/s"
"1/m"

Appendices

MODEL partofbed
PARAMETER
Length
dbed
dp
rho
Ebed
qmax
kappaL
CAP
A
B
C
diff
SBET
lambda

AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS

REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL

#
#
#
#
#
#
#
#
#
#
#
#
#
#

Reactor Length
Reactor diameter
Particle diameter
density particle
External porosity
Maximale belading
Specific Conductivity
Differential Capacitance
interaction coeff. 1
interaction coeff. 2
interaction coeff. 3
diffusivity coefficient
specific surface area
Ohmic loss factor

DISTRIBUTION_DOMAIN
Axial AS (0:Length)
Axial2 AS (0:Length)
VARIABLE
q
J
cB
V
K
kappa
cBinlet
Vzero
u
D
kL
Abed
F

AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS

[m
]
[m
]
[m
]
[kg/m3
]
[m3/m3
]
[mol/m3
]
[1/Ohm.m ]
[F/m2
]
[m3/mol.V2]
[m3/mol.V ]
[m3/mol
]
[m2/s
]
[m2/kg
]
[]

# coordinate for the concentration

# coordinate for the voltage

DISTRIBUTION(Axial) OF loading
# in solid phase
DISTRIBUTION(Axial) OF FLUX
# adsorption flux
DISTRIBUTION(Axial) OF CONCENTRATION # in reactor
DISTRIBUTION(Axial2) OF VOLTAGE # Elect potential
DISTRIBUTION(Axial2) OF ADSCONST
# ads. constant
CONDUCTIVITY
# Effective conductivity
CONCENTRATION
# Concentration in inlet
VOLTAGE
# ext applied potential
speed
# interstitial velocity
dispersion
# Dispersion coefficient
transfer
# mass transfer coeff
specificsurface
# surface area installed
flow
# volume flow

BOUNDARY
# @ x = 0
# At the inlet
cB(0) = cBinlet;
# closed boundary
# cB(0) = cBinlet + D/u * PARTIAL(cB(0),Axial) ; # open boundary
# V(0) = Vzero + lambda*Length* PARTIAL(V(0),Axial2) ; # cocurrent
operation
PARTIAL(V(0),Axial2) = 0 ;
# countercurrent operation
# @ x = Length

# At the outlet

206

Appendices
PARTIAL(cB(Length),Axial) = 0 ;

# no flux condition

V(Length) = Vzero - lambda*Length*PARTIAL(V(Length),Axial2);

# countercurrent
# PARTIAL(V(Length),Axial2) = 0 ; # cocurrent
EQUATION
# Specific interfacial area
Abed = 3*(1 - Ebed) / (dp/2);
# Interstitial liquid speed through column
u = F/(Ebed*3.1415*0.25*dbed^2);
# Dispersion coefficient
D = (2*Ebed/(3-Ebed)) * diff + (u*dp/2) ;
# Mass transfer correlation of Wilson and Geankoplis, 1966
kL = 1.09*diff/(Ebed*dp)*(dp * u /(Ebed*diff))^0.3333 ;
# Component B balance fluid phase
FOR x := 0|+ TO Length|- DO
$cB(x) = (D * PARTIAL(PARTIAL(cB(x),Axial),Axial)) (u * PARTIAL(cB(x),Axial)) - (J(x) * Abed) ;
# Diffusion convection and adsorption parts
END # for loop
# Component B balance solid phase
FOR x := 0 TO length DO
$q(x) = J(x) * Ebed * Abed /(rho*(1 - Ebed)) ;
END # for loop
# Transport Equation to solid phase with LF isotherm
FOR x := 0 TO Length DO
J(x) = kL * (cB(x) - (q(x)/((qmax -q(x) )*K(x) ))^2 ) ;
END # for loop
# Langmuir-Freundlich Adsorption Constant from Electric Field
FOR x := 0 to Length DO
K(x) = (A * V(x)^2) + (B * V(x)) + (C) ;
END # for loop
# Effective conductivity from specific conductivity
kappa = ((2 * Ebed) / (3-Ebed)) * kappaL ;
# Electric potential differential equation
FOR x := 0|+ TO Length|- DO
$V(x) = (Ebed*kappa/(rho*SBET*(1-Ebed)*CAP)) *
PARTIAL(PARTIAL(V(x),axial2),axial2) ;
END # for loop

207

Appendices
END # Model
#####################################
##
##
## MAIN MODEL STARTS HERE
##
##
##
#####################################
MODEL totalbed
PARAMETER
nbeds AS INTEGER : DEFAULT = 4 # number of submodels
UNIT
FOR i := 1 TO nbeds DO
Bed(i) AS partofbed
END # for
EQUATION
FOR i:= 1 TO (nbeds 1) DO
Bed(i).cB(Bed(i).Length) = Bed(i+1).cBinlet ;
Bed(i).F = Bed(i+1).F ;
Bed(i).Vzero = Bed(i+1).Vzero ;
END # for
END # totalbed
#####################################
##
##
## SIMULATION STARTS HERE
##
##
##
#####################################
PROCESS ElecBuis
UNIT Elecbuis AS totalbed
MONITOR
# The variables that will be reported
Elecbuis.B1.cB(2) ;
Elecbuis.B1.cB(26) ;
Elecbuis.B1.cB(50) ;
Elecbuis.B2.cB(26) ;
Elecbuis.B2.cB(50) ;
Elecbuis.B3.cB(26) ;
Elecbuis.B3.cB(50) ;
Elecbuis.B4.cB(26) ;
Elecbuis.B4.cB(50) ;

208

Appendices
SET # Constants
WITHIN ElecBuis DO
FOR i:= 1 TO nbeds DO
Bed(i).Axial := [BFDM, 2, 50]; # Discretization method for mass
Bed(i).Axial2 := [CFDM, 2, 50]; # Discretization method for charge
Bed(i).qmax
:= 10
;
# [mol/kg
] monolayer adsorption
Bed(i).Length := 0.02
;
# [m
] of the bed
Bed(i).dbed
:= 0.7
;
# [m
] diameter of the bed
Bed(i).dp
:= 5E-4
;
# [m
] particle diameter
Bed(i).kappaL := 5.88
;
# [1/Ohm.m ] electric cond.liquid
Bed(i).rho
:= 1146
;
# [kg/m3
] particle density
Bed(i).Ebed
:= 0.45
;
# [m3/m3
] porosity of the bed
Bed(i).CAP
:= 0.15
;
# [F/m2
] diff capacity system
Bed(i).A
:= -0.1
;
# [m3/mol.V2] 2th order infl pot
Bed(i).B
:= 0
;
# [m3/mol.V ] 1th order infl pot
Bed(i).C
:= 0.2
;
# [m3/mol
] 0th order infl pot
Bed(i).diff
:= 9E-10
;
# [m2/s
] diffusion coeff
Bed(i).SBET
:= 1.1E6
;
# [m2/kg
] BET surface area
Bed(i).lambda := 0.666 ;
# [] ext loss parameter
END # for
END # within
ASSIGN # more variables
WITHIN Elecbuis DO
Bed(1).F := 33.33E-5 ;
Bed(1).cBinlet := 0.5 ;
Bed(1).Vzero := 0.0 ;
END # within
INITIAL # Differential Variables
WITHIN ElecBuis DO
FOR i := 1 TO nbeds DO
FOR x := 0|+ TO Bed(i).Length|- DO
Bed(i).cB(x)= 0 ;
Bed(i).V(x) = 0 ;
END # for
FOR x := 0 TO Bed(i).Length DO
Bed(i).q(x) = 0;
END # for
END # for
END # WITHIN
SOLUTIONPARAMETERS
GPLOT := ON ;
REPORTINGINTERVAL := 2000 ;
SCHEDULE
SEQUENCE
CONTINUE FOR 80000

209

Appendices
RESET # Desorption starts here
Elecbuis.B1.F := 33.33e-5 ;
Elecbuis.B1.Vzero := 1.0 ;
Elecbuis.B1.cBinlet := 0.0 ;
END # reset
CONTINUE FOR 320000
END # sequence
END # Process

Program listing 3. A single current collector is used in combination

with a variable differential capacity and a non-empirical influence from
the capacity and the potential on the adsorption equilibrium.
#######################################################
##
##
##
gPROMS model source
##
##
model: Langmuir-Freundlich, ext resist, 1D
##
##
single current collector, variable capacity
##
##
October 2000
##
##
Vincent Fischer
##
##
##
#######################################################
DECLARE
TYPE
VOLTAGE
FLUX
CONDUCTIVITY
CONCENTRATION
loading
ADSCONST
speed
dispersion
flow
transfer
specificsurface
END

=
=
=
=
=
=
=
=
=
=
=

0.5 :
0
:
1.0 :
1
:
1
:
10 :
1e-6:
1e-9:
1e-6:
1e-4:
1E6 :

-1E15
-1E15
-1E15
-1E30
-1E30
-1E15
-1e15
-1e15
-1e15
-1e15
-1e15

:
:
:
:
:
:
:
:
:
:
:

1E15
1E15
1E15
1E15
1E15
1E15
1e15
1e15
1e15
1e15
1e15

# DECLARE SECTION

#####################################
##
##
## MODELS START HERE
##
##
##
#####################################
MODEL Tube
PARAMETER

210

UNIT
UNIT
UNIT
UNIT
UNIT
UNIT
unit
unit
unit
unit
unit

=
=
=
=
=
=
=
=
=
=
=

"V"
"mol/m2.s"
"1/Ohm.m"
"mol/m3"
"mol/kg"
"m3/mol"
"m/s"
"m2/s"
"m3/s"
"m/s"
"1/m"

Appendices
Length
dbed
dp
rho
Ebed
qmax
kappaL
K0
diff
SBET
lambda

AS
AS
AS
AS
AS
AS
AS
AS
AS
AS
AS

REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL
REAL

#
#
#
#
#
#
#
#
#
#
#

Reactor Length
Reactor diameter
Particle diameter
density particle
External porosity
Maximale belading
Specific Conductivity
open circ ads coeff. 3
diffusivity coefficient
specific surface area
Ohmic loss factor

DISTRIBUTION_DOMAIN
Axial AS (0:Length)
Axial2 AS (0:Length)

[m
[m
[m
[kg/m3
[m3/m3
[mol/m3
[1/Ohm.m
[m3/mol
[m2/s
[m2/kg
[-

]
]
]
]
]
]
]
]
]
]
]

# coordinate for the concentration

# coordinate for the voltage

VARIABLE
q
AS DISTRIBUTION(Axial) OF loading
# in solid phase
J
AS DISTRIBUTION(Axial) OF FLUX
# adsorption flux
cB
AS DISTRIBUTION(Axial) OF CONCENTRATION # in reactor
V
AS DISTRIBUTION(Axial2) OF VOLTAGE # Electric pot
K
AS DISTRIBUTION(Axial2) OF ADSCONST # adsorption const
CAP, epsRA, epsRB, CGC, CHA, CHB
AS DISTRIBUTION(Axial2) OF VOLTAGE # diff capacity
kappa
AS CONDUCTIVITY
# Effective conductivity
cBinlet AS CONCENTRATION
# Concentration in inlet
Vzero
AS VOLTAGE
# ext applied potential
u
AS speed
# interstitial velocity
D
AS dispersion
# Dispersion coefficient
kL
AS transfer
# mass transfer coeff
Abed
AS specificsurface
# surface area installed
F
AS flow
# volume flow
BOUNDARY
# @ x = 0
# At the inlet
cB(0) = cBinlet;
# closed boundary
# cB(0) = cBinlet + D/u * PARTIAL(cB(0),Axial) ; # open boundary
# V(0) = Vzero + lambda*Length* PARTIAL(V(0),Axial2) ; # cocurrent
operation
PARTIAL(V(0),Axial2) = 0 ;
# countercurrent operation
# @ x = Length
PARTIAL(cB(Length),Axial) = 0 ;

# At the outlet
# no flux condition

V(Length) = Vzero - lambda*Length*PARTIAL(V(Length),Axial2);

countercurrent
# PARTIAL(V(Length),Axial2) = 0 ;
# cocurrent
EQUATION
# Specific interfacial area
Abed = 3*(1 - Ebed) / (dp/2);
# Interstitial liquid speed through column
u = F/(Ebed*3.1415*0.25*dbed^2);

211

Appendices
# Dispersion coefficient
D = (2*Ebed/(3-Ebed)) * diff + (u*dp/2) ;
# Mass transfer correlation of Wilson and Geankoplis, 1966
kL = 1.09*diff/(Ebed*dp)*(dp * u /(Ebed*diff))^0.3333 ;
# Component B balance fluid phase
FOR x := 0|+ TO Length|- DO
$cB(x) = (D * PARTIAL(PARTIAL(cB(x),Axial),Axial)) (u * PARTIAL(cB(x),Axial)) - (J(x) * Abed) ;
# Diffusion convection and adsorption parts
END # for loop
# Component B balance solid phase
FOR x := 0 TO length DO
$q(x) = J(x) * Ebed * Abed /(rho*(1 - Ebed)) ;
END # for loop
# Transport Equation to solid phase with LF isotherm
FOR x := 0 TO Length DO
J(x) = kL * (cB(x) - (q(x)/((qmax -q(x) )*K(x) ))^2 ) ;
END # for loop
# Langmuir-Freundlich Adsorption Constant from Electric Field
FOR x := 0 to Length DO
K(x) = 0.1+0.5*K0*EXP(- ( (CHA(x)*CGC(x)/(CHA(x)+CGC(x)) )
- (CHB(x)*CGC(x)/(CHB(x)+CGC(x) )) )*V(x)^2*34.092028) ;
END # for loop
# Effective conductivity from specific conductivity
kappa = ((2 * Ebed) / (3-Ebed)) * kappaL ;
# Electric potential differential equation
FOR x := 0|+ TO Length|- DO
$V(x) = (Ebed*kappa/(rho*SBET*(1-Ebed)* CAP(x) )) *
PARTIAL(PARTIAL(V(x),axial2),axial2) ;
END # for loop
# dielectric constants as function of V
# fit function taken to describe the real function
FOR x := 0 TO Length DO
epsRA(x) = 1/(0.036245262 - 0.00031153174*V(x) +
0.023023228*V(x)^2 - 0.0063672479*V(x)^3 );
END # for loop
FOR x := 0 TO Length DO
epsRB(x) = 1/(0.24572379 - 0.0010757573*V(x) +
0.041396901*V(x)^2 - 0.011482038*V(x)^3 );
END # for loop
# Gouy Chapman capacity
FOR x := 0 TO Length DO
CGC(x) = (0.03327*epsRA(x))^0.5*COSH(19.47203*V(x));
END # for loop
# Helmholtz capacities

212

Appendices
FOR x := 0 TO Length DO
CHA(x) = epsRA(x)/80.49156;
END #for loop
FOR x := 0 TO Length DO
CHB(x) = 1/(59.81714/epsRB(x) + 45.40607/epsRA(x));
END # for loop
# total capacity (potential and coverage dependent)
FOR x:= 0 TO Length DO
CAP(x) = 0.5 *( (q(x)/qmax)/(1/CHB(x) + 1/CGC(x))
+ (1 - q(x)/qmax)/(1/CHA(x) + 1/CGC(x)) );
END # for loop
END # Model
#####################################
##
##
## SIMULATION STARTS HERE
##
##
##
#####################################
PROCESS ElecBuis
UNIT Elecbuis AS Tube
MONITOR
# The variables that will be reported
#Elecbuis.V(1) ;
#Elecbuis.V(4) ;
Elecbuis.K(1) ;
Elecbuis.cB(50) ;
Elecbuis.q(50) ;
Elecbuis.V(1) ;
Elecbuis.V(17) ;
Elecbuis.V(33) ;
Elecbuis.V(50) ;
Elecbuis.CAP(25) ;

SET # Constants
WITHIN ElecBuis DO
Axial := [BFDM, 2, 50]; # Discretization method for mass
Axial2 := [CFDM, 2, 50]; # Discretization method for charge
qmax
:= 10
;
# [mol/kg
] monolayer adsorption
Length := 0.02
;
# [m
] of the bed
dbed
:= 0.7
;
# [m
] diameter of the bed
dp
:= 5E-4
;
# [m
] particle diameter
kappaL := 5.88
;
# [1/Ohm.m ] electric conduct liquid
rho
:= 1146
;
# [kg/m3
] particle density
Ebed
:= 0.45
;
# [m3/m3
] porosity of the bed
K0
:= 0.2
;
# [m3/mol
] open circ equil const
diff
:= 9E-10
;
# [m2/s
] diffusion coefficient
SBET
:= 1.1E6
;
# [m2/kg
] BET surface area
lambda := 0.666 ;
# [] external loss parameter
END # within

213

Appendices

ASSIGN # more variables

WITHIN Elecbuis DO
F := 33.33E-6 ;
cBinlet := 5 ;
Vzero := 0.0 ;
END # within
INITIAL # Differential Variables
WITHIN ElecBuis DO
FOR x := 0|+ TO Length|- DO
cB(x) = 0;
END # for
FOR x := 0 TO Length DO
q(x) = 0;
END # for
FOR x := 0|+ TO Length|- DO
V(x) = 0;
END # for
END # WITHIN
SOLUTIONPARAMETERS
GPLOT := ON ;
REPORTINGINTERVAL := 2000 ;
SCHEDULE
SEQUENCE
CONTINUE FOR 80000
RESET # Desorption starts here
Elecbuis.F := 33.33e-6 ;
Elecbuis.Vzero := 1.0 ;
Elecbuis.cBinlet := 0.0 ;
END # reset
CONTINUE FOR 320000
END # sequence
END # Process

214

Appendices

APPENDIX C: ECONOMICS
General data
Packed bed data

Bed diameter
Bed length
Area of the bed
Number of layers
Volume of the bed
Density of the carbon
Porosity of the bed
Size of the surface
Electrical capacity
Amount of carbon

D
L
A = D2
x
V=ALx
p

SBET
C
m = V p (1-)

0.7 m
0.02 m
0.385 m2
10500
0.1543.848 m3
1146 kg/m3
0.45
1100 m2/g
0.15 F/m2
972426 kg

A
meOH
cp,A
cp,steam
cp,meOH
Hg
Hg,meOH
Egas
Kgas
KmeOH

1000 kg/m3
791 kg/m3
4180 J/(kg K)
2000 J/(kg K)
2500 J/(kg K)
2.26 106 J/kg
1.1 106 J/kg
31.65 106 J/m3
0.159 /m3
0.159 k/m3

Physical properties

Density water
Density methanol
Heat capacity water
Heat capacity steam
Heat capacity methanol
Heat of vaporisation
For methanol
Energy of heating gas
Price of heating gas
Price of methanol
Waste water data

Flow
Pollutant in
Pollutant out
cout
Seconds in one year
Amount removed

c
t
M =F t (c-cout)
215

20 L/min
5 mol/m3
0.05 mol/m3
3.15 107 s
5.256 104 mol

Appendices

Isotherm data

LF constant 1
LF constant 2
LF constant 3

n
K
qmax

0.5
0.2 (m3/mol)
10 mol/kg

Max loading

Kc n
q = q max
1 + Kc n

3.09 mol/kg

# regeneration cycles

beds = M /(m q)

34017

Case 0) No regeneration
Investment

44.94 + 1.365(Lx )1.31

Column
=
2.2 10 year
(WEBSCI and WUBO, 1999)

2.075.19 k

Operation costs

Waste carbon
(Staal, 2001)
New carbon
(Staal, 2001)

= 0.455 k/tonne beds m

38.27 k

= 3.182 k/tonne beds m

267.9 k

Case 1) Steam regeneration (300 C, 90% efficiency)

Investment

2 44.94 + 1.365(Lx )1.31

Column
=
2.2 10 year
(WEBSCI and WUBO, 1999)
hr 1.5m
150
tonne 8hr
Steam house
=
2.2 10 year
(WEBSCI and WUBO, 1999)
Initial bed
= 3.182 k/tonne m

216

4.1410.37 k

0.31612.66 k
0.78831.502 k

Appendices

Operation costs

Labour costs
= 0.02742beds
(27 /hr, 4 hr, 2 men per regeneration)
Waste carbon
(0.455 k/tonne beds m)/10
New carbon
(3.182 k/tonne (beds-1)m)/10
Steam needed
msteam = 1.5 m beds
(Suzuki, 1990)
Q1 = msteam cp,A 80 K
Q2 = msteam Hg
Q3 = msteam cp,steam 200 K
Q1 + Q 2 + Q 3
80% efficiency heater
Q tot =
0.80
Q tot K gas
heating costs
=
E gas

74.21.855 k
3.827 k
26.7123.64 k
126.3 tonne

2.376 k

Case 2) Methanol regeneration (66% efficiency)

Investment

2 44.94 + 1.365(Lx )1.31

=
Column
2.2 10 year
(WEBSCI and WUBO, 1999)
33.28 + 1.662(20V )0.7
Ktank =
Storage tank
2.2 10 year
(WEBSCI and WUBO, 1999)
Separation section
= 2Ktank
Initial bed
= 3.182 k/tonne m

4.1410.37 k

1.8325.748 k
3.66511.5 k
0.78831.502 k

Operation costs

Labour costs
0.02742beds
74.21.855 k
(27 /hr, 4 hr, 2 men per regeneration)
Waste carbon
= (0.455 k/tonne beds m)/3 12.76 k
New carbon
=(3.182 k/tonne (beds-1)m)/3 89.078.8 k
Volume methanol
VmeOH = 20 V beds
2672 m3
Methanol costs
= KmeOH 0.1VmeOH
42.5 k
Energy needed
Q1 = VmeOH meOH cp,meOH 45 K
217

Appendices

80% efficiency heater

Heating costs

Q2 = VmeOH meOH Hg,meOH

Q1 + Q 2
Q tot =
0.80
Q tot K gas
=
E gas

16.1 k

Case 3) Electrosorption (50% efficiency)

Investment

5 69.88 + 4.117(Lx )0.289

=
Column
2.2 10 year
(WEBSCI and WUBO, 1999)
Electric circuit
= 100/10
Initial bed
= 3.182 k/tonne m

16.6518.11 k
10 k
0.78831.502 k

Operation costs

Labour costs
0.02742beds
74.21.855 k
(27 /hr, 4 hr, 2 men per regeneration)
Waste carbon
= (0.455 k/tonne beds m)/2 19.14 k
New carbon
=(3.182 k/tonne (beds-1)m)/2 133.6118.2 k
Charging the bed
U1 = (1 V)2 m SBET C
Desorbing phenol
U2 = (1 V)2 m SBET (0.238-0.055)
0.1 10 3 beds
= U1
0.088 k
Charging costs
2.2kWhr 3600
0.1 10 3 beds
=U2
0.053 k
Desorption costs
2.2kWhr 3600

218