in.

{a “.53."...

.3“ "map F
4 up...

33“. .1 J. ,
.3 , . r .;
1.1.11 11...?
Ii :3 :
. x. .1.
11
3h“
”8.35.10.
.Tilv'

m. fix:
0.. Ibis! t!-
lxn nun... ‘oixll
an: 2- i?! it."
1.15. L!
.39.
It

( u.
.5
5mm?
x... ,.
.mm.

1.7.9. ‘\
.ll. .2
!NHI.: It

53...”...

. 77‘ . ‘3.
Au...|:;:|\ 3
.v;
..
Ir‘ 7:».

A I,\vr (L
rt .‘

‘i.

:9. . 2. t
11.11:“;

I 1.32... .
L a 1 :1. 5i .1.
: .3Hr......|:i:ly...'l\

magi.
, {awhékafifigfigéyfik é

. 4. «a...
g... .6. V

43...;5? A

 

THESIS

(V

                                                                                                                           

Adsorption:

1293 01410 7431

This is to certify that the

dissertation entitled

A Study Adapting Cubic Equations of State

h

presented by C * '

Ramkumar Subramanian

has been accepted towards fulfillment

Ph.D.

of the requirements for

degree in Chemical Engineering

 

Date /0//0/?5-

 

Majorvprofessor

 

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

‘—

 

 

LIBRARY
Michigan State
University

 

 

 

PLACE II RETURN BOX to romovo this chockout from your rooord.

TO AVOID FINES rotum on or botoro dot. duo.

 

DATE DUE DATE DUE DATE DUE

4%}? ll:—

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

MSU to An Atfirmottvo MbnlEqunl Opporumtty Intuition
mam

 

ABSORPTION: A STUDY ADAPTING CUBIC EQUATIONS OF STATE

Ramkumar Subramanian

A DISSERTATION

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

DOCTOR OF PHILOSOPHY

Department of Chemical Engineering

1995

ABSTRACT

ABSORPTION: A STUDY ADAPTING CUBIC EQUATIONS OF STATE
By

Ramkumar Subramanian

The Simplified Local Density (SLD) method is a new engineering approach to
model adsorption based on spatial invariance of the chemical potential along with an
equation of state. This work extends previous pure fluid applications to binary mixtures
and more complex adsorbent geometries (slits, pores). Model predictions are compared
to molecular simulations and experimental data. The SLD model can represent
adsorption isotherms of Types I-V at subcritical conditions. Supercritical behavior has
not been classified, however the SLD approach can accurately represent the complex
behavior exhibited at supercritical conditions. In addition, clustering (molecular
charisma) in supercritical fluids is modeled by representing infmite dilution by
calculating the fluid density in the region around a single solute molecule.

Simple engineering models such as Langmuir and Freundlich cannot represent the
variety of experimental adsorption isotherm shapes. On the other hand, molecular
simulations represent the behavior, but are not suitable for routine process design. This
SLD approach bridges the gap by retaining both the essential physics of the adsorption
problem and the simplicity of an equation of state. The SLD approach is demonstrated to
be a powerful engineering tool for prediction of adsorption, and requires fewer

parameters than previous engineering models. In addition, the SLD approach can be

extended to non-ideal gas phase regions without additional parameters. The single
parameter used for each adsorbate is temperature- and pressure-independent, permitting

extrapolations and predictions of adsorption behavior.

Dedicated to my parents Krishnambal Subramanian & K. R. Subramanian, to my sisters
Shyamala Nageswaran & Rama Ganesh, and to my late Kadayam grand father

iv

Acknowledgments

First of all, I would like to thank my advisor, Dr. Carl Lira, for his guidance
throughout my research, and also for his patience during tough times. I would like to
thank Drs. D. J. Miller, R. B. Buxbaum and P. M. Duxbury for serving on my committee
and for making appropriate suggestions to guide my research. I would like to thank Dr.
Bharath Rangarajan for all the discussions that we have had during this research, and for
his valuable suggestions.

I would like to thank Arvind and Shalini Mathur for helping me during times of
need and for their wonderful company throughout my stay here in East Lansing. I would
like to thank my wonderful roommates S. Iyer, and R. Easwar for being the best
companions anyone could ask for. I would like to thank Sanjay Yedur for all the
wonderful times we spent at coffee shops discussing worldly matters! I would like to
thank Mike Bly, Sanjay Padaki, Valerie Adegbite, Satyadev Chilukuri, Kartik Pashupati,
Nishant Rao and other unnamed friends for all their support and friendship during my
stay at Michigan State. I would like to thank the secretaries, Faith, Julie, Beth and Candy
for being so wonderful. I would like to thank my parents Mr. and Mrs. K. R.
Subramanian, and my sisters, Shyamala and Rama for supporting me throughout my
study. Finally, I would like to thank my wife Jayashree for her love and understanding
for the last few months.

Michigan State is a great place to live and to study, and for the last six years I
have had the greatest time of my life. Thank you Michigan State, and Go Spartans!

Table of Contents

LIST OF TABLES

LIST OF FIGURES
NOMENCLATURE (SLD MODEL)
CHAPTER 1: INTRODUCTION

CHAPTER 2: BACKGROUND, LITERATURE REVIEW AND
PRELIMINARY RESULTS

DEFINITION
CLASSIFICATION OF ABSORPTION ISOTHERMS

ADSORPTION MODELS
Henry’s Law
Langmuir Isotherm
BET Isotherm
Toth & UNILAN Equations
Two Dimensional Equations-of State
Potential Theories
Polyani Potential Theory
Dubinin’s Theory
Frenkel-Halsey-Hill Theory
van der Waals Model of Barrer and Robins
Statistical Mechanical Models
Inteng Equation Theory
Density Functional Theory
Molecular Simulations
Mixed-Gas Isotherms
Ideal Adsorbed Solution (IAS) Theory

PRELIMINARY RESULTS - SIMPLIFIED LOCAL DENSITY MODEL OF
RANGARAJAN, LIRA AND SUBRAMANIAN
Model Development
Outline of Algorithm to Solve for Density Profile
Results of SLD Model
Ideal Gas - Hard Wall
Ideal Gas - Attractive Wall
Attractive Fluid - Hard Wall
Attractive Fluid - Attractive Wall
SLD Model - Discussion

vi

viii

ix

xi

35
35
41

45
45
45
46
55

Table of Contents (contd.)

CHAPTER 3: QUANTITATIVE MODELING OF ADSORPTION AND

EXTENSION OF THE SLD APPROACH TO PREDICT CLUSTERING 58
INTRODUCTION 58
MODEL DEVELOPMENT 60

Algorithm for Solving Density Profile 66
RESULTS 67
DISCUSSION 84

CHAPTER 4: ADSORPTION OF PURE GASES IN SLITS AND PORES,

AND ADSORPTION OF BINARY MIXTURES 87
INTRODUCTION 87
MODEL DEVELOPMENT - PURE GAS 89
Slit-like Pores 89
Cylindrical Pores 92
MODEL DEVELOPMENT - BINARY MIXTURES 96
Algorithm to Solve Density Profile 97
RESULTS 98
DISCUSSION 1 l 0
CHAPTER 5: CONCLUSIONS 1 19

APPENDICES 122

REFERENCES 196

vii

Table 2-1:
Table 3-1:

List of Tables

Molecular Properties Used in the SLD Model
Molecular Properties Used in the SLD Model

Tables in Appendix 3 - Properties of Gases

Table 3-1
Table 3-2
Table 3-3
Table 3-4
Table 3-5
Table 3-6
Table 3-7

Supercritical Ethylene
Supercritical Methane
C02

Argon

Propane

Methane

Ethylene

viii

47
68

187
189
190
191
192
193
194

Figure 2-1:
Figure 2-2:
Figure 2-3:
Figure 2-4:
Figure 2-5:
Figure 2-6:
Figure 2-7:

Figure 2-8:

Figure 2-9:

Figure 2-10:
Figure 2-11:
Figure 2-12:
Figure 2-13:
Figure 2-14:

Figure 3-1:
Figure 3-2:
Figure 3-3:
Figure 3-4:
Figure 3-5:
Figure 3-6:

Figure 3-7:

List of Figures

Adsorption of Ethylene on Graphon -Data of Findenegg (1983)
Adsorption of Propane on Graphon - Data of Findenegg (1983)
Adsorption of Krypton on Graphon - Data of Findenegg (1983)
Variety of Behaviors at Subcritical Conditions

Comparison of Molecular Simulation (van Megen and Snook, 1982)
to Experimental Data (Specovius and Findenegg, 1978)

Grand Canonical Ensemble Monte Carlo Simulation of Ethylene
in a Slit (T/Tc = 0.85)

Grand Canonical Ensemble Monte Carlo Simulation of Ethylene
in a Slit (T/T, = 1.2)

Density Functional Theory - Comparison to Grand Canonical
Ensemble Monte Carlo Simulation

F ugacity Vs Pressure for a Pure Fluid

Adsorption of Ethylene on Graphon

Adsorption of Ethylene on Graphon

Density Profile of Ethylene on Graphon at 283.15 K

Weak Adsorption of Ethylene on Graphon

Adsorption of Krypton on Graphon

Adsorption of Ethylene on Graphon

Adsorption of Krypton on Graphon

Adsorption of Propane on Graphon

Adsorption of Argon on Graphon

Adsorption of Methane on Graphon

Density and Excess Number of Solvent Molecules within a
Sphere of Radius R around a Solute

Number of Excess Solvent Molecules within a Sphere of
Radius L around a Solute

10

28

3O

31

32
43
49
50
51
53
54
69
71
72
73
74

76

77

Figure 3-8:
Figure 3-9:
Figure 3-10:

Figure 3-11:

Figure 3-12:

Figure 4-1:
Figure 4-2:
Figure 4-3:

Figure 4-4:
Figure 4-5:
Figure 4-6:
Figure 4-7:
Figure 4-8:
Figure 4-9:
Figure 4-10:
Figure 4-11:
Figure 4-12:

Figure 4-13

Figure 4-14
Figure 4-15:

Clustering of CO2 around Naphthalene

Clustering of CO2 around Naphthalene

Comparison of Cluster Size with Fluorescence Spectroscopy -
CO2 around Naphthalene

Comparison of Cluster Size with Fluorescence Spectroscopy -
CO2 around Pyrene

Comparison of Cluster Size with Fluorescence Spectroscopy -
Ethylene around Pyrene

Slit-Like Pores

Cylindrical Pores

Grand Canonical Ensemble Monte Carlo Simulation of Ethylene
in a Slit (T/T, = 0.85)

Ethylene on Carbon Slits at T/Tc of 0.85

Ethylene on BPL Carbon at 212.7 K

Ethylene on BPL Carbon at 260.2 K

Ethylene on BPL Carbon at 301.4 K

Methane on Activated Carbon at 243.15 K

Adsorption of Ethylene-Methane Mixture on BPL Carbon at
212.7 K and Initial Bulk Ethylene Concentration of 0.74
Adsorption of Ethylene-Methane Mixture on BPL Carbon at
212.7 K and Initial Bulk Ethylene Concentration of 0.74
Adsorption of Ethylene-Methane Mixture on BPL Carbon at
260.2 K and Initial Bulk Ethylene Concentration of 0.235
Adsorption of Propylene-Ethylene Mixture on BPL Carbon at
293.15 K and 1.013 bar

a(r)/ab for a Pore that can hold 2.5 Ethylene Molecules in the
Horizontal Plane

Adsorption of Homogeneous Ethylene on Graphon
Adsorption of Homogeneous Ethylene on BPL Carbon at
212.7 K

Figure Appendix 3-1

Error in the Peng-Robinson Calculation of Saturation Density

78
79

80

81

82

90

93

99

100

102

103

104

106

107

108

109

111

113

' 115

116

195

Nomenclature (SLD Model)

a, b Constants of the van der Waals or Peng-Robinson equation
f fugacity
g(r) radial distribution function
N A Avagadro’s number
pressure

T temperature

2 distance from surface of the wall. Intermolecular distance is z+ofS/2
Z compressibility factor

k Boltzmann’s constant

v molar volume

nex number of molecules

A s Surface Area

 

Greek

1 F“ surface excess
4) fluid-fluid pair potential
[.1 chemical potential
‘1’ fluid-wall potential (physical adsorption), fluid (solvent)-solute potential

(clustering)

p molar density
a molecular diameter

a pair interaction parameter

xi

Subscripts
at! attractive property

rep repulsive property

jf contribution due to fluid-fluid interaction
f fluid property
3 solid or solute property (in clustering)

fir contribution due to fluid-solid interaction (physical adsorption) or fluid
(solvent)-solute interaction (clustering)
b bulk property (uniform fluid)

0 standard state

xii

CHAPTER 1: INTRODUCTION

Physical adsorption is a separation process in which certain components of a fluid
phase are transferred on to the surface of a solid adsorbent [McCabe, Smith & Harriott,
1984]. It is the result of intermolecular forces of attraction between the fluid and solid
molecules. If the attractive force of the solid on the fluid is greater than the intermolecular
forces among the fluid themselves, the fluid will condense upon the surface of the
adsorbent [T reybal, 1983],and may be accompanied by the evolution of heat. In
engineering applications, adsorption is usually used an alternate to extraction. Adsorption
tends to have a smaller capacity but a higher selectivity than extraction. It can also be
more gentle for e.g. proteins may be denatured by organic extraction but may be adsorbed.
Separation processes using adsorption are on the increase while those using extraction are
on the decrease [Belter, Cussler & Wu, 1988].

Physical adsorption is used in gas purification processes such as the removal of
volatile organic compounds from stack gases, as a means of fractionating fluids that are
difficult to separate by other methods, and in adsorbent regenerations using supercritical
fluids [Tan and Lion, 1988, 1990]. Physical adsorption is also of interest in transportation
and storage of fuel and radioactive gases, separation and purification of lower
hydrocarbons, solid-phase extractions fluids, in supercritical extractions and
chromatography [Findenegg, 1983; Barton, 1983; Strubinger and Parcher, 1989], and in

critical point drying [Rangarajan and Lira, 1992]. Understanding the thermodynamics and

structure of the gas- solid interface is essential to the understanding of heterogeneous
catalysis and wetting phenomena [Findenegg, 1983]. While there are a large number of
theoretical and experimental studies of adsorption below the critical temperature, there are
few studies of adsorption near or above the critical temperature of the fluid. Such studies
are especially relevant to the storage of methane at ambient temperatures (T, = 1.57) at
reasonably high densities.

Theoretical approaches to understanding and predicting adsorption range from
simple empirical fits (Freundlich/Tom isotherms) to theoretically-sound methods such as
molecular dynamics (MD) and Monte-Carlo (MC) simulations. Computer simulations
such as the grand canonical ensemble Monte-Carlo semi-quantitatively predict the cusp-
like behavior near the critical point [van Megan and Snook, 1982]. However, such
methods are computationally intensive. Simulations are difficult near the critical point due
to fluctuations, and require a large number of molecules and consequently significant
amounts of supercomputer time. Statistical mechanical theories such as the density
functional theory are also computationally intensive although they are about two orders of
magnitude faster than computer simulations [Gubbins, 1990]. On the other hand, the
traditional empirical and semi-empirical methods which are computationally undemanding
are unable to account for the wide variety of shapes of adsorption isotherms seen near the
critical region. There is an engineering need for a model that bridges the gap between the
simple and sophisticated models, and has predictive capability. Further, the development
of processes utilizing supercritical fluids requires engineering models capable of spanning

large pressure ranges. Experiments with supercritical fluids are more difficult and

expensive than experiments at atmospheric conditions, and a model that can predict
supercritical behavior would be a well-received developmental tool.

The focus of this thesis is to develop the simplified local density (SLD) model, and
demonstrate its application in predicting pure fluid adsorption isotherms over wide
pressure and temperature ranges by comparing model predictions with experimental data
from literature. The SLD model developed in this work superimposes the fluid-solid
interaction potential on the van der Waals and Peng-Robinson equations—of-state. My
intention is to present the concepts necessary to adapt common cubic equations of state
for describing the adsorption phenomenon. The SLD model is intended to bridge the gap
between the computationally intensive but more theoretically sound statistical mechanical
models and the undemanding, empirical methods. The SLD concept should be a useful
engineering supplement to other available models, and is not intended as a replacement for
the statistical mechanical models.

With the above objectives, the thesis has been divided into the following chapters:

1) Introduction
2) Background, Literature Review and Preliminary Results.

3) Quantitative Modeling of Adsorption and Extension of the SLD
Approach to Predict Clustering

4) Adsorption of Pure Gases in Slits and Pores, and Adsorption of Binary
Mixtures

5) Conclusions and Recommendations

In chapter 2, a detailed literature survey of the various theories of adsorption is
given, followed by the development of the van der Waals SLD model. High pressure
experimental adsorption data that show fascinating behavior are also given in this chapter.
The basic assumptions leading to the SLD model, and the model equations are described
in this chapter..

In chapter 3, the van der Waals equation is replaced by the Peng-Robinson
equation-of-state, thus making the model a quantitative predictive model. Model results
for adsorption on flat walls are compared to experimental data. The model is modified to
describe clustering in supercritical fluids, and model predictions are compared to
experimental fluorescence spectra.

The model is extended to describe adsorption in slits and cylindrical pores in
chapter 4. Model predictions for adsorption of pure fluids are compared to both
experimental data and molecular simulations. The model is extended to describe
adsorption of binary mixtures in slits. Model predictions of both the concentration and
amount adsorbed are compared with experimental data.

Finally, the results of this model are summarized in chapter 5, followed by

recommendations for future work.

CHAPTER 2: BACKGROUND, LITERATURE REVIEW AND
PRELIMINARY RESULTS

DEFINITION
While several different methods are used to measure adsorption, adsorption is un-

ambiguously defined by the surface excess, the excess moles per unit area.

It
A

S

(X

I‘CX =

 

L:[p(z)-p.]dz [2-1]

where p (z) is the density at a distance 2 from the wall, and pb the bulk density. 2,, is often
defined as the plane above the surface of the solid where the fluid-solid potential is zero.
Before reviewing the theories for adsorption, some experimental results will be
presented. Reliable experimental measurements of adsorption of fluids on well character-
ized solids are essential in order to test different theories of adsorption. While most com-
mercial adsorbents are porous and heterogeneous, graphitized carbon black (g.c.b.) or
exfoliated graphite or the basal planes of graphite provides an atomically smooth homo-
geneous and non-porous surface to study adsorption. Graphitized carbon black (or Gra-
phon) provides the ideal surface to understand physical adsorption and wetting, and theo-
ries for adsorption on homogeneous surfaces can be extended to heterogeneous and po-
rous materials. Some of the best data available are those of Findenegg and co-workers
who have examined the adsorption of a variety of gases such as ethylene, propane, meth-

5

6
ane, krypton and argon on graphon [Findenegg, 1984; Specovius and Findenegg, 1978,
1980]. Some of their results are presented in Figures 2.1 - 2.3, where one sees the ad-
sorption isotherms of subcritical propane and ethylene, near-critical propane and ethylene,

and supercritical ethylene and krypton, on graphon.

CLASSIFICATION OF ADSORPTION ISOTHERMS

Adsorption isotherms (subcritical) are classified into six types [Sing et al., 1982;
Sing 1983], as shown in Figure 2.4. Type I isotherms are normally seen in microporous
adsorbents, Type IV & V on mesoporous adsorbents and Type II, III and VI on non-
porous adsorbents. Types III & V are similar to Types II & IV respectively except for the
stronger magnitude of the adsorbate-adsorbent interaction in the latter cases. Type VI
isotherms indicate layering, and are observed for the adsorption of Ar/Kr on g.c.b. at low

(liquid N2) temperatures [Bienfait 1980 and references therein].

ADSORPTION MODELS

In this section some of the basic theories for physical adsorption will be briefly re-
viewed. There is an enormous amount of literature available on adsorption and the reader
is referred to several monographs [Young and Crowell, 1962; Brunauer, 1945; de Boer,
1953; Defay and Prigogine, 1966; Nicholson and Parsonage, 1982; Yang, 1987; Lee,

1988].

 

a

 

 

 

 

 

 

 

 

 

43—26315 K
40 .. +273.15 K
+283.15 K
35 __ +293.15K
I +313.15 K
a" i
g 30 ..
8
3: 25 --
III
e 20 --
0
X
I.”
I5 -
IO -
5 -
0 u I I I I I
0 20 4O 60 80 100 120

Pressure (bar)

Figure 2-1: Adsorption of Ethylene on
Graphon - Data of Findenegg (1983)

 

 

 

 

 

 

 

 

30 --
25 --
A
N
20 --
g
o
S.
.v 15 --
In
in
o
o
x
In 10 --
_. . . - __ ~ +273.15K
5 ; .‘ ~—~ ‘ +323.15K
(1,7,; " +343.15K
.. if +363.15 K
0 ' I : : I : : :'

0 5 IO ‘15 20 25 30 35
Pressure (bar)

Figure 2-2: Adsorption of Propane on Graphon -
Data of Findenegg (1983)

40

14

 

 

 

 

 

   

 

 

12 --
IO ‘-
‘MA
:2
O 8 .. ..
E
:3.
V
8 6 --
O
O
X
Ill
4 -_ +253.15K
+273.15K
+298.15K
2 -_ +323.15K
I +348.15K
"I +373.15K
0 I I I
0 4O 80 120

Pressure (bar) .

Figure 2-3: Adsorption of Krypton on Graphon -

Data of Findenegg (1983)

 

 

WomwQanmhv/N «Clogs/N

 

 

 

Fig

 

 

 

 

 

 

 

 

 

j

.0

(D

e.

O

U)

.0

<3:

'5

3

O

E

< Y In
P / P58“

Figure 2-4: Variety of Behaviors at Subcritical
Conditions

11
Henry's Law
In the low relative pressure region, the adsorption is linear, and the amount ad-
sorbed (n) is proportional to the pressure (P). This region is often called the Henry's law

region.

n = k P [2'2]

where k is the Henry’s law constant. Sometimes the pressure at which this law is obeyed
is below experimentally observable pressures. If the adsorption is measured at sufficiently
low pressures discontinuities in the adsorption-pressure curve are sometimes noticed indi-
cating phase transitions in the adsorbed phase as shown by the results of Ross and Clark
[1954], Fisher and McMillan [1957], and more recently by a host of modern techniques by

Thomy and Duvall [Bienfait 1980].

Langmuir Isotherm
One of the most popular empirically-adjusted models for adsorption is the Lang-
muir isotherm, a two parameter equation. This is one of the most popular isotherm due to

its extreme simplicity, and is often expressed as

 

e = = [2-3]

12
where W is the mass adsorbed, and W... the amount adsorbed in a monolayer. The Lang-
muir equation is only derived for monolayer adsorption and therefore predicts only Type I
isotherms. However, it is frequently applied to microporous materials. It is more suitable

for describing chemisorption.

BET Isotherm

The Brunauer, Emmett and Teller (BET) theory extends the Langmuir theory to
multilayers. It assumes that molecules adsorb in stacks or layers and that the uppermost
molecules in the adsorbed stacks are in equilibrium with the vapor. The BET theory treats
the first adsorbed layer separately and assumes that all layers above the first layer are con-

densed, and treats them identically as a liquid . This theory predicts that

I’ll”,

C
P P P
(“EX“TCE)

9-1-1- [2...]
Wm NM

 

where P, is the vapor pressure, N the number of moles adsorbed and C a constant related
to the fraction of surface that is uncovered. The BET equation is often used in deterrnin-

ing pore-size distributions and is often expressed as

1 1 C—lfi

W01; _1) = CW," + CW", P0 [2.5]

 

 

 

 

 

 

Tl

13

In pores where the number of layers are limited to n the BET equation is :

C{l—(n+1)[—I:] +n[£] }
N _ P. P.

N — n+1

.. (£_1){1+(C—1)P_C[_P_) }
P P0 P0

 

 

W
e _ W. _ [2-6]

 

This equation reduces to the previous equation when n is infinity and the Langmuir
isotherm for n = 1. The BET theory and its modifications are very p0pular since they
Show all the different types of subcritical isotherms. However, some of its limitations are
due to the assumptions made in the developing the theory 1) the surface is energetically
homogeneous; 2) there are no lateral adsorbate interactions; 3) adsorption is localized -
adsorbed molecules are fixed on sites and do not laterally move; 4) the heat of adsorption
in the second and higher layers is the same as the heat of condensation. Most of these as-
sumptions are valid for pressures in the range 0.05 < P/P‘ < 0.35, where the high energy
sites are filled but extensive multi-layer condensation has not commenced. The BET the-
ory is limited to subcritical temperatures, and assumes an ideal gas vapor phase, although

corrections can be applied.

Toth & UNILAN Equations

The Toth equation can be written as

l4

mP
n = ———l',7 [2-7]
(b + P’)
This is a three parameter equation m, b, and t, where m is the saturation capacity (mol/kg),

b and tare constants of the Toth equation, n is the specific amount adsorbed, and P is the

pressure. A similar equation called the UNILAN equation can be written

 

n = 31n[c+ Pi] [2-8]
c+Pe

23
where sis a constant. Both these equations reduce to the Langmuir equation (Toth - t =
l, UNILAN - s = 0). Both these equations are empirical equations, but are for heteroge-
neous adsorbents. The constants needed are temperature-dependent and atleast two iso-
therms must be measured experimentally before calculating the temperature coefficient of
adsorption. However, these theories seem to be the best overall for microporous adsorb—

ents, and have been used as the basis for multicomponent adsorption models which will be

discussed later in this chapter.

Two Dimensional Equations-of-State

The BET theory of multi-molecular adsorption is a generalization of Langmuir's
unimolecular theory. Both theories are based on localized physisorption, and both neglect
lateral interactions between fluid molecules, i.e. interactions between molecules in the

Same layer. A different approach treats the adsorbed phase as a two-dimensional fluid,

15
and uses a two-dimensional equation-of-state to represent the adsorbed layer [Hill, 1946,
1947, 1948; de Boer, 1953]. The two-dimensional equations-of-state allow for interac-
tions between fluid molecules in the first adsorbed layer. Neither the Langmuir nor the
BET theory allow for interactions between molecules, either on the surface, or between
molecules in different layers. When a fluid is adsorbed onto a surface it exerts a pressure
opposing the surface tension of a clean surface. If yo is the surface tension of the clean

surface and 7 that of the surface with an adsorbed fluid, then the spreading pressure It is

defined as

”=70“? [2-9]

The surface tension 7, or 1: can be related to the Gibb's surface excess

1‘" = "1391 [2.101

where a is the activity. For a dilute solution, the activity (a) may be replaced by the con-

centration C or pressure P, and a two-dimensional equation of state is used to give 1:. For

an ideal gas the equation is

mlzRT [2-11]

16
where A is the area per mole. This equation implies no forces or interactions between the
fluid molecules, which is what one would expect at low coverage. These two equations
directly lead to the Henry's law isotherm. Such an approach along with activity coeffi-
cients has been used to model adsorption from gas mixtures and from solutions [Myers
and Prausnitz 1965; Radke and Prausnitz 1972].

Several attempts have been made using a two-dimensional equation-of-state to
represent the adsorbed layer which accounts for mobile interactions and includes lateral
interactions in the first layer (The Langmuir and BET models account only for localized
adsorption - i.e. they neglect lateral interactions between fluid molecules) [I-Iill, 1946,
1947, 1948; de Boer, 1953]. The two dimensional van der Waals equation predicts a two-

dimensional critical temperature Tc; which is one half the bulk critical temperature (T c2 =
0.5 T C) and a two-dimensional critical pressure (I), = 0.361 P, ( Vc / N )1”. Here the sub-

script c2 is used to describe the critical constants for the two dimensional state. Similar
results for Tc; are obtained with different systems. A Lennard-Jones Devonshire [1937]
model gives T c2 = 0.53 Tc. An excellent discussion of two dimensional critical tempera-
tures and pressures from the van der Waals equation-of-state is given by de Boer [1953].
Typically Tczl T c is about 0.4 with extremes of 0.36 and 0.56. This lies in-between mean
field theory prediction of 0.5 and a 2-D Ising model prediction of 0.37 [Bienfait 1980].

The two-dimensional van der Waals model leads to the following equation

 

O
e]exp(-k, G) [2-12]

P
—=kOex
P . ,[1

0

17

The constant k, depends on a2 and b2 the two-dimensional van der Waals parameters, and
therefore only the adsorbate. A high value of k. implies strong inter-molecular forces
(fluid-fluid). The constant k; measures the strength of adsorption, and the stronger the
adsorption the smaller the k2.

This equation predicts phase transitions on the surface for weak adsorption forces
and strong intermolecular attraction (i.e. large k, and k2). If the fluid-fluid intermolecular
attraction is weak, no condensation is seen. For strong adsorption forces, a tendency to-
wards saturation (Type I isotherm) is seen, and the adsorbed fluid is either a compressed
supercritical fluid or a condensed two dimensional liquid.

Relative to the BET equation, the Hill de-Boer model predicts lower pressures for
coverage's below 0.8, but over-corrects the errors of the BET equation. This is due to the
fact that lateral interactions are allowed, and not due to the mobility. If lateral interactions
are allowed in a localized first layer, condensation takes place at very low pressures, and
jumps from almost zero to near unity at a critical coverage 9, of 0.5. Such a result could
probably be deduced by percolation on a lattice, since the similarity between the abrupt

jump at 9c of 0.5 and a critical percolation probability of pc of 0.5 does not seem coinci—

dental.

18
Potential theories
Several other popular theories are the Polyani potential theory [Polyani 1932] and
modifications of it such as the Dubinin-Radushkevich [Dubinin 1947, 1975], and the

Frenkel-Halsey-Hill [Halsey 1948; Hill 1950] theory.

(i) Polyani Potential Theory

The Polyani potential theory assumes that a potential field exists at the surface of a
solid that attracts the molecules of the surrounding fluid. These attractive forces are es-
sentially London dispersion forces or van der Waals forces, which vary as 1/r3 (where r
is the distance from the surface). The adsorption potential (8;), at a distance i from the
surface is the work done in bringing a molecule from the gas phase at density Px to that
point where the density is p i° Polyani noticed the similarity with gravitational potential,

and used the equation of hydrostatics to calculate the adsorption.

e, =[p'vdp [2-13]

I

It Was assumed that the adsorbent potential at a given distance (or surface at) i is e,- , re-
gardless of the number and kind of molecules located between the i'th surface and the
SOlid or gas, an assumption which was justified later by London's theory of van der Waals
fOl'ces. Each of the potential surfaces along with the surface of the solid adsorbent en-

clOSes a volume 6;. At the surface of the adsorbent 80 is maximum and (l),- is zero. The

19
gas phase is assumed to be the region in which 8 = 0, and limits the adsorption volume to
d) = ‘l’max- In actuality, the potential never goes to zero but for practical purposes 8 = 0 at
a finite distance and therefore a finite volume from the surface. The theory assumes that
de/dT = 0 implying that do/dT = 0 and therefore that e = f (4) ) for all temperatures for a

given adsorbent. This is used to form a characteristic curve, which is experimentally de-

termined. The amount adsorbed (n) is given by
41....
n= [0 (p—p.)d¢ [2-14]

In using the theory, first a characteristic curve 8 = f (q) ) is created. This curve should
span the entire adsorption space, and is therefore typically evaluated at a moderate tem-
perature. The second step consists of evaluating adsorption isotherms from these charac-
teristic curves. The characteristic curve equations can be simplified if the temperature is
significantly below or above the critical temperature. These equations ignore adsorption
in vapor form, and the compressibility of the liquid layer, as well as assume ideal behavior.
Near Tc these assumptions result in errors, but the different errors fortuitously cancel re-
sulting in fairly accurate predictions. For high temperatures and pressures the Polyani
theory has been modified by using empirical correction factors to account for the com-

pressibility and thermal expansion of the adsorbed liquid [Berenyi 1923; Brunauer 1945].

20
(ii) Dubinin's Theory
An empirical relationship was found between 80 of Polyani theory and the van der

Waals constant 'a' of the fluid :

5,, = NE [215]
which originates from the mixing rules of Berthelot, who proposed that the attractive po-
tential between two molecules at a fixed distance is proportional to ,Ia,a2 where a 1 and a;
are their van der Waals constants. Using this, the potential theory was extended by Du-
binin and co-workers to include adsorption of mixed gases onto a solid. They expressed

the characteristic curve for all gases on the same adsorbent as :

8 = Bf(¢) [2-16]

where B is called the coefficient of affinity and is defined as the ratio of adsorption poten-
tials (B = 8 ISO) of a gas and a reference fluid, on the given adsorbent. They found that the
molar volumes of the adsorbate in liquid state were closely proportional to the coefficient
of affinity. The characteristic curves resembled the positive side of a Gaussian curve
which lead Dubinin and Radushkevich to assume that the adsorption volume V occupied

by the liquid adsorbate (and V, by the reference) followed

V = V0 exp[—K[%) ] - [247}

where K is a constant depending on the shape and size distribution of pores. This leads to

2 2
RT P0
ln( W) — ln(VoP) - K17?) [ln(?)] [2-18]

Subcritical vapor pressure data have been extrapolated into the supercritical region and a

good fit has been obtained using the extrapolated fugacities [Maslan et al., 1953].

(iii) Frenkel-Halsey-Hill Theory

This theory presumes that the adsorbed film consists of a liquid at a density pl

(density of bulk liquid at that temperature) with a thickness of h. The surface excess (fer)

is given by

I‘" = h(pI —p) [2-19]

Where p is the density of the bulk gas. Assuming that the adsorbate and adsorbent interact

tlll‘ough a 12-6 Lennard-Jones potential (and integrating to a 9-3 potential), this leads to

ln(-5] = —-k— [2-20]

where k is a constant.
In a further modification of the Frenkel—Halsey-Hill theory, multiple liquid like
slabs and fugacities have been used to represent the adsorbed layer that has a density p...

for 2 S 2. and p; for z, < z < 2:, where the fugacity is related to adsorption by :

“’35 5"— [2-21]
kT 2,1

5‘
m
b» Ix
w

I:

Here f, is the fugacity of the saturated vapor, z, the distance of the i'th statistical layer of
molecules, AI: is the perturbation energy of a molecule in the first layer at a distance 21.
For a single step density profile if the density of the layer is p; for z < z,- and p for z > 2,,

then 2,. can be determined from experiments and

 

[2-22]

and the parameter m fit from a plot of 1n(z,-) vs. ln(-ln(fo/f)). Such methods with multiple-
step density profiles have been used to describe adsorption generally in the sub-critical

range [Findenegg 1983].

23

van der Waals model of Barrer and Robins

This theory, developed by Barrer and Robins [1951], has been unrecognized for
more than forty years. It uses the fact that the chemical potential is constant at all dis-
tances from a wall, and can be determined by the temperature and pressure far from the
surface. This potential consists of two components - potential due to the interaction of the
fluid molecules themselves, and a potential due to interaction of the fluid and the wall.
Since the potential exerted by the wall and the bulk fluid potential are known the density
profile and surface excess can be calculated. They assumed that the fluid obeyed the van
der Waals equation-of-state and the fluid-wall interaction can be approximated by the 9-3

Lennard-Jones potential. This leads to the following equation

 

 

 

0 exp( 0 _2a0_ SIN )= 9,, exp 0,, _2a0b [2_23]
1-9 1-9 bRT rRT l-Ob l-Bb bRT

where a and b are the van der Waals constants, N and n are the total number of molecules

and moles in the system respectively, 0 = bn/V, where V is the total volume occupied by
closely packed molecules, r is the distance between the fluid molecule and the surface, C is
an empirical constant and the subscript b stands for bulk properties. If v is the molar vol-
ume of the fluid being adsorbed on the surface of the adsorbent then equation 2-23 may be

written as

24

l b 20 CN 1 b 2a
exp - — 3 = exp — [2-24]
v—b v—b vRT rRT vh—b vh—b vhRT

 

 

 

 

This model predicted the different adsorbent isotherms whose magnitude could be
manipulated by changing the empirical constant C. Hill [1952] suggested that this model
can be improved by treating the fluid not as a bulk fluid as Barrer and Robins had done,
but as an inhomogeneous fluid i.e. the chemical potential should not be calculated from the
local density but from the entire density profile. He also suggested replacing the inaccu-

rate van der Waals equation by a better equation.

Statistical Mechanical Models

There are different types of models derived from chemical physics, including the
Lattice Gas model, Integral Equation theory, and Density Functional theory. The Lattice
Gas model is a case of localized adsorption, where the molecules of the solid are large
compared to the molecules of the fluid. This results in the formation of a lattice of ad-
sorption sites that are separated from each other by high potential barriers. This model
simplifies to the Langmuir or BET approaches and predicts several phase transitions on
the surface.

Theoretical interpretations of adsorption include virial expansions in terms of bulk
density at low pressures [Steele 1962]. A van der Waals theory for adsorption has been
proposed by Sullivan [1979]. This approach has been generalized to examine wetting

transitions at the solid-gas interface [Tarazona and Evans 1983], adsorption and wetting

25
transitions in binary fluid mixtures by Telo de Gama and Evans [1983]. The van der
Waals density functional model for inhomogeneous fluids has also been used to study
wetting transitions at surfaces [Kung et al., 1990; Teletzke et al., 1982a, b]. More recent
applications of statistical mechanics have involved adaptation of integral equations devel-
oped for homogeneous fluids to inhomogeneous adsorbed fluids. These include integral
equations for liquid state theories such as the BBGKY hierarchy for non-uniform liquids
[Lee, 1987; Henderson, 1992] and density functional theory [Saam and Ebner, 1978;

Fisher and Methfessel, 1980; Tan and Gubbins, 1990; Kierlik and Rosinberg, 1992].

(i) Integral Equation Theory

Statistical mechanical theories use distribution functions to describe the correlation
between fluid molecules and hence the short range order. The primary cause of this short
range order, or structure, of a uniform fluid are the repulsive forces, with attractive forces
forming the background. In order to get a complete theory of liquid state an expression is

needed for g( r) the radial distribution function [McQuarrie, 1976; Rowlinson and Widom,

1982], or the pair correlation function.

-pu
2 I e ”dr...dr
V N. i I 3 N [245]

gm = N2(N-2)! 2,,

Where, there is a system of N particles in a volume V and at temperature T, ZN is the con-

flgur‘ational integral, r1, r2, r3 are coordinates in space. If the total potential of a N-body

26
system can be assumed to be pairwise additive, then all the thermodynamic properties of

the system can be derived from the radial distribution function. For e. g. the pressure is

related to g( r) by

3:— "Pi” ' 4 2d [226]
kT—p 6kT 0 ru(r)g(r) Itr r -

A number of approximate equations have been used to solve Equation 2-26. These in-
clude the BBGKY hierarchy and Omstein-Zernike equations with the Percus-Yevick and
hypernetted-chain closures. These approaches are reviewed by McQuarrie [1976] and

Henderson [1992], and are widely used in predicting clustering in supercritical fluids.

(ii) Densig Functional Theog

The density functional approach was derived by Saam and Ebner [1978]. It is

based on the fact that there exists a functional 9 of the local density p( r), such that at

equilibrium this functional is minimized.

F = F[f(r)s p(r), C(r‘ ,r2 )] [2'27]

Where f( r) is the local free energy, p( r) is the local density and c( r1, r2) is the direct corre-
lation function. Gubbins and co-workers have used this approach to describe adsorption is

slits and cylindrical pores and have successfully compared their results to molecular simu-

27
lations. Kierlik and Rosinberg [1991, 1992] proposed a simplified version of the free en-
ergy functional for the inhomogeneous hard-sphere fluid mixture that requires four distinct

weight functions and generates a triplet direct correlation function for the one component

fluid.

Molecular Simulations

While there a large number of complicated theories they are not exact. The best
solutions, for a given potential, are obtained by computer simulations where the equations
of motion are solved for the particles of an ensemble. Observable macroscopic quantities
are then calculated by time-averaging the appropriate microscopic equations. Such a cal-
culation is called a molecular dynamics (MD) calculation. In the Monte-Carlo (MC)
scheme a large number of configurations are sampled, and particles moved to find the
most stable state, and averages are calculated to give a certain macroscopic property.

Computer simulations are difficult to perform near the critical point due to the
large density fluctuations. However, these techniques are very powerful, and yield the
best results for a given potential. van Megan and Snook [1982, 1985] have simulated the
adsorption of fluids near the critical region using a grand canonical ensemble Monte-Carlo
method. They compare the results of their simulation with experimentally measured
[Specovius and Findenegg, 1978] isotherms of ethylene on graphite, just above the critical
temperature and find reasonably good qualitative agreement with data (Figure 2.5). The
simulation results show cusp-like maxima which are close in magnitude to experimental

values, but the location of the cusp is off by 25% (~15 bar). Their results cannot be

28

 

 

 

 

 

20 ‘ L0 50 80 ‘ 100
plborl -

Figure 2-5: Comparison of Molecular Simulation (van
Megen and Snook, 1982) to Experimental Data

(Specovius and Findenegg, 1978)
'- - Molecular Simulation (T / r, = 1.03)

-- Experimental Data (T / Tc = 1.02)
- Experimental Data (T / T‘: = 1.04)

29
quantitatively compared to experimental data because their simulation of bulk ethylene
uses a shifted potential and does not represent the experimental bulk properties of ethyl-
ene; e.g. the simulation critical temperature (Tc, ,,-,,,) of ethylene is 226 K, as opposed to an
experimental value of 282.4 K. However, their results have been the basis for comparison
for the density functional theory using the same shifted potential. Figure 2.6 and 2.7 Show
the adsorption of ethylene on a slit like pore at sub and supercritical temperatures. At low
slit widths Type I isotherms are seen. There is also a slight decrease in the surface excess
with increasing pressure. As the slit width increases type IV isotherms are seen. In flat
walls, type H isotherms are seen. Figure 2.8 shows the density functional theory’s predic-
tion of van Megen and Snook’s simulation [Tan and Gubbins, 1990], and it can be seen
that the density functional theory does a very good job of reproducing the results of the
molecular simulation. Snook and van Megen conclude that traditional adsorption theories
will not be able to account for the wide variety of shapes of adsorption isotherms encoun-

tered in their computations.

Mixed-Gas Isotherms

Current models that predict mixed-gas isotherms use pure gas models such as the
Toth or UNILAN equations. One such theory is the ideal-adsorbed-solution theory of
Myers and Prausnitz [1965]. This theory is thermodynamically consistent and exact at the
limit of zero pressure. The success of such theories depends of their ability to predict the
pure gas isotherms accurately. Small errors in the prediction of pure-gas experimental

data may cause large deviations in the mixed-gas predictions. Other sources may be the

30

 

 

 

 

 

 

0 '2 ’4 -6 -8 ,1-0
Figure 2-6: Grand Canonical Ensemble Monte Carlo

Simulation of Ethylene in a Slit (TIT c=0.85)
-h=~, - - - 11:10.0, ----h=5.0, - ----- h=3.5

31

 

 

 

 

 

 

O .1 «2 ..
D
Figure 2-7: Grand Canonical Ensemble Monte Carlo

Simulation of Ethylene in a Slit (T IT e=1 .2)
-h=~,~uh=50,~~uh=35 .

32

 

 

 

 

OUIIIISIIIIILI'4':IIII
0 0.1 0.2 0.3 0.4

. p
Figure 2-8: Density Functional Theory - Comparison to
Grand Canonical Ensemble Monte Carlo Simulation

 

 

33
neglect of surface heterogeneity and adsorbate-adsorbent interactions, both of which may
cause non-ideal behavior [Valenzuela and Myers, 1989]. Later models of this type, did
not neglect these heterogeneities and were called non-ideal adsorbed solution theories.
Since the basis of mixed-gas isotherm modeling is the IAS theory, we shall discuss it in

more detail.

Ideal Adsorbed Solution (IAS) Theogz

For a multicomponent system containing N adsorbates, the IAS equations for a

perfect gas are

Py.=P."x. (i=1,2, ..... ,N) [2-28]

where the adsorption isotherm of pure ith adsorbate is n,"(P,-"), y,- is the mole fraction in the

gas phase and x,- is the mole fraction in the adsorbed phase.

‘I’,°(P,") = ‘11: (P20) = ...... = \I';(P,3) [249]

ix, =1 [2-30]

The function ‘I’,°(P,~") is the integral for the spreading pressure (ILA/RT) for the adsorption

0f Pure ith adsorbate at the same temperature as the mixture.

34

 

M_Pn _»wmm
——_[ FdP—J d [2.311

n
RT 0 0 d ln(n)

The independent variables are T, P, y1,...y~.r, leaving us with 2N unknowns (Pf, x.) and 2N
equations (2-28 - 2-30). After these three equations are solved, the total amount n, and

the amount of ith component are found as

l N x.
— = 2-3— [2—32]
n: i=1 "i

n.- = M.- [2-33]

The general algorithm for solving this theory is given in Valenzuela and Myers
[1 989]. Further development of non-ideal adsorbed solution theories, and the develop-
ment of the fast-IAS theory, where the integrals for the spreading pressure are written as a
series expansion in terms of the central moments of the adsorption energy distribution, are
also given in Valenzuela and Myers [1989]. All these theories require the pure gas iso-

therm parameters (for Toth/UNILAN it is 3), and hence are not predictive models.

35

PRELIMINARY RESULTS - SIMPLIFIED LOCAL DENSITY MODEL OE
RANGARAJAN, LIRA AND SUBRAMANIAN

This model, originally developed by Rangarajan [1992], and Rangarajan, Lira and
Subramanian [1995] forms the basis for the rest of this dissertation, and hence will be dis-
cussed in detail. Some of the predictions of this model will also be presented and com-
pared with experimental data. The model has a basis similar to the Polyani potential the-
ory, where the adsorbent exerts an attractive potential on the adsorbate. The model is also

Similar to the approaches of Barrer and Robins [1951], Hill [1949], Sullivan [1979], and

Kung et al. [1990].

Model Development

The attractive potential between the fluid and solid, at any point z, is assumed to
be independent of temperature and the number of molecules at and around that point. At
equilibrium, the molar chemical potential [.1 must be uniform throughout the system

[Denbigh 1981; Guggenhiem 1967] and is calculated by contributions from fluid-fluid and

fluid-solid interactions.

u = u, = 11,12) + 11.12) [2-341

where the subscript 'b' refers to the bulk, the subscript 'fj“ to the fluid-fluid interactions and

the subscript 'fs' to the fluid-wall interactions. Equation 2-34 requires that at a distance 2

36
from the wall, the sum of the chemical potential due to the fluid-fluid interactions and the
attractive potential exerted by the solid on the fluid remains constant, and equal to the
bulk chemical potential. While such an equation has been written from the principles of
Cherrlical equilibrium, this fundamental result for inhomogeneous systems has been derived
by minimizing the grand potential [Evans 1979; Rowlinson and Widom 1982; Davis and
Scriven 1982]. In using this equation care should be taken to ensure a consistent basis for
chemical potential, whether it be per molecule or mole. If ‘I’( z ) is the potential exerted

by the wall on a fluid molecule , then on a molar basis,

11,, (z) = N.‘P<z) 12-351

where NA is Avogadro's number. Therefore
11,12) = u. - N.‘P(z) [2-361

For a non-ideal fluid :

MT): n.1T)+RT1n[%) 12-371

0

Where fb is the bulk fugacity and f, is 1 bar. For an inhomogeneous fluid,

37

 

3.12)] [”8]

0

pfl(T) = po(T)+ RTln[

Where ff,(z) is the fluid-fluid contribution to the fugacity of the fluid at a position z. The

Equations 2-36 - 2-38 lead to

 

-‘P(z)] [2_39]

We...) ..

It should be noted that ‘I’ is negative when attractive, leading to an increased fluid-fluid

fugacity near the wall. The fluid chemical potential consists of a repulsive contribution

prep and an attractive contribution um, [Vera and Prausnitz, 1972]

 

 

”fl = ”rep + ”at! [2-40]
For a homogeneous fluid, using the van der Waals equation-of-state
RT b
= RT 1n + 2-41
”W l: (v-b) v—b] [ l]

Where v is molar volume and b is the van der Waals constant. The attractive potential is

given by

38

u... = [unfldv [2—421
V V

where V denotes the volume of integration (all space occupied by the fluid), (p(r) the two-
body interaction potential, g( r) the radial distribution function, which is taken to be a con-
stant outside the hard sphere diameter for a van der Waals fluid [McQuarrie, 1976]. In

the case of a homogeneous fluid (v ¢ v(z)) Equation 2-42 simplifies to the common form

“at: = —_- [2‘43]

where ab is the van der Waals constant. In order to evaluate the integral in Equation 2-42

for an inhomogeneous fluid, if the density changes with position more gradually than ¢( r)
does, then one needs to include density as a function of position. Equation 2-42 suggests
that a feasible approximation might be to use a density at the local position for evaluation

of the integral. In other words, the fluid at point z is treated as a homogeneous fluid at a
density of p(z) (such an approximation is the one used by Barrer and Robins [1951]).
Since the two body potential ¢( r) is short ranged (oc 1/r6), and the product of 1])( r) and
p(z) appears in the configurational integral (Equation 2-42), it seems reasonable that most
of the value of the integral (at 2) results from contributions of (p(r) and p(z) near 2. Far

from the given point z, where p(z) changes, the two-body potential (¢(r)) will approach

39

zero. Since the density of the fluid is larger closer to the attractive wall and smaller farther
away from the wall, some of the errors introduced by this approximation cancel, although
the approximation will not be as good near a phase transition. The term 'local' refers to
the fact that all thermodynamic properties at point z are calculated using a single density
value, p(z), a 'local average' density, and are not calculated from gradients in density
about the point z. This approximation is called the local density approximation.

Using an equation-of-state, an expression for fugacity in terms of the molar density
or volume can be derived. The van der Waals equation-of-state leads to the following ex-

pression for fugacity for a homogeneous fluid:

 

 

 

fb=R

T exp[ b 2“”) [2-44]
v-b

V-b-VRT

For an inhomogeneous fluid, the problem is one of calculating the density profile p(z)
which satisfies Equation 2-39. Following Hill [1951], Sullivan [1979], and Kung et a1.

[1 990] we assume the short-range repulsive forces determine pm, by the local density. If
We assume that pm, is also determined by the density of p(z), (local density approximation)

then we can rewrite Equation 2-44 to calculate density as a function of z

40

[ b _ 20(2)] [245]
v(z)-b v(z)RT

 

f” (2) = v(z) — b ex

where a(z) is evaluated from the integral given in Equation 2-42, in which the molar vol-
ume in the denominator is approximated by the local molar volume at 2. In the model of
Barrer and Robins [1951], the van der Waals a is independent of z (a(z) = ab). In the pre-
sent model, we exclude from the integral of Equation 2-42 that portion of space occupied

by the solid adsorbent. This leads to the following expressions for a(z) (see Appendix 1):

 

a(z)=a,, —5—+£-£— for 0.55—z—Sl.5 [2—46]
16 16c” 01,
a(z)=a,, 1- 1 3 for ass—Ls 1.5 [2.47]
[z I] or
g __-__
L “If 2 _

 

 

The van der Waals constant a is related to the critical temperature (Tc) and the critical

Pressure (PC) of the fluid. If the excluded volume does not change, then the b term re-

mains constant. Since

41

 

 

2
. £7211; [2-481
b = 18:; [2-49]

The wall critical temperature is related to the bulk critical temperature as

__a...,, = 0.5 = _____(Tm" / Tc )2 [2-50]
ab Pcwall ch

and combined with the equation for b yields Tm” = TC/Z. This is the mean-field theory

prediction of a reduction of the critical temperature at the surface.

Outline of algorithm to solve for densigr profile

For a fluid at a given temperature and pressure, the bulk fugacity fb can be evalu-

ated. For a given adsorbent if the potential ‘P(z) is known, then Equation 2-39 specifies
the local fugacity. With the assumption of uniform local density, Equation 2—45 can be
801 Ved, at each point z, for v(z) or its reciprocal the density. For a given temperature, bulk
PreSsur-e (fugacity) and adsorbate-adsorbent system, the molar volume (or density) is
Spe<2ified at each point by Equation 2-39 and 245. However, below the critical point ~
thel‘e may be three values of v which will satisfy Equation 2—45, and the question then be-

Comes that of choosing the correct value of v. Since the fugacity rather than pressure is

42
Specified, the Maxwell equal-area construction cannot be used to select the root. The
different roots of Equation 2-45 give the same fugacity, but give different values of pres-

sure P, where P is given by the van der Waals equation :

RT _a(z)
v(z)-b v2

 

P(z) = [2-51]

This difference in pressure is used to select the correct root. A plot of fugacity versus
pressure for a pure fluid is shown in Figure 2-9. Regions AB and BE represent the stable
vapor and liquid phases respectively. Regions BC and DB are the meta-stable parts of the
isotherm, and CD represents the unstable portion of the van der Waals isotherm. It may
be noticed from the curve that for a given fugacity, the stable root is always at the highest
pressure. The roots to the fugacity Equation (2-45) are substituted into the van der Waals
equation and the root which gives the highest pressure is selected.

Equation 2-45 can be re-written in the following forms:

 

RT b 2a(z)
= b —— — 252
Va) + { f(z)iexP[ v(z) — b v(z)RT] [ l

 

 

 

43

Figure 2-9: Fugacity vs Pressure for a Pure Fluid

b

f(z){v(z) - b}] + 2a(z)
RT vRT

[2-53]

 

 

v(z) = b +
ln[

Equations 2-52 and 2-53 are solved by successive substitution. Equation 2-52 is written
to provide an increasing value of v(z), and for an initial guess the value of v(z) is set equal
to RT/f (2). Equation 2-53 is started off with v(z) = 1.1 b, and rapidly converges on the
liquid root. Checks are placed to ensure that all the realistic roots at a point are obtained.
The computer program used here is attached in Appendix 2. At a given temperature, this
program utilizes the critical temperature of the fluid, the fluid-solid interaction parameter,
and the size of the solid and fluid molecules to calculate the surface excess as a function of

pressure, and the density as a function of position.

Results of the SLD Model

The SLD model discussed here, superimposes the attraction of the fluid molecules
to the wall onto the attraction of fluid molecules to one another. The asymmetry of fluid-
fluid interactions near an inert wall is incorporated into the equations for a, the attractive
Constant. Several limiting cases are discussed below and compared with theory, computer

Simulations and experiments.

45

1. Ideal Gas - Hard Wall

For an ideal gas and hard wall, there are neither fluid-fluid intermolecular attrac-
tions, nor fluid-wall attractive forces: Therefore a = 0, b = 0, and ‘I’(z) = 0, f = P, and

Equation 2-39 reduces to P( z) = Pb" When combined with the ideal gas law, this reduces

to p(z) = p1,, and as expected there is no adsorption.

2. Ideal gas - Attractive Wall
Here a = 0, b = 0, f = P. Equation 2-39 exactly reduces to P(z) = Pb exp [-

‘l’(z)/(kT)]. With the ideal gas law this gives p(z) = p), exp [- ‘I’(z)/(kT)]. The surface ex-

cess is
re! = L [P(z) — Pbldz = pb L {exp[——k-T—] - 1} dz [2-54]

Since the integral is constant, and p), proportional to the pressure Pb, this leads to Henry's

Law.

3. Attractive Fluid - Hard Wall

 

For an attractive fluid and hard wall, a > 0, b > 0, and ‘I’(z) = 0, f(z) = fb. The
InOde] predicts a reduced density near the wall, because of attraction of the fluid molecules

near the wall to fluid molecules in the bulk. In fact a vapor-liquid phase transition is seen

46

adjacent to the wall for certain cases when the bulk fluid is liquid. This is similar to the

observations of Abraham and Singh [1978].

4. Atflctive Fluid - Attrzyltive Wall

For this case, a > 0, b > 0, ‘I’(z) < 0. In general, for an attractive fluid near an at-
tractive wall an increase in density is exhibited near the wall depending on the magnitude
of ‘I’(z). If the wall—fluid forces are stronger than the fluid-fluid forces, the fluid Will wet
the wall, otherwise there will be a rarefaction of gases near the wall. To represent the in-

teractions of the adsorbate and adsorbent we have chosen the partially integrated 10-4

potential model [Lee, 1988]

. 4 1
‘I’(z) = 47rpamef,o,3[5:—:f— m— -Zx—4] [2-55]
1’1 I
Where pm," = 0 382 atoms/ A2 , x,- is the intermolecular distance between fluid molecular
Centers and the ith plane of solid molecules, where we have truncated the interactions at
the fourth plane of solid atoms, where the interplanar spacing is 3.35 A. sf, is the well
depth of the fluid-solid potential, and is obtained from the Lorentz—Berthelot mixing rules

Which state that efs = (8,; 8:3)0'5 and of, = 0,5 (of, + 0'”). Table 2-1 summarizes other pa-

r atneters used in this work. Note that Efl‘ is not tabulated because the fluid-fluid contribu-

tion is calculated from Equations 2-46 - 2-48 and not directly from 8,]. Steele [1974],

47

Table 2-1: Molecular properties used in SLD model

 

 

 

 

 

System 6/] (nm) of, (nm) £,..//< (K)
Ethylene-Graphon 0.422 0.38l 450. 23
Krypton-Graphon 0.35 0.345 225

 

 

 

 

 

48
Nicholson and Parsonage [1982], and Lee [1988] offer a further discussion of fluid-solid
potentials.

Figure 2-10 shows adsorption isotherms that were calculated using the 104 po-
tential of Equation 2-47 and the parameters for ethylene and graphon in Table 2-1. The
curves are qualitatively similar to Figure 2-2. A knee is present, the magnitude of which
can be increased by increasing the solid-fluid interaction energy 8,, (Figure 2-10 also
shows calculated adsorption isotherms for condensed phases at subcritical temperatures
above the bulk vapor pressure which are not present in Figure 2-2). Some characteristic
features of adsorption of near-critical fluids are the cusp-Shaped isotherms seen near the
critical point, and the 'cross-over' of adsorption isotherms at different temperatures. Con-
sider the effect of temperature at a fixed pressure below the critical pressure. Below the
critical pressure, adsorption decreases with increasing temperature, but this trend is re-
versed above the crossover region at pressures above the critical pressure. The crossover
of the isotherms is eliminated by plotting the calculated surface excess with respect to cal-
culated bulk density as shown in Figure 2-11. The reason for the cusp-like behavior of a
supercritical isotherm in Figure 2-10 can be understood from the calculated density pro-
files shown in Figure 2-12. Below the critical pressure region, the adsorbed layer in-
creases in thickness faster than the bulk density increases. Above the critical pressure re-
gion, the bulk fluid becomes increasingly incompressible, and the bulk fluid density ap-

Proaches the density of the adsorbed fluid, causing the surface excess to decrease.

49

 

30

25 --

Excess (ll mol/mz)
61' B

—l
O
I

  

 

. ---263.15K
'1 : ------ 273.15K
, .' —283.15K
,' ----293.15K
I ----303.15K

 

 

 

 

 

 

20 40 60 80 100

Pressure (bar)
Figure 2-10: Adsorption of Ethylene on
Graphon

120

Excess (ll mollmz)

25

50

 

20 -~

—I
01
l

d
O
l

 

 

-—-283.15 K
--- 293.15 K
------ 303.15 K

   

 

 

 

.....
'O

l l l l

 

 

0

l
I I I I I

‘ 0.002 0.004 0.006 0.008 0.01

Density (moles/ems)
Figure 2-11: Adsorption of Ethylene on
Graphon

 

0.012

Density (moles/ems)

51

 

 

 

 

 

     

 

 

 

0.02
- - - 38.83—bar
0.02 -- ------ 42.76 bar
—47.26 bar
00] __ ----so.04 bar
' - - - — 54.22 bar
0.01 --
0.01 -— ‘5:
0.01 -- \ ................... _
l ‘,
I :
0.01 —- I g
1 '. a.
I 2‘ -----------
[1m —-— \\ .
\\ ________________________________________
0.00 —— """"""""""
0.00 I : I
0 5 10 15

Zle
Figure 2-12: Density Profile of Ethylene on
Graphon at 283.15 K

 

20

52

Type III isotherms are normally observed when the attraction between the solid
and fluid is weak. Figure 2-13 demonstrates the capability of the SLD model for predic-
tion of Type III isotherms. In this case, the fluid parameters are for ethylene, but the
magnitude of the 8,, has been decreased. The smaller 8,, results in elimination of the knee
and yields a Type III adsorption isotherm. At very low pressures and temperatures, a dis-
continuity in the adsorption isotherms is predicted (not shown here), indicating a phase
transition in the adsorbed phase similar to those discussed for a two-dimensional model
[Ross and Olivier, 1964]. At high pressures, when the bulk fluid exists as a liquid, the
SLD model can predict a negative surface excess when Efs is small. Negative surface ex-
cesses have also been predicted by Sullivan [1979]. Figure 2-14 shows some predicted
adsorption isotherms of krypton on graphon, at temperatures far above the critical tem-
perature. When compared with experimental data shown in Figure 2-3, the predictions are
again qualitatively correct. This model predicts some gas-liquid transitions on the surface,

similar to the transitions seen in the experiments of Thomy [Bienfait, 1980].

53

 

35

 

——203K
---233K
------ 263K

 

 

 

25 ~-

Excess (ll mollmz)
51' B

._.I
O

l

T

 

 

“’—
-—-‘-"—
...-.....................C.........-..

 

Pressure (bar)

Figure 2- 13: Weak Adsorption of Ethylene on
Graphon

 

54

 

 

 

 

 

 

 

 

Figure 2-14: Adsorption of Krypton on
Graphon

14
12 i
10 ‘— II” \\\\
“A [I ........................ \.
g z’/ \K‘;--
0 8 -- x’ ——
E I, ................ \
=- I .o -----
V / o .d’.
m I .'. "/
m 6 "' /.’ ’z"
a) .- .-
0 I. .1
. x 0‘ I
In .' ,o’
4 '1 [O
i ——253K
J ---273K
2 ------ 298K
—348K
- --373K
0 i I I 4 4 I I
0 20 40 60 80 100 120
Pressure (bar)

 

55
SLD Model - Discussion

The SLD model builds on the approaches of the simple theories such as Langmuir,
BET, two-dimensional equations-of—state and the van der Waals approach of Barrer and
Robins. By treating the fluid with a van der Waals equation with a suitably modified a, the
SLD model allows for interactions between adsorbed molecules at various distances from
the wall. All of the models mentioned above (including the SLD) assume an energetically
homogeneous surface. The effect of heterogeneity will be pronounced at extremely low
pressures and coverages, where the high energy sites are unfilled. To represent a hetero-
geneous surface, one could fit the potential to adsorption in the Henry’s law region, and
then work with a pseudohomogeneous surface.

The limitations of the model can be attributed to: (i) the lack of structure in the
fluid; (ii) the use of the local density approximation; and (iii) the use of the van der Waals
equation to describe the fluid properties. Since we are using the van der Waals equation
and a mean field approach, we do not predict any of the fluid structure seen using either
computer simulations or density functional theory [Snook and Henderson, 1978; Vander-
lick et al., 1988; Kierhk and Rosinberg, 1991]. This model cannot describe discrete fluid
structure or be used to study packing near a wall. The use of the local density approxima-
tion also results in the physically unrealistic prediction of abrupt vapor-liquid interfaces.
Despite the fundamental theoretical limitations of the model, since model predictions
mimic experimental trends, this model may be acceptable for engineering calculations. In
order to estimate the errors introduced by the use of the local density approximation,

Equation 2-42 was solved keeping the density inside the integral, leading to an integral

56
equation (IE) [Pyada, 1994]. When compared, solutions to the SLD model and IE ap-
proach showed differences which were dependent on pressure, temperature and magnitude
of the gas-solid interaction potential. For the adsorption of ethylene on graphon, the SLD
model tends to underpredict adsorption (relative to the integral equation (IE)) by about
10-20% between 0.9 < Tr < 1.1, except where the surface excess increases steeply at high
pressures. At pressures near 1 bar, the differences are less than 1%. These comparisons
show that the local density assumption can provide a reasonably approximate solution to
the adsorption problem.

A limitation of the van der Waals-based SLD model is the fact that predicted
magnitudes of adsorption and vapor pressures are incorrect. For accurate prediction, the
equation-of-state must be able to correctly predict the vapor pressure and liquid density.
The inaccuracies of the van der Waals equation in this regard are well documented. In
order to represent the fluid properties any equation-of-state can be used, provided the re-
pulsive and configurational contributions can be separated for use in Equation 2-39. The
selection of the van der Waals equation as the basis of this work was due to the fact that it
is the simplest and most easily adapted equation-of-state with a theoretical basis. The ob-
jective of this work is to demonstrate that the proposed approach provides qualitative
predictions with the van der Waals equation, and lay down the framework required to use
a more accurate equation-of-state. A comparison of Figures 2-2 and 2-10 shows that the
SLD model exhibits poor prediction of the vapor pressure for subcritical isotherms, but -
does better at predicting the pressure of the maxima in the supercritical isotherms. This is

because the van der Waals parameters a and b were obtained from the critical temperature

57
and pressure. The 8/, in Table 2-1 have been selected to provide semiquantitative fit to the
magnitude of experimental adsorption.

Cubic equations are widely used in industry, and a method that adapts them to the
adsorption problem could find widespread use in process calculations. In the next chapter
we will see how changing the fluid equation-of-state significantly improves the quantita-
tive modeling of adsorption. We will also show the extension of this model to predict

Type I, IV and V isotherms, and multicomponent mixtures.

CHAPTER 3: QUANTITATIVE MODELING OF ADSORP-
TION AND EXTENSION OF THE SLD APPROACH TO
PREDICT CLUSTERING

INTRODUCTION

In the previous chapter, we presented an engineering model that adapts the van der
Waals equation-of-state to describe the fascinating behavior seen in the experiments of
Findenegg [1983]. The fluid-solid potential was superimposed on the van der Waals
equation-of-state, and the configurational energy integral in the inhomogeneous fluid
phase is simplified with a local density approximation. While this model does a good job
of qualitatively describing the adsorption behavior, quantitatively it does not predict the
isotherms very well. The primary reason for this shortcoming was that the van der Waals
equation does not predict accurately the vapor pressure and density of the adsorbing fluid.
In this paper we use the Peng-Robinson equation to describe the fluid properties, which
significantly improves the model predictions relative to the experimental results.

An adsorption model may also be adapted to describe clustering in supercritical
fluids [Lee et. a1, 1991]. Clustering is a phenomena that occurs in supercritical fluids
whereby a large number of solvent molecules collapse around each solute molecule
[Debenedetti, 1987; Eckert et. a1, 1986], and the effect is most pronounced near the criti-
cal point. One of the first experimental studies on clustering involved the measurement of
partial molal volumes near the critical point for solvents such as carbon dioxide and ethyl-
ene around solutes such as naphthalene [Eckert et. al., 1986], demonstrating that the par-

tial molal volumes become extremely negative near the critical point. Debenedetti [1987],

58

59

using a fluctuation analysis converted these negative partial molal volumes to number of
solvent molecules that cluster around the solute. Kim and Johnston [1987 a, b] state that
at 300.5 K and 79.8 bar, the partial molal volume of naphthalene in supercritical carbon
dioxide at infinite dilution is -7800 cc/mol, which corresponds to the condensation of 80
solvent molecules around a solute molecule. In a related study, Johnston et. al. [1987]
show that there is a shift in the solvatochromic data of phenol blue in ethylene indicating
an increase in the number of solvent molecules around a solute. Wu et. al. [1990] provide
integral equation calculations on other supercritical systems that also suggest that cluster-
ing does occur in supercritical mixtures.

Brennecke and Eckert [1989] measured the intensity ratio of the fluorescence
spectra to probe the local density of dilute organics in pure supercritical fluids. The in-
tensity ratio for naphthalene in CO; at 4 K above the critical point was much greater than
at 14 K above the critical point, indicating the clustering effect near the critical point.
Such an increase was also demonstrated for systems such as pyrene-CO; and pyrene-
ethylene. Lee et. a1. [1991] have reviewed both the integral equation and spectroscopic
measurements that give evidence of clustering in supercritical fluids.

Clustering is used to explain the unusual behavior in supercritical fluids such as
increased solubilities, synergistic effects of mixed solutes and entrainer effects. The forces
acting on the fluid molecule which result in this behavior are the fluid-fluid and fluid-solute
intermolecular forces. On a molecular level, there is a similarity to the fluid-fluid and
fluid-adsorbent intermolecular forces in physical adsorption on solid surfaces. In cluster-

ing, this leads to an increase in density of the fluid molecules in the region around the sol-

60

ute similar to the increase in density around the adsorbent in physical adsorption, and the
fluid can be modeled in the inhomogeneous system.

In this chapter, the results of this simplified local density (SLD) model using the
Peng-Robinson equation-of-state are given. The basic concepts of the SLD approach
have been discussed in the previous chapter. The modifications made to this model to in-
corporate the Peng-Robinson equation, and model predictions of experimental results for
adsorption on flat walls are discussed in detail in this chapter. The modifications made to
this model to predict the clustering phenomenon are also discussed in detail. The model is
capable of predicting the density profile and number of solvent molecules around a solute
to indicate long-ranged interactions, giving evidence of clustering around a solute near the
critical point. The model is shown to correlate with fluorescence spectroscopy measure-

ments.

MODEL DEVELOPMENT

The approach of the model is the adaptation of a cubic equation-of-state to define
the properties of a fluid that is in the external potential field of an adsorbent. First, we de-
velop the model for flat solid adsorbents (incorporating the Peng-Robinson equation) and
then generalize the results for clustering. A fluid molecule at any position from the solid
interacts with both the solid and fluid molecules. Therefore, the chemical potential of the

fluid may be written as a sum of the fluid-fluid and fluid-solid interactions

Il(z)=llb =ug(2)+u,.(2) [3-11

61

where ttfflz) and 11st z) are the chemical potentials due to the fluid-fluid and fluid-solid

interactions, and z is the distance normal to the surface of the solid, and the subscript b
stands for bulk properties. If the fluid-solid potential is given by 111(2) on a molecular ba-

sis, then on a molar basis

Hf, (Z) = N A‘I’(Z) [3'2]

For a non-ideal bulk fluid, the chemical potential jib of Equation 3-1 is related to

the fugacity as

u. = #0 + R T ln(-£42] [3-3]

0

Details of the development of this model are given in the previous chapter.
Using the Peng-Robinson equation-of-state, can write an expression for the bulk

fllgacity in terms of molar volume.

ln(-1%):(Zb - 1) - ln(Zb - B) -

 

 

A, ln(z, +2414 B] {3.4}

2.828 B z, - 0.414 B

where

62

 

P . . ..
Z, = RV; IS the compressrblllty factor

 

 

a P
Ah R2bT2
3:21:
RT
P= RT ab

vb -b_ vb2 +2bvb -b2

where a, b are the constants of the Peng-Robinson equation [1976].
Rearranging Equations 3-2 - 3-4, the component of the fluid fugacity due to fluid-

fluid interaction is

‘1’( )
f 17 (z) = f), “PL—[TL] [3-5]

Using the PEng-Robinson equation, the fluid-fluid fugacity at any point becomes

_ b _ 0(2) V(z)
lnlffld] — v(z)—b v(z)2 +2bv(z)—b2 RT

- 1,1[M] _ 0(2) 1nl:v(z)+2.414 b]

 

 

13-6]

 

RT 2.828bRT v(z) - 0.414 b

63

where v( z) is the molar volume at any point z. a(z) is calculated using the SLD approach,

as given in the previous chapter, assuming a fluid-fluid pair potential oc f6 and is given by

 

a(z) = a, —5— + ii for 0.5 s —z— s 1.5 [3-7]
16 166, of
1 z
a(z) = ab 1- for 1.5 S — S 00 [3-8]
8(—Z— — 0.5)3 0 x
G

where of is the diameter of the fluid molecule.

In modifying this model to arrive at Equations 3-6 - 3-8 for the Peng-Robinson
equation assumptions are necessary to evaluate the configurational integral. Returning to

the integral from which the equation-of state a parameter is based, the attractive chemical

potential is given by

u... = jurfifldv 13-91

V

Where V denoted the volume of integration, (p(r) is the two-body interaction potential, and

8( r) is the radial distribution function.

64

In the SLD model, the fluid at a point z is treated as a homogeneous fluid at a

density of p(z). Solving Equation 3-9, the attractive pressure term of the van der Waals

equation for a homogeneous fluid becomes -a/v2'. In the case of the Peng-Robinson equa-

a
v2+2bv—b2°

 

tion, the attractive pressure term for a homogeneous fluid is — The Peng-

Robinson equation, as an empirical equation-of-state, has an unknown equivalent form of
the integral given by Equation 3-9, and therefore some assumptions are necessary to ap-
proximate this integral. The term v2 in the denominator of the attractive part of the van
der Waals equation is a result of the volume derivative of Equation 3-9 . The term v2 +
2bv- b2 in the Peng-Robinson equation may be treated as an empirical equivalent of the
term v2 of the van der Waals equation, and can be subject to an analogous SLD approxi-

mation. If the density is uniform up to the surface of the adsorbent, then we know that at

this limit 11qu = rib/2, resulting in a..." = al./2 for the van der Waals equation. The Peng-
Robinson equation should have the same limits. Also, the functionality of the a/ab should
have the same qualitative curvature in both the van der Waals and Peng-Robinson equa-
tions. Since the true form of Equation 3-9 is not known for the Peng-Robinson equation,
as a first approximation, we use the van der Waals ratio of a(z):ab for the Pong-Robinson
SLD model, and retain v2 + 2bv- b2 as the denominator of the attractive pressure term,
Which results in Equations 3-6 - 3-8.
The other pair potential needed for this model is the fluid-solid potential. The ad-
Sol‘bate-adsorbent interactions are given by the integrated 10-4 Lennard-J ones potential '

[Lee, 1988] as

65

0‘}. 1 4 1
‘I’(z)=47rpa,m8fi0'fi EFT-2'27 [3-10]
1 i=1 i

For modeling clustering in supercritical fluids, the approach is essentially similar,
with the exception of the configurational integral calculation. The wall in adsorption onto
a flat surface is replaced by a solute molecule with the solvent molecules clustering around
it. Geometrically, this means replacing a planar surface with a spherical particle. The co-
ordinate system is changed from rectangular geometry to spherical geometry and the fu-
gacity, molar volume and the Peng-Robinson a term become functions of r, the radial dis-

tance from the center of the solute molecule. The corresponding a( r) become

 

 

'8 2 , 2 , 1 ‘
3+;(r-r “30’ ( + ’)3
r r
—a(r) =3 ’ for 0.50, s r-91 s 150, [3-11]
+— —
_ r {7’2 +r2 —277‘” r’2 +r2 +2rr’}_

 

 

P16 26) 1 203, 1
—- +
a(r)_i 3 3 (r—r’f 3 (r+r')’

ab 16 +03{ 1 l }
r r’2+r2—2rr’ r’2+r2+2rr’_

 

 

 

for 1.50, s r--62—‘ s cop-12]

 

 

 

 

Where

_(o,+o,)

I

66

r’2 -l-r2 —O'2
n_ f

r — I
2r

 

where the subscript f stands for fluid, and the subscripts for the solute. Details are given

in the appendix 1.

The solute-solvent interactions are given by the Lennard-Jones potential

0.12 0.6
‘I’(r) =4e,,(-;’2i--r7" [3-13]

The component of the fluid-fluid fugacity due to fluid-fluid interactions, i.e. the

equivalent form of Equation 3-5 for clustering is

 

‘i’( )
fflm = fbexp[- k; ] [344]

Using the Peng-Robinson equation-of-state, the fluid-fluid fugacity can be expressed

Similar to Equation 3-6, except that the molar volume, a and fugacity are functions of r

instead of 2.

Algorithm for solving density profile

To calculate the surface excess or cluster size, the density profile must be known.

The density can be evaluated by first calculating the bulk fugacity and density from a given

67

temperature and pressure. Knowing the fluid-solid potential, the fugacity at any point can
' be calculated by Equation 3-5 (physical adsorption) or 3-14 (clustering). The density is

then calculated iteratively using Equation 6 or by the modified form using v( r) and a( r) for
clustering. Once the density profile is known the surface excess can be calculated by equa-

tion 2-1. Similarly, for clustering, the cluster size (Nu) can be calculated as
N" = f(p(r)— pb)dV [3-15]

The FORTRAN programs for calculating the adsorption isotherms and cluster sizes are

available in appendix 2.

RESULTS
The diameters of the fluid, solid and solute are the Lennard-Jones parameters as
tabulated in Reid et. al [1987] and in Lee et. al. [1991]. The Of; is calculated as the arith-

metic mean of the fluid and solid (or solute in clustering), while the fluid-solid potential

(83) is used as an adjustable parameter to adjust the magnitude of adsorption. The values
Used for all the parameters are summarized in Table 2-1 or 3-1.

Figure 3-1 shows the surface excess of ethylene adsorbing onto graphon at various
Sub— and supercritical temperatures, and over a wide pressure range. This model is capa-
ble of quantitative fits for the adsorption isotherms over this entire pressure and tempera-
ture range, using a single parameter. The model also does a good job of accurately pre-

dicting the pressure at which the peaks occur, showing the characteristic cusp-like behav-

Table 3-1: Molecular properties used in the SLD model

68

 

 

 

System odnm) 0, (nm) Cg (nm) EfS /k (K)

Ethylene - Graphon .422 .34 .381 140
Krypton - Graphon .3655 .34 .35275 95
Propane - Graphon .51 18 .34 .4259 110
Argon - Graphon .3542 .34 .3471 68
Methane - Graphon .3758 .34 .3579 107
Carbon dioxide - Naphthalene .3794 .6199 .49965 800
Carbon dioxide - Pyrene .3794 .714 .5467 850

Ethylene - Pyrene .422 .714 .568 800

 

 

Excess (It mollmz)

69

 

 

 

 

 

 

     
  

 

 

Pressure (bar)

Figure 3-1: Adsorption of Ethylene on
Graphon

 

I20

45
1:1 1:) 263.15 K
40 __ A A 273.15 K
A o 283.15 K
0 293.15 K
35 “r A
g - 313.15 K
o
A o
30 " A °
1:1 A o
l} o
25 "" I I o;
I
20 -- 9! °.-"
4’ 8" 03968
15 - I f" 9”"\
8 03’ o ‘
O ’6’ *‘g —!-‘~
10 - 1,..-"""o \x “x-
’ 90 “ ‘\..~ “‘~-.
°o o -~a::”f~-—_
5 ‘ .7 ° 0 o O” i i
0 u I I I I I
(I 20 40 60 80 HI)

70

ior seen near the fluid critical point. The model predicts the correct temperature-
dependent crossover of adsorption isotherms in the supercritical region. One of the pri-
mary reasons that the Peng-Robinson SLD model does a better job than the van der Waals
SLD model is that the Peng-Robinson equation predicts more accurate density values near
the solid surface, and a more compressible adsorbed layer. A comparison of experimental
saturation densities for ethylene [Angus et. a1, 1974], and both these equations [Pyada,
1994] shows that while the van der Waals equation underpredicts the saturation densities
by 80-90% between 211 and 260 K, the Peng-Robinson equation is within about 6%. van
Megen and Snook [1982] have done Monte-Carlo simulations of Lennard- Jones mole-
cules (that simulate ethylene) on graphon at supercritical temperatures [8] (T IT c = 1.03)
(see Figure 2-5). Their results significantly overpredict the experimental values of Fin-
denegg, and their maxima is around 42 bar, as opposed to the experimental value of 58.7
bar, while the SLD model has a maxima at 57.8 bar.

Figures 3-2 and 3-3 Show the adsorption isotherms of krypton onto graphon at
supercritical temperatures, and propane at subcritical temperatures. The krypton iso-
therms 253 and 273 K crossover at around 140 bar. Once again, the model predicts this
crossover at the correct pressure. This model was also tested on gases like argon and
methane (Figures 3-4 and 3-5). While, the model predictions for argon were good, the
model did not do a very good job of predicting the isotherms for methane, the reason
maybe that the Peng-Robinson equation does not do a good job of predicting the fluid
properties of liquid methane [Prausnitz, 1980]. The PVT behavior of these different gaSes

and the error in predicting density values is shown in appendix 3.

H

 

 

 

 

 

 

 

 

 

14
a D
:1 1:1
12 "' n D
a n
o o 0 o
10 -- ° ° -
01‘ o ’, """"""" ~-‘o
E n o ’z” a
= x . .
O 8 -- l/’ A ........ 0
g. D °/ A .f L--—""‘" __-,_.....-._.
V II 9 ,I”’ " .— "- .-.‘_.._ _.
m 9’ .. . ' ". ‘ ‘—
8 6 ‘T l I 4' .,( _,r '3 ‘ A
O I l.’. ,': ° A ‘ ‘
I '0 ./ .o ‘
4 a- 79' / '.’ ‘ ‘ 13253.15K
..' . ,
"Ii/’3. o ‘ 0 273.15 K
'a'/k . ‘ g A 298.15 K
2 -:'.' ‘ l323.15K
:5"; 0348.15K
t 1373.15K
0 I I I
0 40 80 120 160

Pressure (bar)

Figure 3-2: Adsorption of Krypton on Graphon

72

 

   

 

 

 

 

 

 

 

 

30
i .
I .
I I
25 -- I i
I I
l "O
l .
l 18
"A 20 -— '1 I;
g ’ a
— 1:1 I 8.
O I I
E i ,0
.3- 15 I’ o"
v) I 00/
8 /'
O 99/
X W'
Lu 10 a 'wo’
A“! ,0"? A
“0’99" 1 13273.15 K
{ 0323.15 K
41343.15 K
0363.15 K
20 30

 

Pressure (bar)

Figure 3-3: Adsorption of Propane on Graphon

4O

Excess (u mol/mz)

73

 

 

 

 

 

 

 

 

 

9
8 --
7 .-
6 -—
5 ..
4 --
3 -- I
.' Fl .
"7' 0.; n253.15K
.. ' O
2 . /'/ «1273.15 K
1",” 4298.15K
0 323.15 K
O -- I I I
0 40 80 120

Pressure (bar)

Figure 3-4: Adsorption of Argon on Graphon

I4

74

 

12~

IO-

Excess (It mollmz)

 

 

 

 

a 253.15 K
o 273.15 K
A 298.15 K
0 323.15 K

 

 

 

l I l
I l l

80 l 20
Pressure (bar)

Figure 3-5 : Adsorption of Methane on
Graphon

 

160

75

Figures 3-6 - 3-12 describe the clustering phenomenon. The density profile and
the excess number of solvent molecules within a sphere of radius R around a solute are
plotted in Figure 3-6. As shown in the figure, the density goes through a sharp increase
near the surface of the solute but quickly drops to the bulk value. The excess number of
solvent molecules is slightly negative very near the surface before increasing sharply. Such
a trend has been shown in the integral equation calculations of xenon in neon [Wu et. al,
1990] (Figure 3-7). The figure also shows that the excess number of solvent molecules
increases sharply up to 5 molecular diameters, and then slowly till 9 or 10 molecular di-
ameters, showing the long-ranged effect of this clustering phenomenon. This trend has
also been reported by Wu et. al [1990].

Figures 3-8 & 3-9 shows the results for clustering of C02 molecules around naph-
thalene at 2 K, 4 K and 14 K above its critical point. As seen in Figure 3-8, the cluster
size increases sharply near the critical point, but drops quickly as we move away from the
critical point. The results of the fluctuation calculation of Debenedetti [1987] on the same
system at 308 K are also plotted in Figure 3-8. The model predicts the correct density at
which maximum number of C02 molecules cluster around naphthalene, and does a fair job
of predicting the number of CO; molecules. The model also predicts a slightly broader
density range over which this clustering phenomenon occurs. However, the model does
accurately predict the increase in the number of C02 molecules around a single naphtha-
lene molecule as the critical point is approached, and also the variation of this peak with

temperature and pressure. When this cluster size is plotted against bulk density, cross-

Nex

76

 

80 30000

 

L-25000

l

L

r200m

L

~15000

«10000

 

 

Nex

' - " - P

-10 4. I I 0
0.00 5.00 10.00 15.00 20.00

Rlo

Figure 36: Density and Excess number of
solvent molecules within a Sphere of Radius R
around a Solute

 

 

 

 

 

 

 

(sumowfi) d

77

 

a...
.......
0'
.o
o
o
.o

701

10‘

I [g,2(r) - l] ‘Irtr2 dr

((511)3 0

L

l

 

O

 

 

 

I
H
O

13 Is
L10“

 

 

C
U11

F1 QUre 3-7: Number of Excess Solvent Molecules Within

a Sphere of Radius L Around a Solute
A: Xenon in Neon (Attractive)
F1: Neon in Xeonon (Repulsive)
P: Pure Lennard-Jones

78

 

100

 

 

D - - - 306.2 K
303.33 K

 

ll
a 308.38 K,inferred from expt.

 

 

 

ll

H

I: ------ 318.2K
I

 

ul-
-

 

 

1-

80 100 120 140
Pressure (bar) -

Figure 3-8: Clustering of CO2 around
Naphthalene

160

79

 

 

.' - - - 318.2 K
90 __ ------ 306.2 K

-' ‘ a 308.38 K, inferred from expt.
308.38 K

 

 

 

Nex
8

 

 

 

 

0 : l l n

4000 8000 12000 16000 20000

Density (gmollm3)

Figure 3-9: Clustering of C02 around
Naphthalene

Localized Cluster Size

80

 

 

 

 

 

 

 

 

 

0.38
58 --
4:- 0.36
54 -—
.. 034
w --
A .. 0.32
A
A
AA
46 " -—303.2 K ._ 03
------ 313.2 K
a Data - 308.2 K
A Data - 313.2 K
42 i i i 023
4000 8000 12000 160m 20000
Density (gmollma)

Figure 3-10: Comparison of Cluster Size with

Fluorescence Spectroscopy - CO2 around

Naphthalene

VII I-I

Localized Cluster Size

88

81

1.1

 

 

 

 

 

 

 

 

 

32 j- a
.T '|
76 T
-- 0.9
70 --
64 T “ 0.8
58 uh
-- 0.7
52 --
-.—303.2 K
------ 323.2 K L' 0-6
46 " a Data 303.2 K
. A Data 323.2 K
40 J : : : 0.5
0 5000 100(1) 15000 20000

Density (gmollm3)

Figure 3-1 1: Comparison of Cluster Size with
Fluorescence Spectroscopy - CO, around

Pyrene

Elll-l

Localized Cluster Size

82

 

 

 

 

 

 

 

 

 

72 0.8
68 --
-- 0.75
64 ....
60 -- 0.7
56 -- . " 0-65
D ,5. '
52 -- °
a -- 0.3
43 -- ° a
—234.2 K
------ 303.2 K “ 0-55
44 "" a Data - 234.2 K
a Data - 303.2 K
40 1 ' i 1 i 0.5

0 3000 6000 9000 1 2000 1 5&1)
Density (gmollma)

Figure 3-12: Comparison of Cluster Size with
Fluorescence Spectroscopy - Ethylene around

Pyrene

Ell l-l

83

overs of Figure 3-8 disappear, and the maxima appear at the same bulk density (Figure 3-
9).
The fluorescence intensity ratios are an indication of the total number of molecules

around the solute within a short range, and hence the total amount within the first two

. o. o, .
molecular diameters [—'— S r S -—‘ + 20 f] outsrde the solute was chosen as a measure of

2

the local density. Figures 3-10 - 3-12 show the total number of molecules around a solute
within the first two molecular diameters against bulk density. This is cross-plotted against
the fluorescence signal of Brennecke [1989], Brennecke and Eckert [1989] and Brennecke
et. al. [1990]. The increase in the intensity of the signal near the critical point suggests
that the solvent molecules cluster around the solute. For naphthalene, the stable signal is
the ratio between fluorescence peaks 1 and 4, while for pyrene it is 1 and 3. These figures
show the predicted total number of C02 molecules clustering around naphthalene and py-
rene, and ethylene clustering around pyrene. The model predicts the trends seen in the
fluorescence spectra. For CD; clustering around pyrene and naphthalene at 308.2 K, the
data show a sharp increase at low densities, then becomes rather flat before increasing
again at high densities. Also, at higher temperatures, the flat region is rather small. Once
again this trend is captured by the model. This difference is the trend of the fluorescence
spectra is even more pronounced for ethylene clustering around pyrene. The figures also
show that the intensity of the ratio of the fluorescence spectra and the excess number of
C02 molecules around pyrene is more than around naphthalene. The fluid-solvent inter—
action parameter is higher for pyrene than for naphthalene, and pyrene being a larger

molecule provides a larger surface area for the solvent molecules to collapse upon. Com-

84

paring the C02-pyrene and ethylene-pyrene systems, pyrene attracts a larger number of
C02 molecules than ethylene. The reason may be because of the bigger size of the ethyl-
ene molecule which sterically limits the number of ethylene molecules that can surround
the solute. In all the cases presented here, the model predicts greater clustering near the

critical point which is in agreement with the spectroscopic data.

DISCUSSION

The SLD model uses the fluid-fluid and fluid-solid interaction parameters to model
adsorption over wide pressure and temperature ranges. It characterizes the fluid mole-
cules using cubic equations-of-state. Predictions of the SLD model using the van der
Waals model were shown earlier [Rangarajan et. al, 1995, chapter 2]. Here, we use the
Peng-Robinson equation to describe the fluid, thus significantly improving the predictions
of the fluid properties. By comparing Figures 3-1 - 3-2 with the results seen in chapter 2,
the significant improvement is obvious. However, for any particular system, the equation-
of-state that does the best job for representing its properties may be chosen, and the SLD
approach adapted for that purpose.

The fluid-wall potential is used as an adjustable parameter in this approach. It is
fitted to any one of the different temperatures given in the data, and then the same poten-
tial is retained for the other temperatures. This makes the model a semi—predictive model,
in that, while one of its parameters needs to be fitted to experimental data at a particular
temperature, it then becomes a predictive model to describe the adsorption characteristics

at other temperatures. The geometric mean of the Lennard-J ones potentials tends to un-

85

derpredict the values of the surface excess. We have consistently used a fluid-solid inter-
action parameter that is about 1.8 to 2.2 times the geometric mean.

The SLD model superimposes the Lennard-Jones potential for solute-solvent in-
teractions on the Peng-Robinson equation of state that describes the properties of the su-
percritical solvent. If the solute is attractive, then the density of solvent molecules near
the solute increases resulting in the phenomenon of clustering. This is very similar to the
adsorption of a fluid on to a flat wall, except that the wall-fluid interactions are slightly
different than the solute-solvent interactions. As shown by the results of both these two
cases, there are many similar behaviors between the two cases. The crossover of the ad-
sorption isotherms is also seen in the case of clustering (Figure 3-8). This crossover dis-
appears when the excess is plotted against bulk density (Figure 3-9). The peak in the su-
percritical adsorption isotherms decreases as we move to higher reduced temperatures.
This phenomenon is also seen in clustering, where the number of excess solvent molecules
drops sharply as the temperature is increased from 2 to 4 K above the critical point. One
of the differences is that below the critical pressure, the surface excess of a fluid adsorbing
on a flat wall increases gradually and then drops sharply as seen at 283 and 293 K for eth-
ylene (Figure3-l). However, in clustering, as the fluctuation calculation of partial molal
Volume data indicates, the number of molecules increases sharply as the critical pressure is
approached and decreases rapidly past the critical pressure. This model’s predictions for
Clustering are more in line with the adsorption data on a flat wall i.e., the excess number of
Solvent molecules increases more gradually below the critical pressure, and then drops

Sharply after the critical pressure. This may be an artifact of the SLD assumptions.

[it
he

H0

D0

T10

[11.

if:

fen

10]

ill

86

Wu et. al. [1990] showed that for the attractive mixture of neon clustering around
xenon, there is a negative N“ close to the surface of the xenon molecule (Figure 3-7).
The SLD model also predicts such a negative value close to the surface. At the surface,
the attractive and repulsive parts of the Lennard-J ones potential cancel, and the solvent
molecules are not attracted to the solute, and this negative value is the result of the sol-
vent-solvent interactions. By choosing the correct solute-solvent potential parameter, this
model may be used to describe a repulsive mixture where there is depletion of solvent
molecules e. g. xenon molecules cavitating around neon.

The SLD model serves as a good first approximation of the adsorption problem.
The use of the Peng-Robinson equation helps resolve one of the limitations of the previous
SLD model viz. use of an equation that can better represent the fluid properties, thus sig-
nificantly improving the capabilities of this model.

High-pressure adsorption and clustering in supercritical fluids are phenomena that
are of importance in various industries. A simple engineering model to describe these dif-
ferent phenomena will enhance their applications in industry. Modifications of this model
to predict adsorption in slits and pores, and for describing adsorption of mixtures onto

different surfaces are discussed in the next chapter.

{1

CHAPTER 4 ADSORPTION'OF PURE GASES IN SLITS
AND PORES, AND ADSORPTION OF BINARY MIXTURES

INTRODUCTION

in two previous chapters, we presented a new model, called the Simplified Local
Density (SLD) model, for describing physical adsorption of gases on flat walls using cubic
equations-of-state. In this model, the fluid-solid potential is superimposed on a cubic
equation-of-state, and the configurational energy in the inhomogeneous fluid phase is
simplified with a local density approximation. In the first chapter, we used the van der
Waals equation, which gave qualitative predictions of the experimental data of Findenegg
[1984]. The van der Waals equation was used because it is the simplest fundamentally
sound equation-of-state and provided a basis for qualitative studies. This model predicts
Type II and III isotherms for adsorption on flat walls, shows the characteristic cusp-like
behavior and crossovers seen near the critical point. In the second chapter, we modified
the model to incorporate the Peng-Robinson equation-of-state, which represents the fluid
properties more accurately, including vapor pressure and compressibility, and showed
dramatic improvement in the quantitative representation of experimental measurements of
Ffindenegg.

Most commercial adsorbents are porous, activated carbon being one of them. Ad-
sorption in pores and the phenomenon of pore filling among carbon based adsorbentshas
received considerable attention [Dubinin, 1975]. Gregg and Sing [1982] have classified

pores into three categories, macropores (50 - 100 nm), mesopores (2 - 50 nm) and micro-

87

88

pores (<2 nm). Adsorption properties of microporous solids enjoy wide spectrum of ap-
plications including chromatography and other separation processes, filtration, industrial
effluent cleanup, ion exchange, biological applications etc. [Gubbins, 1990; Pendleton and
Zeetllemoyer, 1984]. In microporous solids, adsorbed fluids show many unusual proper-
ties such as preferential adsorption of certain fluid components, chemisorption at particu-
lar sites, hysterisis effects and a variety of new and unusual phase transitions [Gubbins,

1 990]. Dubinin [1975] says that all existing theories of physical adsorption proceed with
the same physical image in describing adsorption in porous and nonporous adsorbents.

Design of columns for thermal swing adsorption and pressure swing adsorption

require the simultaneous solution of partial differential equations for the material, energy
and momentum balances describing the dynamics of adsorption in columns, in conjunction
with the kinetic and equilibrium pr0perties of the adsorbents [Sircar and Myers, 1985].
The models for multicomponent adsorption isotherms are discussed in chapter 2. There
are plenty of experimental and molecular simulation data available for adsorption of pure
gases and mixtures onto slits and pores at low pressures [van Megen and Snook, 1982,

l 985; Barton et al., 1984; Powles et al., 1988; Valenzuela and Myers, 1989; Tan and
Gubbins, 1990; Kieer and Rosinberg, 1992]. In this paper, we extend the SLD model to
describe adsorption of pure gases and binary mixtures in confined spaces, i.e. slits - two

infinite parallel flat surfaces and cylindrical pores.

89

MODEL DEVELOPMENT - PURE GAS
Slit-like Pores

These can be thought of as two parallel semi-infmite plates separated by a small
distance (Figure 4-l). Any particle in between the walls will be subjected to attractive
forces from the two walls. The distance between the two walls may be quite small (<2
nm), which implies that only particles smaller than this will be adsorbed i.e. exclusion will
be important.

The chemical potential of the fluid in the vicinity of a wall may be written as a sum
of the fluid-fluid and fluid-solid interactions, as described in the previous two chapters.

Recapping from the last chapter, using the Peng-Robinson equation-of-state, we
can write an expression for the bulk fugacity in terms of molar volume [Sandler, 1989].

The fluid-fluid fugacity at any point is an analogous expression

b _ a(z) V(z) _
v(z) — b v(z)2 + 2bv(z) - 1)2 RT

ln[v(z)-b] _ a(z) ln[v(z)+2.4l4 b]
RT 2.828bRT v(z)-0.414 b

 

1“”.ng =
[4'1]

where v(z) is the molar volume at any point z, a and b are the constants of the Peng-
Robinson equation [Sandler, 1989; see previous chapter also]. a(z) is calculated using the
SLD approach after suitable modifications, and details are given in appendix 1. The equa-

tions for a(z) in slits are

90

plane of symmetry

9(2)

p:

 

 

 

 

 

 

 

 

 

 

 

Potential

Energy

Slit-Like Pores

1

Figure 4-

'he

ten

91

 

 

 

 

 

 

 

 

 

 

 

 

“(0:2 _Z_+.5__ 1 3 , for assists [4-2]
ab 8 017 6 L-Z 0’7
3 ———0.5
a(Z)_§ §_ 1 _ l
a -8 3 3 _ 3 ,
b 3[__Z__05] {L 2—05]
or of a
[4-31
for 153-3—3—5—15
o, 617
0(2):; L-z+%_ 1 3 , for BL_1.5_<_O_Z_gal;.—0,5 [4-4]
a 0'
” 1’ 3[—5——0.5J ” f H
0'17 -

 

 

Wher e L+o, is the distance between the centers of the atoms of the two surfaces. Note

that L has to be at least 3 Gfl' in length for the above equations to be used. If L < 3 Of]

the“. the equations will have to be slightly modified and the corresponding equations are

 

0'
“(3:2 —L —l , for 0.SS—z-Sl.5,andifz+0'fl >L-'—fl [4’5]
ab 8 of 2

92

 

 

 

 

' 1
“(0:2 2 _ 6} +—5_
a 8 O' 0' 3 6 ’
b H {Ir—21’ ) [4-6]
0'
for 0.5s-Z—s15 andif “of, 3L...__1?'_.
or 2

a(z)___3 L-Z 0f]

a 8 0' 0' 3
” 1" 3(z-—”-(L—z))

 

 

0” [4-7]

, for 1.53-5—sL———
o 2

 

 

The adsorbate-adsorbent interactions in slits are given by the partially-integrated
10-4 Lennard—Jones potential [Lee, 1987] (Equation 3-10), where the interactions are

summed over the two walls.

Cylindrical Pores
Any molecule in a cylindrical pore is subject to a force from all sides of this cylin-

der (Figure 4-2). The configuration integral in pores is given by

 

)=}-°}i”(‘,‘2—r‘— +22) drd0dz [43]

where r and z are the radial and axial distances. Using the SLD approach, the region in-
side the pore is divided in three parts: 1) at the pore surface; 2) near the pore surface; and

3) far from the pore surface, and the details of these calculations are given in appendix 1.

 

93

Axis of symmetry

&\ R

i
l‘ p=pa

 

 

 

l.’
I

. U‘
‘1 ‘in”. .‘0
"V? ' “1":

y." .t" ‘ 1"
his,» h .t
=5". 3.. ~
! v ‘- "J.
K’u-uz, ,‘ ,. .-a~
. .l 't ’>_
J' "

at;
atavv

~e.u

tr
\

 

1
)’

 

-: al'f ‘15 ’ .2.

—- A v 1
'03 :31 r
‘5: ” "
7“ "ii“?

5.9
4 :
at

 

 

Figure 4-2: Cylindrical Pores

94

The resulting solution leaves us with an integral to be solved for a(rl):ab, and it is the

longest part in the calculation of the amount adsorbed. The results can be summarized as

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

follows:
— 62 _Zz a. 2 2+ 2 q
3 r” 14 cos'l ———fl dz- I 1t— r, Z 3,,
a(fi)=3°r ° 4617 2’: "1’2z(4rI +2)
0.. 7‘ 1 , ... . 1 n ’
—— d0dz+—
L 41:]0 4132c0520+z2 246},r
[4-9]
G
orr-R—-—fl-
f l 2
j” cos-'(v)dz-—j°’ j°°“( 1 d0dz
3 ° 49} (w 222+2)
“(0:303
ab 1t —-l-r 1; - 0’5 -tan"(ofl./w) (19+ 1;
b 4 0 4w3 2w2(w2+o}) 2W3 120}
[4-10]
for R— ”Sr,< —-——-fi—
363 3O
“(4): 1’ “3 —1 —d0, for osn<R—-—”— [4-11]
a, 11: 30'” 160w3 2
2
0'
r12+0'i,r Zz-(R-EL] O’ 2
Where V: , and w=r,c089+\[(R——H) —rlzsin20-,R+
2r,‘/O'§,- 2 2

0J2 is the pore radius, and r, 0, and z are the cylindrical coordinates.

95

If the fixed particle is near or at the wall, then the configurational energy calcula-
tion at each point takes about 45 seconds for the fortran calculation. If the fixed particle is
away from the wall, the calculation takes about one second to compute the configurational
integral. Note that the calculation has to be done only once for a given z/of. Similar to the

30'
slits, if the pore radius is smaller than 3 of; + om then rI > R—ji- in all cases and

Equation 4-11 is not required.

In cylindrical pores the adsorbate-adsorbent interactions are given by the potential

of Tjatopolous et al. [1988]

 

_ 2
‘I’(x,Rfi)—n0;fl EM 1 [4-12]

 

 

 

where x is the distance of closest approach between the fluid molecule and the pore sur-
face, Rf; = R-t-o,/2 is the pore radius, F [ a, B; Y: 8] is the hypergeometric series with pa-

rameters 0t, [3, y, and n is the number density of carbon atoms.

Algorithm for solving the density profile

The algorithm for solving the density profile and calculation of surface excess are

given in the previous chapters. In slits and pores the only difference is the calculation of

96

a/ab, which is given by Equations 4-2 - 4-7 for slits and 4-8 - 4-11 for pores. The amount

adsorbed in slits is calculated as
L—ofl/Z
Amt =j p(z)(S. A. )dz [4-13]
non/2

and in pores as

 

 

Amt = LR-” 27:7Lp(r) dr = IOR-ak r S”: p(r) dr [4-141
R _ S
2

where S.A. is the experimental surface area (for e.g. BPL activated carbon has a surface
area of 988 m2/g) and L is the length of the cylinder. It was assumed that the experimental

surface area was 21tL(R-o,/2).
MODEL DEVELOPMENT - BINARY MIXTURES

In binary mixtures, the van der Waals mixing rules are used to determine the prop-

erties of the mixture [Prausnitz et al., 1986]

M
i=1 j=l

a=izyiyjaij [4-15]

b = 2m 14-161

whe

COP. \

 

Wher

97

a,j =aj, =,/a,,afl(1—k,,) {4.17}

where y,- is the mole fraction of component i, and k,-,- is the binary interaction parameter, a
constant introduced to obtain better agreement in mixture equation-of-state calculations.

Equation 4-1 may be written for the fugacity of each component as

22 y .14.,-
,-

3, .-
lnty—p]=-;<z-I>-m<z-B>-m. .. 7

 

 

Where A, B and z are as discussed in chapter 3.

in =1 [4-19]
i=1

Algorithm to solve density profile

For a binary mixture, the above equations result in solving two simultaneous equations
(fugacity of each component) for both composition and density with position, subject to
the constraint that the sum of the mole fractions of the two components is 1. A technique
developed by Asselineau et al. [1979] for calculating VLE data has been adapted for this
purpose. In this method, Equations 4-17 - 4-19 are written such that the right hand side

equals zero. These equations are called the objective equations, G], G; and 6;.

GI =l'9'—P—1 [4-20]
fl

98

P
(32 =flL—1 [421]
f2
G3 = l-yl "Y2 [4’22]

Partial derivatives with respect to the independent variables (pressure and local composi-

tion) are then written in matrix form, which are then inverted to recalculate the objective

 

functions such that G = \/G,2 + G: + G; is minimized. Further details of this method

may be obtained from Asselineau et al. [1979], and the resulting computer program is at-

tached in the appendix 2.

RESULTS

The diameters of the fluid, solid and solute are the Lennard-Jones parameters as
tabulated in Reid et al. [1987], and are given in chapter 3. The of, is calculated as the

arithmetic mean of the fluid and solid, while the fluid-solid potential parameter (cf/k) is
adjusted to match the magnitude of adsorption.

van Megen and Snook have done Monte-Carlo simulations of Lennard- Jones
molecules (that simulate ethylene) on slit-like pores [1985] (Figure 4-3). Their simulation
results cannot be quantitatively compared to either experimental or SLD model results be-
cause their simulation of bulk ethylene does not represent the experimental bulk properties
0f ethylene; e. g. the simulation critical temperature (Tum) of ethylene is 226 K, as op-

posed to Tc = 282.4 K seen experimentally. Figure 4-4 shows the corresponding model

99

 

 

 

 

 

0 5 -£ is 3 ,1-0

P

Figure 4-3: Grand Canonical Ensemble Monte Carlo

Simulation of Ethylene in a Slit (TIT 5°35)
-h=eo, - - - 11:10.0, ----h=5.0, - ----- 11:3.5

nNEEOE ..€ mwuwvaXmH

100

 

 

 

 

 

 

 

 

30
---h=3.5 , — -
------ h=5 .1
25 "' ----h=10 .1
I
I.
@20 ."
.. I
s ,-
O .
E ,/
$15 " / I
(D ‘,'
'” . .r. ..........................
C3 ----------- .1
g 'O l".’
I.” 10 "' '.,o".’./
”’_~7.:él".. ___________________________
5 '25.."
l
0 02 DA 00 03
PIP“at

Figure 4-4: Ethylene on Carbon Slits at Tl'l'c of
085

101

predictions for adsorption of ethylene on a slit. Their calculations were done at T/T* of
0.95 (where T* = 202 K), yielding a ratio of TIT“...m of 0.85. We used the same ratio of
WTC (resulting temperature was 239 K) for our calculations. To convert the excess in
1.tmol/m2 to the units used by van Megen and Snook, divide by 9.323. ‘b’ is the separation
between the carbon atoms in terms of ethylene molecular diameters. Qualitatively, the
SLD model shows all the trends seen in the Monte-Carlo simulations. At low slit widths
the pore fills quickly, resulting in Type I isotherms. There is also a slight decrease in the
excess as the pressure increases, as reported by van Megen and Snook. As the slit width
is increased Type IV isotherms are seen.

Figure 4-5 shows the adsorption of ethylene in pores of different diameters at low
pressures. These results are compared to the data of Reich et. al. [1982], who have stud-
ied the adsorption of ethylene on BPL activated carbon with a surface area of 988 m2/g.
The value of the ef/k that was used in Figures 3-1 and 44 has been retained in this case.
The model predictions with a pore diameter of 2.45 nm match the experimental data rather
well. However, activated carbon has a pore size distribution with pores ranging from mi-
cro to macropore range, with an average pore size of 3.2 nm, with largest portion of pores
in the 1.5 - 2 nm region. As seen in Figures 4-6 and 4-7, at temperatures of 260.2 and
301.4 K, the model predictions at 2.45 nm are not as accurate. A factor in this tempera-
ture dependency may be due to the fact that in order to truly represent the experimental
data, we need to sum the amount adsorbed over this pore size distribution. At small pore
widths type I isotherms are seen while type IV isotherms are seen at large pore widths.
The amount adsorbed for a pore width of 4.56 nm is actually smaller than the amount ad-

sorbed in a 3.29 nm pore for pressures up to 6 bar at 212. 7 K and for the entire pressure

Amount (mmol/g)

102

 

 

 

 

 

 

 

 

 

 

 

20
- - - - Pore Dia = 4.56 nm
18 __ ----PoreDia=3.29nm _._____
Pore Dia = 2.45 nm I: "
- - - Pore Dia = 1.60 nm __ ........ } .-'- ------- ‘
16 ‘- ------ Slit Width = 1.61 nm ,-’ .
:1 Experimental Data l.’ .t’
/
14 -. // ,
0" I
12 "" I" l I
./ /'
0’ I .
10 "' " 1
. 1' .................
......... S
T u
D
0 1 2 3 4 5 6 7

Pressure (bar)

Figure 4-5: Ethylene on BPL Carbon at 212.7 K

Amount (m mollg)

16

103

 

14-

12-

O
1
r

 

 

------ Pore Dia = 4.58 nm
- - - Pore Dia = 3.29 nm
Pore Dia = 2.45 nm
---- Pore Dia: 1.60 nm
:1 Experimental Data

 

 

 

 

 

 

 

 

Figure 4-6: Ethylene on BPL Carbon at 260.2 K

Pressure (bar)

I

8

12

16

Amount (mmol/g)

104

 

- — - Pore Dia =1.60 nm

9 -.. -— - - Pore Dia = 4.56 nm
-— - — Pore Dia = 3.29 nm
Pore Dia = 2.45 nm
:1 Experimental Data /

 

 

 

 

 

 

 

 

0 4 8 12 16
Pressure (bar)

Figure 4-7: Ethylene on BPL Carbon at 301 .4 K

105

range at 260.2 K. Model predictions of adsorption of methane at high pressure were also
conducted. The results obtained were qualitatitively similar to the data of Barton et. al
[1984] (Figure 4-8).

In Figures 4-9 and 4-10, the model predictions of adsorption of a binary mixture of
ethylene and methane are shown. In this case the diameters of both ethylene and methane
are assumed to be 0.4 nm to simplify calculations near the surface of the pore. The fluid-
wall parameter for ethylene is retained as 140 K, while for methane it is taken as 107 K,
which was calculated by fitting it to the pure gas adsorption data. A binary ethylene-
methane interaction parameter of 0.022 was used [Sandler, 1989]. The composition of
ethylene in the adsorbed phase is plotted as a function of pressure at 212.7 K with a bulk
ethylene composition of 74%. The model predictions for the selectivity shows that the
composition in the adsorbed phase is inversely proportional to the pressure, as seen by the
experimental results of Reich et al. [1982]. Note that the model assumes a slit-like geome-
try to model the system. In Figure 4-5, we showed that a 1.6 nm slit shows trends similar
to that of a 2.45 nm pore for ethylene adsorbing on activated carbon at 212.7 K. We used
the slit-like geometry to represent this binary system to show that the SLD model is capa-
ble of showing the right trends with respect to composition and amount adsorbed in binary
mixtures, and as seen by the figure, the SLD model shows great promise for mixture
modeling. Figure 4-10 shows the model prediction for the amount adsorbed as a function
of pressure. As seen by the figure the SLD model qualitatively predicts the amount ad-
sorbed, and shows all the right trends seen in the figure. 1

Figure 4-11 shows the composition of ethylene in an ethylene-methane mixture at

260.2 K mixture. At this temperature, too, the SLD model predicts the selectivity of eth-

0.3

106

 

0.25 --

.0
M
1

Amount (9 CH4lg carbon)
0 E
L. 01

   
 

 

 

 

 

 

 

 

 

0-05 ' ---,--oPore Radius 1.11 nm
- — - Pore Radius 2.05 nm
Average
0 . 1 1 : i
o 40 80 120 160

Pressure (bar)

Figure 4—8: Methane on Activated Carbon at
243.15 K

 

Composition of Ethylene

107

 

0.985

 

 

 

 

1:1 :1 Experimental Data

0.98 -- —Model Prediction
3,’
a 0.975 ""
.C
0. 1:1
'0
3 0.97 -- a
L
O
3 0935
<
0
5
.E 0.96 --

0.955 --
0.95 1 1 1 1 1 1

 

 

 

0 l 2 3 4 5 6
Pressure (bar)

Figure 4-9: Adsorption of Ethylene - Methane
mixture on BPL carbon at 212.7 K and Initial
Bulk Ethylene Concentration of 0.74

108

 

 

Amount (mmollg)
\1
(II

 

 

 

 

 

 

 

7 __
O
6.5 ~-
0
6 " Eth
— ylene - Amt

5.5 ~~ 0 Data - Ethylene Amt.

’5 1 1 ' fi' 1 1

1 2 3 4 5 6

Pressure (bar)

Figure 4-10: Adsorption of Ethylene - Methane
mixture on BPL carbon at 212.7 K and Initial
Bulk Ethylene Concentration of 0.74

Composition of Ethylene in the

Adsorbed Phase

109

0.8

 

 

— Model Prediction
0.78 ~-

 

 

:1 Experimental Data

 

0.76 --

0.74 1

I

r

0.72 1

0.66 -

I

 

 

 

0.34 , 1 1
0 5 1 0 1 5
Pressure (bar) '

 

Figure 4-1 1: Adsorption of Ethylene - Methane
mixture on BPL Carbon at 260.2 K,’and Initial Bulk
Ethylene Concentration of 0.235

|||||| Ill-ll
.........

1
11
H
8
\f
11
D1
3

110

ylene very well. Figure 4-12 shows the adsorption of two members of a homologous se-
ries propylene and ethylene adsorbing on activated carbon. This activated carbon is called
Nuxit-AL, which is made by a Hungarian manufacturer and its characteristics could not be
obtained readily, it was approximated to have the same physical properties as BPL carbon.
Since ethylene and propylene are similar molecules, the binary interaction parameter was
assumed to be zero. The experiments conducted by Szepesy and Illes [1963] were done
by maintaining the pressure constant and varying the concentration of the bulk propylene.
Here we wanted to see the model predictions over a wide range of concentrations. As
seen in Figure 4-12, the model can predict the concentration in the adsorbed phase to a
very high degree of accuracy. Once again, modeling the carbon as a pore, and incorporat-

ing the pore-size distribution data into the model will enhance the model predictions.

DISCUSSION

As discussed in the two previous chapters, the limitations of the SLD model can be
attributed to: 1) the local density approximation; 2) 1ack of structure in the model; 3)
choice of the fluid equation-of-state to describe the fluid properties. In the previous
chapter, we showed how modifying the fluid equation to incorporate the Peng-Robinson
equation dramatically improved the model predictions for adsorption in flat surfaces.
Here, we modify the SLD approach to predict adsorption in slits and pores. Model pre-
dictions are compared to both molecular simulations and experimental data, and as seen by
the results the SLD approach can quantitatively predict adsorption in activated carbon by

approximating the BPL activated carbon to be a 2.45 nm pore at 212.7 K. However, as

111

 

 

 

 

 

 

 

 

 

0.9 :1
0.8 -- n

.E

‘a’

2 0.7 ~-

3 3 .

2 2

D. n, 0.6 ~-

h

0 8

C .D

g 3 0.5 --

3 11

'E <

3 0.4 -_

g :1

U 0 3 __ —Model Prediction

' :1 Experimental Data
0.2 1 1 1 1
0 0.1 0.2 0.3 0.4

Concentration of Propylene in Bulk Phase

Figure 4-12 : Adsorption of Propylene-
Ethylene mixture on BPL Carbon at 293.15 K
and 1.013 bar

SU

CV

HO
Car
bui
C165
din

are

112

the temperature is increased the model predictions get worse. It is hoped that by suitably
incorporating the pore-size distribution in the model i.e. by calculating the total amount
adsorbed as the sum of the amount adsorbed over different pore sizes, the model predic-
tions will be quantitative. Another important feature is the fact that the fluid-solid interac-
tion parameter for pores was the same as that of flat walls. This suggests that the fluid-
solid interaction parameter, while not a theoretically determined parameter, can be found
by fitting the data at one temperature and for one surface. Once such a table can be writ-
ten for a number of fluids, this model can be used as a predictive model over different
surfaces. It also strengthens the universality of the SLD approach.

Another limitation of the SLD model is the lack of structure in the model. How-
ever, at pore sizes smaller than 3 molecular diameters, Figure 4-13 shows that the ratio
a(r):a,, is symmetric about the center of the pore but the maximum is off center. As a re-

sult, the calculated maximum density would be where a(r)/ab is maximum. The trends can

be attributed to the fact that at pores smaller than 3 6f, there is a packing effect of the
molecules such that if the particle is at the center of the pore then in the horizontal plane
no other particles can be present, while if the particle is off-center, then another particle
can exist in the same horizontal plane. This suggests that there is some inherent structure
built in this model due to the sizes of the molecules which may lead to exclusion of parti-
cles. This ‘structure’ is not seen in slit-like pores. Slits are infinitely long in two-
dimensions, while pores are infinitely long only in one-dimension. Hence, more molecules

are in the vicinity of a given molecule in a slit, and the packing effect is less important.

113

 

 

 

 

 

 

 

0.2 d-

0.15 ~-
{1
\
E
m

0.1 —~

0.05 --

— Pore Dia = 1.395 nm
-0.5 0 0.5
RIO”

Figure 4-13: a(r)/a, for a Pore that can hold 2.5
ethylene molecules in the horizontal plane

 

 

114

The final limitation of the SLD model may be attributed to the local density ap-
proximation. Figures 4-14 and 4-15 show the adsorption of ethylene on graphon and on
activated carbon at 212.7 K, by treating ethylene as a homogeneous fluid (a(z)=ab.), simi-
lar to the Barrer and Robins van der Waals model. By comparing these figures with Fig-
ures 3-1 and 4-5, it is obvious that the volume occupied by the adsorbent must be ex-
cluded. The local density approximation makes the SLD model tractable for routine calcu-
lations. It can also solve the entire density profile over wide pressure ranges within a few
seconds, and can give quantitative predictions for the surface excess and amount ad-
sorbed. Hence, the local density approximation, while a limitation, significantly improves
the model predictions over homogenous models. 1

There is also an incorrect limit in a/ab in the model for slits and pores as a slit/pore
width of one molecular diameter is approached. To discuss this inconsistency let us dis-

cuss slits. In such cases, for slits, the SLD model predicts that

_"..=.3_[_P__1]=o atEIL=1 [4—23]

This requires that as the wall separation of the slit approaches 6”, aka, approaches zero.
The physical inference of this result is that either the force of attraction between the parti-
cles is zero or that particles are not present. However, particles can be present in a single
Plane when the slit is one molecular diameter in width and they certainly attract each.
other. This suggests that this SLD methodology will not work in its present form for slit

SiZes approaching 0'5. The problem arises due to the application of mean field theory that

115

 

 

 

 

 

 

 

 

 

 

45
263.15K
—— 273.15K
— -—293.15K
- - - - 313.15K
A
"E
=
3
:3.
V
m
rn
8 gag
fi,‘ 4 \
° ,-—o\---~-
-" 00°” \° ‘~‘
0 ow \\ \\‘~
\\
O
60 80 100 120

Pressure (bar)

Figure 4-14: Adsorption of Homogeneous
Ethylene on Graphon

116

 

b 0'! 0

Amount (mmollg)
00

 

 

 

Du

 

— Model Prediction
:1 Experimental Data

 

 

 

 

Pressure (bar)

Figure 4-15: Adsorption of Homogeneous
Ethylene on BPL Carbon at 212.7 K

 

117

neglects particle size, at an interface where particle size is most important. Mean field
theory usually assumes that fluid particles are present in all space between the walls of a
slit. At an interface this is incorrect; the reason can be visualized by considering a flat
wall. In flat walls, we begin integrating from z=0.5 of, because half the volume of this
particle will not be in contact with any fluid particles. This means that no fluid centers can
exist between 0 and 0.5 0,7, which is consistent with particles having finite size. As we
move the fixed particle away from the wall more and more of the particle surface will be in

contact with other particles, and it can be surrounded by other particles at z 2 20' f. In a

confined slit the SLD approach is applied for each wall, excluding particles inside 0.5 03

from each wall. When the width is just one molecular diameter, the limits overlap, thus

the entire volume is excluded.

One way to correct this inconsistency is to refitg/ab such that the integration be-
gins from 2:0 6])“, as opposed to 2:05 of}. This would neglect particle sizes for all but the
fixed particle. The problem with this solution is that ‘I’,, includes particle size, and at
Z<0.5 6,; ‘19, becomes large and positive, and hence excludes particles, thus introducing a
mathematical inconsistency.

Another approach would be to calculate the configurational integral for a slit of
one molecular width as a base reference, and add to this integral as slit widths are in-
creased. Presently, however, no clear and simple solution to this problem exists, because

the problem of finite particle size is simply transferred from the first layer to each subse-

quent layer.

118

In order to investigate the behavior of a/ab in confined spaces approaching off, we
can easily consider packing of disks in a two-dimensional channel. Solving for the con-
figurational energy of closest packing in a two-dimensional channel for disks of a finite
size by summing a l/r6 potential, the configurational energy for a channel of one molecular
diameter width is approximately 2% lower than the configurational energy of a channel 1.5
molecular diameters wide. Based on this two-dimensional result, it seems reasonable to
assume that a/ab is constant for slit widths below 1.5 off. We chose 1.5 molecular diame-
ters because 1 molecular diameter represents total constraint in packing, while 2 molecular
diameters represents total freedom in packing. Extension of the two-dimensional summa-
tion technique to three-dimensions is most easily and realistically achieved by conducting
molecular simulations in small slits, and determining dependence of configurational energy
on slit size to further refine a/ab. Molecular simulations can be conducted also for pores

to see if the approximation a/ab = constant is followed for pore widths < 1.5 ofl.

CHAPTER 5. CONCLUSIONS

There are many different models that describe adsorption ranging from the simple
empirical fits to the theoretically sound molecular simulations. The simple are widely used
because of their simplicity. Models based on statistical mechanics and molecular
simulations are theoretically sound but significant calculation time. A model that can
retain the simplicity associated with the empirical models and contains the basic physics of
the adsorption problem can be very useful. This is the basis for the development of the
SLD model.

The SLD model is a simple model based on spatial invariance of the chemical
potential, along with a cubic equation of state to describe the fluid properties. The SLD
approach can model a variety of adsorption isotherms - from the low pressure region to
the supercritical region. The model provides information about the density profile of the
adsorbed fluid. In the subcritical region the SLD model predicts Type II and Type III
isotherms on flat walls. It also predicts phase transitions at the surface (under certain
conditions), as well as an abrupt transition between an adsorbed gas layer and fluid layer
(again under certain conditions). Negative adsorption is sometimes predicted, similar to
the predictions of Sullivan [1979]. In the supercritical region the model shows the
experimentally observed cusp-like behavior near the critical point as well as the crosSovers

seen at higher pressures.

119

120

Quantitative modeling, extension to Types 1, IV and V isotherms, prediction of
clustering in supercritical fluids, and adsorption of multicomponent mixtures are
discussed in chapters 3 and 4. By using the Peng-Robinson equation, the SLD model can
quantitatively predict the adsorption of pure gases on flat walls. By modifying the solid-
fluid potential to be two body fluid-fluid potential, the SLD model can predict clustering
in supercritical fluids. The local density approximation calculations were modified to
account for this change in geometry. Comparisons to experimental fluorescence spectra
showed that the SLD model can predict the phenomenon of clustering. In slits and pores,
the SLD model predicts types I, IV and V isotherms, and can also quantitatively model
adsorption. Using the pore-size distributions will improve the model predictions.

By approximating activated carbon to be a slit, adsorption of binary mixtures were
carried out. Calculations were done for two mixtures - ethylene & methane and
propylene & ethylene. The ethylene-methane mixture calculations were done at two
different temperatures - 212.7 and 260.2 K. These calculations show that the SLD model
can qualitiatively predict the selectivity of a binary mixture. The SLD model
qualitatively predicted the amount adsorbed. By modeling the activated carbon as a pore
and suitably incorporating the pore-size distributions, the model predictions can be
enhanced.

The SLD model serves as a good approximation and can provide rapid surveys of
adsorption behavior to guide more detailed and time-consuming simulations. This. model
is not intended to replace any of the more complete theories such as Monte-Carlo
simulations, molecular dynamics, or density functional theory, but instead provides a

simple and approximate solution to the adsorption problem, The model may be viewed

121

as a compromise or a bridge between the two-dimensional equations-of-state models,
Frenkel-Halsey-Hill theory, and the more rigorous density functional or integral equation

approaches.

APPENDICES

The van der Waals equation and configurational energy calculations - Flat Walls

If we consider a canonical ensemble, the partition function Q is given by [Vera and Prausnitz, 1972]

Q._1_V~(m,-_2_)”fl"
N! ’ 21:1' A3
where A = h/(21r ka)”. The Helmholtz energy A and the pressure are related to the partition

function:
F . - 91 . u .332
3" r.» 3" rN

31an N arm

61’ 61’

where the molecular rotational, vibrational, electronic partition function, qme, is assumed to have no

volume dependence. The van der Waals free volume V}, for a system with N molecules, is given by

 

 

V, . V- .11. 1.
NA
The repulsive part of the pressure is
P"? - km 1 RT
V - Lb v b
NA

123

 

 

If <1) is the sum of two-body interactions between an arbitrarily selected central molecule and all
other molecules around it, and (b( r) is the two-body interaction potential, and g( r) the radial

distribution function, then,

d>(r) g(r) 4urzb

VIZ

° ' 1
For a van der Waals fluid, following McQuarrie [1976], g( r) = O, for r s o, and g( r) = constant, for r
> o, where o is the hard-core radius of the fluid molecule. To simplify notation, the subscript fl is
omitted from 0 throughout the appendix.

If we assume that the molecules interact with an infinite hard-core repulsive potential and an

inverse-sixth attractive potential, then for r > o,

6

Mr) - -e,°—6
f

where 6,, is the maximum energy of attraction between a pair of fluid molecules. In the bulk fluid,

this leads to

a” - 21: 61.03 Ali/3

 

124

Nonhomogeneous Systems. Consider a fluid molecule of diameter 0 at location 2 in the vicinity of
a wall, where 0.50 s z s 1.50. For convenience of integration we define a cylindrical coordinate
system with the origin at an arbitrary molecular center at z, and let y be a dummy variable to denote

axial distance and let r denote radial distance. The configurational energy can be calculated from

the two integrals:

With the local density approximation: p02) z p(z)
1b - -1/2e

5
NA _.—
,o u p(z)“ 6)

Since a is proportional to <I>, for 0.50 s z s 1.50, we find

 

3_(Q_ em _(5 6 2)

am «puma» 16 ' TE :

For 1.50 s z < co, the integration is performed in three parts:

11111 the

125

With the local density approximation p(y) = p(z) these can be integrated to give:

 

126

and

27 RITE T

a
M" 64 PC TM

 

 

NI.—

Note that atz= 0/2,a=abu,k/2, and at z - a, a - am

Note that a and its derivative with respect to z are continuous functions at z = 1.50.

127

Calculation of Configurational Integral for Clustering

This calculation is similar to the configurational integral calculation for adsorption on a flat wall and
further details are given in the appendix of a previous publication [Rangarajan et. al., 1995].

If O is the sum of two-body interactions between an arbitrarily selected central molecule and
all other molecules around it and (p(r) is the two-body interaction potential, and g( r) is the radial

distribution function, then

O = I: %¢(r)g(r)4m'2dr

For a van der Waals type fluid, following McQuarrie [1976], g( r) = O for r S 0, and g( r) =
constant, for r > 0' , where 0 is the hard-core radius of the fluid molecule. Note that the subscript f
has been dropped from 0 for the appendix to simplify the equations.

If we assume that the molecules interact with an infinite hard-core repulsive potential and an

inverse sixth attractive potential, then for r > 0' ,

O.6

¢(r)=_efl-r_6-

where e 1,. is the maximum energy of attraction between a pair of fluid molecules. In the bulk fluid,

this leads to

_ ~47refl. 0'3N,l
”’ 3

 

p

128

27: 61,. O'BNj
ab = 3

 

The relation between the partition function, the attractive and repulsive part of the pressure, a

and O are given in previously [Rangarajan et. al., 1995; Vera and Prausnitz, 1972]].

Nonhomogeneous Systems. Consider a fluid (solvent)molecule of diameter 0 at location r in the

 

vicinity of a solute molecule of diameter 0, (typically larger than 0'), where 0.50 S r — O; S 1.50_

Using cylindrical coordinates with the origin at the center of the solute molecule, we integrate across
all space with r denoting radial distance from the center of the solute molecule, x be the radial
distance from the center of an arbitrary molecule (parallel to the intermolecular vector) and y be a
dummy variable denoting axial distance (perpendicular to the intermolecular vector). The

configurational integral can be calculated from five integrals:

_ 6 °° °° 270"!de
(D1 - ‘ Ear 0' NALL P0006337);
~ ~ 2nydydr
__ 6

(DZ _ 617 0' NA110-10 p(r,y) (r2 +y2)3
_ 6 r—'" °'° ZWdydx
‘1’: “"517 0 ”Ala 1WP(xJ)——(x2 +y2)3

r+r' 0° ZWdydx
(D4 = - EH 06NA £_r.jf——r.2_(,_x)z p(x,y) (x2 + yz )3

Znydydx

O=—e 0'6N .. .. x, ———-
5 fl AJIr+r'J0 p( y) (x2 + y2 )3

where

_ (0+03)

I

129

r’2 + r2 —(1'2
2r’

II

 

O=O, +O2+O3+O4+O5

With the local density approximation: p(x, y) = p( r) and since a is proportional to O, we find that

 

 

 

a(r)_¢(r)
ab (Db
andleadsto
3 8 2 , a3 1 1 2 3 1
ar =—a —+— r-r +—- — +—O'
() 16 b(3 O'( ) r(r’2+r2—2rr” r’2+r2+2rr’) 3 ((r+r’)3))

Similarly, for 1.50’ S r — 2.24 < oo , the integration is performed as before with an additional

integral. As the fluid molecule under consideration is moved to a distance greater than 1.50 the

additional integral accounts for the cohesive energy in the region between 0' I, S x S r —0. This leads

to a six-part integral resulting in

3 16 203 l 03 1 1 203 1
a(r): ab( - I3+ ( I2 2 I ’2 2 I)+ I3)
16 3 3 (r-r) r r +r -2rr r +r +2rr 3 (r+r)

 

 

 

These modified a( r) are then used in the calculation of the fluid fugacity.

130

Calculation of Configurational Integral for Slits

Definitions: 2 is the distance from the surface of one of the walls. L is the distance
between the two wall surfaces. r is the radial distance. x is a dummy variable that denotes

axial diatance.

Case 1. Particle is near the wall: 92-5 2 s 3%

The integral is carried out by splitting the regions as follows:

Region 1. x ranges from 0 to z-0/2

2-2 21v
(D1: 4€JO6NAJO 2 Efip(x)mdrdx

1
ch =—2 3an [£-—]
1 8J0 A (Z) O 2

Note that the steps needed to carry out the integration are given in detail in the integral for

flat walls.

Region 2. x ranges from O to 0'

<1), = -4eflo6N,j:j;;27p(x)

131

( 22.:tr2)3 rdx
x r

(I)2 = 46170 3N,lrtp(z)

Region 3. x ranges from 0 to L-0/2-z

(D3 = ‘48306NA LL-olz-zjo p(x)

O

(I)3 =-—481,0‘5

2
(D3 = -38fi6 3NA1tp(Z)

—28fl63NA1tp(z) C:

 

 

 

 

~ 22nr23lrdx
(x +r)
L-0I2—zp(x)
NARI —de
6 x
r -
1
1-
(5-1-5f
_ 0 2 0 -
+<I>2+<I>3
5 1
—+3- L 1 3
aim-£1
_ 0 2 0 _

 

 

132

Note that for particle near the other wall, the result is a mirror image of the above result

i.e. replace 2 with L - 2.

Case 11. For Large Slits, there is also a region 1.5 S S — 1.5

Cl|1~

i
0

In this case, the axial distance is split into 4 parts

Region 1.x going from 0 to z-0/2
6 Z-O/Z 0-
<1>1 = -4e is N, j [0 p(x)———3drdx

2m
(x2 + r2)

2
(bl = _EEJO3NA1tp(z) 1"

 

 

Region 2. x going from O to 0 - the result is same as case I - region 2.

Region 3. x going from O to 0 - the result is same as case I - region 2.

Region 4. x going from 0 to L-0/2-z - the result is same as case I - region 3.

<I>=<I>,+<I>2+<I>3+<I>4

133

 

(I) = -28fl.03NArrp(z) -Z—+——

 

 

Case III. In small slits i.e less than three fluid molecular diameters big, all three regions
may not be present. In such cases, case I - region 2’s upper limt will become the same as
case I - region 3’s upper limit i.e. region 2’s upper limit will become L—0/2-z. Region 1

will be the same as case I - region 1.

. L-0/2—z . 21V
O2 =—4SJG6NAJO JWP(X)WH1X

. L 1
.2 z -..,.a~,.....[;-_-£]
O = O, + O2
3 L
(I) = -28 50’ NARP(Z)[; - I]

In the limit that L = 0, the integral reduces to a two-dimensional integral. Here the ‘

solution is given by

134

«1 ,21trdr
4),, “we”, p T

where p’ is given by mol/area.

 

 

For slits smaller than 2 0 in width the density is a constant. To convert density in

moles/volume to moles/area multiply by the slit width L. i.e.

However, noting that in a slit there are two walls, i.e. the surface area is twice that of a flat

Wall. the equation has to be modified to become

USing this equation we get

 

Note that at the limit L: 2 0

a;

ab 16

However, the derivative is not constant at L: 2 0.

Q|r~
ll
colt.»
II
A

136

Calculation of Configurational Integral for Cylindrical Pores

Definitions: The distance between the fixed particle and any particle is given by a. R is

the distance from the center to the pore surface. r is the radial distance, and r, is the

distance from the center of the pore and the fixed particle. 6 is the angle formed between
the center of the pore, the fixed particle and the projection of any particle on the
horizontal plane. 2 is the distance between any particle and the horizontal plane in which

the fixed particle is located.

a=\/r2 +22

The cylinder makes a circle in the horizontal plane. The locus of the circle is given by the

equation

where x and y are the x and y coordinates from the center of the circle (r is the radius of

this circle).

x = r1 — rcos(9), y = rsin(0)

Solving the above two equations we get

137

 

r = rl cos(9)+\/(R—(—;-)2 -r,2 sin2(9)

where 0' is the fluid diameter. Note that the subscript ff has been dropped from 0 to

simplify notation.

Now, the configurational energy is given as
_ 6 r _ r
<1>_—4e,,N,o Jilfidrdedz — [SH-[firdedz

where B = -4eflNA0°.
I. Particle away from the wall

Region 1. z varies from 0 to °<

2 Fr J4. cosl)+"(R—%)2-r.2 sinze rdrdedz
0 0

B: o (r2+zz)3

138

 

 

 

when

~ 1 1
=J;J:< 4Z4 - 2 _2 F9 dz
4l:’i cosO+\/[R-%) ""12 sinZBJ +22
. 1

 

 

 

2
Let b = r1 c059 +\/(R—%) --r,2 sin29

Solving this equations yields

 

O 1: 1.11 11 0 _tan"(0/b) 9

F 1203 '71" o 4123 -2b2(b2+02) 2b3

Region 2. z = 0 to 0

_ = cg]: J: qucose+,’(R-—)’-rfisin’0 _r___drd0dz

(r +2“)

Just as before solving the above equation yields

 

2: 1: _1r 0 +tan"(0/b)
B 403 4o 2b2(b2+02) 2b3

139

Hence, for particle away from the wall i.e. O S rl S R — 30' / 2, O is the sum of the two.

II. Particle nearthe wall i.e. R-30 /2Srl SR—o /2

Region 1. For 2 = 0 to °<, the integral remains the same as in 1.

Region 2. For z = o to 0

In this case, 6 cannot go from O to 1t, but can only go from O to

O’
r12 +02 -22 _(R__)2
cos"I 2 , because of exclusion due to the presence of the fixed
ZnJoz—z

 

particle.

 

Note that at z=‘/02 —(R-%-rl)2, c050 =—1, 6 =1:

The integral in the r direction is the same as before. This yields

 

d) _ l w -l l cos"w 1
F— 4641.0 COS (WMZ-ZKJO (b2 +Zz)2 dedz

140

r,2 +02 -z2 -(R-‘—;-)2

where W:
2rl ’62 -22

 

Hence for particles near the wall i.e. R- 30 /2 S rI S R—0 l2

 

2=L'l£[ 1t 0' tan"(0/b) 0+—1—-J‘:cos"(w)dz-%I:J’:os-

13 121:3 4 m—2b2(b2+02)- 2b3 46‘

III. Particle at the wall

In this case 6 varies only from O to 1t.

Region 1. z varies from 0 to °<

Here,

 

rl c089 +\/(R-%)2 - r,2 sin2 9 = 2rl cosO

Hence, the integral becomes

1

(b2 +22)2

141

’2 Moose _r__drd9dz
B :41]: 1" (r +z 23)

 

O ~ 12 l 12 1
-- = —-d9dz - l dad
[3 L1: 424 0": 4[4r,2 c0529 +22]2 2

 

 

2_ 1r - 1t(2r,2+zz)
.10 223(4’12+ZZ)3/2 Z

[V2403-

Region 2. For 2 = O to 0

Once again 6 cannot go from O to 1t, but, can only go from O to

0'
2 2 2 2

 

 

 

cos"1 2 J 2 2 2i , because of exclusion due to the presence of the fixed
r1 0 — 2
particle.
112:”1 +6 -z2—(R——) 62-12
21', JO 2 - z 2’1

Doing the calculations as before,

 

 

ch 1 _ ch-z2 1 - cos-'w 1
«B— = J: COS '(T)dz -Z‘[" L d6dz

4 [4r]2 cos2 6 + 22]2

142

 

 

 

1r ‘1. 1r(2r,2+zz) d -l _.0;J
“40

O
B _ 240.3 °223(4r,2 +zz)3/2

         

—:Z)d —_'10 J11“ “[4rl cos 1‘9+z 1d9dz

Note that in all cases 2— = _a_ = (b

O, a, —l6efl03NA1tp/3

 

0000000000 00000000000

A ndix2- rams

FLAT WALL - Pure Gas
Program to calculate surface excess or adsorption isotherms
The program uses Peng-Robinson equation of state, with a modified
a to account for exclusion. The strength of the adsorption
potential and its dependence on position is given by Psi.
The fluid-fluid potential uff = 11 total - u ads.
From this potential the various fluid properties are calculated.

The variables used are defined below.
All units used are SI unless otherwise stated.

COMMON R,TC,PC,OMEGA

Define Psi as a statement function

Use a 10-4 potential where PSI is the negative of the
intermolecular potential.

SIGFW is sigma fluid-wall, in Angstroms.

EPSFW is the fluid-wall potential in K

ALPHA is the ratio of the spacing of the graphite basal planes
(3.35 A) to SIGFW.

The number density of C atoms in graphite is 0.382 atoms/A42.
The equation used is that suggested by Lee (1.5).

PSI(ETA)=4.0*3.1415926*O.382*SIGFW*SIGFW*EPSFW*(-0.2

OD

OOOOOOO

l/ETA* *10+O.5/ETA**4+0.5/(ETA+ALPHA)**4+0.5/(ETA+2.0*
*ALPHA)* *4-1-0.5/(ETA+3.0*ALPHA)**4+0.5/(ETA+4.0*ALPHA)
**#4)

- WRITE(*,*)'Temp Plim Epsilon Sigff Sigww Tc Pc omega'
READ(*,*) T,PLIM,EPSFW,SIGFF,SIGWW,TC,PC,OMEGA

WRITE(*,*) 'EPSILON‘

READ(*,*) EPSFW

T=212.7

PLIM=7.E5

SIGFF=4.22

SIGWW=3.4

TC=282.4

PC=50.4e5

OMEGA=0.089

OPEN(UNIT=7,FILE='denpro.dat',STATUS='UNKNOWN')
OPEN(UNIT=8,FILE='adsorp.dat’,STATUS='UNKNOWN‘)

R = Gas constant J/K/mol

Tc = Critical Temperature K

Pc = Critical Pressure Nlm2
OMEGA = Acentric Factor

SIGMA = Sigma fluid-wall in m.

T = Temperature K

FLIM = Bulk fugacity limit.

143

0000000

010

M0

144

EPSFW = Fluid-wall potential well depth in K

SIGFW = Sigma fluid-wall in Angstroms

ALPHA = Spacing of graphite planes / Sigfw
SIGFF = Sigma fluid-fluid in Angstroms
SlGWW = Sigma wall-wall in Angstroms

AB = Bulk Pong-Robinson 'a'

BB = Bulk Peng-Robinson 'b'

R=8.3 14
POMEG=0.37464+1 .54226‘OMEGA-0.26992*(OMEGA "2)

TR = T/TC
SIGFW=(SlGWW+SlGFF)/2.
SlGMA=SlGFW*IE-10
ALPHA=3.35/SIGFW
ALPHA = Spacing of graphite planes / Sigfw
SlGFF = Sigma fluid-fluid in Angsuoms

Calculate a bulk and b
AB=0.45724*(R‘TC'(1+FOMEG*(1-SQRT(TR))))“2/PC
BB=0.0‘T780"R‘TC/PC

Write parameters to files

PCB=PC/l.OES
IF(TR.GE.1) GO TO 2
CALL FSAT(AB,BB,T,FUGS,PS,VV,VL)

DELP=PLIMI40.
Loop for bulk fugacity
PB=0.0ES
1:0

FOR EACH BULK FUGACITY VALUE
PB=PB+DELP

PB=PLHW
B=BB
A=AB

First calculate bulk density (DENB) and bulk pressure (BP)
CALL BVCAL(A,B,T,PB,PS,V.FB)
CALL PV(A,B,T,V,P)

DENB=1 .OIV
BP=P

145

C WRITE(7,*) PRESSURE. P
C WRITE(7,")'BULK DENSITY', DENB

DELETA=0.1
ETA=.9
EXCESS=0.0

K=0

ETA=ETA+DELETA

0‘

Calculate local fugacity
Use the following formula to get the local
fugacity (F), for a given position ETA.

F=FB*EXP(PSI(ETA)/l‘)

0000

ALNFB=ALOG(FB)
ALNF=ALNFB+PSI(ETA)/l"
PSIV=PSI(ETA)/l'

C Calculate local 11
CALL ACALCCETA,AB,A.SlGWW,SIGFW,SIGFF)

C Using local parameters calculate V and local density DENL
CALL VCALC(A,B,T,ALNF,BP,V)
DENL=1.0/V
EXCESS=SIGMA*(DENL-DENB)*DELETA*l.E6+EXCESS
C DENLG is the density in gmoles/cc
DENLG=DENU1 .0136
WRITE (7,102)ETA.DENLG
WRWBO*) 'A'O.O'DAO.O'O'B'D',"BO'Q'DETA
K=K+l
IF(K.LT.200) GOTO 6
9O CONTINUE
BP=P/l .0E5
DENBG=DEN B/ 1.136

WR1TE(6,*) BP, ',', EXCESS, ',', DENBG
WRITE(8,*) BP, ',', EXCESS, ',', DENBG

J=J+l
[F (J .LT.40) GOTO 5

146

100 CONTINUE

C WRITE(*,*)'Press l to continue'

C READ(*,*) INT

C [F (INTEQJ) GOTO 566
CONTINUE

102 FORMAT(1X,F8.2,2X.F12.5)

END

C?Itfiifififit.‘fittfitltfififitfitt.¥¥¥¥fittli***.**tt*#*ti¥ll¢¥¥$*********¥

C234567
C SUBROUTINES

Cttttttttttttttittttttttttttttttttt*ttttutti!tttttttttttttttttttttt
SUBROUTINE VCALC(A,B,T,ALNF,BP,V)
COMMON R.TC.PC.OMEGA
IV1=1
1V3=1
DPDV1=1
DPDV3=1

P1=O
P3=0

V=0

C 1F(BP.GT .PS) GOTO 16
CALL VICALC(A,B,T,ALNF,V1,IV1)
IFGVI.BQ.1)GOTO 10
CALL DPDV(A,B,T,V1,IDPDV1)
IFODPDVIJSQJ) GOTO 10

CALL PV(A.B,T,V1,PI)

10 CONTINUE
CALL V3CALC(A.B,T,ALNF,V3,IV3)

1F (IV3.EQ.1) GOTO 14

147

CALL DPDV(A,B,T,V3,IDPDV3)
IF(IDPDV3 .EQ.0) THEN

CALL PV(A,B,T,V3 .P3)
ELSE

CONTINUE

ENDIF

l4 CONTINUE
IF (P1.GT.P3) THEN
=Vl
ELSE
=V3
ENDIF
CONTINUE
C WRITE(2,")P1,P3
GOTO 15
C 16 CALL V3CALC(A,B,T,ALNF,V3,IV3)
C IF(IV3.EQ.O) GOTO 21

C 16 CALL VICALC(A,B,T,ALNF,V1,IV1)
C V=Vl

C GOTO 15

C 21 CALL PV(A,B,T,V3,P3)
C V=V3

15 RETURN

END

 

J‘AAJALJ.‘.LJ-AAAA-lAAALJJAJAAAAAJJJAJJJJAJJJJAAJA¢44al_JJ_JJJJJ
vsTTrTTVTVTTTTTTTTTTCTTTTTVTTVVTTTTTTTTTTTWTvvvveuvlvvvlt*****

148
SUBROUTINE DPDV(A,B,T,V,1DPDV)
COMMON R
Input A,B.T,V and Output = IDPDV
THIS SUBROUTINE EVALUATES THE DERIVATIVE dP/dv

AND RETURNS A VALUE OF 1 FOR IDPDV IF THE SLOPE
IS POSITIVE. i.e. THE ROOT IS UNSTABLE.

000 0

IDPDV=0
DP=-R"T/(V-B)/(V-B)+2.0"A*(V+B)/((V*V+2*B*V-B *B)**2)
IF (DP.GT.0) IDPDV=1

RETURN
END

C.t.‘..*.***..*UUtfittfitfittfiltlfififitlt*Ililt##3##*¥*******¥*******¥

SUBROUTINE VICALC(A,B,T,ALNF,V1,1V1)
COMMON R,TC,PC,OMEGA

IV1=O

C This subroutine uses a successive substituion method
C to get the small root Vl.
C Inputs A,B,T,ALNF and Output V1, IV]

C Starting guess for calculating v
Vl=l.05"B
C Iterate 60 times

D020 I = 1,60
ARG=(V1-B)/R/T

IF((V1.LE.B).OR.(ARG.LE.0.0)) THEN
IV 1 =1
C WRITE(9,")'INV1CAL V1: ',V1,' 3 = ',B,'ARG = ',ARG
C I,'I =',I,' ALNF,ALNF
GOTO 30
ELSE
DENO=ALNF+ALOG(ARG)+A*V1/R/r/(V1 '"2+2"'B *V l -B .”2)
DENO=DENO+AI2.82/Rfr/B’ALOG((V1+2.414*B)/(V1-0.4 14*B))

VI=B+BIDENO

ENDIF
C WRI'I'E(6,"')'V1,B',V1,B
C STOP

149

20 CONTINUE

ALNFC=B/(vr-B)-A*v1/Rrr/(VI*vr+2*B*v1-B*B)
ALNFC=ALNFC+ALOG(R'T/(Vl-B))
ALNFC=ALNFC+A/2.82/B/RII“ALOG((Vl-O.414*B)/(V1+2.414*B))

ERROR=EXP(ALNFC-ALNF)-l .0
IF (ABS(ERROR).GE.0.001) 1V1 =1

30 CONTINUE

C

WRITE(9,“)'VI ',ALNFCALNFN l
RETURN
END

COIOOOOQIOOOOIOIOOOttfittttfitllfiltllfifitttttttttt***¥¢**#****#****13*

C234567

00

SUBROUTINE V3CALC(A,B,T,ALNF,V3,IV3)
COMMON R,TC.PC.OMEGA
IV3=O

STARTING GUESS FOR V3 = RT/F

ALNV3=ALOG(R"I')-ALNF
V3=EXP(ALNV 3)
ITERATE TILL DONE OR 40 TIMES

DO 201:1,40

IF (V 3.LE.B) THEN
WRITE(9,")'V3 LE 8', V3,‘ B ',B
IV3=I
GOTO 30
ENDIF
ARG=B/(V3-B)-A*V3/R/T/(V3"2+2*B‘V3-B“2)
IF(ARG.GT.40.) THEN
D3=1./V 3
WRITE(9.‘)'V3 = ',V3,‘ ARG LARGE = ',ARG,'D3 =',D3
lV3=1
GOTO 30
ELSE
DENO=A/2.824/Rfl'/B I"ALOG((V3-0.4 14*B)/(V3+2.4 14 *B))
ALNV33=ALOG(R‘T)-ALNF+ARG+DENO
V3=B+EXP(ALNV3B)
ENDIF

20 CONTINUE

ALNFC=B/(V3-B)vA"V3/R/r/(V3‘V3+2*B *V3-B *B)

ALNFC=ALNFC+ALOG(R*T/(V3-B))

150

ALNFC=ALNFC+A/2.824 I"ALOG((V3-0.4 l4 I"B)/(V3~l-2.4 l4*B))/B/R/I‘

C warmer-y 1N v3. 1n fc, v3 ',ALNFC,V3
ERROR=EXP(ALNFC-ALNF)-l.0
IF(ABS(ERROR).GE.0.001) 1v3=1

30 CONTINUE

C WRITE(4,"')'V3',ALNFC,ALNF,V3
RETURN
END

Ctitfitfitfifitfitfiitttitttltlfiltitifitilttittttittttittii*1!**********¥

C234567

SUBROUTINE PV(A,B,T,V,P)

COMMON R
C This subroutine returns the value of P given a,b,T,v
C using the van der Waals equation of state.

P=R*T/(V-B)-A/(V*V+2"B*V-B‘8)

RETURN
END

C.ititttfittllfitlttltlfifilltlfitt.ttfifiliiltttttlttilttttlttt

C234567
SUBROUTINE ARANGE(R1,RZ.R3)

C PROGRAM TO PUT 3 NUMBERS IN DESCENDING ORDER
DO 20 J=1,3

IF(R2.GT.R1) THEN
TEMP=R 1
R1 =R2
R2=TEMP

ENDIF

1F(R3.GT.R2) THEN
TEMERZ
R2=R3
R3=TEMP

ENDIF

20 CONTINUE

RETURN
END

Ct.tit.*lttttittit...ttitfittlfiififilfitltiltttltttfittttttltit*********************

C234567
SUBROUTINE ACALC(EI'A,AB,A,SIGWW,SIGFW,SIGFF)

C This subroutine calculates the van der Waals a term
C after taking into account the effect of exclusion.

151

AB = Value of a in the bulk '

ETA = reduced distance from the center of wall

The main program sends a reduced distance ETA from the
center of the wall molecule, which is the basis for

the integrated 9-3 potential.

However the integrations for configurational energy have
been done from the edge of the wall molecule.

Therefore it is necessary to translate the coordinate.

BETA is the distance from the edge of the wall

in reduced units

0000000000

BETA =( ETA - (0.5 *SIGWW/SIGFW))*(SIGFW/S IGFF)

IF (BETA.LE.1.5) T1 IEN
A=AB*(5.0+6.0*BETA)/ 16.0
ELSE
A=AB"'(1 .0-1 .0/8.0/(BETA-0.5)“3)
EN DIF

RETURN
END

(:t******ttttttttttttttttttttttitittttfitittittit$ttttttttttt

&7l 2 3 4 5 6 72

O.‘ttttlttIttitfitttlltfitlfiO...fit.*ttttt'*tlfififitfititlttfiifit****t***

" SUBROUTINE CUBIC "

I‘********¢$tttttfillltltttltttit.##1##.#I*************************

" THIS SUBROUTINE FINDS THE ROOTS OF A CUBIC EQUATION OF THE
"' FORM X“3 + A2‘X“2 + AI‘X + A0 = 0 ANALYTICALLY. *

*ttfififitttttttttfi*tlt*‘fittttttt*lfiitliitttt*ttlttttt******t********

"' VARIABLES
SUBROUTINE CUBIC(A2,A l,A0,Rl .RZ,RB,C1,C2,C3,IFLAG)
DOUBLE PRECISION CHECK,DAO,DA1.DA2,PI,P2,Q,R,SS1,552
COMPLEX ES 1 ,ES2,S l .8221 ,ZZZ3,CCHECK
DAO = DBLE(AO)
DAI = DBLE(AI)
DA2 = DBLE(A2)
Q = DA1/3.D00 - DA2*DA2/9.DOO
R = (DAI *DA2-3.D00*DAO)/6.D00 - (DA2/3.D00)**3
CHECK = Q"3 + R‘R
IF(CHECK.GT.0.0) THEN
IFLAG = 1
P1 = R+DSQRT(CHECK)
P2 = R-DSQRT(CHECK)
IF(PI .LT.0.0) THEN
$81 = -DEXP((DLOG(-1.D00*PI))/3.D00)
ELSE
$81 = DEXP((DLOG(PI))/3.DOO)
ENDIF
IF(PZ.LT.0.0) THEN
SS2 = -DEXP((DLOG(-I.D00*P2))/3.DOO)

00000000

 

152

ELSE
552 = DEXP((DLOG(P2))/3.DOO)
ENDIF
R1 = $51 + $52 -DA2/3.D00
R2 = -(SSI+552)-DA2/3.D00
R3 = R2
C l = 0.0
C2 = (SQRT(3.))"(551-552)/2.D00
C3 = - C2
ELSE IF (CHECK.LT.0.0) THEN
IFLAG = 3
RR = 1."‘R ‘V
RECK =1 .*CHECI(
CCHECK = CMPLXCRECK,0.0)
E51 = CLOG(RR+CSQRT(CCHECK))/3.
E52 = CLOG(RR-CSQRT(CCHECK))/3.
51 = CEXP(ES 1)
52 = CEXP(ESZ)
21 = (51-1-52) - A2/3
22 = ~(51+52)/2 - A2/3 +(CMPLX(0.0,3“.5))"(S1-52)/2
23 = -(51+52)/2 - A2/3 - (CMPLX(0.0,3".5))"(S1-52)/2
R1 = REAL(ZI)
R2 = REAL(Z2)
R3 = REAL(Z3)
C1 = 0.0
C2 = C1
C3 = C l
ELSE
[FLAG = 2
1F(R.LT.0.0) THEN
551 = -DEXP((DLOG(- l .DOO‘R))/3.D00)
ELSE IF(R.EQ.0.0) THEN
551 = 0.0
ELSE
551 = DEXP((DLOG(R))/3.D00)
ENDIF
552 = 551
R1=551+SS2-DA2/3.D00
R2 = -(551+SS2)/2-DA2/3.D00
R3 = R2
C I = 0.0
C2 = C1
C3 = C2
ENDIF
RETURN
END

 

C
C

153

PROGRAM FOR CLUSTERING

ittfififiiifittttfitfititltttttt*t.‘ttttlttttttttttttttt*ttttlfllttttttttfil

C234567

0000000000 0000000000000

00000000000000

Program to calculate surface excess or adsorption isotherms

The program uses Peng-Robinson equation of state, with a modified
a to account for exclusion. The strength of the adsorption

potential and its dependence on position is given by Psi.

The fluid-fluid potential uff = u tOtal - u ads.

From this potential the various fluid properties are calculated.

The variables used are defined below.
All units used are 81 unless otherwise stated.

COMMON R.TC.PC.OMEGA

Define Psi as a statement function

Use a 10-4 potential where P81 is the negative of the
intermolecular potential.

SIGFW is sigma fluid-wall, in Angstroms.

EPSFW is the fluid-wall potential in K

ALPHA is the ratio of the spacing of the graphite basal planes
(3.35 A) to SIGFW.

The number density of C atoms in graphite is 0.382 atoms/AAZ.
The equation used is that suggested by Lee (1.5).

P51(ETA)=4.‘EPSGS‘(1./ETA“6ol JETA"12)

R = Gas constant I/K/mol

Tc = Critical Temperature K

Pc = Critical Pressure N/m2

OMEGA = Acentric Factor

SIGMA = Sigma fluid-wall in m.

T = Temperature K

FLIM = Bulk fugacity limit.

EPSGS = Fluid-wall potential well depth in K
SIGFW = Sigma fluid-wall in Angstroms
ALPHA = Spacing of graphite planes / Sigfw
SIGFF = Sigma fluid-fluid in Angstroms
SIGWW = Sigma wall-wall in Angstroms
AB = Bulk Peng-Robinson 'a'

BB = Bulk Peng-Robinson 'b'

WRITE(*,")' Plim ,T,EPSILON'
READ("',") PLIM,T.EPSGS

R=8.3 14
TC=282.4
PC=50.4E5
OMEGA=0.089
SIGFF=4.22

 

154

SIGWW=7. 14
SIGFW=(SIGWW+S IGFF)/2.
SIGMA=SIGFW"1E-10

Write parameters to files
OPEN(UNIT=7,FILE='denpro.dat',STATUS='UNKNOWN')
OPENCUNIT=8,FILk‘adsorp.dat’,STATUS='UNKNOWN')

WRITE(7,*) ' Ethylene at ',T,‘ K with an epsilon of ',EPSGS
WRITE(8,"') ' Ethylene at ‘,T,‘ K with an epsilon of ',EPSGS
WRITE(8,"‘) 'P (Pa) Excess (micro mol/m"2) Den (mol/cc)'
WRITE(8,*) 'Clus (# molecules solvent/molecule solute)’
WRITE(8,997)

WRITE(8,")

WRITE(8,998)

997 FORMAT(7X,'P',9X,'EXCES5',8X,'CLUS',9X.'DEN')
998 FORMAT(7X,'0,',10X,'0,',13X,'0,',13X,'0')

v.0

0H0f5

TR=T/T C
FOMEG=0.37464+1 .54226*OMEGA'0.26992*(OMEGA"2)

Calculate a bulk and b
AB=0.45724‘(R‘TC’(I+FOMEG*(I-SQRT(TR))))"2/PC
BB=0.07780’R‘TCIPC

PCB=PCI 1 DE
IFCI'R.GE.1) GOTO 2
CALL FSAT(AB,BB.T,FUGS.PS.W,VL)

DELP=PLIMI40.

Loop for bulk fugacity
PB=0.0E5
I=0

INN!EAKELBULKIRKENCHWFVALUE
PB=PB+DELP

PB=PLHH

B=BB

A=AB

First calculate bulk density (DENB) and bulk pressure (BP)
CALL BVCAL(A,B,T,PB,PS,V,FB)
CALL PV(A,B,T,V,P)
DENBN=1.0/V
In this method, we calc. the bulk density from the program
in which eta starts at 20.9 and comes inward. Sometimes do check
DENBN and DENB and make sure
WRITE?!) 'BULK DENSITY = ?'
READ("',"') DENB
BP=P

WRITE(7,"')'P = ', P,’ Pa',‘ DEN CALC = ',DENBN,’ gmol/m3'
WRITE(7,')'DENB = ',DENB,’ gmol/m3',’ FB = ',FB,' Pa'

 

155

WRITE(7,") 'DEN (gmol/m3) PMV (cc/gmol)’

WRITE(7,“) 'CLUS (it excess solvent moles/moles. solute)‘
WRITE(7,“) 'AMOUNT (# total solvent molec/molec. solute)‘
WRITE(7,"')

WRITE(7,999)

999 FORMAT(7X,'ETA',8X.'DEN',10X,'PMV',9X,'CLUS',7X,'AMOUNT)

C

O

0000

C

Start iterating in ETA
DELETA=0.01
ETA=.99
EXCESS=0.0
CLUS=0.0
K=0
ETA=ETA+DELETA

Calculate local fugacity

Use the following formula to get the local
fugacity (F). for a given position ETA.
F=FB*EXP(PSI(ETA)/I')

ALNFB=ALOG(FB)
ALNF=ALNFB+PSI(ETA)/T
PSIV=PS 1(ET A)/T

Calculate local a
CALL ACALC(ETA.AB,A,SIGWW.SIGFW,SIGFF)

IF(TR.GE.1.) GO TO 20
CALL FSAT(A,B,T.FUGS,PS,VV,VL)

Using local parameters calculate V and local density DENL

20 CALL VCALC(A,B,T,ALNF,BP,PS.V)

00000

C

DENL=1 .O/V
IF((DENUDENB.GT.0.999).AND.(DENL/DENB.LT.1.001)) THEN

DENL=DENB '
ENDIF
CLUS=SIGMA"3’(DENL-DENB)“4.'3.1416‘ETA‘ETA‘DELETA*6.023E23+CLUS

IF(K .EQ.0) THEN
CLUS=CLUS-DENB*4./3*3.I416*SIGMA"3*6.023E23
write(",") eta,clus

ENDIF

AMT is the excess amount (moles) of solvent moles in the system

TV is the vol of the system at any given eta
EV is the equivalent volume in the pure system
TW is the vol of the solvent molecules in the actual system

provided they are present at bulk density

AMT=4."'3.1416"SIGMA"3‘ETA‘ETA*(DENL-DENB)*DELETA+AMT
TW=4.I3.*3. 1416‘SIGMA"3"(ET A"3-1 .)
TV=4./3."3.1416*(ETA’SIGMA)"3

EV=A MT/DENB+TW

PMV=TV-EV

PM V=PMV *6.023E29

AMOUNT is the total # of solvent molecules in the system

156

C TN is the total it of molecules in the pure system at any eta

C EXCESS is the excess # of solvent molecules in the system
AMOUNT =4.*3.14 l6*5IGMA**3‘ETA'ETA*DENL‘DELETA*6.023E23+AMOUNT

TN=DENB"6.023E23‘TV
EXCESS=AMOUNT-TN

C DENLG is the density in gmoles/cc
DENLG=DENUI .OE6

IF((K.EQ.1 14).OR.(K.EQ. 190)) THEN
WRITE(".*) ET A,AMOUNT

EN DIF
WRITE (7,")ETA,',',DENL,',',PMV,',',CLUS,',',AMOU NT

=K+1
IF(K.LT.2000) GOTO 6

90 CONTINUE

BP=P/l .OES
DEN BG=DENBI I .56

WRITE(8.’)BP.',',EXCESS,',',CLUS,',',DENBG
WRITE(6,*) BPEXCESS,CLUS.PMV

I=I+I
C IF (J.LT.40) GOTO S

100 CONTINUE

102 FORMAT(1X.F8.2,2X.F 12.5.2X,FIO.5)
END

C***¢*****fitttttttttfiIt.*tlttifittfittlttittttfl***¥************tI“!!!

C234567

C SUBROUTINES
Cttaste-teammat-tantrum*tttstseasesetetsasteeeeatatttatetetsetta:tum

C Use subroutine VCALC as per fiat wall program

C Use subroutine DPDV as per flat wall program

Use subroutine VICALC as per fiat wall program
Use subroutine V3CALC as per fiat wall program

Use subroutine PV as per fiat wall program
Use submarine ARANGE as per fiat wall program
SUBROUTINE ACALC(ETA,AB,A,SIGWW,SIGFW,SIGFF)

This subroutine calculates the van der Waals a term
after taking into account the effect of exclusion.

AB = Value of a in the bulk

C

C

C

C

C

C

C

C

C ETA = reduced distance from the center of wall

 

157

The main program sends a reduced distance ETA from the
center of the wall molecule, which is the basis for

the integrated 9-3 potential.

However the integrations for configurational energy have
been done from the edge of the wall molecule.

Therefore it is necessary to translate the coordinate.

BETA is the distance from the edge of the wall

in reduced units

BETA = (ET A*(SIGFF+SIGWW)-SIGWW)/(2.'SIGFF)

00000000

SIGMA1=SIGFF

SIGMA2=SIGWW
R=(SIGMA1+SIGMA2)/Z.
Z=ETA*(SIGMAI+SIGMA2)/2.
R1=(R*R-SIGMA1*SIGMA1+Z"Z)/(2"Z)

IF (BETA.LE.I .5) THEN
Al=8./3.+2."'(Z-RI)/SIGMAI
A2=SIGMAI *rs/(zrmemz-z-zrm *2»
A3=(SIGMAI"3)/(Z*(R*R+Z*Z+2.*R*Z))
A4=2.*SIGMAI"3/(3."(Z+R)"*3)
A=AB"(AI+A2-A3+A4)"3116.

ELSE
A l=16J3.-2."SIGMAI"3/(3.*(Z-R)“3)
A2=2.*SIGMA1”3/(3.*(Z+R)"3)
A3=SIGMA1 I""‘3/(Z"‘(R"'R+Z"Z-2."Z"'R))
A4=SIGMA1“3/(Z*(R"'R+Z"Z+2.*Z"R))
A=AB"(AI+A2+A3-A4)"3./16.

EN DIF

RETURN
END

C Use subroutine CUBIC as per fiat wall program

C
C

158

SLITS - PURE GAS

Ilfifitttttttttttltlil*i*****$***¥**#******************INN!!!*********

C234567

nonooonoonoon

0000000000 0000000000000

Program to calculate surface excess or adsorption isotherms

The program uses Peng-Robinson equation of state, with a modified
a to account for exclusion. The strength of the adsorption

potential and its dependence on position is given by Psi.

The fluid-fluid potential uff = u total - u ads.

From this potential the various fluid properties are calculated.

The variables used are defined below.
All units used are SI unless otherwise stated.

COMMON R,TC,PC,OMEGA

Define Psi as a statement function
Use a 10-4 potential where P81 is the negative of the
intermolecular potential.

SIGFW is sigma fluid-wall, in Angstroms.

EPSFW is the fluid-wall potential in K

ALPHA is the ratio of the spacing of the graphite basal planes
(3.35 A) to SIGFW.

The number density of C atoms in graphite is 0.382 atoms/AAZ.
The equation used is that suggested by Lee (1.5).

P511(ETA)=4.0"'3.1415926*0.382*SIGFW‘SIGFW'EPSFW*(-O.2
I/EI‘A“ lO+0.5/ETA"4+0.5/(ETA+ALPHA)" l"4«l'0.5/(ETA‘1-2.0"I
*ALPHA)"4+0.5/(ETA+3.0"ALPHA)"*4+0.5/(ETA+4.0‘ALPHA)

ttt4)

PSI2(XI)=4.0*3.I415926‘0.382*SIGFW*5IGFW‘EPSFW*(-O.2
l/XI"10+0.5/X1"4+0.5/(X1+ALPHA)"4+0.5/(X1+2.0"'
‘ALPHA)'*4+0.5/(XI+3.0*ALPHA)**4+0.5/(X1+4.0‘ALPHA)

*tt4)
PSI(ET A)=PSI 1 (ET A)+PSI2(XI)

R = Gas constant J/K/mol

Tc = Critical Temperature K

Pc = Critical Pressure N/m2

OMEGA = Acentric Factor

SIGMA = Sigma fluid-wall in m.

T = Temperature K

FLIM = Bulk fugacity limit.

EPSFW = Fluid-wall potential well depth in K
SIGFW = Sigma fluid-wall in Angstroms
ALPHA = Spacing of graphite planes / Sigfw
SIGFF = Sigma fluid-fluid in Angstroms
SIGWW = Sigma wall-wall in Angstroms
AB = Bulk Peng-Robinson 'a'

159

C BB = Bulk Pang-Robinson 'b'
C ALPHA = Spacing of graphite planes / Sigfw
C SIGFF = Sigma fluid-fluid in Angstroms

WRITE(‘,*)' Plim ,TEPSILON, WIDTH'
READ(",*) PLIM,T,EPSFW.ZET A

R=8.3 l4
TC=282.4
=50.4E5
OMEGA=0.089
SIGFF=4.22

SIGWW=3.4
SIGFW=(SIGWW+S 1GFF)/2.

SIGMA=SIGFW*1E-10
ALPHA=3.35/51GFW

C Write parameters to files
OPEN(UNIT=7,FILE='denpro.dat',STATUS='UNKNOWN')

OPEN (U NIT =8,FILE='adsorp.dat'.5TATUS='UNK NOWN' )

WRITE(7,"') ' Ethylene at ',T,’ K with an epsilon of ',EPSFW
WRITE(8,") ' Ethylene at ',T,’ K with an epsilon of ',EPSFW
WRITE(7,*) ' with slit width of ',ZETA .' ethylene mol dia'
WRITE(8,"') ' with slit width of ‘,ZETA ,' ethylene mol dia'

WRITE(8,"‘)
WRITE(8,*)'P bar Exc um/m2 Amt mmol/g Amtl mmol/g Den gmol/cc'

Wimsf.) 0’" '9 09" '90," '9 0

TR=TIT C
POMEG=Q37464+I .54226*O MEGA 41.26992‘(O MEGA l”2)

C Calculate a bulk and b
AB=0.45724*(R'TC*(I+FOMEG‘(I-SQRT(TR))))"2/PC

BB=0.07780*R"'TC/PC

PCB=PCII .OES
IF(TR.GE.1) GOTO 2
CALL FSAT(AB,BB,T,FUGS.PS,VV,VL)

2 DELP=PLIMI40.

C Loop for bulk fugacity
PB=0.0ES
1:0

C FOR EACH BULK FUGACITY VALUE
C PB=PLIM
3 PB=PB+DELP

B=BB

A=AB

C First calculate bulk density (DENB) and bulk pressure (BP)
CALL BVCAL(A,B,T,PB,PS,V,FB)

 

0000

C

C
C

160

CALL PV(A,B,T,V,P)
DENB=1.0/V

BP=P

WRITE(7,"‘)'P = ', P,‘ bar'

WRITE(7.")'DENB = ',DENB,‘ gmol/m3',‘ FB = ',FB.‘ bar'

WRITE(7.*)
WRITE(‘7,")' ETA',', ',' DEN gmol/m3',', ','Log(Fugacity)'

BET A=0.5-(ZET A-1 .)l 100.
DELET A=(ZE'TA-1.)/ 100.
EXCESS=0.0
AMOUNT=0.0
AMOUNT1=0.0

L20
Start iterating in ETA

BET A=BETA+DELETA
ETA=(BETA*2.*SIGFF+SIGWW)/(51GFF+SIGWW)

XI=(ZETA“SIGFF+SIGWW/‘2.-BETA*SIGFF)/SIGFW

Calculate local fugacity
Use the following formula to get the local

fugacity (F), for a given position ETA.

F=FB*EXP(PSI(E'TA)/I')
ALNFB=ALOG(FB)
ALNF=ALNFB+PSI(ETA)/I‘
PSIV=PSI(ETA)/I'

Calculate local a
CALL ACALC(BETA,AB,A.ZETA)

IF(TR.GE.1.) GO TO 20
CALL FSAT(A,B,T.FUGS,PS,VV,VL)

Using local parameters calculate V and local density DENL

20 CALL VCALC(A.B,T,ALNF,BP,PS,V)

DEN I: 1 .ON
EXCESS=5 IGMA‘(DENL-DENB)*DELETA" 1 .E6/2.+EXCESS

Amount is den(gmol/m"3)"sigma(m)‘deleta'S.A.(m"2/g)

to convert Amount to mmol/g multiply by 1000
AMOUNT=DENL‘SIGMA*DELETA*988."'I.E3+AMOUNT

AMOUNT1=DENL‘S IGMA‘DELETA*988.I2."1 .E3+AMOUNT1

DENLG is the density in gmoles/cc

DENLG=DENU1 .0156
WRITE (7,"')BETA,',',DENLG,','.ALNF

L=L+I
IF(L.LT.101) GOTO 5

BP=P/l .OES
DEN BG=DENBI1 .E6

 

161

WRITE(8,*)BP,',',EXCESS.','.AMOUNT,',',AMOUNTI ,‘,'.DEN BG
WRITE(6,") BP,EXCESS .AMOUNT,AMOUNT1

I=I+l
IF(I.LT.40) GOTO 3

102 FORMAT(1X,F8.2,2X.F 12.5,2X,F10.S)
END

Cflt$ttititttltttltittt##3##tttttttttttttfi*******tilt**************

C234567

C SUBROUTINES
Ctttttttttttttttttttttttttttttttitttttttttttrtt*IrrlrtrtrrtntrhkhtrrlrrlrIt:Irwin":

C Use subroutine VCALC as per flat wall program
Use subroutine DPDV as per flat wall program

Use subroutine VICALC as per fiat wall program
Use subroutine V3CALC as per fiat wall program

C

C

C

C Use subroutine PV as per fiat wall program
C

Use subroutine ARANGE as per fiat wall program

C234567
SUBROUTINE ACALC(BETA,AB,A,ZET A)

This subroutine calculates the van der Waals a term
after taking into account the effect of exclusion.
AB = Value of a in the bulk

ETA = reduced distance from the center of wall
The main program sends a reduced distance ETA from the

center of the wall molecule, which is the basis for

the integrated 9-3 potential.
However the integrations for configurational energy have
been done from the edge of the wall molecule.
Therefore it is necessary to translate the coordinate.
BETA is the distance from the edge of the wall

in reduced units

SA=ZETA-0.5-BETA

IF (BET A.LE.1.5) THEN
A I=BETA+0.5
A2: I .- l .ISA"3
A3=1./3."'A2
A=AB*(A1+A3)*3./8.
END IF

IF(BETA.GT.(ZETA-1.5)) THEN

AI=SA+I.
A2=l .- 1 J(BETA-0.5)"3

 

162

A3=A2I3.
A=AB‘(A1+A3)*3./8.
END IF

IF((BETA.GT .1 .5).AND.(BETA.LE.(ZETA-l .5))) THEN
AI=7./3.-1./3./(BETA-0.5)""“3
A2=1.-1./SA"3
A3=A2/3.
A=AB*(A1+A3)*3./8.
END IF

RETURN
END

C Use subroutine CUBIC as per flat wall program

163

C PORE -PURE GAS

C *ltttfitttfitifitt***t********¥*t*¥**************#*******************

C234567
Program to calculate surface excess or adsorption isotherms

The program uses Pcng-Robinson equation of state, with a modified
a to account for exclusion. The strength of the adsorption

potential and its dependence on position is given by Psi.

The fluid-fluid potential uff = u total - 11 ads.

From this potential the various fluid properties are calculated.

The variables used are defined below.
All units used are SI unless otherwise stated.

0000000000000

DIMENSION H90(30000),H91(30000),H92(30000),H93(30000),H94(30000)
DIMENSION H3(30000),H31(30000),H32(30000),H33(30000),H34(30000)
DIMENSION AC(30000),RHO(30000),D15T(30000)

REAL TC.PC,R,PLIM,AB.A,BB.B.T.OMEGA,DELP

COMMON R.TC.PC.OMEGA.AC

PSI(ETA)=0.382*3.14159’3.141S9'SIGFW'SIGFW’EPSFW*(3.‘SIGFW"*4“
I(H3(K)"(RFW/(RFW*RFW-ETA‘ETA'SIGMA*SIGMA))"4+H31(K)*((RFW+ALPHA)
*/((RFW -1-ALPHA)"l I*2-ETA *ETA *5lGMA‘SIGMA))**4+H32(K)*((RFW+
*2.‘ALPHA)/((RFW+2.*ALPHA)**2-ETA'ETA‘SIGMA'SIGMA))**4+H33(K)
"((RFW+3.‘ALPHA)/((RFW+3.*ALPHA)**2-E'TA‘ETA *SIGMA*SIGMA))**4+
*H34(K)‘((RFW+4."'ALPHA)/((RFW+4.‘ALPHA)*‘2-ETA‘ETA‘SIGMA*SIGMA))
“‘4)-63./32."'SIGFW“10"(H90(K)"'(RFW/(RFW*RFW—ETA"'ETA*SIGMA*SIGMA))
"*10+H91(K)’((RFW+ALPHA)/((RFW+ALPHA)**2-ETA*ETA*SIGMA*SIGMA))

"' "" 10+H92(K)"((RFW +2. I"Al..PHA)/((RFW +2. l"ALPHA)"'"‘2-E’TA"'E1‘A"‘SIGMA"‘
'SIGMA))*" 10+H93(K)*((RFW +3.‘ALPHA)/((RFW +3.*ALPHA)"2-ETA*ET A
“SIGMA‘SIGMA))"10+H94(K)*((RFW+4.*ALPHA)/(CRFW+4.*ALPHA)**2-
I"ETA"EI‘A"SIGMA"SIGMA))""" 10))

DO 10 1=1.sr.1
OPEN(UNIT=2,FILE='hyperS-9.dat',STATUS='UNKNOWN')
READ(2,") ET,F0,Fl ,F2,F3,F4
C WRITE(*,*)ETAI

J=ET*1000

H90(J)=F0

H9I(I)=Fl

H92(I)=F2

H93(J)=F3

H94(I):F4
C WHIP—('a') J.Er.H90(J).H9l(J).H92(J).H93(J)
10 CONTINUE

DO 15 M=1,8I,I
OPEN(UNIT=4.FlLE='hypcr5-3.dat'.STATUS='UNKNOWN')
READ(4,“) ET Al ,GO,G 1 ,GZ,G3.G4
J=ETA1"'1000

C

164

H31 (I)=G 1
H32(I)=G2
H33(I)=G3
H3(I)=GO
H34(J)=G4
WRITE("‘."‘) J.ETAI .H3(J).H31(I).H32(J).H33(J)

15 CONTINUE

0000000000000000 0o~0

DO 6 MA=I.81.1
OPEN(UNIT=3,FILE='ethdia5.dat'.STATU5='UNKNOWN')
READ(3."') ETAIANEW

=ETA1 "' 1000
AC(J)=ANEW
WRITE(‘.*) JETAIACU)

CONTINUE

write(*.*) AC(O)

R = Gas constant I/K/mol

Tc = Critical Temperature K

Pc = Critical Pressure N/m2

OMEGA = Acentric Factor

SIGMA = Sigma fluid-wall in m.

T = Temperature K -

FLIM = Bulk fugacity limit.

EPSFW = Fluid-wall potential well depth in K
SIGFW = Sigma fluid-wall in Angstroms
ALPHA = Spacing of graphite planes / Sigfw
SIGFF = Sigma fluid-fluid in Angstroms
SIGWW = Sigma wall-wall in Angstroms
AB = Bulk Peng-Robinson 'a'

BB = Bulk Pong-Robinson 'b'

ALPHA = Spacing of graphite planes / Sigfw
SIGFF = Sigma fluid-fluid in Angstroms

WRITE("'."')' Plim ,T,EPSILON'
READ("',"') PLIM.T.EPSFW

R=8.3 14

TC=282.4

PC=50.4E5

OMEGA=0.089

SIG MA=4.22

5 IGWW=3.4

ALPHA=3.35

RFW= 1 2.25
SIGFW=(SIGWW+S IGMA)/2.
SIGGS=SIGMA*IE-10

Write parameters to files
OPEN(UNTT=7,FTLE='denpm.dat'.STATUS='UNKNOWN')
OPEN(UNIT=8,FILE='adsorp.dat'.STATUS='UNKNOWN')

WRITE(7.") ' Ethylene at ',T.‘ K with an epsilon of ',EPSFW
WRITE(8,"') ' Ethylene at ',T.‘ K with an epsilon of ',EPSFW

00.10

M

0000

165

WRITE(8,*)' RFW is ',RFW. ' angstroms '. 'R is RFW - 1.7'

WRITE(8.*)

WRI’I'E(8.")' Pbar',', ‘,'Excess um/m2'.'. ','Den molcs/cc'
WRITE(8.") 0.', '. 0.'. '. 0

TR=T/T C

FOMEG=O.37464+I .54226’OMEGA-0.26992*(O MEGA "2)

Calculate a bulk and b
AB=0.45724"(R*TC*(1+FOMEG*(I-SQRT(TR))))**2/PC
BB=0.07780*R*TC/PC

PCB=PC/1.0E5
IF(TR.GE. 1) GOTO 2
CALL FSAT(AB,BB.T.FUGS,PS,VV,VL)

DELkPLIM/40.

Loop for bulk fugacity
PB=0.0E5
I=0

FOR EACH BULK FUGACITY VALUE
PB=PB+DELP

PB=PLIM
B=BB
A=AB

First calculate bulk density (DENB) and bulk pressure (BP)
CALL BVCAL(A,B,T,PB,PS,V,FB)
CALL PV(A,B.T,V,P)
DENB=1.0/V
BP=P
WRITE(7,"')'P = '. P.‘ bar'
WRITE(7.‘)'DENB = ',DENB.‘ gmol/m3'.’ FB = ',FB.‘ bar'
WRITE(7,")
WRITE(7.')' ETA',', DEN gmol/m3'.'. ‘,'Log(Fugacity)'

Start iterating in ETA
DELETA=0.025
E‘I‘A=2.025
EXCESS=0.0
AMOUNT=0.0
AMOUNTI=0.0
KDEL=25
K=2025
L=0
ETA=ETA~DELETA
K=K-KDEL

Calculate local fugacity

Use the following formula to get the local
fugacity (F). for a given position ETA.
F=FB"'EXP(PSI(ETA)/I')

 

0000

C

166

ALNFB=ALOG(FB)
ALNF=ALNFB+PSI(ETA)/T
PSIV=PSI(ETA)/T

Calculate local a
WRITE(".*)K,H3(K).H3I(K).H32(K).H33(K)
WRITE(*,*)K.H90(K).H91(K),H92(K).H93(K).PS1(ETA)
WRITE(*.*) K. AC(K)

AC(O)=0.9262953
CALL ACALC(ETA,K,AB.A)

IF(TR.GE.1.) GO TO 20
CALL FSAT(A.B.T.FUGS.PS,VV,VL)

Using local parameters calculate V and local density DENL

20 CALL VCALC(A.B.T.ALNF.BP.PS.V)

C

C
C
C
C

DENL=1 .0/V

 

Using Rectangular Rule: .

EXCES S=SIGGS *(DENL-DENB)"DELETA* l .E6+EXCESS
Amount is den*r*dr*2 Pi L; 2 Pi L = 988/Pore Radius.
Amt.=denl*siggs'eta*siggs‘deleta‘988/(2.5'siggs) [mmol/g]

AMOUNT =SIGGS*DENL'ETA‘DELE’TA"988/2.5* l .e3+AMOUNT

 

RHO(L)=DENL
DIST(L)=ETA

DENLG is the density in gmoles/cc
DENLGzDENUI .0E6
.WRITE (7.*)ETA.’.'.DENL.'.‘.ALNF

L=L+l
IF(L.LT.81) GOTO 5

CALL SIMP(RHO.DIST.EXI.AMT1.L,DENB.SIGGS.A3)

BP=P/ I .0E5

DENBG=DENB/1.E6

WRITE(8."') BP, ',', EXI . '.'.AMTI .'.',DENBG
WRITE(6.*) BP, ',', EXCESS. ',',AMOUNT,'.',DENBG
WRITE(6.*) BP, ',',EX1.',‘. AMT] .A3

J=I+l

IF (J.LT.40) GOTO 3

102 FORMAT(1X.F8.2.2X.F12.5.2X.F10.5)

END

CIit.Ililtlttttlttttittfitttttttlfitfilttt*ttIII!!!******************

C234567

C

SUBROUTINES

Ct...fittitl##1##tit.tit.*tfitttlttttitltttifi*tt*tfitttit#**********

C

THIS IS A SIMPSON'S RULE SUBROUTINE

 

167

SUBROUTINE SIMP(RHO.DIST,EX1.AMT1.L.DENB,SIGGS.A3)
DIMENSION RHO(30000),DIST(30000)

E1=0
A1=0
E2=0
A2=0

DO 250 I=0.L-1
A3=A3+SIGGS I"RHO(I)"'DIST(I)“0.025 I"988./2.5"' 1 .E3
E3=E3+SIGGS ‘(RHO(I)-DENB)*0.025*1 .E6

250 CONTINUE

DO 200 I=1,L-2,2
E1=E1+(RHO(I)-DEN B)
A 1=A1+RHO(I)*(DIST(I)"SIGGS)

200 CONTINUE

DO 225 1=2.L-3.2
E2=EZ+(RHO(1)-DENB)
A2=A2+RHO(I)‘(DIST(I)’SIGGS)

225 CONTINUE

000000

0000000

EXI=RHO(O)-DENB+RHO(L-1)—DENB+4.*EI+2.*EZ
EXI=0.025*SIGGSI3.*EXI'I.E6
AMTI=0.025/3.*(DIST(L.I)*SIGGS ‘RHO(L-l)+4.*Al+2.*A2)
AMT I=988./2.5"AMT1 *1.E3
RETURN
END

Use subroutine VCALC as per fiat wall program

Use subroutine DPDV as per fiat wall program

Use subroutine VICALC as per flat wall program

Use subroutine V3CALC as per fiat wall program

Use subroutine PV as per fiat wall program

Use subroutine ARANGE as per fiat wall program

SUBROUTINE ACALCCETA.K .AB.A)

This subroutine calculates the van der Waals a term

after taking into account the effect of exclusion.

AB = Value of a in the bulk

ETA = reduced distance from the center of wall

The main program sends a reduced distance ETA from the
center of the wall molecule, which is the basis for

the integrated 9-3 potential.

DIMENSION AC(30000)

168

COMMON R.TC.PC.OMEGA.AC
A=AB‘AC(K)

IF (AC(K).EQ.0.0) THEN
K 1=INT(K/lOO)
K2=KI '100
AK2=FLOAT(K/100.-KI)
K3=(K 1+ I)"'100 ,
AC(K) = AC(K2)+(AC(K3)-AC(K2))*AK2
WRITE(*.*) K.K1.K2.AK2.K3.AC(K)
A=AB’AC(K)

ENDIF

RETURN
END

C Use subroutine CUBIC as per fiat wall program

 

000000000000 0

g.

169

PROGRAM FOR BINARY MIXTURES on SLITS

Define Psi as a statement function

Use a 10-4 potential where PSI is the negative of the

intermolecular potential.

SIGGSI and SIGGSZ are sigma fluid-wall for components 1 & 2. in Angstroms.
EPSGSI & 2 are the fluid-wall potential in K

EXCI & 2 are the excess in micro mol/m"2

AMTI & 2 are the amount adsorbed in mmol/g

YI & 2 are the compositions of the two components in the adsorbed phase
PI & 2 are the fugacity coefficients for the two components

Fl & 2 are the fugacities

p is the pressure in Mpa

TC, PC and W represent the critical temperature. pressure and omega’s

This program is to calculate the adsorption in a mixture
REAL K12
CHARACTER RFILE*80

PSI1(ETA)=4.0*3.1415926*O.382*SIGGSI*SIGGS1‘EPSGSI*(-O.2
1/ETA" l O-l-O.5/E'I‘A"'"‘4+0.5/(E'I'A+ALPHA)"'"'4+-0.5/(ETA-l-2.0"l
‘ALPHA)"4+O.5/(ETA+3.0*ALPHA)"4+O.5/(ETA+4.0‘ALPHA)
Oit4)

P512(ETA)=4.0*3.1415926‘0.382*51G652‘SIGGS2*EP5GS2*(-O.2
l/ETA“10+0.5/ETA"*4+0.5/(ETA+ALPHA)‘*4+0.5/(ETA+2.0"'
*ALPHA)"4+O.5/(ETA+3.0’ALPHA)‘*4+0.5/(ETA+4.0*ALPHA)
tfit4)

PSI3(XI)=4.0*3. 1415926‘0.382*51GG5 1 "SIGGSI I"EPS GS 1 I"(~02
1/X1**l 0-t~0.5/XI"'"‘4¢l-t'0.S/(XI+ALPHA)"“"4+0.5/(XI+2.O"l
"ALPHA)"4+0.5/(XI+3.0*ALPHA)“4+0.5/(X1+4.0‘ALPHA)
$.14)

PSI4(XI)=4.0*3. 14 l 5926*0.382*SIGGSZ‘SIGGSZ*EPSG52*(-O.2
l/XI”10+0.5/XI"4+0.5/(XI+ALPHA)"4+0.5/(Xl+2.0"
‘ALPHA)"4+0.5/(Xl+3.0*ALPHA)**4+O.5/(Xl+4.0*ALPHA)
OOI4)

PSIA(ET A)=PSI 1 (ET A)+PSI3(XI)
PSIBfETA)=P512(ETA)+PSI4(XI)

SIGGS 1=4.0
SIGGSZ=4.0
ALPHA:3.35/SIGGS l

R = 8.314
TC2=282.4

PC2=5 .04
W2=0.089

 

170

TC1=364.9
PC1=5.49
W1=0.13
K12=0

PC2 = 1000000‘PC2
PCI = IOOOOOO‘PCI
2 PRINT“. 'ENTER T(K). P(MPA),EPSI, EPSZ, WIDTH'
READ“. T.P.EPSGS I .EPSGSZZETA
P=1000000*P
3 PRINT“, 'ENTER Yl'
READ", Y1
EXC1=0.0
EXC2=0.0
AMT1=0.0
AMT2=0.0
TYPE'(X.A.$)'.'ENTER THE FILE NAME FOR RES ULTS :'
ACCEPT'(A)'.RFILE
OPEN(UN1T=2.NA ME=RFILE.TYPE='UNKNOWN')

Y2=1~Y1
YIIN=YI
Y21N=Y2

DDYI = 1E-5
DDY2 = IE-S

TR1=T/TC1
TR2=TITC2

" SET DELTA P FOR PARTIAL DERIVATIVES
C DDP=.01‘P

" CALCULATE COMPRESSIBILITY FACTOR (Z)
“ CALCULATE A1.A2. THEN A

CALL ACALC(R.WI.TC1.TR1.PC1.A11)

CALL ACALC(R.W2.TC2.TR2.PC2,A22)

CALL BCALC(BI,TC1.PC1 .R)

CALL BCALC(B2,TC2.PC2.R)

A12 = SQRT(A1 1‘A22)’( I-K12)

CALL AMIX(YIIN.Y21N,A11.A22,A12.A)

CALL ASTAR(A1 1.A22.A12.A.ASI.AS2.A512.AS.P.R,T)
CALL BMIX(Y 11N.Y21N.BI.B2.B)

CALL BSTAR(B 1 .BZ,B.BS I ,BSZ.BS,P,R.T)

CALL SCUB1C(AS.BS.R1.R2.R3.IFLAG)
IF(IFLAG.EQ. 1) z=R1
IF(IFLAGEQ.2) THEN
z.-.-Rt
IF(R2.GT.RI) Z=R2
END IF
IF(IFLAG.EQ.3) THEN

 

171

CALL ARANGECRI ,RZ.R3)
Z=Rl
END IF

CALL PH(YIIN,Y21N.A51.A512.A52.A5.BSI.BSZ.BS.Z.P1.P2)
C WRITE("',*) 'A & B star ',AS,BS.'FUG’.P1,P2

DENB=PfZlRfT
PBULK=PI 1000000

WRITE(2.*) 'T ='.T.' K'.'P ='. PBULK.‘ MPa'

WRITE(2."') 'BULK DEN ='. DENB,‘ gmol/m"3'.' Width = ',ZETA
WRIT'E(2,"')'Epsilon 1 & 2 ='. EPSGSI .EPSGSZ

WRITE(2,*) 'inital mol. frac. Y1 = ',Yl ,' Y2 = '. Y2

FIBULK=YIIN*P1*P
F28ULK=Y21N*P2*P
PETA=P

C ------------- Start iterating in ETA

 

BET A=ZETA/2.+(ZET A/2.-O.5)/ 100.
DELETA=(ZETA/2.-0.5)/100.

L20

WRITE (2,102)

5 BETA=BETA-DELETA
ET A=(BETA‘2J'SIGFF+SIGWW)/(SIGFF+SIGWW)
XI=(ZETA*SIGFF+SIGWW/2.-BETA"SIGFF)/5IGGS1

"' SET DELTA P FOR PARTIAL DERIVATIVES
DDP=.01*PETA

GOLD=IE5
K=0
H = 0

CALL AZCALC(BETA.A1 1.A1.ZETA)
CALL AZCALC(BETA,A22.A2.ZETA)
A12 = SQRT(A1 l"A2)"'(1-K 12)

F1 = FIBULK‘EXPO’SIAG'TTAVT)
F2 = FZBULK*EXP(PSIB(ETA)/T)
P=Fl/P1+F2/P2

C WRITE("',") 'F1.F2.P'.F1.F2,P

10 CALL AMIX(Y1.Y2.A1,A2.A12,A)

C WRITE(",") 'A='. A
CALL ASTAR(A1.A2.A12.A.ASI.A52.A512.A5.P.R,T)

CALL BM1X(Y 1.Y2.BI.BZ.B)
C WRITE("'.*) 'B=', B

172

CALL BSTAR(B I ,BZ.B,BS l .BSZ.BS.P.R.T)
CALL SCUBIC(AS.BS.RI ,R2,R3.IFLAG)
IF(IFLAG.EQ.1) THEN

ZI=RI

CALL PH(Y1,Y2.ASI.A512,A52.AS.BS 1.852.851] ,P1,P2)
CALCULATE OBJECTIVE FUNCTIONS

G1 = Y1*P1"‘P/Fl-1

02 = Y2'P2’P/F2-l

G3 = l - Y1 - Y2

G0=SQRT(GI"'GI+ G2*G2 + G3*G3)

NROOT=1

GOTO 25
ENDIF

IF(IFLAG.EQ.2) THEN
2] =R 1
Z2=R2
IF(R2.GT.R 1) THEN
ZI =R2
22=R 1
ENDIF
EN DIF

IF(IFI..AG.EQ.3) THEN
CALL ARANGE(R1.R2.R3)
Zl =R 1
22=R3

EN DIF

CALL PH(Y1.Y2,ASI.A512.A52.AS.BS1352.35.21.PIA,P2A)

CALCULATE OBJECTIVE FUNCTIONS

GIA = Y1*P1A*P/Fl-l

G2A = Y2*P2A*P/F2-1

G3A = l - Y1 - Y2

GOA=SQRT(GIA*GI A+ GZA*G2A+ G3A‘G3A)
WRITE("'.")'GIA.G2A'.G1A.GZA

GO=GOA

GI=GIA

GZ=G2A

G3=GBA

2:21

P1=P1A

P2=P2A

NROOT =1

IF(ZZLTBS) THEN
PRINT". 'SKIPPING 22'

GOTO 25
EN DIF
CALL PH(Y 1.Y2,A51.AS12,A52.AS.BSI.BS2.BS,ZZ,PIB,PZB)

CALCULATE OBJECTIVE FUNCTIONS
618 = Y1‘P18‘P/F1-1

173

623 = Y2’P28‘P/F2-1
038 = 1 - Y1 - Y2
GOB=SQRT(GIB*GIB + G2B‘G2B + G3B*G3B)

IF (GOA.GT.GOB) THEN
GO=GOB
G 1:0 1 B
GZ=G2B
G3=G3B
2:22
P I =P I B
P2=PZB I
NROOT=2

ENDIF

25 CONTINUE
PRINT*,'GO=',GOA,GOB ,'NROOT='.NROOT,'IFG='.IFLAG

 

F1CALC=Y1*P1"‘P
F2CALC=Y2"‘P2"P
YT=Y I +Y2
WRITE(","') F1CALC.F1,P
WRITE(*,*) F2CALC,F2,H
PRINT", 'ETA,Y I ',ETA,Y2
IF (K.LT.5.AND.GOLD.LT.GO) THEN
K=K+I
Y I=Y I -DYI
Y2=Y2~DY2
P=P-DP
DY] = DY 1/2
DY2 = DY2/2
DP = D?”
GOTO 100
END IF
K=0
GOLD=GO
9999 FORMAT(4GI$.5.IX,I3)

000

C IF(GO.LT.2E-2) THEN
lF(GO.LT.1E-3) THEN
C WRITE (*,"‘) 'Bubblc Pressure Calculation'
C WRITE (‘,101)
C WRITE (',103) X1.X2,T

DEN=P/Z/8.3 14/1‘

PETA=P

EXC1=SIGGS 1 *(Y I *DEN-Y I IN‘DENB)*DELE'TA* 1 .E-4+EXCI
AMTI=Y1*DEN‘SIGGSI‘DEETA‘988.*1.E-7+AM'T1

EXC2=S 1668 1 *(Y2‘DEN-Y2IN‘DENB)*DELETA *1 .E-4+EXC2
AMT2=Y2"DEN"S 1068 1 ‘DELE’TA‘988J‘I .E-7-I-AM'T2

WRITE (‘,104) Y1,DEN,P,BETA,H
WRITE (2,104) Y1 ,DEN,P,BETA,H
101 FORMAT( I 0X,'X I',lOX,'X2',10X,'T—K')

174

102 FORMAT(4X,'Y I ',7X,'DEN'.9X,'PRESS3,9X,'ETA'.8X,'# ITER')
104 FORMAT(F8.5,2X,ZG12.5.2X.GI2.5.2X,F5.l)

GOTO 200

ENDIF

YY1= Y1 + DDYI
YY2 = Y2 + DDY2
PP=P+ DDP

"' CALCULATE DERIVATIVES WRT P
CALL ASTAR(A1,A2,A12,A,ASI ,ASZAS 12,AS,PP,R,T)
CALL BSTAR(B 1 ,BZ,B,BS 1.BSZ,BS ,PP,R,'I')
CALL SCUBIC(AS ,BS,R 1 ,R2,R3,IFLAG)

IF(IFLAG.EQ.1) THEN
Z=R1
CALL PH(Y 1,Y2.ASI ,ASI2,A52,AS.BSI ,B82,BS .Z,PN1,PN2)
GOTO 30

ENDIF

 

IF(IFLAG.EQ.2) THEN
B IG=R1
S MALL=R2
IF(RI .LT.R2) THEN
BIG=R2
S MALL=R1
ENDIF
R 1 :8 IO
R3=S M ALL
EN DIF
IF(IFI..AG.EQ.3) THEN
CALL ARANGE('RI ,RZ,R3)
ENDIF

IF(NROOTEQ. 1) THEN

Z=R 1

CALL PH(Y 1,Y2,ASI ,AS 12,ASZ,AS,BS I ,BSZ,BS ,Z,PN 1,PN2)
ELSE

Z=R3
CALL PH(Y 1,Y2,ASI,ASI2,ASZ,AS,BSI,BSZ.BSZ.PN1.PN2)
ENDIF

30 D! DP = (Y1*PN1*PP-Y1‘Pl"P)/DDP/Fl
D2DP = (Y 2*PN2*PP-Y2*P2*P)/DDP/F2
DBDP = 0

" CALCULATE DERIVATIVES WRT Y1

CALL AMIX(YY1,Y2,AI,A2,A12,A)

CALL BMIX(Y Y1 ,Y2,B I ,BZ,BV)

CALL AS'TAR(AI,A2,A12,A.AS1.AS2.ASI2,AS,P,R,T)
CALL BSTAR(B 1 ,B2,B.BS 1 ,BS2.BS,P,R,T)

CALL SCUBIC(AS ,BS,RI ,R2,R3,IFLAG)

175

IF(IFLAG.EQ.1) THEN
Z=R1
CALL PH(YYI .Y2.As 1 .As 12,A82,AS,BS 1 352.351.9111 .PN2)
0010 40

ENDIF

IF(IFLAG.EQ.2) THEN
BIG=R1
SMALL=R2
IF(RI.L'T.R2) THEN
BIG=R2
SMALL=R1
ENDIF
RI=BIG
R3=SMALL
EN DIF
IF(IFLAGEQ.3) THEN
CALL ARANGE(R1.R2.R3)
EN DIF

IF(NROOTEQ. 1) THEN

Z=Rl

CALL 911(le ,Y2,AS 1.As 12,A32,As,as1,Bsz,Bs,zJ>N1.1>N2)
ELSE

Z=R3

CALL PH(YYI ,Y2.AS 1 ,As 12,A52.As 351,352.35 .Z,PN1,PN2)
ENDIF

40 D1 D1 = (PNI l"YYI *P-PI‘YI *P)/DDYI/F1
D201 = (PN2'Y2‘P-P2‘Y2‘PVDDY IIF'Z
D3D1 = -1

"' CALCULATE DERIVATIVES WRT Y2

CALL AMIX(Y 1,YY2,A1,A2.A12,A)

CALL BMIX(Y1,YY2,BI,BZ,B)

CALL ASTAR(AI.A2.A12.A.AS1.A52,ASI2,AS.P,R.T)
CALL BSTAR(B 1 ,BZ,B.BS I ,BSZ,BS,P.R.T)

CALL SCUBIC(AS ,BS,R 1 ,RZ,R3,IFLAG)

IF(IF'LAGEQJ) THEN
Z=R1
CALL PH(Y I ,YY2.AS 1,AS 12.A82,AS,BSI ,BSZ,BS,Z,PN1.PN2)
GOTO 50

ENDIF

IF(IFLAG.EQ.2) THEN

BIG=R I

SMALLzRZ

IF(R1.LT.R2) THEN
BIG=R2
SMALL=RI

ENDIF

R 1=BIG

 

176

R3=SMALL
ENDIF

IF(IFLAGEQ.3) THEN
CALL ARANGE(R1,RZ,R3)
EN DIF

IF(NROOTEQ. 1) THEN

Z=R1

CALL PH(Y 1,YY2,AS 1 ,AS l2.A52.AS ,BS 1 ,BSZ,BS Z.PN1,PN2)
ELSE

Z=R3

CALL PH(Y 1 .YY2,ASI ,A812.A82.AS ,BS 1 ,B82,BS .Z,PN1.PN2)
ENDIF

 

50 D1D2 = (PNI *YI l"P-PI I"Y1 'P)/DDY2/FI
D2D2 = (PN2‘YY2‘P-P2’Y2'PVDDY2/F2
D3D2 = -1

* CALCULATE INCREMENTS
CALL INV(D 1 DP,DZDP.D3DP,DI D1 ,DZDI,D3D1,D1D2,D2D2.D3D2,DET)

DY1 = -GI"DID1 - G2‘D1D2 - G3‘D1DP
DY2 = -Gl"D2D1 - GZ‘D2D2 - G3‘D2DP
DP = -Gl"‘D3D1 - 62*D3D2 - G3‘DBDP
PRINT‘, ETA,Y1,DY1,P,DP
IF(ABS(DY1).GT.YI) THEN
PRIN'T*,'WARNING, LARGE DY',DYI
DYI=0.5*Y1*DY1/ABS(DY1)
DY2=-DY1
C P=9.E7
C DP=0
END IF
100 Y1 = Y1 + DY1
Y2 = Y2 + DY2
P = P «1» DP
IF(Y1.LT.0)Y1=0
IF(Y1.GT.1)Y1=1.
IF(Y2.LT.0)Y2=0
IF(Y2.GT. 1)Y2= I .

H=H+1
IF (H.LT.200) GO TO 10
PRINT", 'ITERATIONS EXCEEDED'
C DDYI=1E-3
C DDY2=1E-3
C DDP=.01"P
C GOTO 10
200 CONTINUE
L==L+l
IF(L.LT.101) GOTO 5
WRITE(",") EXC1,EXC2,DENB,AMT1,AMT2

177

WRITE(2,"') EXC1,',',EXC2,',',DENB
WRITE(2,"') AMT 1.',',AMT2

PR1NT*,'ENTER 1 FOR SAME T, DIFFERENT x1'
PRINT", 'ENTER 2 FOR NEW T AND NEW x1'
PRINT‘, ENTER 0 TO QUIT

READ", IDB

IF (IDB.EQ.1) GO To 3

IF (IDBEQ.2) GO To 2

END

.**ilfitfifitttfififitltil*Ititltlifi*ltt***********************¥*************

" SUBROUTINES BELOW 1!!!

*

SUBROUTINE ACALC(R,W,TC,TR,PC.AN)

FW = 0.37464 + 1.54226‘W - .26992 " W**2

AN = 0.45724'R“2"TC"2"(1+FW*(1-SQRT(TR)))"2/PC
RETURN
END

SUBROUTINE AMIX(MI ,M2.A1.A2.A12,AA)
REAL AA, A1, A2, A12, M1, M2

AA = M1"2"A1 + 2*M1‘M2‘A12 + M2"2 * A2
RETURN
END

SUBROUTINE ASTAR(AI ,A2.A12,A,AS 1 ,AS2.AS l2.AS .P,R,T)
AS] = AI‘P/(R"2*T”2)
A52 = A2'P/(R"2*T”2)
A812 == A12’P/(R"2*T“2)
AS = A‘P/(R"2"T"2)
RETURN
END

 

SUBROUTINE BCALC(B,TC.PC,R)
B = 0.07780‘R'TC/PC

RETURN

END

SUBROUTINE BMIX(M1,M2,B 1 ,B2,BB)
REAL BB, Bl. B2, M1, M2
BB = M1’BI+ M2*B2
RETURN
END

SUBROUTINE BSTAR(BI,132,B.Bs1 ,BSZ,BS,P,R,T)
1351 = Bl‘P/(R‘T)
1332 = BZ‘P/(R‘T)
as = 3*9/(Rm

RETURN

END

 

SUBROUTINE TM2(BS,T2)
T2 = BS - 1
WRITE (‘3') T2

178

RETURN
END
SUBROUTINE TM1(AS,BS,T1)
T1 = AS - BS"(2+3"BS)
WRITE (‘,‘) T1
RETURN
END
SUBROUTINE TMO(AS,BS,T0)
'11) = BS“(BS"'*2 + BS - AS)
WRITE(*,“‘) TO
RETURN
END

 

SUBROUTINE PH(MI ,M2,AS I ,A812A82,AS,BS 1 ,Bsz,Bs,z,I>1,I>2)
REAL M1,M2
Q = M1’ASI+ M2*A512
CALL PHICLC(BS l ,AS.BS.Z.Q.P1)
Q = Ml‘ASlZ + M2-Asz
CALL PHICLCCBSZ,AS,BS.Z.Q.P2)
RETURN
END

 

SUBROUTINE PHICLC(B,AS,BS,Z.Q.PHI)
TERMI = (BlBS)*(Z-1)

TERMIA = (2*Q/AS) - (BIBS)
TERMIB = 2"Q/AS
TERM2 = (2*Q/AS - B/BS)"'(AS/(2*SQRT(2.0)"BS))

TERM3 = I/(Z-BS)

TERM4 = ((Z-(SQRT(2.0)-1)*BS)/(Z+(SQRT(2.0)+ 1 )"‘BS ))
WRTTE(*.") BS, TERMIA, TERMIB
WRITE(*,") TERM2, TERM3, TERM4

PHI = TERM3*(TERM4“TERM2) I"EXP('I'ERM I)
RETURN
END

 

SUBROUTINE SCUBIC(ASBS.RI.RZ.R3.IFLAG)
CALL TM2(BS,T2)
CALL TM1(AS,BS,T1)
CALL TMO(AS,BS,'m)
WRITE (*,"') T2,TI ,TO
CALL CUBICCI‘2,TI,T0,RI,R2.R3,C1,C2,C3,IF'LAG)
WRITE(",*) IFLAG
WRITE("',") 'R1 R2 R3'
WRITE("',") R1,R2,R3
RETURN
END

179

 

C
SUBROUTINE ARANGE(R1,RZ,R3)

C PROGRAM TO PUT 3 NUMBERS IN DESCENDING ORDER
DO 20 1:13 ’

IF(R2.GT.R 1) THEN
TEMP=R I
R I =R2
R2=TEMP

ENDIF

IF(R3.GT.R2) THEN
TEMkR2
R2=R3
R3=TEMP

ENDIF

 

20 CONTINUE

RETURN
END
C

 

SUBROUTINE AZCALC(BETA.AB,A.ZET A)
SA=ZETA.0.5-BETA

IF (BETA.LE.I.5) THEN
A 1=BETA+0.5
A2=I.-I./SA*"'3
A3=I./3."'A2
A=AB *(AI+A3)"3./8.

END IF

IF(BETA.GT.(ZETA-1.5)) THEN
AI =SA+1.
A2=1.-1J(BETA-0.5)"3
A3=A2/3.
A=AB“(A1+A3)*318.

END IF

IF((BETA.GT .1 .5).AND.(BETA.LE.(ZETA-I .5))) THEN
A1=7./3.-1 l3./(BETA—0.5)"3
A2=1.-1./SA“'*3
A3=A2/3.
A=AB*(A1+A3)*3./8.
END IF

RETURN
END

 

C

180

SUBROUTINE INV(DIDP,D2DP,D3DP,D1D1,D2D1,D3D1,DID2,D2D2,D3D2,DET)
C011 = D2D2’D3DP-D3D2‘D2DP
C012 = -(D2D1*D3DP-D2DP*D3DI)
C013 = D2D1*D3DZ-D3D1"D2D2

C021 = -(D1D2*D3DP-D3D2*D1DP)
C022 = D1DI*D3DP-D3D1"'D1DP
C023 = -(D1 DI *D3D2-D3DI *D1 D2)

C031 = D1D2’D2DP-D2D2’DI DP
C032 = -(D1D1*D2DP-D2D1 I"DI DP)
C033 = D1D1’D2D2-D2D1'D1D2

DET: D1DI'C011+ D1D2‘C012 + DIDP‘COI3
WRITE("',*) 'DET'
WRITE("',*) DET

DIDI = C01 1/DET

DID2 = C021/DET

D1DP= CO31/DET

DZDI = C012/DET
D2D2 = C022/DET
D2DP = CO32/DET

D3Dl = COI3/DET
D3D2 = COB/DET
D3DP = C033/DET
RETURN
END

Use Subroutinc CUBIC as per other

 

 

C
C
C

15

17

181

ACALIDOC
This is a program to calculate the value of the vanderwaals 'a'
in a cylindrical pore for T] varies from R-3 sigma/2 to R-sigma/‘Z

IMPLICIT DOUBLE PRECISION(A-Z)
INTEGER I,J,I(

WRITE(*,"‘) 'ETAI',' SIGMA'.' PORE SIZE IN NUMBER OF SIGMA'

READ("'.*) ET A1,SIGMA.P
R1=SIGMA *ETAI
=SIGMA'P

WRITE(*,*) R,Rl
CASE 1 rl varies from R-3 sigma/2 to R-Sigma/Z

A 1 =0.0
P1=3.1415926535897932

z1=(SIOMA*SIOMA-cR-SIGMA/2.-R1)"2)
WRITE(*,") 21
ZI=DSQRT(ZI)
WRITE(*,*) 21
DO 10 1:0.1000
bl'Zl/IOOO.
WRITE(*,*) z
B=Rl *RI+SIOMA*SIGMA-z*z
B I=(R-SIGMA/2.)*(R-SIGMA/2.)
32=SIGMA*SIGMA.z*z
WRITE(*,*) 132
IF(BZ.EQ.0.0) THEN
WRITE(*,*) '32 EQ 0'
GO TO 15
ENDIF

B3=DSQRT(B2)
A=(B-B I )flJRI/B3
WRITE(*,"‘) A
IF(A.GT.1.0) THEN
WRI'I'E(","') 'A.GT.1.0 ',Z,A
A=1 .0
ENDIF
IF(A.LE.-1.0) THEN
WRITE(*,*) Z,A
=—I .0
ENDIF
F=DACOS(A)
GO TO 17
=PI

A1=F’DSQRT(SIGMA*SIGMA-(R-SIGMA/2.-R1)**2)/1000.+A1

10 CONTINUE

A1=Al/4./(SIGMA**4)

 

182

WRITE(*,"') 'A1= ',AI

DO 20 J=0.10(X)
ZI =DSQRT(S IGMA‘SIGMA-(R-SIG MHZ-R I )**2)
Z=J*ZI/l 000.
D0 30 K=0,1000
B=RI *R 1+SIGMA"'S IGMA-Z‘Z
B 1=(R-SIGMA/2.)*(R-SIGMA/2.)
BkSIGMA‘SIGMA-l‘z
IF(BZ.EQ.0.0) GO TO 25
B3=DSQRT(BZ)
A=(B-B ”fl/R 1/33
IF(A.GT.I .0) THEN
WRITE(*,"') 'AA GT 1.0 'Z.A
A=I .0
ENDIF
IF(A.LT.-I.0) THEN
WRITE(*,*) Z, A
=-1.0
ENDIF
GO TO 27
25 =-I.0
27 THETA=DACOS (A)*K/I (XX).
C WRITE(*.*) Z, A,THETA

C=Rl *DCOSCT HET A)
C1=(R-SIGMAI2.)"(R-SIGMA/2.)

C2=R1 I"R1"DSIN('THE’I‘A)"‘DSIN('T1-IE'I‘A)
C3=DSQRT(C 1 -C2)

C4=C+C3

F I =C4 *C4+Z"Z

F=I JF III: I

AA=F"'ZI/1000."DACOS(A)/I000.+AA

30 CONTINUE
20 CONTINUE

AA=AA/4.
WRITE(*,*) 'AA= ',AA

=DSQRT(SIGMA*SIGMA-(R-SIGMA/Z-R1)"'*2)

DO SO K=0.IOOO
THETA=P1"K/1000.
C=R1 ‘DCOS (THET A)
C 1=(R-SIGMAf2.)"(R-SIGMA/2.)
C2=Rl"Rl"DSIN('1'I-IEI‘A)*DSIN(TI-1ETA)
C3=DSQRT(CI -C2)
C4=C+C3
C5=C4*C4

F1=SIGMA/2JC5/(C5+SIGMA"SIGMA)
F2=Yl2JC5/(C5+Y"Y)

 

183

F3=DATAN(SIGMA/C4)/2./CS/C4
F4=DATAN(Y/C4)/2./C5/C4
F=F1-F2+F3-F4
AB=F*PI/1000.+AB

50 CONTINUE

AB=PI*(SIGMA-Y)/4./SIGMA**4-AB/4.
WRITE(*,‘) 'AB = ',AB

C PART 21.e. 2 going from sigma to infinity

C

D0401=0,1000

THETA=I"'Pl/1000.
B=R l *DCOSCI‘HETA)
B 1=(R-SIGMA/2.)*(R-SIGMA/2.)
Bth’Rl *DSINCI'HETA)’DSIN(THETA)
83=DSQRT(Bl-82)
A=B+B3
F1=PI/4./NA/A
F2=SIGMNZJAIAKA *A+SIGMA’SIGMA)
F3=DATAN(SIGMA/A)/2JAIA/A
=Fl-F2-F3
WRITE(*,*) Fl ,F2.F3

A2=F‘PI/1000.+A2

40 CONTINUE

A2=A2/4.
A3=PI/1ZJSIGMA/SIGMNSIGMA-A2

WRITE(*,*) A2,A3
A4=A1-AA+AB+A3
A5=A4*3.*SIGMA‘SIGMA‘SIGMA/PI
WRITE("'."‘) 'A/ABULK= '. A5.A4

END

_
‘

C
C
C

C
C
C

184

ACALZJ:
This is a program to calculate the value of the vandcrwaals 'a’
in a cylindrical pore rl varies from 0 to R-3sigma/2

IMPLICIT DOUBLE PRECISION(A-Z)
INTEGER I

WRITE(*,*) 'ETAI'.’ SIGMA'.’ PORE SIZE IN NUMBER OF SIGMA'
READ(*,*) ET AI,SIGMA.P

R1=SIGMA"'ETA1

R=SIGMA‘P

WRITE(*,*) R,RI
CASE 2 rl varies from 0 to R-3sigma/Z

A1=0.0
PI=3.1415926535897932

D0 10 1:0.1000
THE'TA=I"'PI/1000.
B=R1 *DCOSCTHET A)
B 1=(R-S IGMA/Z.)‘(R-SIGMA/2.)
B2=R 1 *RI *DSIN(T HETA)‘DSIN(THET A)
B3=DSQRT(B 1-B2)
A=B+B3
F1=A*A"A
F=1./FI
WRTTE("',") THETA .A,F
F3=2.*A *A‘A
F=Fl +F2/F3
. A1 =F‘PI/1000.+A I

10 CONTINUE

A2=1./3.-SIGMA*SIGMA*SIGMA/16.*A1
A3=A2*3.
WRITE("‘,"') 'A/ABULK ', A3.A1,A2

END

 

 

 

C
C
C

185

ACAL3.F
This is a program to calculate the value of the vandcrwaals 'a'
in a cylindrical pore r1 = R-sigma/Z

REAL ET A1 ,SIGMA,P.R,R1 .A,A 1 .A2

WRTTE(*."') 'E'TAI'. ' SIGMA',’ PORE SIZE IN NUMBER OF SIGMA'
READ("',*) ETA1,SIGMA,P
R1=SIGMA"E'TA1

=SIGMA‘P

WRITE(*."‘) R,Rl
CASE 3 T] = R-sigma/Z

A1=0.0
PI=3. 1415927

DO 101:0.100
Z=I"'SlGMA/IOO.
B=SIGMA‘SIGMA-Z‘Z
Bl=SQRT(B)
THETA=ACOS(Bl/2./Rl)
F=TI-1ETA/4./SIGMA"4

WRI'I'E(*.‘) Z.A.F

A1=F"SIGMAI100.+A1

10 CONTINUE

WRITE("',"‘) 'A1= ',Al

DO 20 1:1,100
Z=J*SIGMA/100.
B=SlGMA’SlGMA-Z‘Z
Bl=SQRT(B)
THETA=ACOS(Bl/2./Rl)

C=2."'RI*R1+Z"'Z
Cl=4."Rl*Rl+Z"‘Z
C2=CI M1 .5
C3=ClC2flflfl
write(*,"') C,CI ,C2,C3
IF(Z.EQ.SIGMA) GO TO 25
C4=Z*TAN(THETA)/SQRT(CI)
C5=ATAN(C4)
GO TO 27

C5=PI/2.

C6=C3"'C5

write(*,*) C3.C5.C6

D1=R1*R1*SIN(2‘THETA)
D2=Z"'Z*CI*(2.*R1*RI+Z*Z+2.*RI I"R1 *COS(2."‘THETA))

I86

D3=Dl/D2
C WRITE("',"') D1,D2,D3

F=C6~D3

A2=F‘SIGMA/100.+A2
C WRITE(*,"') Z, A,THET A

20 CONTINUE
A2=A2/4.

AA=A 1-A2
WRTTE("',"') 'A2= ',A2,‘ AA: ',AA

C PART 2 i.e. 2 going from sigma to infinity

DO 40 1:0,2000
Z=SIGMA+I*SIGMA/ 100.
B=R1*RI+Z*Z
BI=4.“R1*R1+Z*Z
82:81 “1.5

FEB/82m
A3=F*SIGMA/100.+A3

40 CONTINUE
A3=A3*PI/8.
A4=PII24JS TOMA/S IGMA/SIGMA
A5=A4-A3

WRITE(*,"‘) A3,A4,A5

A6=AA+A5
A=A6*SIGMA‘SIGMA‘SIGMA'BJPI
WRITE(*,*) 'A/ABULK: ', A6.A

END

 

 

 

Appendix 3 - Properties of Gases

 

'Table: Appendix 3 -1 Supercritical fihylene

 

40 F
P(Psia)
500
600
700
800
900
1000
1 100
1200
1300
1400
1500

V(f13llb)

0.26035
0.18181
0.04754
0.04511
0.04368
0.04267
0.04187
0.04122
0.04067
0.04019
0.03977

P(MPa)

3.447367
4.1 3684
4. 82631 3
5.51 5787
6.20526
6.894733
7. 584207
8.27368
8. 9631 53
9.652627
1 0.3421

V(cm3lmol V (PR)

455.9669
318.415
83.2597
79.0039

76.49945

74.73058

73.32949
72.191 1

71 .2785
70.3872

69.65163

443.6109J
307.7956
93.08558
85.36646

81 .0538
78.03798I
75.72067
73.84304
72.26802
70.91416
69.72904

 

160 F
P(Psia)
500
600
700
800
900
1 000
1 100
1200
1 300
1400
1 500

 

V(ft3llb)

0.29183
0.21835
0.16006
0.10129
0.05804
0.05119
0.04821
0.04637
0.04506
0.04404
0.04322

P(MPa)

3.447367
4. 1 3684
4.82631 3
5. 51 5787
6.20526
6.894733
7.584207
8.27368
8.9631 53
9.652627
10.3421

V(cm3/mol v (PR)

51 1.0997
382.4097
280.3229
177.3954

101.649
89.65217
84.43312
81.21061
78.91633
77.12994
75.69383

498.076
371.1377 ‘
271.7131

175.101
109.1489!

95.9797
89.38151

85.037
81.82833
79.3
77.225

 

 

187

 

188

 

 

Table: Append"ix_3_-1 SLupereriticT-I Ethylene

 

80 F
P(Psia)
500
600
700
800
900
1 000
1 100
1 200
1 300
1400
1 500

V(ft3llb)
0.31926
0.24619
0.19172
0.14789
0.1 1013
0.08024
0.06464
0.05726
0.05314
0.05049

0.0486

100 F
P(Psia)
500
600
700
800
900
1 000
1 100
1 200
1 300
1400
1 500

V(ft3/lb)
0.34426
0.27016
0.21611
0.17436
0.14072
0.1 1315
0.09162

0.079677
0.06737
0.06131

0.0572

 

P(MPa)

3.447367
4. 13684
4. 82631 3
5.51 5787
6.20526
6.894733
7.584207
8.27368
8.9631 53
9.652627
1 0.3421

P(MPa)

3.447367
4. 1 3684
4.82631 3
5.51 5787
6.20526
6.894733
7.584207
8.27368
8.9631 53
9.652627
1 0.3421

V(cm3lmol V (PR)

559.1 395
431 . 1676
335. 7709
259.0088
1 92.8774
140.5292

1 13.208
1 00.2829
93. 06733
88.42622
85. 1 161 5

V(cm3/mol V (PR)

602.9236
473.1477
376.4666
305.3673
246.4515
196.1665
160.4597
139.5432
117.9892
107.3759
100.1779

545.9039T
419.54
326.25
252.5661
191.167
143.56
117.688
104.796
97.089I
91.807
87.86828.

589.77
461.53
368.65
297.77
241.692
196.747
161.9475
137.3667
121.2451
1 10.5723
103.106

 

 

189

 

TT‘abIe: Appendix 32 Supercritical Methane

 

0 F
P(Psia)
400
500
600
700
800
900
1 000
1 100
1200
1400
1 600
1 800

120 F
P(Psia)
400
500
600
700
800
900
1 000
1 100
1200
1400
1600
1800

 

V(fl3/|b)

0.70345
0.54934
0.44659
0.37323
0.31829
0.27569

0.2418
0.21431

0.1917
0.15722
0.13293
0.11561

V(ft3llb)

0.9341
0.74051
0.61 165
0.51978
0.45106
0.39777

0.3553

0.3207
0.29202
0.24737
0.21439
0.18923

P(M Pa)
2.757893
3.447367

4.13684
4.82631 3
5.515787

6.20526
6.894733
7.584207

8.27368
9.652627
1 1.03157
12.41052

P(MPa)

2.757893
3.447367

4.13684
4.826313
5.51 5787

6.20526
6.894733
7.584207

8.27368
9.652627
1 1.03157
12.41052

702.642

548.709

446.077
372.8013
317.9244
275.3733
241 .523
214.0638
191 .4798
157.0394
132.7773
1 15.4772

933.0271
739.6594
610.9475

519.183
450.5419
397.3131
354.8919
320.3316
291.6846

247.0859 .

214.1438
189.0126

V(cm3lmol v (PR)

602.89
542.64
503.81
476.417
455.814
439.614
426.457
415.5051
406.21 1
391.2046
379.533
370.137

V(cm3/mol V (PR)

868.8
738.296
657.42
602.41 8
562.4746
532.0424
508
488.47
472.25
446.74
427.5

41 2.38

1

 

 

 

190

 

mi. .23 83.2- 8686 235° 835... 5.35 0: «38 3% 8m
53. 83. ER? $58... 82.3 3.355 934.8 2.5 $.68 83 new
22 9:” 59.3- 983 $236 «323 «82.3. 8.8 ~33 «3 8m
Sod so.” 883- $385 8°86 3836 «2.3.9. 8.3 9.8..” 3:. 8.3.
new 8E 3.8%.? «$83 333° $28.0 88.8 $6.. 8.8... «2.» EN
53 B? 583 6233 388.0 898.6 383.3. e.g. 8&8 33 8a
83 23 «833 3386 858.0 9.886 9.5.33. 8.; .38 a: 8w

«is. was. 6:6 9. ea 6.8 5.2820 5230 66.2 x

.62“. >03“. 33. 633 P £22. £966 use: mOe<> Ema L.

3 BERN 2820

32 was. $3 "on.

293 v. Neon no»

 

 

 

80 a-” 52.2.2 3:3.

 

 

 

191

«Nomad 505 Ed rumored momumdm
aovuwmd 333.0 3:. 3.0 mam—him
wwnommd mvmm _.o.o «cam 36 Em 30. Pm
womommd mmm..~o.o Nwmomod vmmoflv

nonnmmd omvmmod «mammod ammodv

33.006 mmmowod mmmmwod 32.95

932.6 888.0 «8.4.8.0 $2.13
98 0.8
t 2.28 Emcee
9.6.8
3R6
~83

nmdm
omdm
05.3
mode
ms. _.v
mean
2.5

4023.20 4023.20
050...

”08.. $06
mad: vomé
Nmémw mmmé
8.8.. mom...”
«Yew word
mmdwm mmod
mmfimo rum...
mOm<>

was. wmwé

x mdmw

._.<wa

momooé.

m6:
a:
mg
m3
0!.
on?
our

"(Ow—20

 

 

Em; 1n xicona< 8381...

 

 

192

000.0 000.0 0000.0 P- 5 00000 0005000 N000000 080.wa ‘va 0.000 05.0 0 P .000
.00... 00% 5503.? 3053.0 .0 .0000 0000000 0000.005 0.02. 50.050 000.0 0 _. .00
30... 30... 0000.30 50050.0 0000 50.0 .000 «0.0 0050.50 00.50 0.00:. .05.? 0 _. .000
50010 PONv0 00033.0 00005.0 0000 3.0 V000 50.0 00 500.00 05.05 0.05Nv 35v0 0 _. .050

was. was. 6:6 .5 96 6.8 6.2520 40.2320 662 v.
40:... >03“. .9 £23 £28 050... «92> Ewe h
58.3. 3mm? .20sz
$26 was. 9.? "on.
3mm... x 38 1.9.

 

 

 

2.2.2.. m-» 5.52:? 333.

 

193

 

 

 

 

00.0
500.0
500.0
000.0
505...
1050.?
00000

was.
40:“.

..00.0
500.0
000.0
000.0
505. ..
050. F
5000.0

was.
>03".

0005.0?
53.0..-

5005?.5-
00000.0-
00055.0
030.50
0050? ...0

.25 .x.

0000000
000.000
0 90050.0
00030.0
00 .000

33000
0000050

h...

0000.00
5000 3.0
0000 3.0
5505 30
05.0 .00
550.000
00.0000

axm
36:2.

00.. .. .00
3.030
00:00
0050.00
00.00
0000000
0V0v000

0.8
£23

0000 _. .05
50000.05
00500.00
0003.50
50VVO. .0
303.5?
..000.vv

0'
5500.0
2.0.0

00.00
3.00
0.00
00.00
00.00
006v
V0. 3.

402.020 40.20.20
050...

0100.. 000 .v
v0.03. ..00.v
0100.. 000.0
50.000 F0 .0
00.000 0V0 .0
P0 .0 F0 000. ..
50.050 0.0. _.
mr...)—
mOn.<> .500
0000 .. .e
0&2 V006.
v. 0.00..

00..

00 _.
00 _.
00 ..
05 ..
00 _.
00 _.

v.

...

"(Gm—20

non.

no...

 

 

2.2.8:. 0.0 xwucona.‘ no.0:

 

 

194

 

 

 

Fable: Appendix 3-7 Ethylene

 

TC= 282.4 K

PC= MPa
OMEGA=

T density density Tr
K calc exp

269 0.012189 0.012938 0.95255
255 0.014574 0.014641 0.902975
241 0.016372 0.015908 0.853399
225 0.018031 0.017078 0.796742
211 0.019248 0.017965 0.747167

 

°/o error

-5.78915
-0.45762
2.916771
5.580279
7.141664

 

 

_A_ - Am.-

 

% Error in Density

195

 

 

 

 

 

 

 

 

 

 

 

15
10 o-
5 4r-
0 --
-5 -—
— Ethylene ‘::\ ‘\
-10 -- --— Methane ‘~\ \
------ Propane \
- . _ . Argon A“ \
- - - - C02 ‘
45 ’r % i i [L
0.7 0.75 0.8 0.85 0.9 0.95
T

Figure Appendix 3-1: Error in the Peng-
Robinson Calculation of Saturation Density

LIST OF REFERENCES

 

 

LIST OF REFERENCES

Abraham, F. F .; Singh, Y. J. Chem. Phys. 1978, 68, 4767.
Angus, 8.; Armstrong, B.; de Reuck, K. M.; Featherstone, W.; Gibson, M. R. Eds.,
‘ International Union of Pure and Applied Chemistry; Butterworths: London, 1974,
39.
Asselineau, L.; Bogdanic G.; Vidal, J ., Fluid Phase Equil. 1979, 3, 273.
Barrer, R. M.; Robins, A. 3., Trans. Faraday Soc. 1951, 47, 773.

Barton, S. S.; Dacey, J. R.; Quinn, D. F. In Fundamentals of Adsorption, Myers, A. L.;
Belfort, G. Eds., Engineering Foundation: New York, 1983, 65.

Berenyi, L. Z. Phys. Chem, 1923, 105, 55.

Bienfait, M. In Phase Transitions in Surface Films, Dash, J. G.; Ruvalds, J. Eds., Plenum
Press: New York, 1980, 29.

Brennecke, J. F.; Eckert, C. A. In Supercritical Fluid Science and Technology; Johnston,
K. P.; Penninger, M. L. Eds., ACS Symposium Series 406, American Chemical
Society: Washington DC, 1989,14.

Brennecke, J. F. Ph. D. Thesis, University of Notre Dame, 1989.

Brennecke, J. F.; Tomasko, D. L.: Peshkin, J .; Eckert, C. A. Ind. Eng. Chem. Res. 1990,
29, 1682.

Brunauer, S. The Adsorption of Gases and Vapors Vol. 1, Princeton University. Press,
Princeton: New Jersey, 1945.

Davis, H. T.; Scriven, L. E. In Advances in Chemical Physics, Prigogine, 1., Rice, S. A.,
Wiley: New York, 1982, 49, 357.

196

 

 

197

Debenedetti, P. G. Chem. Eng. Sci. 1987, 42, 2203.
de Boer, J. H. The Dynamical Character of Adsorption, Clarendon Press: Oxford, 1953.
Defay, L.; Prigogine, 1., Surface Tension and Adsorption, John Wiley, New York, 1966.

Denbigh, K. The Principles of Chemical Equilibrium, 4th ed., Cambridge University
Press: Cambridge, 1981.

Dubinin, M. M. Proc. Acad. Sci. USSR. 1947, 55, 331.
Dubinin, M. M. Progr. Surf Membrane Sci. 1975, 9, 1.
Eckert, C. A., Ziger, D. H.; Johnston, K. P.; Ellison, T. K. J. Phys. Chem. 1986, 90, 2798.
Evans, R. Advances in Physics 1979, 28, 143.

Findenegg, G. H. In Fundamentals of Adsorption, Myers, A. L., Belfort, G. Eds.,
Engineering Foundation: New York, 1983, 207.

Fisher, B. B., McMillan, W. G. J. Am. Chem. Soc., 1957, 79, 2969. '
Fisher, J .; Methfessel, M. Phys. Rev. A. 1980, 22, 2836.

Gregg, S. J .; Sing, K. S. W. Adsorption, Surface Area and Porosity, Academic Press:
London, 1982.

Guggenheim, Thermodynamics, 5th ed., North-Holland Publishing: Amsterdam, 1967.

Henderson, D. Fundamentals of Inhomogeneous Fluids; Marcel Dekker: New York,
1992.

Hill, T. L. J. Chem. Phys. 1946, I4, 441.
Hill, T. L. J. Chem. Phys. 1947, 15, 767.
Hill, T. L. J. Chem. Phys. 1948, 16, 181.
Hill, T. L. J. Phys. Chem. 1950, 54, 1186.
Hill, T. L. J. Chem. Phys. 1951, I9, 261.

Hill, T. L. In Advances in Catalysis, Frankenburg, W. G.; Rideal, E. K.; Komarewsky,
V. I. Eds., Academic Press: New York, 1952, 4, 211.

 

 

198

Johnston, K. P.: Kim, S.; Combs, J. In Supercritical Fluid Science and Technology;
Johnston, K. P.; Penninger, M. L. Eds. ACS Symp. Ser. 406, American Chemical
Society: Washington DC, 1989.

Kierlik E.; Rosinberg, M. L. Phys. Rev. A. 1990, 42, 3382.

Kierlik E.; Rosinberg, M. L. Phys. Rev. A. 1991, 44, 5025.

Kim, S.; Johnston, K. P. Ind. Eng. Chem. Res. 1987a, 26, 1206.

Kim, S.; Johnston, K. P. AIChE J. 1987b, 33, 1603.

Kung, W. C.; Scriven, L. E., Davis, H. T. Chem. Phys. 1990, 149, 141.

Lee, L. L. Molecular Thermodynamics of Non-Ideal Fluids; Butterworths: Stoneham,
MA, 1988.

Lee, L. L.; Debenedetti, P. G.; Cochran, H. D. In Supercritical Fluid Technology:
Reviews in Modern Theory and Applications, Bruno, T. J .; Ely, J. F. Eds., CRC
Press: Boca Raton, F1, 1991, 193.

Lennard-Jones, J. B.; Devonshire A. F. Proc. Royal Soc. A 1937, I63, 132.

Maslan. F. D.; Altman, M.; Abereth, E. R J. Phys. Chem. 1953, 57, 106.

McQuarrie, D. A. Statistical Mechanics; Harper & Row: New York, 1976.

Myers, A. L.; Prausnitz, J. M. AIChE J. 1965, 11, 121.

Nicholson, D.; Parsonage, N. G. Computer Simulation and Statistical Mechanics of
Adsorption, Academic Press: New York, 1982.

Peng, D. Y.; Robinson, D. B. Ind. Eng. Chem. Fund. 1976, 15, 59.
Polyani, M. Trans. Faraday Soc. 1932, 28, 316.
Powles, J. G.; Rickvayzen, G.; Williams M. L., Mol. Phys. 1988, 64, 33.

Prausnitz, J. M. Proceedings of the 2nd International Conference, EFCE, DECHEMA:
Berlin, 1980. ‘

Prausnitz, J. M., Liechtenthaler, R. N., de szedo, E. G. Molecular Thermodynamics of
Fluid Phase Equilibria, Prentice-Hall, Englewood-Clifi‘s: New Jersey, 1986.

Pyada, H. Masters Thesis, Michigan State University, 1994.

 

199

Radke C. J., Prausnitz, J. M. AIChE J. 1972, 18, 761.

Rangarajan, B.; Lira, C. T. In Better Ceramics Through Chemistry V, Hampden-Smith,
M. J ., Kemperer, W. G., Brinker, C. J. Eds. , Materials Research Society: Pittsburgh,
PA, 1992, 559.

Rangarajan, B.; Lira, C. T.; Subramanian, R. AIChE J. 1995, 41, 838.

Reich, R., Ziegler, W. T., Rogers, K. A. Ind. Eng. Chem. Proc. Des. Dev.l980, 19, 336.

Reid, R. C.; Prausnitz, J. M. Poling, B. The Properties of Gases and Liquids, 4th ed,
McGraw-Hill: New York, 1987.

Ross, S., Clark H., J. Am. Chem. Soc. 1954, 76, 4291.
Ross, S.; Olivier, J. P. On Physical Adsorption, Interscience Publishers: New York, 1964.

Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity, Clarendon Press: Oxford,
1982.

Saam, W. F.; Ebner, C. Phys. Rev. A. 1978, 17, 1768.

Sandler, S. L. Chemical and Engineering Thermoaynamics, 2nd ed., Wiley: New York,
1989.

Snook, I. K.; Henderson, D. J. Chem. Phys. 1978, 68, 2134.
Specovius J .; Findenegg, G. H. Ber. Bunsenges Phys. Chem. 1978, 82, 174.
Specovius J .; Findenegg, G. H. Ber. Bunsenges Phys. Chem. 1980, 84, 690.

Steele, W. A. In The Solid-Gas Interface Vol 1, Flood, E. A. Ed., Marcel Dekker: New
York, 1967, 307.

Steele, W. A. The Interaction of Gases with Solid Surfaces, Pergamon Press: Oxford,
1974.

Strubinger, J. R.; Parcher, J. F. Anal. Chem. 1989, 61, 951.
Subramanian, R.; Pyada H.; Lira, C. T. Ind. Eng. Chem. Res. (in press).
Sullivan, D. E. Phys Rev. B. 1979, 20, 3991.

Tan, C. -S.; Liou, D. -C. Ind. Eng. Chem. Res. 1988, 27, 988.

200

Tan, C. -S.; Liou, D. -C. Ind. Eng. Chem. Res. 1990, 29, 1412.

Tan, 2.; Gubbins, K. E. J. Phys. Chem. 1990, 94, 6061.

Tarazona, P.; Evans, R. M01. Phys. 1983, 48, 799.

Teletzke, G. F.; Scriven, L. E.; Davis, H. T. J. Colloid Interface Sci. 1982a, 87, 550.
Teletzke, G. F.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1982b, 77, 5794.

Telo de Gama, M. M.; Evans, R. Mol. Phys. 1983, 48, 687.

Tjatopolous, G. J .; Feke, D. L.; Adin Mann, Jr., J. J. Phys. Chem. 1988, 192, 4006.

Valenzuela, D. P., Myers, A. L. Adsorption Data Handbook, Prentice-Hall, Englewood-
Cliffs: New Jersey, 1989.

Vanderlick, T. K.; Scriven, L. B.; Davis, H. T. J. Chem. Phys. 1989, 90, 2422.
van Megen, B.; Snook, I. K. Mol. Phys. 1982, 45, 629.

van Megen, B.; Snook, I. K. Mol. Phys. 1985, 54, 741.

Vera, J. H.; Prausnitz, J. M. Chem. Eng. J. 1972, 3, 1.

Wu, R. -S.; Lee, L. L.; Cochran, H. D. Ind. Eng. Chem. Res. 1990, 29, 977.

Yang, R. T. Gas Separation by Adsorption Processes, Butterworths, Stoneham, MA,
1987.

Young, D. M.; Crowell, A. D. Physical Adsorption of Gases, Butterworths, London,
1962.

"llllllllllllllll