.o i.
... _ .23." .
.5 ..

an

.73”., 71 . , u

2:. 2..

. . mar. .
a? a

.u...

. .. t
um .2

.3. .1
5:43.?!
:9.

 

3:. a
.n .
O. ‘.._Al_.~.

3153...!!!

 

it (ISL?
A1429?

.1. a.

,.i..hrx(u. g7.“

.. :05.
b

2.:
2.?

 

v\

) mltAbbhodgu.

o¢¢UArou¢a0 (‘1
”Ail.

 

 

:ltln n , [elf-5!; .
.1..- rr. . Y‘RI...’ I ‘1
itito. IX; 1 n‘filnin‘
5.530 \I‘)’ 50. Ill
‘ if! .

Altitnvk
.li b.I ._
.7.)

5
37... flank... W15:

IAOODI‘ . 1.

 

..r x .9}.
16.... o n .21. .Y ..

. yin...‘au~1|?..— : 1‘ V..o|
V ........J..., , . as? 23 :323 an 1.1..
v\ ‘. -. .. .. - i: . . . ‘ . . 35.9.3"? ‘ $52.1... .15... 5%-... .3!.;

, 2.25.53: 7. ‘ . .

4.6 3,.

| I.
II

 

llllllllllllll“HUI!!!H)llHlHllllllfllfllllHill

3017718119

LIBRARY

Michigan State
University

 

 

 

 

This is to certify that the

thesis entitled

Sfiiuéu es O'QC? as Mdgcxr 1mm of]
AC‘kH/OU‘PA Caf‘éah ama Ohfrrs F/avor

ficcovery From C146”)! P
presentedby

Aaf‘on D So UK l 6

has been accepted towards fulfillment
of the requirements for

Ma$1€r o {NSC Pncedegreein Chem C 11 Eng nei’r 03

(Made

V
Major professor

Date fill/‘7?

0-7639 MS U is an Afimaliw Action/Equal Opportunity Institution

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINE return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

we chIRCJDaanpGS—p.“

STUDIES OF GAS ADSORPTION ON ACTIVATED CARBON AND CHERRY
FLAVOR RECOVERY FROM CHERRY PITS

By

Aaron D. Soule

A THESIS
Submitted to
Michigan State University

in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Chemical Engineering

1998

ABSTRACT

STUDIES OF GAS ABSORPTION ON ACTIVATED CARBON AND CHERRY
FLAVOR RECOVERY F ROM CHERRY PIT S

By
Aaron D. Soule

This thesis considers two thrusts. The first thrust considers the application of
adsorption, as well as some other separations processes, to the recovery of natural cherry
flavoring fi'om cherry pits. Some benchmark work has been established relating to the
filtration, extraction, and distillation methods presented here. Each experimental
apparatus is presented in detail. Some adsorption data is also presented pertaining to
benzaldehyde and benzyl alcohol on synthetic resins. It was found that benzaldehyde
adsorbs more strongly in each case with the two aromatics showing signs of competitive
adsorption.

The second thrust of this thesis involves mathematical modeling of gas adsorption
on activated carbon. Both mixtures and pure systems will be discussed. The Simplified
Local Density (SLD) model treats the carbon surface area as slits with some effective slit
width. Fluid-solid interaction parameters are used to calculate density profiles across the
slit, thus, yielding an adsorbed amount upon integration over the slit width. The ESD
(Elliott-Suresh-Donohue) equation of state is used to calculate local compositions,
pressures, and fugacities. This method of adsorption prediction works well for gases
such as methane and ethylene. It also works well for ethylene-methane and ethane-
methane binary mixtures. However, for a COz-toluene fit, the model can only

demonstrate qualitative properties at this point in time.

I dedicate this finished product to my family for their unending support and love.

no.

ACKNOWLEDGMENTS

I would like to thank Dr. Carl Lira for his patience and support as I have struggled to
produce this finished product from my education. I would also like to thank Dr. Charles
Petty for giving me a second chance to complete the transport phenomena final exam
when my energy and motivation were failing me.

Regards to the friends who have come and gone over the years. You will see me again.

Special thanks to Cassandra Smith for her contributions in chapter 3.

Credit goes to Greg Donath and Ryoko Yamasaki for assisting in the adsorption runs and

issues related to the GC.

TABLE OF CONTENTS

List of Tables .................................................................................. vi
List of Figures ................................................................................. viii
Nomenclature .................................................................................. xi
Chapter 1: Introduction ...................................................................... 1
Chapter 2: Recovery of Benzaldehyde and Benzyl Alcohol from Cherry Pits ...... 3
Chapter 3: Adsorption Modeling With the BSD Equation of State .................... 40
Chapter 4: Carbon Dioxide: A Study of Supercritical Fluid Adsorption ............ 68
Chapter 5: Adsorption of Binary Mixtures ................................................ 75
Appendix A: Optimizing the Adjustable Parameters .................................... 96
Appendix B: Pure Adsorption Programs and Subroutines .............................. 99
Appendix C: Binary Mixture Adsorption Program ...................................... I47

LIST OF TABLES

2-1: Results of Amygdalin Breakdown in Hydrolyzate ................................. 5
2-2: Results of Mandelonitrile Breakdown in Hydrolyzate ............................. 6
2-3: Kernel and Shell Contributions to Hydrolyzate ..................................... 7
2-4: Effect of Pre-Warming Pits to 30°C before Hydrolysis ............................ 7
2-5: Comparison of Water-to-Pit Mass Ratios in Hydrolysis ........................... 8
2-6: GC Analysis of a Sample Regenerated from XAD-4 .............................. 25
2-7: GC Analysis of a Sample Regenerated from SP-850 ............................... 26
2-8A: Regeneration Samples from Run in Figure 2-12 .................................. 28
2-8B: Benzaldehyde Mass Balance for C02 Regeneration .............................. 29
2-8C: Benzyl Alcohol Mass Balance for C02 Regeneration ........................... 29
2-9: GC Profile of Original Regenerate ................................... ' ................ 30
2-10: Fraction A—The First 0.5 mL Boiled Over ........................................ 31
2-11: Fraction B—A Second Distilled Sample of 1 mL ................................. 31
2-12: Fraction C—The Remaining 1.5 mL of Bottoms .................................. 32

2-13: Mass Distribution of Benzaldehyde and Benzyl Alcohol During Boil-Over...32

2-14: Sample Masses as Percentages of Original Feed ................................... 32
2—15: 100 mL Regenerate Sample Before Boil-Over .................................... 33
2-16: Regenerate Sample After Complete Boil-Over .................................... 33
2472 GC Profile of a 100 mL Regenerate Sample after Boil-Over ..................... 35
2-18: Ethyl Acetate Phase of First Extraction .............................................. 35
2-19: Water Phase of First Extraction ....................................................... 35
2-20: Ethyl Acetate Layer After Second Extraction ....................................... 36

vi

2-21: Water Layer after Second Extraction ................................................. 36

2-22: Ethyl Acetate Solution to be Distilled ................................................ 37
2-23: GC of 1 mL Bottoms .................................................................... 37
2-24: GC of 0.2 mL Bottoms .................................................................. 37
3-1: Pure Component ESD Parameters ...................................................... 49
3-2: Comparison of Experimental Saturation Molar Volumes to the Bulk Peng-Robinson
and Bulk ESD Predictions ..................................................................... 50
4-1: Parameter Sensitivity to ACK Carbon Surface Area ................................. 69
5-1: Pure Component Fluid-Fluid Diameters and ESD Parameters ...................... 78
5-2: Mixing Parameters and Effective Mixture Diameters ................................ 79
5-3: Parameters used in Mixture Calculations ............................................... 82

vii

LIST OF FIGURES

2-1: Filtration Schematic ................................................................... 11
2-2: Adsorption Apparatus ................................................................ 15
2-3: Adsorption Test on 100 mL XAD-4 Resin at 15 mL/min Flow Rate .......... 16
2—4: Adsorption Test on 60 mL SP-850 at 15 mL/min Hydrolyzate Flow .......... 17
2-5: Adsorption Test on 60 mL XAD-4 with Flow at 15 mL/min .................... 18
2-6: Adsorption Test on 60 mL of SP-850 Reused Once at a Flow Rate of 15 mL/min
............................................................................................ 19
2-7: Adsorption Test on 55 mL of Reused XAD-4 with a Flow Rate of 15 mIJmin
............................................................................................ 20
2-8: Adsorption Test on 200 mL XAD-4 with a Flow Rate of 30 mI/min .......... 21
2-9: Adsorption Test on 200 mL SP-850 with a Flow Rate of 30 mL/min .......... 22

2-10: Benzaldehyde and Benzyl Alcohol Inlet Concentrations as Functions of Time
Corresponding to the Run in Figure 2-3 ................................................... 23

2-11: Benzaldehyde and Benzyl Alcohol Inlet Concentrations as Functions of Time
Corresponding to the Run in Figure 2-4 ................................................... 24

2-12: Adsorption Run on 100 mL XAD-4 Bed at 10 mIJmin Hydrolyzate Flow. . .27

3-1: Model of a Slit-Shaped Pore Showing the Variables used to Define Distances in the

Manuscript ...................................................................................... 47
3-2A: Adsorption of Ethylene on BPL Activated Carbon ............................... SI
3-2B: Adsorption of Ethylene on Columbia Grade L Carbon ........................... 52
3-3A: Adsorption of Ethane on BPL Carbon .............................................. 53
3-3B: Adsorption of Ethane on Columbia Grade L Carbon .............................. 54
3-4: Adsorption of Butane on Columbia Grade L Carbon ................................ 55
3-5: Adsorption of Propane on Columbia Grade L Carbon .............................. 56
3-6: Adsorption of Methane on BPL Carbon ............................................... 57

viii

3-7A: Adsorption of Propylene on Columbia Grade L Carbon ........................... 58

3-7B: Adsorption of Propylene on BPL Carbon ............................................ 59
3-7C: Adsorption of Propylene on Black Pearls I ........................................... 60
3-8: Adsorption of Nitrogen on Columbia Grade L Carbon ................................ 61
3-9: Adsorption of Acetylene on Columbia Grade L Carbon .............................. 62
3-10: Adsorption of Carbon Monoxide on Columbia Grade L Carbon ................... 63
3-11: Correlation of Fluid-Solid and Fluid-Fluid Interaction Potentials on Columbia
Grade L Carbon ................................................................................... 64
3-12: Correlation of Fluid—Solid and Fluid-Fluid Interaction Potentials on BPL Carbtg;
4-1: Carbon Dioxide Adsorption on Degussa IV Activated Carbon ...................... 70
4-2: Carbon Dioxide Adsorption on Columbia Grade L Carbon ........................... 71
4-3: Carbon Dioxide Adsorption on BPL Carbon .......................................... I .72
4—4: Carbon Dioxide Adsorption on ACK Carbon ........................................... 73

5-1: Methane-Ethane Adsorption on BPL Carbon at 301.4 K with a Bulk Ethane Mole
Fraction of 0. 733 .................................................................................. 84

5-2: Methane-Ethane Adsorption on BPL Carbon at 301.4 K and a Bulk Ethane Mole
Fraction of 0.501 .................................................................................. 85

5—3: Methane-Ethane Adsorption on BPL Carbon at 301.4 K and a Bulk Ethane Mole
Fraction of 0.255 .................................................................................. 86

5-4: Methane-Ethane Adsorption on BPL Carbon at 212.7 K and a Bulk Ethane Mole
Fraction of 0.733 .................................................................................. 87

5-5: Methane—Ethane Adsorption on BPL Carbon at 212.7 K and a Bulk Ethane Mole
Fraction of 0.501 .................................................................................. 88

5-6: Methane-Ethane Adsorption on BPL Carbon at 212.7 K and a Bulk Ethane Mole
Fraction of 0.255 .................................................................................. 89

5-7: Methane-Ethane Adsorption on BPL Carbon at 260.2 K and a Bulk Ethane Mole
Fraction of 0.255 ............................................................................... 90

5-8: Methane-Ethylene Adsorption on BPL Carbon at 212.7 K with a Bulk Ethylene
Mole Fraction of 0.74 ......................................................................... 9]

5-9: Adsorbed Phase Mole Fractions of Ethylene-Methane on BPL Carbon at 212.7 K
With a Bulk Ethylene Mole Fraction of 0.74 .............................................. 92

5-10: Pressure (Density) Effects of COz-Toluene Adsorption on Degussa WSIV. . .93

NOMENCLATURE

A - surface area per unit weight of adsorbent

a - Pen g- Robinson attractive energy parameter

b - ESD size parameter

c - shape factor for BSD repulsive term

Eta - distance from first wall in slit model/fluid-solid diameter
f - fugacity

H - distance from carbon center to carbon center across slit
k - Boltzmann’s constant

P - pressure

R - gas constant

T - temperature

Xi - distance from second wall in slit/fluid-solid diameter
Y - ESD attractive energy parameter

2 - direction of finite length

2 - compressibility factor

Greek

or - spacing between carbon planes
8 - interaction parameter

1] - reduced density

p - density

¢ - fugacity coefficient

6 - molecular diameter

1'“ - excess

‘1’ - fluid-solid potential

Superscripts

attr - attractive term
LCL - local property
rep - repulsive term

Subscripts

A - component A in binary mixture

B — component B in binary mixture

atoms - refers to atom density

bulk - property of bulk fluid

calc - fugacity defined by the equation of state
ff - fluid-fluid property

fs - fluid—solid property

xi

Chapter 1:
INTRODUCTION AND BACKGROUND

Adsorption is defined as the attraction of fluid particles onto the surface of a solid
with the formation of an “adsorbed phase”. It is a useful separation process which often
requires no addition of heat. This thesis looks at some experiments involving this
problem, as well as the utilization of other separations techniques, in the area of food
processing. The significance of this work lies in the fact that natural cherry flavoring
(benzaldehyde) is much more valuable than artificial flavoring. Even though the natural
and artificial flavors are identical in taste, “natural flavoring” labeling is important to the
food industry. Furthermore, cherry pits are an immense waste product which cherry
companies must deal with. If a process could be developed to recover some of the
remaining flavor from the cherry pits, it might prove economically beneficial to such
companies. There is also an environmental issue with cherry pits in cyanide release. The
thesis work of Cassandra Smith has suggested that air stripping may be the solution to
making any liquid waste products from this process meet the environmental regulations.

Adsorption modeling is also dealt with in this thesis. Unlike in the case of cherry
flavoring which involves liquid systems, the models presented here will be for gas
systems interacting with activated carbon. Some of the primary gases of interest are
methane, ethane, propane, ethylene, propylene, carbon dioxide, toluene, and nitrogen.
The ESD equation of state will be combined with the Simplified Local Density model
first developed in the dissertation of Bharath Rangarajan. The SLD model was first used
with the van der Waals (Rangarajan, 1992) and Peng-Robinson (Subramanian, 1995)

equations. While these equations did give some promising results, the ESD equation of

state does better in many cases of pure component adsorption (e. g. methane, ethylene).
Furthermore, it does a much better job in calculating mixture adsorption than does the
Peng-Robinson equation as will be seen with various fits of ethane-methane and
ethylene-methane. However, a fit to toluene-carbon dioxide adsorption data shows that
the ESD still only gives qualitative predictions in some cases. Currently, the mixture
model assumes uniform particle size, which is an unrealistic assumption for some cases.

Despite some its weaknesses, the ESD equation of state still serves as a good

predictor of bulk fluid conditions. It also gives excellent engineering estimates for
adsorption in many systems. Eventually, it is hoped that this equation will be adopted for
use in process simulators and other software, both for bulk and adsorption calculations.
Further work still needs to be done with the SLD model. First, it must be modified to
deal with the different particle sizes present in mixtures. Secondly, it must be modified
to handle mixtures with more than two components. Once code is developed for these

applications, then it might more easily be adopted for commercial use.

REFERENCES

Rangarajan, Bharath; “Shrinkage Characterization in the Production of Aerogels and a
Model for Physical Adsorption of Fluids on Solids—From Vacuum To
Supercritical Conditions”; Doctoral Dissertation, Michigan State University,
1992.

Smith, Cassandra A.; “Adsorption of Water on Carbon: A Study Using the BSD
Equation of State”; Master’s Thesis, Michigan State University, 1997.

Subramanian, Ramkumar; “Adsorption: A Study Adapting Cubic Equations of State”;
Doctoral Dissertation, Michigan State University, 1995.

Chapter 2:
RECOVERY OF BEN ZALDEHYDE AND
BENZYL ALCOHOL FROM CHERRY PITS

INTRODUCTION

A Phase I USDA grant was provided for developing a process to recover natural
flavoring (benzaldehyde and benzyl alcohol) from cherry pits. This work was performed
in conjunction with Natura, Inc (Lansing, MI). It is motivated by the fact that market
prices for natural flavors are much higher than for artificial flavors. While being
identical in taste quality, artificial and natural flavor components can be distinguished by
their deuterium distributions (Hagedorn, 1992).

This work required several experimental steps of process development. It was
necessary to develop a hydrolysis procedure which would provide the greatest
concentration and quantity of benzaldehyde while minimizing its oxidation to benzoic
acid. Filtration was also required in order to eliminate plugging in the bed and tubing.
Adsorption, the heart of the process, involved a large amount of experimental work
aimed toward optimizing several variables. Success was measured by the characteristic
breakthrough curves. Finally, several efforts have given limited success in the

regeneration and purification of the aromatics.

HYDROLYSIS

Based on past experimental work on this project, the same hydrolysis conditions
were agreed upon throughout all of the adsorption runs performed. The pits were ground

up and then added to water preheated to 50°C. For each solution, a 5:] water-to-pit mass

ratio was used. A discussion of this choice will be given later on. The water and pits
were mixed together for 1 hour and then. put through the filtration process described in
the next section.

Due to differences in the pit yields from previous work, various tests were
performed analyzing optimal hydrolysis conditions as well as the chemical processes
occurring during hydrolysis. Tests were done with water-to-pit ratios, amygdalin
conversion, mandelonitrile conversion, kernel and shell contributions, and pre-heating of
the pits.

Based on information from two different papers (Li et al., 1992; Zheng and
Poulton, 1995), the pits are expected to contain quantities of amygdalin and
mandelonitrile. These compounds are known to break down in water producing both
cyanide and benzaldehyde. A 0.5 gram quantity of amygdalin (molecular mass: 457.42) .
was added to 400 mL of hydrolyzate at 50°C and found to mostly break down into
benzaldehyde (molecular mass: 106.12) after about 15 minutes of mixing. Table 2-1
gives the gas chromatography data of this experiment. It shows the hydrolyzate before
and after the amygdalin addition. The lower and higher retention times correspond to
benzaldehyde and benzyl alcohol (molecular mass: 108.1) respectively. The data shows
a net gain of 9.64E-4 moles of benzaldehyde. Since 10.93E-4 moles of amygdalin were
added, this experiment shows that 0.88 moles of benzaldehyde were created per mole of
amygdalin added. One mole of amygdalin stoichiometrically breaks down into one mole
of benzaldehyde when completely hydrolyzed. This experiment is not far from
confirming this result. Note the 10% reported drop in benzyl alcohol concentration

during this addition. This is likely to be due to the normal variations seen on this piece of

equipment. Apparently, amygdalin plays little or no direct role in the presence of benzyl

alcohol.
Table 2-1: Results of Amygdalin Breakdown in Hydrolyzate

H l
2.76 0.4592 25.4 9.57
4.22 1 alcohol 0.9324 51.5 19.06

+Am dam
2.72 deh 5.0862 281.0 105.92
4.20 1 alcohol 0.8345 46.1 17.06

A similar test with similar results was also done for mandelonitrile (molecular
mass: 133.15). Mandelonitrile was more difficult to quantitatively measure due to its
tendency to stick to the sides of glassware; thus, a larger margin of error should be
allowed for this experiment. Approximately 0.5 mL of mandelonitrile was added to 300
mL of hydrolyzate. The results are shown in Table 2-2. The specific gravity- of
mandelonitrile is 0.95 grams per milliliter. Thus, a total of 356.7E-5 moles were added to
the hydrolyzate. The experiment shows a corresponding 297E-5 mole increase in
benzaldehyde. This means that 0.83 moles of benzaldehyde are formed per mole of
mandelonitrile added to the hydrolyzate. Once again, the stoichiometric ratio should be
1:1 when completely hydrolyzed. Also, the observed benzyl alcohol concentration
actually decreases by about 7% in this experiment—probably due to neutral GC
variations. As in the case of amygdalin, it is clear that mandelonitrile plays little or no

direct role in the presence of benzyl alcohol.

  

Table 2-2: Results of Mandelonitrile Breakdown in Hydrolyzate

    

H drol
2.77 0.8310 45.9 12.98
4.23 1 alcohol 1.0317 57.0 15.82

H + Mandelonitrile
2.72 19.8320 1096 310
4.19 1 alcohol 0.9733 53.8 14.93

The pits were separated into their constituent kernels and shells in another test.
The kernel accounts for about 80% of the total pit mass, and the shell makes up about
20%. Table 2-3 shows the results of this test on pits taken from the same bucket. First, a
standard hydrolysis (5:1 water-to-pit mass ratio) was performed. A hydrolysis was then
performed on a sample of kernels only at a 5:1 water-to-kemel ratio. Finally, a
hydrolysis was performed on a sample of shells at a 5:1 water-to-shell ratio. The results
show that a sample of kernels gives a 47% higher yield of benzaldehyde and a 117 %
higher yield of benzyl alcohol than an equivalent mass of pits. On the other hand, a
sample of shells gives 14% and 40% of the respective yields of benzaldehyde and benzyl
alcohol in an equivalent mass of pits. Thus, the kernels make the larger contribution to
the hydrolyzate concentrations of these chemicals. Note that the pits yield a 5:6 ratio of
benzaldehyde to benzyl alcohol while the kernels yield a 1:2 ratio and the shells yield a
1:4 ratio. One would expect the hydrolysis with pits to yield concentration ratios
somewhere between the kernel and shell hydrolyses. It seems that the presence of both
shells and kernels might be facilitating some chemical interconversion between
benzaldehyde and benzyl alcohol. More experiments should be done to confirm this

hypothesis. The reasons for such an interconversion are unknown at this time.

Table 2-3: Kernel and Shell Contributions to Hydrolyzate

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area P_a_rts Per Million
Standard Hydrolysis
2.79 (benzaldehyde) 0.9030 49.9
4.26 (benzyl alcohol) 1.0900 60.2
Kernel Hydrolysis
2.75 (benzaldehyde) 1.3247 73.2
4.22 (benzyl alcohol) 2.3679 130.8
Shell Hydrolysis
2.79 (benzaldehyde) O. 1248 6.9
4.24 (benzyl alcohol) 0.4358 24.1

 

 

 

 

One other test involved pre-heating the pits to 30°C overnight in an attempt to
curtail possible side effects from mold and residual fruit moisture on the surface of the
shells. These results are shown in Table 2—4. While the benzaldehyde concentration
increased slightly (within the normal range of data variation), the benzyl alcohol

concentration went down by about 30%.

Table 2-4: Effect of Pre-Warrnin g Pits to 30°C before Hydrolysis

No Warmin .
2.80 0.6608 36.5

4.27 1.3391 74.0

Warrnin
2.80 0.7404 40.9
4.26 0.9189 50.8

 

Another test involved adjusting water-to-pit mass ratios. These results are shown
in Table 2-5. Ratios of 2:1, 3:1, 4:1, and 5:1 were compared for benzaldehyde and benzyl
alcohol yield. 100 grams of pits were used in each experiment with the water quantity

being adjusted to each of the ratios. The 5:1 ratio has the largest benzaldehyde yield,

while the 3:1 ratio has the largest benzyl alcohol yield. The 2:] ratio shows a large
drop-off in the concentration of each component as compared to the 3:1 ratio. It is
suspected that the large quantity of solids in the 2:1 ratio begin to cause mass transfer
limitations since most of the water absorbs into the ground pits. Due to the higher

benzaldehyde yield, the 5:1 ratio has been chosen for the majority of the work.

Table 2-5: Comparison of Water-to-Pit Mass Ratios in Hydrolysis

 

 

 

 

 

 

 

 

 

 

 

 

 

, Retention Time (min.) Peak Area Parts Per Million Yield (g)

_5_:l

2.77 0.8310 45.9 0.0230
4.23 1.0317 57.0 0.0285
4:_1

2.75 0.7160 39.6 0.0158
4.20 l .2307 68.0 0.0272
3:;

2.67 0.6745 37.3 0.01 12
4.20 2.0465 1 13.1 0.0339
2.74 0.2436 13.5 0.0027
4.20 1.2416 68.6 0.0137

 

 

 

 

 

 

Based upon the above results, it is natural to consult the literature for further
guidance regarding the origin and interconversion of benzaldehyde and benzyl alcohol.
Reilly and coworkers (1986) discuss amygdalin conversion in peaches. Amygdalin first
breaks down into prunasin which in turn breaks down into mandelonitrile. In a paper by
Swain and Poulton (1994), black cherry seeds contain larger quantities of amygdalin in
earlier stages of deveIOpment than in later stages of development. Conversely, higher
amounts of prunasin are found in the later stages of seed development. Also, the paper
identifies enzymes associated with the stages of amygdalin breakdown. Amygdalin
hydrolase catalyzes the breakdown of amygdalin to prunasin while GT-H

(UDPGzprunasin glucosyltransferase) catalyzes the reverse reaction. Prunasin hydrolase

catalyzes the breakdown of prunasin to mandelonitrile while GT-I
(UDPszandelonitrile glucosyltransferase) catalyzes the reverse reaction.
Mandelonitrile lyase catalyzes the conversion of mandelonitrile into benzaldehyde and
cyanide. Kawabe and Morita (1994) describe the role that a particular breed of fungus
plays in interconverting benzaldehyde and benzyl alcohol. The paper also describes a
fungus-facilitated mechanism for benzaldehyde and benzyl alcohol formation by
identifying such precursors as L-phenylalanine, t-cinnarnic acid, and 3-phenylpyruvic
acid. A paper by Lamer and coworkers (1996) identifies another fungus responsible for
benzaldehyde and benzyl alcohol interconversion and formulates another mechanism for
their production. It lists styrene, l-phenyl ethanone, and phenylacetaldehyde as other

precursors.

FILTRATION

The hydrolysis procedure yields a solution containing a wide variety of particle
sizes. This made efficient filtration a challenging task—especially for large quantities of
hydrolyzate. A couple of stages were necessary in removing particulates.

The first stage involved removing the cherry pit shells and kernels through a
screen with a sieve size of 180 microns. This process was very quick, and it effectively
removed all of the large particles. The hydrolyzate was poured into a pan with the screen
built onto the bottom. The pits were emptied out of the pan as necessary.

The second stage was more complicated. It was performed in a metal cylinder
with a diameter of about 3 inches and a length of about 15 inches. The following items
composed the filter (bottom to top): 4 grams of glass wool, 32 grams of diatomaceous
earth (normally used in swimming pool filters), 130 grams of sand, and 50 grams of
cherry pits for the purpose of spreading liquid flow over the entire diameter of the filter.
Figure 2-1 shows a schematic of this filter. A pump pushed the hydrolyzate through the
filter from a feed tank. While this filter effectively removed all noticeable particulates
from the hydrolyzate, the filter cake caused it to plug very quickly. The solution to this
problem involved mixing diatomaceous earth into the hydrolyzate feed. It was found that
maintaining 10 grams of diatomaceous earth per liter of hydrolyzate worked quite well
throughout the duration of the filtration. While the diatomaceous earth increased the
volumetric buildup of the filter cake, it made the filter cake sufficiently porous so that

plugging was not as great of an issue.

10

 

__ Hydrolyzate
Feed Mixed With
Diatomaceous
Earth

 

 

 

 

Pressure TYalw
Vent

___LL ( Pump

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

“° -------...-. Extra Space
--....... -- ............... Pits
........... Sand
"""""""""""" Diatomaceous Earth
___ _ ..................... Glass W001
1
I Dispenser

 

Figure 2-1: Filtration Schematic

ADSORPTION/BREAKTHROUGH

Figure 2-2 shows the apparatus used for conducting adsorption runs. The feed
resided in a 20—L glass jar which in turn was placed in an ice bath. The chilled condition
served to slow down bacterial growth in the hydrolyzate. PharMed tubing carried fluid
from the bottom of the jar to the Masterflex pump. The tubing then ran from the pump to
the top of the adsorption bed. The adsorption bed itself consists of a glass tube
measuring 3 centimeters in diameter. Silicon stoppers were placed on the top and bottom
openings in order to hold the contents. Glass wool was placed underneath the adsorbent
in order to prevent spillage into the tubing. A clamped tube was placed at the top of the
column in order to relieve the system of air when necessary. The beds ran completely
filled with hydrolyzate. A three-way valve was placed at the bottom of the column for
the purpose of convenient sarnplin g. '

These sets of experiments were aimed toward determining the optimal operating
conditions for the adsorption of benzaldehyde and benzyl alcohol. The primary variables
were the resin (adsorbent) bed size and the flow rate. Screening procedures in previous
work have identified XAD-4 (from Rohm and Haas) and SP-850 (from Mitsubishi
Chemical), as the best known resins for benzaldehyde adsorption. All of the adsorption
runs performed with new resin at a flow rate less than or equal to 15 bed volumes per
hour (Figures 2—3 through 2-5) show the breakthrough of benzaldehyde to be slow and
linear. Benzaldehyde reaches 20% breakthrough in a 60 mL bed of XAD-4 after the
elution of about 100 bed volumes of hydrolyzate at 15 bed volumes per hour (Figure 2-5).
A similar bed of SP-850 at the same flow rate reaches 20% breakthrough at an elution of

200 bed volumes (Figure 2—4). Thus, SP-850 seems to have a higher capacity for

benzaldehyde. Note that a bed of 100 mL XAD-4 with a flow rate of 9 bed volumes per
hour reaches 20% breakthrough at an elution of 300 bed volumes of hydrolyzate (Figure
2-3). This implies that the slower flow rate relative to the bed volume gives better
adsorption performance. For these same runs, benzyl alcohol seems to reach its
maximum at or above 100% breakthrough despite the more significant scatter.

Figures 2-6 and 2-7 show runs where the resins were reused. They had already
been put through one cycle of adsorption and desorption before undergoing the runs
shown. The regeneration was done with liquid carbon dioxide pressurized to about 1300
psig in the bed. In each case, about a SOC-liter equivalent of carbon dioxide at room
temperature and pressure was passed through the bed in order to remove any adsorbates.
This work will be discussed further in the next section. The reused resin actually seems
to show better performance. Figure 2-6 shows a 20% breakthrough for SP-850 at 600
bed volumes of effluent. Figure 2-7 shows a bed of XAD-4 similar to that in Figure 2-5.
The bed of reused resin once again performs better—20% breakthrough at over 300 bed
volumes of effluent as opposed to the 100 bed volumes of effluent shown in Figure 2-5.
Furthermore, benzyl alcohol seems to take on a more linear breakthrough curve with the
reused resin. In each case, the benzyl alcohol curve extends well beyond 200%
breakthrough. This seems to imply that either the resin is not being completely
regenerated before the second use or else benzyl alcohol is being formed by a chemical
reaction. More investigation is needed in this matter.

Figures 2-8 and 2-9 examine the effects of doubling the superficial velocity. In
both cases, comparatively larger bed volumes and flow rates were tested. As in Figure 2-

3, the flow rates are only 9 bed volumes per hour (30 mIJmin). However, the superficial

l3

velocities are doubled from 2.12 cm/min to 4.24 cm/min. The XAD-4 bed in Figure 2-8
shows a performance similar to its counterpart in Figure 2-3—20% breakthrough at an
elution of nearly 300 bed volumes of hydrolyzate. The SP-850 bed in Figure 2-9 shows a
similar performance despite some irregularities in benzaldehyde data early in the run. It
exhibits 20% breakthrough at about 300 bed volumes of hydrolyzate elution. Thus, in
this range of flows, superficial velocity does not appear to have a large effect on
benzaldehyde adsorption. It is interesting to note that benzyl alcohol seems to exhibit a
maxima in Figure 2-9.

It is important to note that feed concentration is not kept perfectly constant
throughout each of these runs. Figures 2-10 and 2-11 show typical fluctuations for an
adsorption run. During each of these runs, the feed tank was being refilled every few
hours. This accounts for the various spikes in the graphs. Note that the benzaldehyde
and benzyl alcohol concentrations do not change greatly relative to one another. The
benzaldehyde concentration hovers near 45 ppm in both Figures 2-10 and 2-11, while the
benzyl alcohol concentration hovers near 15 ppm in each case. This may appear to
conflict with some of the data presented in the Hydrolysis section, but the hydrolyzate

concentrations do vary with different batches of pits.

l4

 

Rubber Stopper &
J Pump

Chilled, closed Air Relief Tube
hydrolyzate feed

jar "l

 

 

 

 

 

 

 

 

 

 

 

 

/
Silicon
Stopper
Glass 1
Ice Bath Adsorption
Bed Excess
Liquid-Filled
Space
M
Adsorbent
WV
Glass Wool
Silicon Stopper
Waste Stream (to drain) Sample Stream
3-way valve

Figure 2-2: Adsorption Apparatus

15

C(out)/C(in)

 

1.8

 

 

 

 

1.6 4~ I
1.4 ~-
I
1.2 Tl— . I
I I
1 -~ I
I I I
0.8 w I
I
._ I I
0'6 y = 0.0009x - 0.0294
2 -
0,4 -_ R — 0.7971 .
02 + M
O
0 .. L—
-O.2 1% t 4 l l i i
0 50 100 150 200 250 300 350 400
.Bed Volumes Of Hydrolyzate

 

O C(out)/C(in) Benzaldehyde——
I C(out)/C(in) Benzyl Alcohol

 

Figure 2-3: Adsorption Test on 100 mL XAD-4 resin, 15 mL/min
flowrate (9 bed volumes per hour, superficial velocity: 2.12
cm/min).

16

C(out)/C(in)

 

1.6

 

 

 

 

I
1.4 ~-
I
I
I I
1.2 ‘” . I I
1 .1.—
0.8 --
I
0.6 -- y=0.0011x -00174
R2 = 0.9442 9 ’
0 100 200 300 400 , 500

Bed Volumes

9 C(out)/C(in) Benzaldehyde
I C(out)/C(in) Benzyl Alcohol
—Linear (C(out)/C(in) Benzaldehyde)____
Figure 2-4: Adsorption test on 60 mL SP-850 at 15 mL/min
hydrolyzate flow (15 bed volumes per hour, superficial velocity:
2.12 cm/min).

 

 

C(out)/C(in)

 

1.6

1.4 ~~ I

1.2 -- I

0.8 --

y = 0.001 1x + 0.078

0.6 4- R2 = 0.9544 ’

    

 

0 100 200 300 400 500 600
' Bed Volumes

Ar-
-0-

I
T

 

o C(out)/C(in) Benzaldehyde— ’
I C(out)/C(in) Benzyl Alcohol
—Linear (C(out)/C(in) Benzaldehyde) --_.

 

Figure 2-5: Adsorption test on 60 mL XAD-4 with flow at 15 mL/min
(15 bed volumes per hour, superficial velocity: 2.12 cm/min).

 

 

3.5

 

 

 

3 5.
y = 0.0055x — 0.2062
R2 = 0.9516 .
2.5 -- I '
I
g 2 -- a
Q
'5
O
:5 1.5 --
1 __
I1
I
0 5 y = 0.0002x + 0.0155
' a2 = 0.515
O O .
0 J ———- .L r .1 .L i :
0 100 200 300 400 500 600 700

Bed Volumes
_ I“. C(guFI—Cfin) Benzaldehyde
I C(out)/C(in) Benzyl Alcohol
—Linear (C(out)/C(in) Benzaldehyde)
___:-_-_-Linear (C(out)/C(in) Benzyl Alcohol) __ .

 

 

Figure 2-6: Adsorption test on 60 mL of SP-850 reused once with a flow
rate of 15 mL/min (15 bed volumes per hour, superficial velocity: 2.12
cmlmin).

C(out)/C(in)

 

    
    

 

 

 

 

2.5
I
2 ,_
y = 0.0049x - 0.0231
1.5 -- R2 = 0.9539

1 4

y = 0.00071: + 0.0045

0.5 -- R2 =g.6026
o 09L.
0
0 0
° 9
0 i ‘r t : :
0 100 200 300 400 500

Bed Volumes

 

O C(out)/C(in)Benzaldehyde
‘ I C(out)/C(in) Benzyl Alcohol
—Linear (C(out)/C(in) Benzaldehyde)
—Linear (C(out)/C(in) Benzyl Alcohol)

 

Figure 2-7: Adsorption testing on 55 mL of reused XAD-4 with a
flow rate of 15 mL/min (16.4 bed volumes per hour, superficial
velocity: 2.12 cm/min).

20

1.8

 

 

 

 

 

1.6 + I
I I
1.4 -~
1.2 4’ I I
E
V 1 _0_
g
§ 08 -- u ' '
v ' I I
0
0.6 +
0,4 .. y = 0.0005x + 0.053
R2 = 0.6507 .
0.2 "' . W
W .
O 1 t .4. % t :
O 50 100 150 200 250 300 350

Bed Volumes

 

0 CI out /Cin ngBenzaldehyde
I CI out /Cin IBenzy lAlcohoI
—-Linear ( (out)/C(in) Benzaldehyde)

 

Figure 2-8: Adsorption test on 200 mL XAD-4 with a flow rate of
30 mL/min (9 bed volumes per hour, superficial velocity: 4.24
cmlmin).

21

C(outycan)

 

 

 

 

2.5
21-
I
I
1.5 -- '
I I I
I
I
1 ..
I
0.5 -- '
I
e
I
e e
. O O Q . . .
0 + : + .L .L :
0' 50 100 150 200 250 300 350

Bed Volumes
1 e C(out)/C(in)—Benzaldehyde_
I C(out)/C(in) Benzyl Alcohol

 

 

 

Figure 2-9: Adsorption test on 200 mL SP—850 with a flow
rate of 30 mL/min (9 bed volumes per hour, superficial
velocity: 4.24 cm/min).

22

 

70

 

 

«A,

 

 

 

 

 

 

 

E

o.

3. 40

e

o

E

5

0 30

e

a V

0
20 I
10 /
0 . . . . .

O 500 1000 1500 2000 2500 3000
Minutes Elapsed

 

+ Benzaldehyde ppm
—I— Benzyl Alcohol ppm

 

 

 

Figure 2-1 0: Benzaldehyde and Benzyl Alcohol inlet Concentrations
as Functions of Time Corresponding to the Run in Figure 2-3

23

Concentration (ppm)

70

60

10

 

 

 

 

 

 

 

 

 

 

0 200 400 600 800 1000 1 200 1 400 1600 1 800

Time (min)

+ Benzaldehyde ppm L
-: Penzv' Alcohglppm,

 

 

Figure 2-11: Benzaldehyde and Benzyl Alcohol inlet Concentrations
as Functions of Time Corresponding to the Run in Figure 2-4

24

2000

REGENERATION

After adsorption, the resin is regenerated with carbon dioxide pressurized to 1300
psig in order to obtain the desired components. While carbon dioxide doesn’t have the
best solubility properties, it is a relatively inexpensive gas. Furthermore, its liquid and
gaseous properties can both be utilized at room temperature without any heating—only
pressurization and depressurization. In its pressurized state, the liquid carbon dioxide
passes through the bed and removes adsorbates from the pores of the resin. It is then
depressurized into its gaseous state where it relinquishes its solutes. This is where the
product is collected in this stage. Care must be taken to prevent the lines from freezing
during depressurization. The resulting product is an aqueous solution many times more
concentrated with the desired components than the hydrolyzate. The product contains
water because of the moisture remaining in the bed after adsorption. Also, this product is
a brown color. Gas chromatography reveals the presence of various unknowns as shown
in Tables 2-6 and 2—7. Benzaldehyde is at 2.94 minutes on Table 2-6 and 2.92 minutes on

Table 2-7. Benzyl alcohol is at 4.41 minutes on Table 2-6 and 4.40 minutes on

 

 

 

 

 

 

 

Table 2—7.

Table 2-6: GC Analysis of a Samle Regenerated from XAD-4
Retention Time (min) Peak Area Parts Per Milliog
2.94 (benzaldehyde) 124.91 690]

4.41 (benzyl alcohol) 51.63 2852

5.38 3.33

7.10 8.16

13.79 2.29

 

 

 

 

Table 2-7: GC Analysis of a Sample Regenerated from SP-850

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area Parts Per Million
0.18 2.31

0.20 3.92

2.92 (benzaldehyde) 1 12.38 6209

4.40 (benzyl alcohol) 40.09 2215

5.37 3.95

7.06 4.70

8.42 1.49

11.43 1.06

13.69 1.59

 

 

 

 

The focus will now be shifted toward analyzing the efficiency of the regeneration
step. Figure 2-12 shows the benzaldehyde breakthrough curve for the sample under
consideration. This curve is needed for calculating the total amount of benzaldehyde
adsorbed onto the resin. A linear regression will represent the portion of the data
breaking through. The first step is to find the area under the breakthrough curve over the
entire flow range of the run. The second step is to calculate the average C(out)/C(in)
ratio over the entire range of the run, which turns out to be 0.023. This means that, on
average, 97.7% of the inlet benzaldehyde concentration is adsorbed over the course of the
entire run. For this experiment, the inlet benzaldehyde concentration averaged about 44
ppm over the 120 bed volumes fed (12 liters). Thus, 43 ppm of that amount was
adsorbed throughout the run on average. This gives a total of 0.52 grams of

benzaldehyde adsorbed.

26

 

0.09

0.08 -

0.07 ~

   

p
3

y = 0.001x - 0.0462
R2 = 0.6896

9
E"

C(out)IC(in) Benzaldehyde
E

p
9

Average=0.023

 

0.02 i

0.01 ~

 

 

 

 

0 20 40 60 80 100 120
Elapsed Bed Volumes

 

e C(out)IC(in)
Benzaldehyde
—Linear (C(out)IC(in)
Benzaldehyde)

 

 

 

Figure 2-12: Adsorption Run on 100 mL XAD-4 Bed at 10
miJmin hydrolyzate flow.

27

140

This sample was regenerated with the equivalent of 224 liters of carbon dioxide
at room temperature and pressure at a flow rate between 2 and 3 liters per minute
measured after depressurization. Three samples were taken as shown in Table 2-8A.
Sample A was taken from the first 60 liters of carbon dioxide passed through the system.
This is when most of the water is flushed out of the resin bed; hence, the comparatively
large volume. Sample B was taken from the remaining 160 liters passed. This sample
has a smaller volume but is much more concentrated. Sample C was taken from washing
the sample chamber and tubing with ethanol. Judging by the significant concentration of
the sample, this portion of the experiment cannot be ignored. The total benzaldehyde
mass recovered from the three samples was 0.1712 grams. This means that 33% of the
adsorbed benzaldehyde was recovered by regeneration.

Table 2-8A:

es From the Run in F' 2-12

       

   

A 28.5 658 0.0188

B 6.8 9111 0.0620
C 11 8221 0.0904

It is not certain whether the amount of carbon dioxide used maximizes the yield.
Further work contributed by Greg Donath and Ryoko Yamasaki gives more information
regarding possible inefficiencies in this stage. Table 2-8B gives a mass balance on
benzaldehyde adsorption and regeneration from 55 mL of resin. Adsorption is
represented with a positive number while desorption is represented with a negative
number. The amount adsorbed was calculated in a manner similar to the previous system
shown. After adsorption, the resin was rinsed and stored in water until the regeneration
could be done. The benzaldehyde content of this water was analyzed. After the C02

regeneration was completed, the regeneration equipment was rinsed with propanol in

28

order to recover any residual desorbed material. The resin was then soaked in 100 mL
propanol in an attempt to recover anything remaining on the resin. The data shows that

only 40% of the benzaldehyde adsorbed on the resin can be accounted for by current

 

 

 

 

 

 

 

methods.
Table 2-8B: Benzaldehyde Mass Balance for CO2 Regeneration
fl Masfigg) % of Adsorbed Mass
Adsorption 2.69E-4 100
Water Wash/Storage -3.20E-6 -1.2
C02 Regeneration -7.97E-5 -29.6
System Propanol Wash -2.08E-5 -7.7
Resin Wash with Propanol -2.97E-6 -l.1
Mass Unaccounted For 1.62E-4 60.3

 

 

 

 

 

A similar balance was done for benzyl alcohol. The same experimental steps
were analyzed and are shown in Table 2-8C. What is interesting about this experiment is
that 20% more benzyl alcohol is reported as recovered than the amount originally

adsorbed. This apparent phenomenon should be examined more deeply in future work.

Table 2-8C: Benzyl Alcohol Mass Balance for CO2 Regeneration

 

 

 

 

 

 

 

 

 

 

Step Mass (g) % of Adsorbed Mass
Adsorption 5.66E-5 100

Water Wash/Storage -2.29E-5 405

C02 Regeneration -3.19E-5 -56.4

System Propanol Wash -6.09E-6 -10.8

Resin Wash with Propanol -6.44E-6 -11.4

Mass Unaccounted For - l .07E-5 -l9.l

 

29

 

PRODUCT PURIFICATION EXPERIMENTS

Different efforts have been made to further concentrate and isolate benzaldehyde
and benzyl alcohol from the regenerated sample. The GC data shown in the previous
section implies the presence of some heavier substances (Tables 2-6 and 2-7). The most
logical first course of action would be to do an initial distillation, separating the lighter
components from the heavier components. This experiment was tried with the 3 mL
sample of regenerate shown in Table 2-9. This sample is not related to the regeneration
data previously shown. Each component’s peak area is quantitated as a percentage of the
total area of all peaks recorded on the GC.

Table 2-9: GC Profile of Original Regenerate

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area % Area of All Peaks Parts Per Million
0.15 2.07 1.32

0.17 4.80 3.06

0.47 6.28 4.00

0.54 1.93 1.23

2.94 (benzaldehyde) 77.38 49.30 4275

4.42 (benzyl alcohol) 46.33 29.52 2560

5.40 2.35 1.50

7.12 3.30 2.10

 

 

 

 

 

The above sample was separated into three different fractions via a distillation
column with a boiler and water-cooled condenser. Before distilling, argon was put
through the system in order to remove the air. Oxygen has a tendency to oxidize
benzaldehyde to benzoic acid. Table 2-10 shows the composition of Fraction A. This
consists of the first 0.5 mL boiled and condensed. Note that this fraction is 61% more
concentrated with benzaldehyde than the original. It is also 81% more concentrated with

benzyl alcohol than the original sample.

30

Table 2—10: Fraction A—The First 0.5 mL Boiled Over

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area % Area of All Peaks Parts Per Million
0.17 4.53 1.75

0.47 6.95 2.69

0.54 9.16 3.55

2.94 (benzaldehyde) 124.91 48.36 6901

4.42 (benzyl alcohol) 83.68 32.40 4623

7.11 13.72 5.31

 

Table 2-11 shows a second 1 mL sample (Fraction B) boiled from the same feed

pot. It contains 29% and 133% of the respective concentrations of benzaldehyde and

benzyl alcohol in the original sample. Table 2-12 shows the composition of the

remaining 1.5 mL in the feed pot (Fraction C). It contains 0.6% and 36.4% of the

respective concentrations of benzaldehyde and benzyl alcohol in the original sample.

Note the high concentration of the component at 7.1 minutes. While some of it is being

carried over into the distillate products, it is still 16 times more concentrated in the

remaining bottoms than in the original regenerate. This indicates that this product is

being separated away from benzaldehyde and benzyl alcohol.

Table 2-11: Fraction B—A Second Distilled Sample of 1 mL

 

 

 

 

 

 

 

 

 

 

 

 

 

‘ Retention Time (min) Peak Area % Area of All 133$ Parts Per Million
0.17 1.61 1.62
0.19 1.27 1.28
2.93 (benzaldehyde) 22.76 23.00 1257
4.42 (benzyl alcohol) 61.77 62.42 3413
5.39 1.10 1.11
7.13 1.20 1.21

 

31

 

Table 2-12: Fraction C—The Remaining 1.5 mL of Bottoms

 

 

 

 

 

 

 

 

 

 

‘ Retention Time (min) Peak Area % Area of All Peaks fins Per Million
0.17 3.19 6.02
0.20 2.12 3.99
0.25 1.26 2.37
2.60 1.06 2.00
2.94 (benzaldehyde) 0.53 0.99 29
4.42 (benzyl alcohol) 16.88 31.82 933
5.40 4.48 8.45
7.13 18.26 34.41

 

 

 

 

 

 

One other issue to consider is the efficiency of this distillation step. Table 2-13
summarizes the amounts of benzaldehyde and benzyl alcohol in each sample. Table 2-14
gives each of these masses as percentages of that in the original sample. Note how
benzyl alcohol has more of a tendency to remain in the bottoms than benzaldehyde. This
is due to benzyl alcohol having a boiling point of 205.2°C while benzaldehyde has a
boiling point of 179.5°C (Lamer, 1996). In this boil-over step, about 74% of benzyl .
alcohol and 37% of benzaldehyde were recovered in a distillate (sum of Fractions A and
B) half the volume of the original sample. Some benzaldehyde and benzyl alcohol are
lost from this system due to possible side reactions such as oxidation to benzoic acid.
This might be explained by the peak at 5.4 minutes which begins to show up in Fractions
B and C (see Tables 2-1 1 and 2-12).

Table 2-13: Mass Distribution of Benzaldehyde and Benzyl Alcohol During Boil-Over

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Sample Benzaldehyde Massgg) Benzyl Alcohol Mass (g)
Original 0.0128 0.00768

A 0.00345 0.00231

B 0.00126 0.00341

C 4.35E-5 0.00140

Table 2-14: Sample Masses as Percentages of Original Feed

Sample Benzaldehyde Benzyl Alcohol

A 27.0 30.1

B 9.84 44.4

C 0.340 18.3

 

 

 

 

 

32

Part of the problem with the above experiment could have resulted from
benzaldehyde and benzyl alcohol sticking to glassware. To test this, a much larger
sample (100 mL) was completely boiled over. This sample consisted of regenerate
combined from many different adsorption runs. A large feed size would dwarf any
residual effects of glassware. Tables 2-15 and 2-16 give the respective GC profiles
before and after boil-over.

Table 2-15: 100 mL Regenerate Sample Before Boil-Over

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area
0.18 0.744

2.97 (benzaldehyde) 49.5

4.44 (benzyl alcohol) 34.6

5.42 4.37

7.15 3.02

8.55 1.70

10.96 0.354
13.82 - 1.62

Table 2-16: Regenerate Sample After Complete Boil-Over (Distillate)
Retention Time (min) Peak Area
0.17 2.12

2.96 (benzaldehyde) 46.7

4.45 (benzyl alcohol) 38.1

5.43 3.16

7.17 2.1 l

 

 

 

Note how the benzaldehyde peak area decreases by about 6% after boil-over. The
benzyl alcohol concentration increases by 10% after boil-over. This fluctuation is within
the range of normal variation for this GC. Nevertheless, the losses of benzaldehyde and
benzyl alcohol from the system are significantly less than in the previous experiment. It
should also be noted that the distillate has no color unlike the original regenerate which is

brown.

33

While boil-over helps to separate out some of the heavier impurities, it is still
desirable to isolate benzaldehyde and benzyl alcohol in more concentrated and pure
forms. Since they both form azeotropes with water, distillation techniques are limited.
Some experiments were consequently done in an attempt to extract these components into
a more concentrated form. Ethyl acetate has shown some promising results.

An experiment was done with a sample of regenerate after boil-over. Table 2-17
gives the GC profile of this sample. 10 mL of ethyl acetate was added to this 100 mL
sample. The mixture was shaken and allowed to separate into two liquid phases. The
phase boundary was somewhat unclear because a viscous white emulsion was present at
the boundary. The recovered ethyl acetate phase was about 2.5 mL in volume. The rest
of it was dissolved into the water phase. Table 2-18 gives a GC profile of the ethyl
acetate layer, and Table 2-19 gives a GC profile of the water layer. Ethyl acetate itself is
present in the GC at 0.5 minutes. Notice how benzaldehyde (3 min) and benzyl alcohol
(4.4 min) are respectively 14 and 11 times more concentrated in the ethyl acetate phase
than in the original sample. Furthermore, the ethyl acetate phase is 64 and 15 times more
concentrated with these two respective compounds than the water phase. However, note
how some of the other components also prefer the ethyl acetate phase (e. g. 5.4 min).
Also, the ethyl acetate phase takes on a brown color similar to that observed in the

regenerate prior to boiling.

34

Table 2-17: GC Profile of a 100 mL Regenerate Sample After Boil-Over

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area
0. 19 2.63

0.22 1.06

0.57 1.07

2.98 (benzaldehyde) 1 12

4.45 (benzyl alcohol) 40.7

5.43 3.03

7.16 2.40

10.98 1 . 10

 

Table 2-18: Ethyl Acetate Phase of First Extraction

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area
0.48 3638

3.08 (benzaldehyde) 1594

4.47 (benzyl alcohol) 464

5.43 41.5

7.15 14.6

10.97 1.06

11.53 3.40
Table 2-19: Water Phase of First Extraction '
Retention Time (min) Peak Area
0.49 551

2.95 (benzaldehyde) 25.1

4.43 (benzyl alcohol) 30.5

5.41 2.09

 

 

 

 

 

This is encouraging in light of the fact that 2.5 mL of ethyl acetate extracted 78%
of the benzaldehyde from the 100 mL water phase. However, only 25% of the benzyl
alcohol was removed from this water phase. These results show some logic in that
benzyl alcohol is a more polar molecule, thus, preferring water as its solvent.

Since fair amounts of benzaldehyde and benzyl alcohol remained in the water
phase after the first extraction, this 100 mL water phase was extracted again with 5 mL of
ethyl acetate. This time, none of the ethyl acetate was lost to the water phase as it was
already saturated. Tables 2-20 and 2-21 show the respective ethyl acetate and water

layers after extraction. Note how the ethyl acetate phase is 47 and 13 times more

35

concentrated in benzaldehyde and benzyl alcohol respectively. This extraction removed
51% and 25% of these two respective compounds from the water phase as can be seen by

the differences in Tables 2-19 and 2-21.

Table 2-20: Ethyl Acetate Layer after Second Extraction

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area
0.48 1272

0.53 2519

3.05 (benzaldehyde) 575

4.49 (benzyl alcohol) 299

5.46 22.3
Table 2-21: Water Layer after Second Extraction
Retention Time (min) Peak Area
0.49 596

2.94 (benzaldehyde) 12.2

4.43 (benzyl alcohol) 22.9

5.42 1.42

 

 

 

Since ethyl acetate is more volatile than the aromatic compounds of interest, a ‘
distillation was done on a mixed sample of the ethyl acetate layers from Tables 2-18 and
2-20. A 4.5 mL sample was used for this test, and the GC profile is given in Table 2-22.
Once again, ethyl acetate accounts for the peaks near 0.5 minutes. This experiment
involved boiling off the ethyl acetate while monitoring the benzaldehyde and benzyl
alcohol concentrations in the bottoms. Readings were taken when 1 mL (Table 2—23) and
0.2 mL (Table 2-24) of solution were left in the feed pot. Note that, in Table 2-24, the
benzaldehyde peak area was split between 3.13 and 3.21 minutes, probably due to a bad
injection. When the 4.5 mL sample is distilled down to 1 mL, benzaldehyde becomes
about 5.6 times more concentrated and benzyl alcohol becomes 5.93 times more
concentrated. These numbers would seem to imply that a chemical reaction is producing

more of each of these chemicals since we would expect the solution to be no greater than

36

4.5 times more concentrated with each of these chemicals if no reaction were occurring.
However, we are not completely certain whether the GC area is a linear function of
concentration in these concentration ranges. When the 1 mL bottoms is distilled down to
0.2 mL, benzaldehyde and benzyl alcohol show 1.4-fold and 1.6-fold peak area increases

respectively. It appears that the aromatics begin boiling off at about this point.

Table 2-22: Ethyl Acetate Solution to be Distilled

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area
0.48 1627

0.51 2203

3.05 (benzaldehyde) 1084

4.46 (benzyl alcohol) 387

5.42 34.0

7.13 12.4

11.5 2.94
Table 2-23: GC of 1 mL Bottoms
Retention Time (min) Peak Area
0.48 2158

3.19 (benzaldehyde) 6043

4.52 (benzyl alcohol) 2295

5.45 241

7.11 100

 

 

Table 2-24: GC of 0.2 mL Bottoms

 

 

 

 

 

 

 

 

 

Retention Time (min) Peak Area
0.48 597

3.13 (benzaldehyde) 4670

3.21 (benzaldehyde) 3702

4.53 (benzyl alcohol) 3707

5.30 , 247

5.44 348

7.06 300

 

 

 

 

 

It should be noted that the product is a dark brown color at this stage. As was
previously noted, this color first reappeared after the ethyl acetate extractions. This

seems to imply that some of the heavier components have not been completely separated

37

out in this sequence of extraction and distillation steps. Nevertheless, Table 2-24

represents the best product that the current techniques are able to yield.

38

REFERENCES

Hagedom, M.L.; “Differentiation of Natural and Synthetic Benzaldehydes by H Nuclear
Magnetic Resonance”; J. Agric. Food Chem, 1992, 40, 634.

Kawabe, T.; Morita, H.; “Production of Benzaldehyde and Benzyl Alcohol by the
Mushroom Polyporus Tuberaster K2606”; J. Agric. Food Chem, 1994,42, 2556.

Lamer, T.; Spinnler, H.E.; Souchon, I.; Voilley, A.; “Extraction of Benzaldehyde from
Fermentation Broth by Pervaporation”; Process Biochemistry, 1996, 31, 533.

Li, C.P.; Swain, E.; Poulton, J.; “Prunus Serotina Amygdalin Hydrolase and Prunasin
Hydrolase”; Plant Physiology, 1992, 100, 282.

Reilly, C.C.; Edwards, J .H.; Okie, W.R.; “Isolation and Characterization of Endogenous
Inhibitors of Nitrate Reductase of Peaches”; Journal of Plant Nutrition, 1986, 9, 1335.

Swain, E.; Poulton, J .; “Utilization of Amygdalin during Seedling Development of
Prunus Serotina”; Plant Physiology, 1994, 106, 437.

Zheng, L.; Poulton, J .E.; “Temporal and Spatial Expression of Amygdalin Hydrolase and
Mandelonitrile Lyase in Black Cherry Seeds”; Plant Physiology, 1995, 109, 31.

39

Chapter 3:
ABSORPTION MODELING WITH THE ESD EQUATION
OF STATE '

By
Aaron D. Soule, Cassandra A. Smith, and Carl T. Lira

ABSTRACT

The simple local density approach (SLD) is used to extend the ESD (Elliott etal.,
1990) equation of state to the modeling of gas adsorption on activated carbon, providing
significant improvement in quantitative modeling compared to the SLD approach using
the Peng-Robinson equation of state (Chen etal., 1997) or the van der Waals equation
(Rangarajan et al., 1995). Compared to the Peng-Robinson and van der Waals equations,
the ESD more accurately represents the contributions of attractive and repulsive forces,
and therefore provides increased accuracy when the attractive term is modified for
adsorption. Isotherms are represented for nine fluids on three different activated carbons
using temperature-independent parameters over temperature ranges of up to 167 K. The
fluid-solid interaction energies are shown to correlate with the fluid Lennard-J ones

parameter.

4O

INTRODUCTION

Adsorption is a topic which has been treated by a variety of models. Monte Carlo
simulations and molecular dynamics are computationally intensive methods for
calculating adsorption. Engineers often desire efficient methods for obtaining good
approximations. The best methods incorporate as few adjustable parameters as possible.
The Langmuir, Toth, and Freundlich models are easy to fit, but require temperature-
dependent parameters. The simplified local density (SLD) approach is an engineering
method that can be used with any equation of state and offers predictive capability with
only two temperature-independent adjustable parameters. This paper focuses on results
obtained with the ESD equation of state.

Previous studies have focused on adsorption modeling with the simplified local
density approach applied to the van der Waals (Rangarajan et al., 1995) and Peng-
Robinson (Chen et al., 1997) equations of state. These equations, while showing some
strengths, have characteristics which limit their ability to model fluid properties. For
example, the van der Waals equation is the most simple cubic equation and offers only
qualitative prediction. The Peng-Robinson, on the other hand, represents adsorption
accurately on flat surfaces. It is quite effective in modeling the adsorption of supercritical
fluids such as ethylene. It is capable of predicting the isotherm crossovers found in
experimental data. In porous materials, some success has been obtained (Chen et al.,
1997), but the general application to porous materials cannot generally fit both the
Henry’s law region and the adsorbent capacity. Further, when fitting data, it is difficult

to maintain good fits with the same parameters over large temperature ranges.

41

The ESD is also a cubic equation of state, however, it’s theoretical grounding is
superior in that the repulsive term is constructed to match computer simulations of
spheres and chains using a scalar shape factor to account for deviations from spherical
geometry. The attractive term consists of an expression for the spherical square-well
potential coupled with a shape factor correction. The improved subdivision of repulsive
and attractive forces is important for the SLD approach, which leads to the greater

accuracy of adsorption modeling.

42

SIMPLIFIED LOCAL DENSITY MODEL USING ESD

The Elliott, Suresh, Donahue (ESD) equation of state (Elliott et al., 1990) consists
of repulsive and attractive terms which are weighted differently than those in the Peng-

Robinson equation of state. The ESD equation takes the form Z=1+Z'°”+Z'm, where

 

z"? =—1j:’;n (3-1)
1+ 1.7745 < 711’ >

Here, c is a shape factor for the repulsive term, q is a shape factor for the attractive term,
I] is the reduced density (n=bp), b is the component’s size parameter, p is the density, Z
is compressibility, and Y is a temperature-dependent attractive energy parameter.
Although the ESD equation also can represent associating fluids, none of the components
presented in this paper have associative characteristics, so the terms are omitted from this

paper. The equation can also be represented in terms of fugacity as follows.

(3-3)

 

 

4 4m] 9511 9.59107 V
1n =——c 1-19 + 1+17745<Y > — —
f” 19 1n( 1)) (1—1972) 17745“ '7) (l+l.T745<Yn>) RT

V is the molar volume, T is temperature, R is the ideal gas constant, and fare is fugacity.
In the slit-shaped pores used in modeling, the fluid-solid interaction potential was
modeled using the same 10-4 potential as in previous work (Chen et al., 1997),

incorporating five carbon layers in the form of:

 

 

 

 

' 0.2 _ 0.5 _ 0.5 1
15m10 Eta" (Eta+a)‘
‘I’ 2 =47: maze > 3-4
‘0 p ” ”t 05 _ 0.5 _ 0.5 ( )
. (Eta+2a)‘ (Eta+3cz)4 (Eta+4a)4,

43

where or is the plane spacing between the solid particles (3.35A/ of, ), 0;, is the average of
the fluid and solid molecular diameters [07, = ( O'fi‘ + cs.) / 2], z is the particle position in
the slit relative to the carbon surface, Eta=(z+o,.)/of. is the dimensionless distance from
the carbon centers in the first plane, and pm,” represents the number of carbon-plane
atoms per square Angstrom (0.382 atoms/A2). The fluid-solid potential in relation to the
second wall, ‘I’2(z), can be calculated by replacing Eta in equation 3-4 with Xi, which is
the distance from the second wall (angstroms) divided by the fluid-solid diameter (see
Figure 3-1). The total potential is expressed as
‘11—; = ‘1’, + ‘15 (3-5)
The thermodynamic constraints of the adsorbing fluid fugacity can be estimated

by using equations 3-6 through 3-8 below (Chen et al., 1997).

 

.ubuuz = ”j(z)+flfi(z) (3'6)
y,(D=p°(n+RT1n[fj(f)] (3-7)
.ubua (T) = #0”) + RT1“[!}EOL] (3'8)

The local chemical potential due to fluid-fluid interactions, 113', is calculated
assuming that the local fluid-fluid fugacity can be estimated using the fluid-solid
potential, ‘1’1, and the bulk fluid fugacity, funk. In these equations pig”, is the bulk
chemical potential, P and 11° are the standard state fiigacity and chemical potential
respectively, and pf, is the fluid-solid contribution to the chemical potential. Note that pa

and fa (fugacity) are functions of 2 (position).

f, (2)4... exp[1‘%—T(’—’] (3-9)

For previous work with the Peng-Robinson equation, algebraic expressions were
developed for the Peng-Robinson attractive equation of state parameter, 8(Z)/abuik (Chen
et al., 1997). The same expressions are used to calculate the ESD Y(z)/Yum, as for the
Peng-Robinson a(z)/ab..1k based on these previous derivations. Thus, in slits, the local
fluid-fluid chemical potential (fugacity) is predicted based on Y(z) which can be
determined from Yum and the published algebraic expression for Y(z)/Yum.

In previous work, fluid closer to the wall than Eta=0.5 or Xi=0.5 (half the
diameter of a fluid particle) was ignored (Chen et al., 1997), assuming that no molecular
centers could be contained in this area. This manuscript assumes that the density cutoff
should be the point where the local fugacity, fly, is one tenth of a percent of the bulk
fugacity, fun. This yields a more realistic density profile near the wall. From Equation
3-8, one can calculate Y(z)/k as approximately 2500 K when fff(Z)/fb=o.001 at 373 K.
This value is used for each temperature calculation, which makes the cutoff distance
dependent on the slit width only. This is merely an approximation, but the density below
this region is always low in all calculations we have checked. For the density
calculations in the region Eta<0.5 or Xi<0.5, we use Y(z)/Ym=0.5.

The local density is obtained at each 2 by using equations (3-3), (3-5), and (3-9).
The difference between the local and bulk densities is integrated in the correct geometric

form over the slit distance using the modified Simpson’s rule to yield excess, 1'“.

r“ = Aj[p(z)—p..t1dz (3-10)

45

The variable A is the surface area per unit weight of adsorbent (e. g., square meters per
gram). In the case of adsorption in a slit with homogeneous parallel walls, this
integration over the entire slit width is divided by two since two walls contribute to the
surface area of a slit.

For each fit discussed in this paper, the value of A was taken from the cited
reference and not used as an adjustable parameter in the model. The values of as are
tabulated Lennard-J ones diameters of each fluid, and 0.. is the reported diameter of
carbon (Reid et al., 1987). In calculating adsorption in slits, two adjustable temperature-
independent parameters were fitted: eg/k (fluid-solid interaction potential in Kelvin) and
H (slit width in angstroms). The parameters were fit to optimize the simultaneous
representation of all data in a given figure rather than optimization of individual

isotherms.

 

 

 

Figure 3-1: Model of a slit-shaped pore showing the variables used to define distances in
the manuscript; Eta = (z + 0..)/or., Xi = (H - Eta’o..)/of.

47

RESULTS AND DISCUSSION

Several sets of pure component adsorption data have been successfully fitted with
the ESD version of the simplified local density model. Table 3-1 lists the pure
component ESD parameters, all of which are obtained from bulk fluid properties. Figure
3-2 shows the adsorption of ethylene on BPL carbon over 167 K. Figure 3-3 shows
ethane adsorption which is also a good fit over 167 K. Other fits include butane over 110
K (Figure 3-4), propane over 167 K (Figure 3-5), methane over 89 K (Figure 3-6),
propylene over 139 K (Figure 3-7) and nitrogen over 111 K (Figure 3-8). It is important
to note that all these fits were performed by simultaneously optimizing all of the
isotherms in a given graph through adjustment of the two parameters. In the cases of
butane and propylene, it seems as though the model has trouble fitting low pressure data.
Better fits seem to be generated when the pressure ranges are wider. Nitrogen shows
inaccurate predictions of temperature dependence. The fit of propane is weak in that, at
the highest temperature, the knee region is overpredicted. While butane, propylene, and
propane deviate significantly from the sphericity assumption of the model, it is not
known why the nitrogen data are not fit more precisely. Figures 3-9 and 3-10 show that

the model can also represent acetylene and carbon monoxide.

48

Table 3-1: Pure Component ESD Parameters

 

 

 

 

 

 

 

 

 

 

 

Component c q gfl/k (K) b (cmT/mole)
acetylene 1.6808 2.2967 190.510 13.053
butane 1.7025 2.338 260.583 29.039
carbon monoxide 1.2367 1.4509 103.784 10.171
ethane 1.3552 1.6765 220.449 16.716
ethylene 1.305 1.581 210.275 15.013
methane 1.0382 1.0728 178.082 10.863
nitrogen 1.1433 1.273 106.155 9.907
ropane 1.5481 2.0441 241.433 22.921
propylene 1.5142 1.9794 241.896 20.890

 

 

 

 

 

 

There is a correlation between the pure fluid Lennard-J ones parameters and fluid-
solid interaction energy parameters. Figures 3-11 and 3-12 graphically demonstrate this
relationship for Columbia Grade L and BPL carbons respectively. The fluid-solid
parameters do show some variation between different activated carbons for the same
component. This effect is more pronounced for some components than for others, but in
the case of this study, the parameters do not vary more than 10 K for a single component.

When examining the ESD equation of state, it is also important to consider the
validity of the bulk properties. Since the adsorbed phases have liquid-like densities, the
representation of liquid molar volumes is important. Table 3-2 gives some saturation
bulk property comparisons between the Peng-Robinson and ESD equations near each
component’s critical temperature. This particular sample of bulk properties was chosen
to give the reader a general perspective on the performance of the ESD as compared to
experimental data and the Peng-Robinson equation of state. Since temperature and
pressure are specified, molar volume is the dependent variable used for comparison.
While the ESD appears to predict gas volumes better than the Peng-Robinson, the ESD is

weaker in predicting liquid properties. As can be seen in the data, the ESD volumes vary

49

from experimental volumes by 20% or more for liquids. Gas volumes vary by up to 10
%. The ability of the ESD to more accurately model adsorption over wide temperature
ranges is attributed to the superior representation of the individual contributions of

attractive and repulsive forces, since the bulk liquid properties are predicted with about

the same or less accuracy.

Table 3-2: Comparison of experimental saturation molar volumes (Starling, 1973)
to the bulk Peng-Robinson and bulk ESD predictions. All data points are near the
critical temperature. '

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Peng-
Tempcrature Reduced Pressure E Robinson ESD
(K) Temp. (MPa) cm Ignol cm’Igmol cm’lgrpol
Methane (L) 188.7 0.99 4.32 70.12 83.30 91.84
(G) 162.51 153.74 160.48
Ethane (L) 302.6 0.99 4.62 102.44 121.00 134.00
(G) 235.98 217.18 231.33
Propane (L) 363.7 0.98 3.81 143.94 155.28 175.33
(C!) 368.35 356.55 382.21
N-Butane (L) 422 0.99 3.62 213.63 219.00 247.39
(G) 371.86 379.11 414.71
Ethylene (L) 277.6 0.98 4.51 86.88 101.23 113.41
(G) 240.70 231.54 245.65
Propylene (L) 360.9 0.99 4.29 131.62 152.18 172.04
(G) 300.22 298.49 321.69
Nitrggen (L) 124.8 0.99 3.17 63.80 74.36 81.28
(G) 146.31 140.28 147.06

 

'50

 

 

8 1* :
212.70K . .

 

301.40 K

Adsorption (mmol/g)

1

-ii—
~41-

 

di-

 

 

 

0 r r 1 1
I I 1 T I

0 0.2 0.4 0.6 0.8 1 1.2 1.4
Pressure (MPa)

Figure 3-2A: Adsorption of ethylene on BPL activated carbon (988 square meters
per gram) where H=l3.7 angstroms and efilk=103 K. Data of Reich, et al., 1980.

51

1.6

1.8

 

310.92 K
0

394.20 K

Adsorption (mmol/g)

 

O
-4--
—1
—1i-
_4
qt-

477.59 K

  
       
  

 

 

0 0.2 0.4 0.6 0.8 1
Pressure (MPa)

Figure 3-ZB: Adsorption of ethylene on Columbia Grade L carbon (1152 sq.
meters per gram) where H=13.7 angstroms and enlk=104 K. Data ofRay
andBox, 1950.

52

4L

1.4

1.6

Adsorption (nunollg)

 

 

 

 

    
 

 

 

 

 

6 q

5 a

4 .

3

2

1 A.

0 i i r e i e i i
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6

Pressure (MPa)

Figure 3-3A: Adsorption of ethane on BPL activated carbon (988 sq. meters

per gram) where H=14.2 angstroms and edk=102 K. Data of Ray and Box,
1950.

53

1.8

 

310.92 K

394.20 K

Adsorption (mmol/g)

0

477.59 K

 

0 0.2 0.4 0.6 0.8 1 1.2 1.4
Pressure (MPa)

Figure 3-38: Adsorption of ethane on Columbia Grade L carbon
(1152 sq. meters per gram) where H=14.3 angstroms and eg/k=104
K. Data of Ray and Box, 1950.

54

     
   
    
  

 

 

366.48 K

422.00 K

Adsorption (mmol/g)

   
 

 

I

1

 

 

O 0.02

T

0.04

f

0.06

'1?-

0.08

Pressure (MPa)

Figure 3-4: Adsorption of butane on Columbia Grade L carbon (1152 sq.
meters per gram) where H=14.1 angstroms and stk=l60 K. Data of
Ray

and Box, 1950.

55

0.1

Adsorption (mmol/g)

 

  

477.59 K

   

 

 

1 L r I l 1 1

 

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

Pressure (MPa)
Figure 3-5: Adsorption of propane on Columbia Grade L activated carbon

(1152 sq. meters per gram) where 11:15.6 angstroms and 5,, =114 K.
Data of Ray and Box, 1950.

56

Adsorption (mmol/g)

 

 

 

 

0 0.5 1 1.5 2 2.5 3 3.5 4

Pressure (MPa)

Figure 3-6: Adsorption of methane on BPL activated carbon (988 sq.

meters per gram) where H=12.1 angstroms and efi/k=73 K. Data
from Reich, et. al., 1980.

57

 

4.5

Adsorption (mmol/g)

 

4.5 -

 

 

366.48 K

 

I 1

—u.

 

 

T I

0.02 0.04 0.06 0.08 0.1
Pressure (MPa)

4L—

Figure 3-7A: Adsorption of propylene on Columbia Grade L carbon (1 152 sq.

meters per gram) where H=14.3 angstroms and ef,/k=126 K. Data of Ray and
Box, 1950.

58

0.12

Adsorption (mmol/g)

 

      
 

 

 

 

5 a-
303.15 K
4 4-
3 1 323.15 K
2 «1-
1o
1 41o
0 41 i i i i
0 0.02 0.04 0.06 0.08 0.1

Pressure (MPa)

Figure 3-7B: Adsorption of propylene on BPL activated carbon (1050-1150 sq.

meters per gram) where H=14.3 angstroms and egik=132 K. Data ofLaukhuf,
ct al., 1969.

59

0.12

Adsorption (mmol/g)

 

2.5 *~

1.5 --

0.5 +

   
 

298.15 K

 

 

 

J I L 1 I I
I I I I I

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09

Pressure (MPa)
Figure 3-7C: Adsorption of propylene on Black Pearls 1 carbon (705 sq.

meters per gram) where H=14.0 angstroms and 8;, =115 K. Data of Lewis,
et al., 1950.

60

0.1

Adsorption (mmol/g)

 

2.5

 

     

 

 

2 a-
1.5 ~-
310.92 K

1 .1.
0.5 ~-

0 1 1 1 1 1 1 1

0 0.2 0.4 0.6 0.8 1 1.2 1.4
Pressure (MPa)

Figure 3-8: Adsorption of nitrogen on Columbia Grade L carbon (1152 sq. meters per
gram) where H=12.6 angstroms and ef,/k=60 K. Data of Ray and Box, 1950.

61

1.6

 

2.5 ‘1'

Adsorption (mmol/g)
Cl

 

O

 
 
  

394.20 K

 

 

.11-
,_

 

 

0 i i i
0 0.02 0.04 0.06 0.08 0. 1
Pressure (MPa)

Figure 3-9: Adsorption of acetylene on Columbia Grade L carbon (1152 sq.
meters per gram) where H=14.5 angstroms and sfslk=103 K. Data of Ray and

Box, 1950.

62

0.12

 

3.5

 

 

 

3 at O
2.5 s.
30
_. 2_ 0
g 2
§ 310.92 K
.S
g 15 ._ 366.48 K .
“:3
<
O
1 -_
477.59 K
O
0.5 +
O
0 i 1 1 i 1 i i 1 i
0 0.2 0.4 0.6 0.8 i 1.2 1.4 1.6 1.8
Pressure (MPa)

Figure 3-10: Adsorption of carbon monoxide on Columbia Grade L carbon (1152 sq.
meters per gram) where H=15.8 angstroms and efslk=70 K. Data of Ray and Box,
1950.

63

£er

 

180

 

      
 
   

 
 

 

 

160 «r O
utane
140 -~
propylene
e
120 1" propane
ethan
100 " acetylene
80 «~
0
60 22 . carbon monoxide
nitrogen
40 -»
20 a»
0 1 1 1 1 41
0 100 200 300 400 500
Err/k

Figure 3-11: Correlation of fluid-solid and fluid-fluid interaction potentials. Fluid-
fluid parameters from Reid, Prausnitz, and Poling, 1987. Columbia Grade L carbon.

 

8171: (K)

 

140

120 -_

100 «r

80 «r

60‘

204

Propylene

 
  

Ethylene

   

Methane

 

 

 

di-

50

«i-

100

I I
I

—11—

l
I

150 200 250 300
811’ K (K)

Figure 3-12: Correlation of fluid-fluid and fluid-solid interaction
potentials. Fluid-fluid parameters from Reid, Prausnitz, and
Poling, 1987. BPL carbon.

65

350

SUMMARY

This paper has shown that the ESD equation of state can be adapted for successful
adsorption modeling. While nitrogen and the higher molecular weight components
presented in this paper can’t be fitted as well as others, it is important to note that good
engineering approximations can still be made with any of these compounds.

Our goal is to fine-tune the approach to make it amenable for use in process
simulation software with multicomponent systems. Further work in this area would
involve looking at the adsorption of mixtures--for example, the prediction and modeling
of azeotropic adsorption. We would like to extend the ESD to modeling of zeolites.
Also, more work needs to be done with hydrogen-bonding fluids such as water, as our
efforts in applying the ESD in this particular area are not yet quantitative (Smith, 1997).
Supercritical fluid adsorption is another subject of immediate interest since adsorption

isotherms exhibit crossovers that can be represented with the SLD approach.

ACKNOWLEDGEMENTS

We express appreciation to J. Richard Elliott for sharing the ESD code which was
adapted for modeling and also acknowledge the MSU Crop and Bioprocessin g Center for

support of Cassandra Smith and AaronSoule.

66

LITERATURE CITED

Chen, J.; Lira, C.T.; Orth, M.; Subramanian, R.; Tan, C.; Wong, D. “Adsorption and
Desorption of Carbon Dioxide onto and from Activated Carbon at High Pressures”;
Ind. Eng. Chem. Res. 1997, 36, 2808-2815.

Elliott, J.R.; Suresh, S.J.; Donohue, M.D., Ind. Eng. Chem. Res., 1990, 29, 1476.

Elliott, J.R.; Lira, C.T. Introductory Chemical Engineering Thermodynamics, First
Edition; Prentice Hall, 1998.

Laukhuf, W.L.S. and Plank, CA. J. Chem.Eng.Data. 1969, 14, 48.

Lewis, W.K.; Gilliland, E.R.; Chertow, B.; Hoffman, W.H. J.Am.Chem.Soc.,
1950, 72, 1153.

Rangarajan, B.; Lira, C.T.; Subramanian, R. AIChe Journal, 1995,41, 838-845.
Ray, GC. and Box, E.O. Ind.Eng. Chem. 1950, 42, 1315.
Reich, R.; Ziegler, W.T.; Rogers, K.A. Ind.Eng.Chem.Process.Des.Dev., 1980, 19, 336.

Reid, R.C.; Prausnitz, J.M.; Poling, B.E. The Properties of Gases and Liquids, Fourth
Edition; McGraw-Hill: New York, 1987, 734.

Smith, Cassandra A.; “Adsorption of Water on Carbon: A Study Using the ESD
Equation of State”; Master’s Thesis, Michigan State University, 1997.

Starling, K.E. “Fluid Thermodynamic Properties for Light Petroleum Substances”, Gulf
Publishing Company, Houston, TX, 1973.

67

CHAPTER 4:
CARBON DIOXIDE: A STUDY OF SUPERCRITICAL
FLUID ABSORPTION

Carbon dioxide is a fluid which has received much attention in adsorption studies
as will be shown in this chapter. Chapter 2 is a typical example of its usefulness.
Furthermore, the numerous methods of producing this gas make it cheap, and thus,
favorable for use in chemical processes. A variety of experimental adsorption data is
available in the literature, much of which at its supercritical conditions (above 304 K).
Some of the different data will be examined and fitted with the SLD model in this
chapter. The carbons presented here exhibit a wide range of surface areas depending on
type: 983 to 1699 meters squared per gram. It can be seen that the fluid-solid interaction
parameters vary widely with different carbons. They range from 77 K on DeGussa IV to
109 K on ACK carbon. The carbon dioxide Lennard-Jones diameter is 3.941 angstroms,
taken from Reid, Prausnitz, and Poling. The fluid-solid interaction parameters and the
slit width are optimized in each of these data sets.

Four different data samples have been fitted in this chapter at temperatures
ranging from 212.7 K to 394.2 K. Figure 4-1 shows carbon dioxide adsorbed on
DeGussa IV activated carbon at pressures up to 15 MPa. The isotherms at 284 K and
300 K exhibit adsorption maxima at 6.5 and 5 MPa respectively. The 324 K isotherm
exhibits a maximum at 4 MPa; however, there is little data for analyzing its behavior
beyond this point. Note that the peaks shown at the two higher temperatures were not
fitted. Figure 4-2 shows carbon dioxide adsorption on Columbia Grade L carbon. The
isotherms are more straight because the data set does not go significantly beyond the

Henry’s Law region. The optimized parameters give a good fit to this data, but these

68

same parameters might not give a good fit at higher pressures. The parameters are only
valid for the pressure range in which they are fitted. Figure 4-3 shows a data set on BPL
carbon which goes beyond the Henry’s Law region at 260.2 K and 301.4 K but not far
enough to reach any maxima. Figure 4-4 shows a set which goes up to 15 MPa on ACK
carbon. Note that, near 8 MPa, the isotherm at 313.2 K crosses over the isotherms at
333.2 K and 353.2 K. The 333.2 K and 353.2 K isotherms cross over one another at
about 10.5 MPa. The SLD model agrees quite well with each of these crossover points in
the data.

One special issue of concern in these fits is carbon surface area sensitivity. Each
surface area is estimated with some margin of error. Changing the surface area does have
an effect on the optimal parameters for a carbon-gas system. A quantitative analysis was
performed on the ACK carbon system reported by Ozawa as shown in Table 4-1. The
reported surface area (983 m2/ g) was adjusted 10% in both the positive and negative
directions. This test shows H and 1:2/k to increase with a decrease in surface area.
Conversely, they both decrease with an increase in surface area. Note that eg/k changes
by about the same amount in both the positive and negative directions while H tends to
change most significantly with the decrease in surface area. Studies such as this should
be considered depending upon the degree of uncertainty in a reported surface area.

Table 4-1: Parameter Sensitivity to ACK Carbon Surface Area

 

 

 

 

 

Area (ml/g) H (angstroms) ¢eg/k (K)
885 16. 1 1 18

983 15 .4 109
1081 15 .2 99

 

 

 

69

 

Adsorption (mmol/g)

 

 

 

 

70 1 9
0
60 r 324 K 3
50 ~- : g 300 K
40 O O
8 O
30 -- s .
20 -- . °
° 0
0
10 -~ ,o ‘ 284 K
0
0 k i i 1 1 i i 41 .
o 2 4 6 s 10 12 14 16
Presmre (MPa)

Figure 4-1: Carbon dioxide adsorption on DeGussa IV Activated
Carbon (1699 meters squared per gram) where H=15.5 A and
efi/k=77 K. Data of Chen, et al., 1997.

70

Adsorption (mmol/g)

 

1.8

1.6 r

1.4 1»

     
    

1.2 ‘-
310.92 K

0.8 4*-

0.6 --

0.4 2
0

394.20 K
0.2 --

 

 

«L-

O 0.02 0.04 0.06 0.08 0.1
Pressure (MPa)

Figure 4—2: Carbon dioxide adsorption on Columbia Grade L Carbon
(1152 meters squared per gram) where H=11.7 A and cfjk=86 K.
Data of Ray and Box, 1950.

71

 

0.12

Adsorption (mmng)

12

 

    
 

212.70 K
260.20 K

 

 

   

 

 

8 -
. e e
6 301.40 K
4
2
0 . 1 1 1 1 1 1 1 1
0 0.5 l 1.5 2 2.5 3 3.5 4

Pressure (MPa)

Figure 4-3: Carbon dioxide adsorption on BPL carbon (988 meters
squared per gram) at H=14.2 angstroms and eg=99 K. Data of Reich, et
al., 1980.

72

4.5

Adsorption (mmoilg)

10

 

 

 

 

 

 

 

 

 

 

 

 

 

 

o 5 f is g {0
Pressure (MPa)

16

 

 

 

0 353.2 K (exp)

I 33.2 K (exp)

A 313.2 K (exp)
11111 353.2 K (calc)
- -- - 333.2 K (calc)
313.2 K (calc)

 

Figure 4-4: Adsorption of Carbon Dioxide on ACK Carbon (983 meters
squared per gram) at H=15.4 and a,.lk=109. Data of Ozawa etal., 1974.

73

 

REFERENCES

Chen, 1.; Lira, C.T.; Orth, M.; Subramanian, R.; Tan, 0; Wong, D. “Adsorption and
Desorption of Carbon Dioxide onto and from Activated Carbon at High
Pressures”; Ind Eng. Chem. Res. 1997, 36, 2808-2815.

Ozawa, S.; Kusumi, S.; Ogino, Y.; From Proceedings of the Fourth International
Conference on High Pressure. Kyoto, 1974.

Ray, GC. and Box, E.0.; IndEngChem. 1950, 42, 1315.

Reich, R.; Ziegler, W.T.; Rogers, K.A. IndEng.Chem.Process.Des.Dev., 1980, 19,
336.

Reid, R.C.; Prausnitz, J.M.; Poling, B.E. The Properties of Gases and Liquids, Fourth
Edition; McGraw-Hill: New York, 1987, 734.

74

Chapter 5:
ABSORPTION OF BINARY MIXTURES

INTRODUCTION

In chapter 3, the SLD model was applied to adsorption of pure gases on activated
carbon. This chapter attempts to extend that theory to mixtures modeled by the ESD
equation of state. This particular work was started by Ramkumar Subramanian and was
based on Pong-Robinson modeling. This current work takes the structure of his original
FORTRAN program and converts it to one using ESD fluid properties. While the basic
concepts do not change much from the modeling of pure fluids, there are mixing rules
which must be obeyed in order to obtain correct firgacity values. This model uses three
adjustable parameters as opposed to two for a pure system. The interactions of each
component must be accounted for. Also, fugacity is solved for by satisfying a set of three

objective functions at each local position.

DISCUSSION

Elliott and Lira (1999) have developed all of the ESD mixing rules for binary
systems along with a corresponding expression for the firgacity coefficient. In a binary
system, each component is associated with a series of four tabulated ESD parameters:
(e/k)a, c (dimensionless), q (dimensionless), and b. The expression for the fugacity
coefi'rent of component A, 11111., is as follows. This equation assumes that no association
occurs.

111(¢,,)=T1+T2+T3+T4+T5 (54)

where

75

Tl= —icA Ina—197)) (5-2)

 

 

 

19
4cbp
T2=fil—9; (5'3)
95 ln(l+l.7745<Yn> " <qu>
T3 - .Y-b , b. — Yb 5-4
17745 <Yb> Ex] ‘4 ”1+ 1"”) <Yb> ‘ A ( )
T4=_9.5<qu> YAbAp (5_5)
<Yb> 1+1.7745<Y77>
and
T5=—ln(Z) (5-6)

The variables are defined as follows: p is density, Y is the ESD attractive parameter, and
n=bp. Z is the compressibility solved from the cubic form of the ESD equation of state.

This cubic form is given as follows.

Z3+AA*Z2+BB*Z+CC=O (5-7)
where

AA=I7745*<Y*B>—l—L9*B (5-8)
BB = —1.9*1.7745*B*< 1H3 > —1.7745*< we >+1.9*B—4*c*B+9.5*<q*Y*B > (5-9)
and

CC =1.9*1.7745*B*< Y‘B > -4*1.7745*c*B*< r13 > -95*l.9*B*<q*Y*B > (5-10)

For the above equations, B=(b*P)/(R*T). P is pressure, R is the ideal gas constant, and T
is temperature. Furthermore, the following rules must be applied in calculating the
mixture parameters from the pure component parameters. The variable x represents mole

fraction.

13
b=2x,b, (5—11)
i=A

76

c: Ergo, (5'12)

8
q =inqr (5-13)
1:11
K = exp(£,- /kT)-1.0617 (5-14)
YAB = CXpU‘AB lkT)‘ 1.0617 (5'15)
3 a
<qu >=p<qu >=pZZx,x,Y,,(b,qj +b,q,)/2 (5-16)
i=Aj=A
8A8 =‘/£A£B (l-kAB) (5'17)
3
<Yn>=p<Y*b >=pr,b,-Y; (5-18)

i=A
Note that equations 5-14 and 5-15 calculate Y values for the bulk phase only. In order to
get Y values at a local position, expressions have been developed for Y/Ym in a paper
by Chen and coworkers (1997). In that paper, they are listed as slam. because the ratio
expressions of the attractive parameter were originally derived using the Pong-Robinson
equation. It is assumed that they also apply to the ESD equation of state.

Table 5-1 lists the parameters associated with the components under discussion in
this manuscript along with their Lennard-Jones fluid diameters (Reid et al., 1987). Also,
each mixture system has a mixture parameter, k;,-,, which has been optimized for the best
EOS fit to bulk data at a relevant temperature (Table 5-2). The data of Miller and
coworkers (197 7) were used for the methane-ethane and methane-ethylene systems. The
data of Ng and Robinson (1978) were used for the toluene-carbon dioxide system. The
toluene fluid-fluid diameter is set equal to that of benzene for this study (Reid et al.,

1987) since no value is readily available in the literature.

77

Table 5-1: Pure component fluid-fluid diameters and ESD parameters

 

 

 

 

 

 

 

m nent gg (e/k )g (K) b (cmj/mol) g 9
Carbon Dioxide 3.941 178.269 10.534 1.8321 2.585
Ethane 4.443 220.429 16.716 1.3552 1.6765
Ethylene 4.163 210.275 15.013 1.305 1.581
Methane 3.758 178.082 10.863 1.0382 1.0728
Toluene 5.349 332.752 36.227 1.9707 2.849

 

 

 

 

 

 

 

One large assumption in this particular model is particle size homogeneity. Each
mixture system uses an effective particle size, 01,-, which is simply the average particle
size of the two pure components. These values are also listed in Table 5-2. This
assumption allows for the WY 11111: values to be equivalent for each component. Ifthe
particle sizes were different, the WY 11111 expressions for each component would be
dependent on different definitions of the local position, Eta (in fluid-solid diameters).
Also, separate exclusion regions for each component would have to be defined at the
edges of the slit (regions where a particle is unable to realistically be present). At each
end of the slit, these exclusion regions are equivalent to half of a fluid-solid diameter.
Thus, separate particle sizes would require different regions of adsorption integration for
each component. Since the work in this manuscript is of a preliminary nature, firture
work will involve formulating a more detailed mechanism for incorporating the
individual particle size of each component.

In calculating wall-particle potentials, the diameter of the carbon itself is also
important. This value, on, is 3.4 Angstroms (Subramanian, 1995). This gives rise to
another important value, the fluid-solid diameter (on), which is also tabulated in
Table 5-2. This value is simply the average of the carbon diameter and the effective

mixture fluid-fluid diameter.

78

Table 5-2: Mixing Parameters and Effective Mixture Diameters

 

 

 

 

 

 

 

 

 

Mixture g3 92;. kg

Ethane-Methane 4.1 3.75 0.0094
Ethylene-Methane 4.0 3.70 0.0209
C02-Toluene 4.6 4.00 0. 1058

 

The algorithm for calculating adsorption is quite similar to that presented in
Chapter 3. However, in this case, fugacity expressions are being matched for two
different components (A and B) instead of just one at any local position or the bulk
phase. Also, the Y/thuc expressions cited in the paper by Chen and coworkers (as a/amk)
must be used to solve for Y at each local position. These expressions are functions of
local position and fluid particle size. As in Chapter 3, Eta represents the position in the
slit in terms of fluid-solid diameters. The following equations define the 10-4 Lennard-

Jones potentials relative to the first wall for each component:

 

 

 

 

 

02 _ 05 _ 05

\111(2) = 47mm “3.16M 251010 051%“ -(Eta 3:)4 - 05 (5-18)
(Eta + 2a)4 (Eta + 3a)4 (Eta +46!)4
0.2 _ 0.5 _ 0.5

‘1’13(z) = 4npmoiflsfij Eta” 0530‘ (Eta $56!)“ 05 (5-19)

 

- (Eta +2a)‘ — (Eta +3a)‘ — (Eta +4a)‘
The above two equations are identical due to the assumption of uniform particle diameter
although each component will still maintain its fluid-solid interaction parameter (e2).
Potentials relative to the second wall (‘1’21L and W23) can be obtained by replacing Eta, the
distance from the left wall to the local position, with Xi, the distance from the right wall
to the local position of interest. Total potentials for each component can then be obtained

from the following two equations:

79

‘l’TA = ‘PlA + “’24. (5'20)
‘PTB = “’18 ‘l' W213 (5'21)
Using these potentials, we can then proceed to define local fugacities relative to the bulk

fugacities for each component as shown in the following equations:

fry-.16?) = fa.“ exP[———-ly]:‘ (2)] (5-22)
T
ruck/.1... 41:14:] (5-23)

Fugacity is also defined as the product of a fugacity coefficient (41) and local pressure (P).
Thus, a series of three objective functions will now be introduced which can be solved for

local composition (x11 and x3) and local pressure (P):

61 = XA*¢A*P / fmA - l = 0 (5-24)
G2 = X31031? / fa]; — 1 = 0 (5-25)
G3 =1—xA—x3 =0 (5-26)

The fugacity coefficients are calculated by the ESD equation of state and are functions of
composition and pressure. The coefficient expression and its pertinent mixing rules have

been calculated by Elliott and Lira as shown in Equations 5-1 through 5-6. The following

expressions define the integration necessary for producing adsorption values:

1‘? = AIIJ‘AMZ) " xAbulkpbulk 10'? (5‘27)

I“? = AI[pr(z) " xB,bulkpbqu l‘t (5'28)

80

A is the surface area in meters squared per gram of activated carbon, 2 represents the
local position in the slit, and p is density. The units of adsorption are millimoles

adsorbed per gram of activated carbon.

RESULTS

As previously noted, data sets for three difi’erent mixtures have been fitted using
this theory. While the results are significantly better than those of Subramanian (1995),
which incorporates the Peng-Robinson equation, there is still much room for
improvement. Three adjustable parameters went into each fit: the fluid-solid interaction
potentials for each component «ea/k» and (an/k» in Kelvin) and the carbon slit width
(H in Angstroms). The component interaction potentials were all obtained from the pure
data fits in chapter 3 with the exception of toluene. Due to a lack of literature data, the
toluene fluid-fluid diameter is estimated to be 5.3 angstroms while its fluid-solid
interaction potential is estimated to be 250 K. Since the ethylene-methane and ethane-
methane data sets were obtained from the same carbon (BPL), an average slit width of 14
angstroms was used based on results from the pure systems in Chapter 3. As mentioned
earlier, one weakness in the SLD theory is that different compounds do not always
predict the exact same slit width for a single carbon.

While these parameters (except for toluene) have been optimized to pure
component data, they have not been re-optimized to mixture data. This is a project which
might be useful to pursue in the future; however, one aim of this particular work is to
extrapolate pure component properties into mixture characterization. The values of the

adjustable parameters are listed in Table 5-3. Note that there are different values for the

81

same component within different mixtures. This is due to the fact that the fluid-fluid
potential for each component was fit using the effective particle diameter of the mixture

rather than the tabulated diameter as in chapter 3.

Table 5-3: Parameters Used In Mixture Calculations
e/k H

Methane 4. 1 66 14
Ethane 4. 1 14

Methane 4.0 68 14
4.0 108 14

Toluene 4.6 250 20
4.6 73 20

 

Figures 5-1 through 5-7 contain methane-ethane fits at 301.4 K, 260.2 K, and
212.7 K on BPL carbon with a surface area of 988 m2/gram. Note that the methane
adsorption is consistently under-predicted. However, the model seems to give an
accurate fit for ethane in most cases. Ethane clearly is the more strongly adsorbed
component. This behavior does not seem to change much over the 90 K temperature
range. The homogeneous particle size assumption is likely to be a major cause of the
error in methane prediction.

Figure 5-8 contains a methane-ethylene fit at 212.7 K on the same BPL carbon.
Once again, methane is under-predicted, but the model fits well to ethylene.
Consequently, the under-prediction of methane results in an over-prediction of ethylene’s
mole fraction in the adsorbed phase as shown in Figure 5-9.

Figure 5-10 shows toluene—C02 isotherms at 308 K, 318 K, and 328 K. With k1,-

equal to 0.1058, it is clear that this is a non-ideal mixture. The SLD model (lines with

82

symbols) seems to consistently over-predict the data of Tan and Liou (lines without
symbols). The predictions seem good at low pressures but then diverge at higher
pressures. There seems to be some slope agreement between the data and the predicted
values over the 20 K temperature range. Part of this error could be due to inaccuracies in
the estimated toluene fluid-solid interaction parameter. Also, there is a comparatively
large difference in the Lennard-Jones particle diameters. Carbon dioxide is 3.941

angstroms while toluene is estimated at 5.3 angstroms.

83

Mlllimoles adsorbed per gram of carbon

 

 

 

 

 

6
5 , I
ethane

4 _

I
3 ..
2 ..
1 _

methane °
0

e

O l l I I
0 0.5 1 1.5 2 2.5
Pressure (MPa)

Figure 5-1: Adsorption of Methane-Ethane on BPL Activated Carbon at
301.4 K with Bulk Ethane Mole Fraction of 0.733 (Data of Reich, et al.,

1980).

84

Adsorption (mmoilg carbon)

 

4.5 -

ethane

9°
01
1

(a)
1

1°
0'1
1

N
1

1.5 -

1 . methane

V

 

 

 

0 1 1
0 0.5 1
Pressure (MPa)

Figure 5-2: Adsorption of Methane-Ethane on BPL Carbon at
301.4 K and Bulk Ethane Fraction of 0.501 (Data of Reich, et
al.,1980).

85

1.5

Adsorption (mmollg carbon)

 

 

 

 

3.5 -
l
3 .1
I
Ethane
2.5 -
2 -
I
1.5 -
s
1 r . Methane
e
0.5 -
e
o I I
0 0.5 1 1.5
Pressure (MPa)

Figure 5-3: Adsorption of Methane-Ethane on BPL Carbon at
301.4 K and a Bulk Ethane Fraction of 0.255 (Data of Reich,
et al., 1980).

86

Adsorption (mmollg)

 

 

 

 

 

7 - I
/
I
Ethane
5 4
5 -
4 .
3 ..
2 -1
1 _
Methane
0 . ‘L 1 ° . O
0 0.1 0.2 0.3 0.4
Pressure (MPa)

Figure 5-4: Adsorption of Methane-Ethane on BPL Carbon at 212.7 K
and a Bulk Ethane Mole Fraction of 0.733 (Data of Reich, at al., 1980).

87

Adsorption (mmollg)

 

 

 

 

Methane

 

 

0.2

I

0.4

Pressure (MPa)

0.6

0.8

Figure 5-5: Adsorption of Methane-Ethane on BPL Carbon at 212.7 K
and a Bulk Ethane Mole Fraction of 0.501 (Data of Reich, et al., 1980).

88

 

Adsorption (mmoilg)

 

Ethane

Methane
e

’../"’"‘"f

 

I I I I

0.2 0.4 0.6 0.8
Pressure (MPa)

Figure 5-6: Adsorption of Methane-Ethane on BPL Carbon at 212.7 K
and a Bulk Ethane Mole Fraction of 0.255 (Data of Reich, et al., 1980).

89

 

Adsorption (mmoilg carbon)

 

 

 

 

5
4.5 -
4 ..
Ethane
3.5 -
3 .1 I
2.5 -
2 _
1.5 -
0
1 -
O
0.5 - ./
Methane
0 1 I I I
0 0.2 0.4 0.6 0.8 1
Pressure (MPa)

Figure 5-7: Adsorption of Methane-Ethane on BPL Carbon at 260.2 K
and a Bulk Ethane Mole Fraction of 0.255 (Data of Reich, et al., 1980).

Adsorption (mmol/g)

 

 

o 9—"

ethylene

methane

_L

 

 

0 0.2

l I

0.4 0.6

Pressure (MPa)

Figure 5-8: Methane-Ethylene Adsorption on BPL Carbon at 212.7 K and Initial Bul

Ethylene Concentration of 0.74. Data of Reich, et al., 1980.

91

0.8

 

Ethylene Fraction in Adsorbed Phase

 

0.9 1-

0.8 ..

0.7 -~

0.6 -~

0.5 --

0.4 ~-

0.3 ~-

0.2 --

0.1 --

 

 

 

 

0 0.2 0.4 0.6 0.8

Pressure (M Pa)

Figure 5-9: Adsorption of Ethylene-Methane Mixture on BPL
carbon at 212.7 K and Initial Bulk Ethylene Concentration of
0.74 (Data of Reich, et al., 1980).

92

1.4

 

1.2 «

P
on

Toluene Adsorption (mmoilg)
O
'01

 

0.4 - /

 

 

 

318 K
0'2 ' 308 K
328 K
0 T . . . . . .
7 a 9 10 11 12 13 14 15

Pressure (MPa)

Figure 5—10: Pressure (density) effects of toluene adsorption on Degussa WSIV.
Cm = 1 mmol/L with data points at

pbulk = 6.3 moi/L, 8.4 moi/L, and 10.8 moi/L for each isotherm. Data of Tan and
Liou, 1990.

93

REFERENCES

Chen, 1.; Wong, D.; Tan, C.; Subramanian, R.; Lira, C.T.; Orth, M.; Industrial and
Engineering Chemistry Research, 1997, 36, 2808.

Elliott, J.R.; Lira, C.T. Introductory Chemical Engineering Thermodynamics; Prentice
Hall, 1998.

Miller, RC; Kidney, A.J.; Hiza, M.J.; JChethennodynamics, 1977, 9, 167.

Ng, Heng-Joo; Robinson, Donald B.; Journal of Chemical and Engineering Data, 1978,
23 , 325.

Reich, R; Ziegler, W.T.; Rogers, K.A. Ind.Eng.Chem.Process.Des.Dev., 1980, 19, 336.

Reid, R.C.; Prausnitz, J.M.; Poling, BE The Properties of Gases and Liquids, Fourth
Edition; McGraw-Hill: New York, 1987, 734.

Subramanian, Ramkumar; Doctoral Dissertation, Michigan State University, 1995.

Tan, 0; Liou, D.; IndEng.Chem.Res., 1990, 29, 1412.

94

APPENDICES .

95

 

Appendix A:
OPTIMIZIN G THE ADJUSTABLE PARAMETERS

The SLITS program, as discussed in the master’s thesis of Cassandra A. Smith,
has two adjustable parameters, 8er (fluid-solid interaction potential) and H (slit width).
This thesis proposes a routine for optimally fitting these adjustable parameters to a data
set. Before the advent of this OPTIMIZATION program, the adjustable parameters were
fitted manually by trial and error. All of the fitted graphs in chapter 3 were generated
from the OPTIMIZATION program in conjunction with the SLITS program. This

chapter will discuss how the program was conceived and how it works.
INPUT

In the user interface, the program asks for several pieces of information. It needs
the ID number corresponding to the pure compound being analyzed. It also needs the
surface area of the activated carbon on which the adsorption occurs. Approximate values
for the two adjustable parameters also need to be supplied. An approximate value for
8ka can be estimated from the correlation shown in chapter 3. H is also not difficult to
guess; based on all of the fits done, H usually falls between 12 and 16 angstroms. Also,
H is a property of the carbon surface; it has no dependence (according to the SLD model)
on the type of component being adsorbed. The program usually has no trouble
optimizing these initial guesses. The user is also asked to select which of the two
parameters to optimize first. This reasoning will be discussed in the next section.

A data file must also be supplied to the program. This data consists of the
experimental isotherms. More than one isotherm can be entered and optimized at once.

Basically, at a given temperature, there will be a list of pressures with the corresponding

96

molar adsorption per gram of carbon. The data file is formatted in such a way that the
program will accept any number of data points and any number of isotherms. The only
requirement is informing the program how many points/isotherms it needs to read in.

This is also built into the data file.

RATIONALE / METHOD

The program reads in the pressures from the input file and runs them through a
modified version of the SLITS program in order to generate a calculation for moles
adsorbed. Once this is done for every data point, the modeled adsorption values are
compared to the experimental adsorption values. A mean square error (2(1" army) is
then calculated. This is essentially the factor that determines how good the fit is.

The optimization is done one parameter at a time. The user chooses which one to
optimize first. For example, if H is chosen , 8ka will be held constant while H is being
adjusted to produce the minimum mean square error. First, H is coarsely adjusted by
increments of l angstrom until the best value is found. Then H is more finely adjusted by
increments of 0.1 angstrom until an optimal point is reached. One could choose to
optimize at even smaller orders of magnitude, but those don’t have much of an effect on
the appearance of the fit. After H is successfiilly fitted, this new H value is held constant
while 8ka is being adjusted to minimize the mean square error of the data. Sfi/k is first
adjusted by increments of 10, then by increments of l and 0.1. After this is complete, the
program determines whether the new parameter values are different from the initial
guesses. If so, the procedure repeats itself until the fitted parameter values no longer

change.

97

OUTPUT

Once the optimal parameters are obtained, the user is given the option to plot a
fitted curve to the data. For the fitted curve, the number of moles adsorbed is calculated
for small increments of pressure up to the highest pressure given in each set of isotherm
data. All of the experimental data and fitted data is then exported into an output file

which can be easily imported into Excel and plotted.

98

Appendix B:
PURE ABSORPTION PROGRAMS AND SUBROUTIN ES

Optimization
$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

Program: OPTIMIZATION

!
l
!
!
i This program will optimize the parameters H and epsilon/K (EPSOVERK)
! for adsorption in a slit using the ESD equation of state.

! For use with a single component.

1

1

$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

**#************¥*****************************************##*#*****

Section 1

Variable Definitions and Formats

!
!
!
!
! **********#****************t**************************************

PROGRAM OPTIMIZATION

DOUBLE PRECISION ERROROLD, ERRORNEw, CHANGEI, CHANGEz, x1, x2,
X3

DOUBLE PRECISION ERROR, ERRORCOUNT, SIGFF

DOUBLE PRECISION EXCESS(30, 100), EXCESSCA(30, 100), PRESSURE(30, 100)
DOUBLE PRECISION TEMP(30)

DOUBLE PRECISION ENTRYNUMBER2(10)

DOUBLE PRECISION PRESSUREFIT(100),ADSORP(100)

DOUBLE PRECISION EPSOVERK3(200), H3(200), ERRORARRAY(200)
INTEGER Q, R T, B, Response, ID, L, LMNO

INTEGER ENTRYNUMBERI, COUNT], COUNT2, MEMORY, ARCHIVE

DOUBLE PRECISION H, H2, EPSOVERK, EPSOVERKZ, AREA, P

1 OPEN (UNIT=60,FILE='CHECK.TXT')
1 ERRORCOUNT = o

ERROROLD = 999999

Q=O

P=10

R=1

T=O

99

! P is the optimization interval; Q determines whether E/k or H will

! be optimized; R determines whether the parameters should be changed
! in the positive or negative direction (1 pos., -1 neg); T counts

! how many times R changes during optimization, and thus, determines
! when the parameters are optimized

WRITE(*,*) "Please prepare INDATADAT with your adsorption data,"
WRITE(*,*) "or you can specify your input and output file names"
WRITE(*,*) "within the code’s OPEN statements."

WRITE(*,*)
WRITE(*,*) "Press 0 to quit or 1 to continue."
READ(*,*) Response
IF (Response = 0) then
STOP
ENDIF

WRITE(*,*) "Enter ID number of pure component."
READ(*,*) II)

WRITE(*,*) "Enter the fluid particle diameter in angstroms."
READ(*,*) SIGFF

WRITE(*,*)'Enter surface area of activated carbon.‘
READ(*,*) AREA

5 FORMAT (P97, 5x, F8.4)

*******##****#***#*********************#***************¥**********

Section 2

l

l

!

l A generic data file will be read from; the data sets of interest
! should be pasted into this file from different files
! Read in experimental data
! Each data set is assigned a temperature

! Pressure (MPa) is in the first column; surface excess (mmol/g)
! is in the second column

1

1

2 arrays will be used to collect the pressure and excess data
**********************#********#**t*******#****************#******

OPEN(UNTT =3 l, FILE='INDATA.DAT')

! The number of temperature profiles will be read in.

100

 

READ(31,*) ENTRYNUMBER1

! Each set of adsorption data for the different temperatures
! will now be read in.

D0 COUNT1 = 1,ENTRYNUMBER1
READ(31,*) 1 Blank line

! The temperature of the specified set will be read in to
! 2 decimal places

READ(31,*) TEMP(COUNT1)

! The number of pressure/excess entries will now be read in
READ(31,*) ENTRYNUMBER2(COUNT1)

! The individual pressure/excess entries will now be read in
! Pressure-4 decimal places, in MPa

! Surface Excess-4 decimal places, in mmol/g

READ(31,*)
READ(31,*)

DO COUNTZ = 1,ENTRYNUMBER2(COUNT1)
READ(31,*) PRES SURE(COUNT1,COUNT2), EXCESS(COUNT1,COUNT2)
END DO
READ(31,*)
END DO
! ***#************************t**********************#****#*****#***
! Section 3
l
! User will provide guess values for H and Elk (EPSOVERK)
! A modified version of the Slit program is called to calculate
! surface excess at the pressures specified by the above data
! file-> the calculated surface excess values will be stored
1
1

in a third array

******t***#******************#*#***********#***********#**#*******

101

WRITE(*,"') "TEMP(1)=",TEMP(1)

WRITE(*,*) "PRESSURE(1,1)=",PRESSURE(1,1 )
WRITE(*,*) "PRESSURE(Z,1)=",PRESSURE(2,1)
WRITE(*,*) "ENTRYNUMBER1=",ENTRYNUMBER1
WRITE(*,*) "ENTRYNUMBER2=",ENTRYNUMBER2

WRITE(*,*) "Enter guess values for H and E/k respectively."
READ(*,*) I-I, EPSOVERK

WRITE(*,*) "Which parameter do you wish to optimize first?"
WRITE(*,*) "E/k (0) or H (1)? "

RE1“\13("‘,"‘) Q

IF(Q=1)THEN
P=l

ELSE

P=10

ENDIF

H2 = H
EPSOVERKZ = EPSOVERK
MEMORY = l

100 IF (MEMORY > 190) THEN
MEMORY = 1 ! resets array after most of the spaces are firll
ENDIF

ERRORCOUNT = o
H3(MEMORY) = H
EPSOVERK3(MEMORY) = EPSOVERK

IF (MEMORY > 1) THEN
DO ARCHIVE = 1,(MEMORY - 1)
IF (H = H3(ARCHIVE) .AND. EPSOVERK = EPSOVERK3(ARCHIVE))
THEN
ERRORCOUNT = ERRORARRAY(ARCHIVE)
ENDIF
END DO
ENDIF

IF (ERRORCOUNT = 0) THEN
DO COUNT1 = 1,ENTRYNUMBER1
Do COUNT2 = 1,ENTRYNUMBER2(COUNT1)

x1 = TEMP(COUNT1)
X2 = PRESSURE(COUNT1,COUNT2)

102

CALL MODIFIEDSLIT(X1,X2,H,EPSOVERK,X3,ID,AREA, SIGFF)
EXCESSCA(COUNT1,COUNT2) = X3

ERROR = (EXCESSCA(COUNT1,COUNT2) -

EXCESS(COUNT1,COUNT2))**2

ERRORCOUNT = ERRORCOUNT + ERROR
END DO
END DO

ENDIF
ERRORARRAY(MEMORY) = ERRORCOUNT

MEMORY = MEMORY + 1

***********************#*************#*****#*******#***************

Section 4

The experimental surface excess values will be compared to the
calculated surface excess values, and thus, a value representing
error will be calculated (ERRORCOMP).

LOOP 4A

A loop will be set up to adjust E/k at constant H in order to
minimize ERRORNEW. The guess value of E/k will either be
increased or decreased by increments of 10, 1, and .1 respectively
until error is minimized. The new set of parameters will be
provided to the modified slit program described in section 3.

LOOP 4B

A loop will be set up to adjust H at constant E/k for the
minimization of ERRORNEW. The guess value of H will either be
increased or decreased by increments of 1 and .1 until error is
minimized. The new set of parameters will be provided to the
modified slit program described in section 3 for each case.

LOOP 4C

The two parameters will be adjusted as many times as necessary by
the two loops described above until a minimum error value is
converged upon.

103

 

| *********#****************#****$******************************#****

ERRORNEW = ERRORCOUNT

! B is a marker cataloging the error change as a function of the
! parameter changes

IF (ERRORNEW > ERROROLD) THEN
B=0

ENDIF

IF (ERRORNEW < ERROROLD) THEN
B=1

ENDIF

IF (T > 1) THEN
B=2

ENDIF

150 IF (Q = 0) THEN
! Minimization of EA: at constant H

WRITE(*,*) "Optimizing E/k"
WRITE(*,*) "ERRORNEW=",ERRORNEW
WRITE(*,*) "EPSOVERK=",EPSOVERK
WRITE(*,*) "H=",H
WRITE(*,*) "B=",B
WRITE(*,*) "EXCESSCA(I,1)=",EXCESSCA(1,1)

2001F(B=1)THEN

WRITE(*,*) "Adjusting E/k"
IF (R =—- 1) THEN
EPSOVERK = EPSOVERK + P
ENDIF
IF (R = -1) THEN
EPSOVERK = EPSOVERK - P
ENDIF

ERROROLD = ERRORNEW
GOTO 100

ENDIF

[F (B = 0) THEN
WRITE(*,*) "Changing direction for E/k"

104

R=-R

B= 1

T=T+ 1
GOTO 200

ENDIF

[F (B = 2) THEN
WRITE(*,*) "Changing E/k increment"
P=P/ 10
B=1
T=0
IF (P < 1) THEN
WRITE(*,*) "E/k optimization completed"
P=1
Q=l
GOTO 150
ENDIF
GOTO 200

ENDIF
ENDIF
IF (Q =1)TI-IEN

! Minimization of H at constant E/k
WRITE(*,*) "Optimizing H"
WRITE(*,*) "ERRORNEW=",ERRORNEW
WRITE(*,*) "EPSOVERK=",EPSOVERK
WRITE(*,*) "H=",H
WRITE(*,*) "B=",B
WRITE(*,*) "EXCESSCA( 1 , 1)=",EXCESSCA( l , 1)

6001F(B=1)THEN

WRITE(*,*) "Adjusting H"
IF (R = 1) THEN
H = H + P
ENDHi
IF (R = -1) THEN
H = H - P
ENDIF

ERROROLD = ERRORNEW
GOTO 100

105

ENDIF
IF (B = 0) THEN
WRITE(*,*) "Changing direction for H"
R=-R
B=1
T=T+l
GOTO 600

ENDIF

IF (B = 2) THEN
WRITE(*,*) "Changing H increment"
P=P/10
WRITE(*,*) "P = ",P
B=1
T=0
IF (P < 0.09) THEN
WRITE(*,*) "H optimization complete"
P=10
Q=0
GOTO 1000
ENDIF
GOTO 600

ENDIF

ENDIF

! The difference in the last two values for each parameter will

1 now be tabulated as a fraction of their current values

1000 CHANGEI = ABS(HZ - H)/H

CHANGE2 = ABS(EPSOVERK2 - EPSOVERK)/EPSOVERK

WRITE(*,*) "CHANGE1=",CHANGE1
WRITE(*,*) "CHANGE2=",CHANGE2

IF (CHANGEI < 0.001 .AND. CHANGE2 < 0.001) THEN
GOTO 1510
ENDIF

H2 = H
EPSOVERKZ = EPSOVERK

106

GOTO 200

**********#***************************#*****#*#****#*#***********#
Section 5

l
1
!
! The experimental and fitted data points will now be plotted
! Side-by-side on the same output file.

1

!

***************************************************#**************

1500 FORMAT(A,F6.2)

1510 WRITE(*,*)
WRITE(*,1500) " EPSOVERK = ",EPSOVERK
WRITE(*,1500) " H = ",H
WRITE(*,*)

1600 Response = 1
WRITE(*,*) "Generate complete fitted curve with data? (Y =1 /N=2)"
READ("',"‘) Response

IF (Response = 2) THEN
OPEN(UN1T=32, FILE='OUTDATA.DAT',STATUS='REPLACE')

2000 WRITE(32,*) "Isotherm Output"
WRITE(32,*)
WRITE(32,*) "The first column is pressure in MPa."
WRITE(32,*) "The second column is experimental adsorption (mmol/g)."
WRITE(32,*) "The third column is fitted adsorption (mmol/g)."
WRITE(32,*) "Each section of data is preceded by temperature (K)."
2010 FORMAT(A,F8.2,A)
WRITE(32,2010) "Adsorbent surface area is ",AREA, &
" square meters per gram."
2020 FORMAT(A,F6.2)
WRITE(32,2020) " EPSOVERK = ",EPSOVERK
WRITE(32,2020) " H = ",H
WRITE(32,*)
WRITE(32,*)

2100 FORMAT(lOX,F10.6,6X,F10.4,6X,F10.4)

107

 

DO COUNT1 = 1,ENTRYNUMBER1

WRITE(32,*)
WRITE(32,'(F6.2)') TEMP(COUNT1)
WRITE(32,*)
DO COUNT2 = 1,ENTRYNUMBER2(COUNT1)
x1 = PRESSURE(COUNT1,COUNT2)
x2 = EXCESS(COUNT1,COUNT2)
X3 = EXCESSCA(COUNT1,COUNT2)
WRITE(32,2100) x1,x2,x3
END DO

END DO
CLOSE (32)
ELSE
OPEN(UNIT=33, FILE='OUTDATA2DAT’,STATUS='REPLACE’)

WRITE(*,*)
WRITE(*,*) "Curve is being generated in OUTDATAZDAT"
WRITE(*,*)
WRITE(33,*) "Isotherm Output"
WRITE(33,*)
WRITE(33,*) "The first column is pressure in MPa."
WRITE(33,*) "The second column is experimental adsorption (mmol/g)."
WRITE(33,*) "The third column is fitted adsorption (mmol/g)."
WRITE(33,*) "Each section of data is preceded by temperature (K)."
2210 FORMAT(A,F8.2,A)
WRITE(33,2210) "Adsorbent surface area is ",AREA &
" square meters per gram."
2220 FORMAT(A,F 6.2)
WRITE(33,2220) " EPSOVERK = ",EPSOVERK
WRITE(33,2220) " H = ",H
WRITE(33,*)
WRITE(33,*)

2230 FORMAT( l OX,F10.6,22X,F10.4)
2240 FORMAT(IOX,F10.6,6XF10.4)

DO COUNT1 = 1,ENTRYNUMBER1
WRITE(33,*)

WRITE(33,'(F6.2)') TEMP(COUNT1)
WRITE(33,*)

108

x1 = TEMP(COUNT1)
x2 = PRES SURE(COUNT1,ENTRYNUMBER2(COUNT1 ))

CALL
SLITCURVES(X1,X2,H,EPSOVERK,ID,AREA,PRESSUREFIT,ADSORP,LMNO,SIG
FF)

DO L=1,LMNO ! LWO is the number of pressure increments
! used to plot the curve
WRITE(33,223 0) PRESSUREFITG.),ADSORP(L)
END D0

D0 COUNT2 = 1,ENTRYNUMBER2(COUNT1)
WRITE(33,2240) PRESSURE(COUNT1,COUNT2),
EXCESS(COUNT1,COUNT2)
END DO

ENDDO

CLOSE(33)
ENDIF
CLOSE (31)
Response = 1
WRITE(*,*)
WRITE(*,*) "Fit is complete."
WRITE(*,"‘)
WRITE(*,*) "Do you wish to fit another data set (Y=1/N=2)?"
READ(*,*) Response
IF (Response == 1) THEN
GOTO l
ENDIF
! CLOSE(60)

5000 END

109

Modified SLITS

SUBROUTINE MODIFIEDSLIT(T,BLKP,H,EPSFS,AMTS,ID,AREA, SIGFF)

1**#**#****#*##*******¥******#******#*****###*********#*#*##*¥t****

! TEMP: Temperature in Kelvin

! PRESSURE: Pressure in MPa

! H: Slit width in angstroms

! EPSF S: F luid-solid attractive parameter
1

1

1

#**************************#******************#****#*****¥**#*****

! THIS PROGRAM CALCULATES EXCESS IN A SLIT FOR MIXTURES
USING

! THE ESD EQUATION OF STATE. A SQUARE WELL POTENTIAL HAS
BEEN ADDED.

1

! Z IS FROM CARBON SURFACE. THEREFORE, ZOSFF=O AT THE
SURFACE OF WALL.

! ZLCL IS DISTANCE FROM CARBON CENTER IN FLUID-SOLID
DIAMETERS.

! ZLCL IS EQUIVALENT TO ETA IN PREVIOUSLY WRITTEN PROGRAMS.

IMPLICIT DOUBLE PRECISION(A-H,K,O-Z)
DOUBLE PRECISION, INTENT(1N) :: T,BLKP,H,EPSFS,AREA
INTEGER,INTENT(IN) ::ID
DOUBLE PRECISIONJNTENT(OUT) :: AMT S
CHARACTERH ANSWER
CHARACTER"? KIND
CHARACTER*20 NAME
PARAMETER (NMX=10)
DIMENSION
TC(N1\4X),PC(NMX),ACEN(NMX),ZFEED(NMX),S(NMX),NAME(NMX)
COMMON/ESD/KCSTAR(N1VD(),DH(NMX),C(NMX),Q(NI\D(),VX(N1\D(),
COMMON/ESDNLK(NI\D(),ND(N1\D(),EOKP(NMX)
COMMON/ETA/ETAL,ETAV,ZL,ZV
COMMON/EQN/KIND,APPX
COWON/CONSTK/KIJWMXNMX),INITIAL
COMMON/FEED/ZFEED
COMMON/KVALUES/ S
COMMON/N AME/N AME
COMMON/LOC AL/DELTAZ
COMMON FKEEP
DIMENSION LCLDENS(200),FUGBULK(NMX),PSI(NMX),FUGLCL(NMX)

DIMENSION FUG(NIVD(),FUGCCALC(N1\D(),XA(N1\D(),FKEEP(1\11\D()

110

COMMON/POSITION/ZLCL
DOUBLE PRECISION LCLDENS,EXCE,LOSFF,RHOD,FKEEP
DOUBLE PRECISION FUG,RHO,PB,SIGFF,FUGLCL,ZOSFF,SIGFS,SIGSS
1 COMMON/INTERACTIONS/SIGFF, SIGFS,SIGSS
DIMENSION X(NI\D(),IER(12),Y(NMX),PSIl(N1\D(),PSIZ(N1\D(),RHOD(200)

1 NC = 1
INITIAL=O

1100 CALL GETCRIT(NC,ID,NAME,TC,PC,ACEN)
IF(ID.EQ.0)GO TO 1100
CALL OETESD(NC,ID,EOKP,KCSTARDHC,Q,VX,ND)
DO 15 I=2,NC
DO 10 J=1,I-l
CALL ESTACT(ID,I,J,EOKP,VX,KU)

10 CONTINUE
15 CONTINUE

30 CONTINUE

MAXIT=11 l l 1
! LOCAL PRESSURE AUTOMATICALLY INITIALIZED FROM BULK
PRESSURE
AN SWER=1
INIT=0

1000 CONTINUE

RGAS=8.3 1434
1

1 x=1 Is THE BULK COMPOSITION (PURE COMPONENT)
1
X(1) = 1

! SQUARE WELL POTENTIAL DISABLED

111

12

71

29

WELLP = 0
WELLD = 0

EXCE=O
AMT S=0
ADSORB=0

SIGFF=0.
DO 12 I=1,NC
CALL SIGMAFIND(ID,I,SIGMA,NC)

SIGFF=(SIGMA(I)+SIGFF)
CONTINUE

SIGS S=3 .4
SIGFS=O.5*(SIGFF+SIGSS)
ALPHASP=3 .3 5/SIGFS

LMNO=65 .
DELP=BLKPILMNO
PB=0

IFLAG=0

ZLCL=0

FACTOR=O. 1

J=1
LOSFF=(H-SIGS S)/SIGFF

13 Pressure increment loop deleted

PB=BLKP ! CHANGED TO BULK PRES SURE

BULK=1

LIQ=0

FSQD=0.

CALL FUGI(TC,RGAS,T,PB,X,NC,LIQ,FUGCCALC,ZVB,IER,RHOB,BULK,

XA,FSQ,ALPHA,FSQD,LOSFF,ZOSFF)

Do 29 I=1,NC
FKEEP(I)=FUGCCALC(I)
ZKEEP=ZVB
RHOKEEP=RHOB

LIQ=1

112

203

17

CALL FUGI(TC,RGAS,T,PB,X,NC,LIQ,FUGCCALC,ZVB,IER,RHOB,BULK,
XA,FSQ,ALPHA,FSQD,LOSFF,ZOSFF)

DO 203 1=1,NC
IF(FUGCCALC(I).GT.FKEEP(I))THEN
FUGCCALCO)=FKEEP(I)
ZVB=ZKEEP

RHOB=RHOKEEP

ENDIF

CONTINUE

PLCL=PB

DO 17 1=1,NC
FUGBULK(I)=X(I)*PB *FUGCCALC(I)

LAST MUST BE AN ODD NUMBER. IT IS THE NUMBER OF POINTS

GENERATED.

5

LAST=199.

EXCE=0

AMT S=0
ADSORB=0

DELTAZ=(LOSFF-2.*FACTOR)/(LAST-1)
ZOSFF=FACTOR~DELTAZ

DO 77 ISTEP=1,LAST

ZOSFF=ZOSFF+DELTAZ
ZLCL=(ZOSFF*SIGFF+0.S‘SIGSS)/SIGFS
XI=(LOSFF*SIGFF+SIGSS/2.-ZOSFF*SIGFF)/SIGFS

ITMAX=MAXIT
DO 32 1=1,NC

PSIl(I)=4.0*3.1415926*0382"SIGFS‘SIGFS*EPSFS*(-0.2 &

/ZLCL**10+0.5/ZLCL**4+0.5/(ZLCL+ALPHASP)**4+0.5/(ZLCL+2.0* &

ALPHASP)* *4+0.5/(ZLCL+3 .0*ALPHASP)"'*4+0. 5/(ZLCL+4.0*ALPHASP) &
t t 4)

PSIZ(I)=4.0*3.1415926*0.382*SIGFS* SIGFS*EPSFS*(-0.2 &

IXI"10+O.5/X[**4+0.S/(XI+ALPHASP)"4+0.5/(XI+2.0* &
ALPHASP)**4+0.5/(XI+3.0*ALPHASP)**4+0.5/(XI+4.0*ALPHASP) &
##4)

113

32

50

PSI(I)=PSIl(I)+PSIZ(I)

SQWL=0.

DISTl=(ZLCL-SIGSS/SIGFS/2.)*SIGFS
DIST2=(XI-SIGSS/SIGFS/2.)*SIGFS

SQUARE WELL POTENTIAL FOR A SINGLE COMPONENT ONLY
IF(DIST1 .LE.(WELLD*SIGFF))SQWL=WELLP
IF(DIST2.LE.(WELLD*SIGFF))SQWL=WELLP

FUGLCL(I)=FUGBULK(I)*DEXP((PSI(I)+SQWL)/T)
CONTINUE ! LN(i)=mu

IF(J.EQ.1)THEN
DO 501=1,NC

IF(PSI(I).LT.-2500.)THEN

p=0

GOTO S

ELSE

p=P+1

endif
CONTINUE

iflpeq. 1)then

F ACTOR=ZOSFF

DELTAZ=(LOSFF-2*FACTOR)/(LAST-l)
p=p+1
ENDIF
ENDIF
=2
1 WRITE(60,*) 'ISTEP = ',ISTEP,' in MODIFIEDSLIT'
1 WRITE(60,*) 'PLCL = ',PLCL,' in MODIFIEDSLIT'
1 WRITE(60,*) 'BLKP = ',BLKP,' in MODIFIEDSLIT'
CALL

BUBPL(TC,RGAS,T,NC,Y,IER,PLCL,FUG,RHO,XA,FUGLCL,LOSFF,ZOSFF)

6 continue

11

DO 11 I=1,12

IF (IER(I).NE.0)IFLAG=1
IF (IFLAGEQ.1)WRITE(6,*)'IER',(IER(I),I=1,6)
IF(IER(2).EQ.1)WRITE(*,*)'ERROR ALL VAPOR'

114

IF(IER(3).EQ.1)WRITE(*,*)'ERROR ALL LIQUID'
IF(IER(4).EQ.4)WRITE(*,*)'ERROR IN FUGI-NEG LOG CALCD'
IF(IER(4).EQ.5)WRITE(*,*)'ERROR IN FUGI-FUGACITY OVERFLOWS'
IF(IER(S).EQ.1)WRITE(*,*)'ERROR VAPOR AND LIQUID ROOTS CLOSE'
IF(IER(6).EQ.1)WRITE(*,*)'ERROR VLE ITERATION NO CNVRG',ITMAX
IF(IER(7).EQ.1)WRITE(*,*)'ERROR VLE ITERATION FAILED To IMPROVE
1F(IER(8).EQ. 1)WRITE(*,*)'ERROR 1N TC,PC,OR X,Y'
IF(IER(9).EQ.1)WRITE(*,*)'ERROR P SPECIFIED < 0'
IF(IER(10).EQ.1)WRITE(*,*)'ERROR T SPECIFIED IS UNREASONABLE'
IF(IER(11).EQ.1)WRITE(*,*)'ERROR MORE THAN 10 COMPONENTS

IF(IFLAG.EQ.1)PAUSE

IF(IFLAG.EQ.1)STOP

LCLDENS IN MOL/CM3
44 LCLDENS(ISTEP)=RHO-RHOB
RHOD(ISTEP)=RHO

IF(ISTEP.EQ.(LAST))THEN
DO 99 I=l,(LAST-2),2
POSIT1=I
POSIT2=I+1
POSI'I‘3=I+2
1 EXCESS IN MICROMOLES/M"2
EXCE=SIGFS*DELTAZ* l0**2*(LCLDENS(POSIT1)+4.*LCLDENS(POSIT2)+
&
LCLDENS(POSIT3))/6.+EXCE

ADSORB=SIGFS*DELTAZ*10**2"‘(RHOD(POSIT1)+4.*RHOD(POSIT2)+
RI-IOD(POSIT3))/6.+ADSORB
1 AMT DIVIDED BY 2 IN SLIT INTEGRATIONS
99 CONTINUE
1 AREA IN MA2/G, AMT IN MMOL/G
AMTS = EXCE*AREA/1000.
ENDIF

77 CONTINUE

End of original pressure increment loop

82 CONTINUE

115

! NOTE: FOR REPEAT BP CALCULATIONS IT IS ASSUMED THAT
BOOTSTRAPPING

! IS DESIRED. OTHERWISE, THE USER SHOULD ANSWER 'N' TO

! THE REPEAT QUESTION BELOW AND GO BACK THROUGH THE MAIN
! PROGRAM BEFORE PROCEEDING.

! THIS REPEAT STATEMENT HAS BEEN CHANGED TO START THE
PROGRAM

! OVER THE ABOVE NOTE IS INVALID.

1 INIT=1

END

116

SLITS Curve Generation

SUBROUTDIE
SLITCURVES(T,BLKP,H,EPSFS,ID,AREA,PRESSURE,ADSORP,LMNO,SIGFF)

1 THIS PROGRAM GENERATES A FITTED PRESSURE CURVE FOR THE
OUTPUT
1 FILE

I************¢**********#*****#*¥*#*¥*#***************#**#$***¢****

1 TEMP: Temperature in Kelvin

1 PRESSURE: Pressure in MPa

1 H: Slit width in angstroms

1 EPSFS: Fluid-solid attractive parameter
1

1

1

*************#***************iii************************#*********

1 THIS PROGRAM CALCULATES EXCESS IN A SLIT FOR MIXTURES
USING

1 THE ESD EQUATION OF STATE. A SQUARE WELL POTENTIAL HAS
BEEN ADDED.

1

1 Z IS FROM CARBON SURFACE. THEREFORE, ZOSFF=O AT THE
SURFACE OF WALL.

1 ZLCL IS DISTANCE FROM CARBON CENTER IN FLUID-SOLID
DIAMETERS.

1 ZLCL IS EQUIVALENT TO ETA 1N PREVIOUSLY WRITTEN PROGRAMS.

IMPLICIT DOUBLE PRECISION(A-H,K,O-Z)
DOUBLE PRECISION, INTENTGN) :: T,BLKP,H,EPSFS,AREA
INTEGER,INTENT(IN) ::ID
DOUBLE PRECISION AMT S .
DOUBLE PRECISION, INTENT(OUT) :: PRESSURE(IOO),ADSORP(100)
INTEGER, INTENT(OUT) :: LMNO
CHARACTER*1 ANSWER
CHARACTER*2 KIND
CHARACTER*20 NAME
PARAMETER (NMX=10)
DIMENSION
TCCNMX),PC(NMX),ACEN(NMX),ZFEED(NMX),S(NMX),NAME(NMX)

COMMON/ESD/KCSTARCNWQ,DH(N1\D(),C(NI\D(),Q(1\II\D(),VX(N1\D(),
COWON/ESD/VLKmWQ,ND(1\II\D(),EOKP(N1\/IX)

117

COMMON/ETA/ETAL,ETAV,ZL,ZV
COMMON/EQN/KIND,APPX
COMMON/CONSTK/KIJ(N1\D(,NMX),INITIAL
COMMON/FEED/ZFEED
COMMON/KVALUES/S
COWON/NAME/NAME
COMMON/LOCAL/DELTAZ
COMMON FKEEP
DIMENSION LCLDENS(200),FUGBULK(NMX),PSI(NMX),FUGLCL(NMX)
DIMENSION FUG(N1\D(),FUGCCALC(NI\D(),XA(N1\D(),FKEEP(N1\D()
COMMON/POSITION/ZLCL
DOUBLE PRECISION LCLDENS,EXCE,LOSFF,RHOD,FKEEP
DOUBLE PRECISION FUG,RHO,PB,SIGFF,FUGLCL,ZOSFF,SIGFS,SIGSS
1 COMMON/INTERACTIONS/SIGFF, SIGFS,SIGSS

DIMENSION X(NMX),IER(12),Y(NMX),PSI1(NMX),PSIZ(N1\D(),RHOD(200)
OPEN(UNIT=55,FILE='OU'I'DATADAT')
1 NC = 1
INITIAL=0
1100 CALL GETCRIT(NC,ID,NAME,TC,PC,ACEN)
IF(ID.EQ.0)GO TO 1100
CALL GETESD(NC,ID,EOKP,KCSTARDI-LC,Q,VX,ND)
' DO 15 I=2,NC
DO 10 J=1,I-1
CALL ESTACT(ID,I,J,EOKP,VX,KIJ)

10 CONTINUE
15 CONTINUE

30 CONTINUE

MAXIT=1 11 1 l
1 LOCAL PRESSURE AUTOMATICALLY INITIALIZED FROM BULK
PRESSURE
AN SWER=1
INIT =0

1000 CONTINUE

118

RGAS=8.31434
1

1 X=1 IS THE BULK COMPOSITION (PURE COMPONENT)
1

X(1) = 1

1 SQUARE WELL POTENTIAL DISABLED
WELLP = 0
WELLD = 0

EXCE=0
AMT S=0
ADSORB=0

1 SIGFF=0.
1 DO 12 1=1,NC
1 CALL SIGMAFIND(ID,I,SIGMA,NC)

1 SIGFF=(SIGMA(I)+SIGFF)

12 CONTINUE

1 SIGFF=SIGFF/NC
SIGSS=3.4
SIGFS=O.5*(SIGFF+SIGSS)
ALPHASP=3.35/SIGFS

LMNO=65. 1 IF YOU ADIU ST THE INCREMENTS, REMEMBER
1 TO ADJUST THE PRESSURE AND
ADSORPTION

1 ARRAY DIMENSIONS AS WELL IF
NECESSARY

DELP=BLKPILMNO
PB=0

IFLAG=O

ZLCL=0
FACTOR=O. 1

J=1
LOSFF=(H-SIGSS)/SIGFF

119

DO

71

29

L=1,LMNO 1 BEGINNING OF PRESSURE CURVE LOOP
PB=L*DELP
PRESSURE(L) = PB 1 ARRAY OF PRESSURES

BULK=1
LIQ=0

FSQD=0.

CALL FUGI(TC,RGAS,T,PB,X,NC,LIQ,FUGCCALC,ZVB,IER,RHOB,BULK,

XA,FSQ,ALPHA,FSQD,LOSFF,ZOSFF)

DO 29 1=1,NC
FKEEP(I)=FUGCCALC(I)
ZKEEP=ZVB
RHOKEEP=RHOB

LIQ=1

&

203

17

CALL FUGI(TC,RGAS,T,PB,X,NC,LIQ,FUGCCALC,ZVB,IERRHOB,BULK,
XA,FSQ,ALPHA,FSQD,LOSFF,ZOSFF)

DO 203 1=1,NC
IF(FUOCCALca).GT.FxEEP(I))THEN
FUGCCALC(I)=FKEEP(I)
ZVB=ZKEEP

RHOB=RHOKEEP

ENDIF

CONTINUE

PLCL=PB

DO 17 1=1,NC
FUGBULK(I)=X(I)*PB*FUGCCALC(I)

LAST MUST BE AN ODD NUMBER. IT IS THE NUMBER OF POINTS

OENERA .
LAST=199.

EXCE=O

AMT S=0
ADSORB=0

DELTAZ=(LOSFF-2. *FACTOR)/(LAST-1)
ZO SFF=FACTOR-DELTAZ

120

5

32

50

DO 77 ISTEP=1,LAST
ZOSFF=ZOSFF+DELTAZ
ZLCL=(ZOSFF*SIGFF+0.5*SIGSS)/SIGF S
XI=(LOSFF*SIGFF+SIGSS/2.-ZOSFF*SIGFF)/SIGFS

ITMAX=MAXIT
DO 32 1=1,NC

PSI](I)=4.0*3.1415926*0382*SIGFS*SIGFS*EPSFS*(-0.2 &
/ZLCL**10+0.S/ZLCL**4+0.5/(ZLCL+ALPHASP)**4+0.5/(ZLCL+2.0* &

ALPHASP)**4+0.5/(ZLCL+3.0*ALPHASP)**4+0.S/(ZLCL+4.0*ALPHASP) &

##4)

Ps12(I)=4.0*3.1415926*0382*SIGFS*SIGFS*EPSFS*(-0.2 &
/XI**10+0.5/XI“*4+0.5/(XI+ALPHASP)**4+0.5/(XI+2.0* &
ALPHASP)**4+0.5/(XI+3.0*ALPHASP)**4+0.5/()G+4.0*ALPHASP) &
##4)

PSI(I)=PSII(I)+Ps12(I)

SQWL=0.
DISTl=(ZLCL-SIGSS/SIGFS/2.)*SIGFS
DIST2=(XI-SIGSS/SIGFS/2.)*SIGFS
SQUARE WELL POTENTIAL FOR A SINGLE COMPONENT ONLY
IF(DISTl .LE.(WELLD*SIGFF))SQWL=WELLP
IF(DIST2.LE.(WELLD*SIGFF))SQWL=WELLP

FUGLCL(I)=FUGBULK(I)*DEXP((PSI(I)+SQWL)/T)

CONTINUE 1 LN(f)=mu
IF(J.EQ.1)THEN
DO 501=1,NC

IF(PSI(I).LT.-2500.)TI-IEN
1 WRITE(*,*) 'NO GOOD'
p=0
GOTO 5
ELSE
P=p+1
endif

CONTINUE
if(p.eq. l)then
FACTOR=ZOSFF

DELTAZ=(LOSFF-2*FACTOR)/(LAST-l)

121

p=p+1
ENDIF
ENDIF
J=2
CALL
BUBPL(TC,RGAS,T,NC,Y,IER,PLCL,FUG,RHO,XA,FUGLCL,LOSFF,ZOSFF)

6 continue

DO 11 I=1,12
11 IF (IER(I).NE.0)IFLAG=1

IF (IFLAGEQ.1)WRITE(6,*)'IER',(IER(I),I=1,6)
IF(IER(2).EQ.1)WRITE(*,*)'ERROR ALL VAPOR'
IF(IER(B).EQ.1)WRITE(*,*)'ERROR ALL LIQUID'
IF(IER(4).EQ.4)WRITE(*,*)'ERROR IN FUGI-NEG LOG CALCD'
IF(IER(4).EQ.5)WRITE(*,*)'ERROR IN FUGI-FUGACITY OVERFLOWS'
IF(IER(5).EQ.1)WRITE(*,*)'ERROR VAPOR AND LIQUID ROOTS CLOSE'
IF(IER(6).EQ.1)WRITE(*,*)'ERROR VLE ITERATION NO CNVRG',ITMAX
IF(IER(7).EQ.1)WRITE(*,*)'ERROR VLE ITERATION FAILED TO IMPROVE
IF(IER(8).EQ.1)WRITE(*,*)'ERROR IN TC,PC,OR X,Y'
IF(IER(9).EQ.1)WRITE(*,*)'ERROR P SPECIFIED < 0'
IF(IER(10).EQ.1)WRITE(*,*)'ERROR T SPECIFIED IS UNREASONABLE'
IF(IER(11).EQ.1)WRITE(*,*)'ERROR MORE THAN 10 COMPONENTS

IF(IFLAG.EQ.1)PAUSE

IF(IFLAG.EQ.1)STOP

LCLDENS IN MOL/CM3
44 LCLDENS(ISTEP)=RHO-RHOB
RHOD(ISTEP)=RHO

IF(ISTEP.EQ.(LAST))THEN
DO 99 I=1,(LAST-2),2
POSIT1=I
POSIT2=I+1
POSIT3=I+2
! EXCESS IN MICROMOLES/MAZ
EXCE=SIGFS *DELTAZ* 10* *2*(LCLDENS(POSIT1)+4. *LCLDENS(POSIT2)+
&
LCLDENS(POSIT3))/6.+EXCE

ADSORB=SIGFS"‘DELTAZ*10**2*(RHOD(POSIT1)+4.*RI-IOD(POSIT2)+
&

122

RHOD(POSIT3))/6.+ADSORB
1 AMT DIVIDED BY 2 IN SLIT INTEGRATIONS
99 CONTINUE
1 AREAIN M"2/G, AMT INMMOL/G
AMTS = EXCE*AREA/1000.
ENDIF

77 CONTINUE

ADSORP(L) = AMT S 1 ADSORPTION ARRAY

END DO 1 END OF PRESSURE CURVE LOOP

82 CONTINUE

1 NOTE: FOR REPEAT BP CALCULATIONS IT IS ASSUMED THAT
BOOTSTRAPPING

1 IS DESIRED. OTHERWISE, THE USER SHOULD ANSWER 'N' TO

1 THE REPEAT QUESTION BELOW AND GO BACK THROUGH THE MAIN
1 PROGRAM BEFORE PROCEEDING.

1 THIS REPEAT STATEMENT HAS BEEN CHANGED TO START THE
PROGRAM

1 OVER THE ABOVE NOTE IS INVALID.

1 INIT=1

CLOSE(S 5)

END

123

Arranger

C##*#**#***********#*

SUBROUTINE ARRANGE(R1,R2,R3)
DOUBLE PRECISION R1,R2,R3
C PROGRAM TO PUT 3 NUMBERS IN DESCENDING ORDER

DO 20 J=1,3

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

ENDIF

IF(R3.GT.R2) THEN
TEMP=R2
R2=R3
R3=TEMP

ENDIF

20 CONTINUE

RETURN
END

124

Bubble Point Calculation

SUBROUTINE BUBPL(TC,RGAS,T,NC,Y,IER,
1PLCL,FUG,RHO,XA,FUGLCL,LOSFF,ZOSFF)

REVISION DATE: FEB 93 (FOR ESD COMPATIBILITY)
REVISION DATE: JAN 92 SJS-FOR -VE PRESSURES
REVISION DATE: SEPTEMBER 5, 1985
PROGRAMMED BY: IR. ELLIOTT, JR (JAN. 1983)

PURPOSE: CALCULATE BUBBLE POINT PRESSURE OF LIQUID BASED
ON TEMPERATURE AND LIQUID COMPOSITION.

ARGUMENT S :

INPUT:
TC() VECTOR CRITICAL TEMPERATURES OF THE COMPONENTS
PCO VECTOR CRITICAL PRES SURES OF THE COMPONENTS
ACENO VECTOR ACENTRIC FACTORS OF THE COMPONENTS
IDO VECTOR STANDARD ID NUMBERS OF THE COMPONENTS
RGAS GAS CONSTANT (EG. 8.31434 CC-MPA/(GMOL-K))
T ABSOLUTE TEMPERATURE
X() VECTOR MOLE FRACTIONS IN THE LIQUID PHASE
NC NUMBER OF COMPONENTS
INIT PARAMETER FOR SPECIFICATION OF WHETHER THE INITIAL
SS

‘51

IS PROVIDED BY THE USER OR SHOULD BE CALCULATED.
INIT = 0 INITIAL GUESS FOR P CALCULATED BY PSTART
INIT = 1 INITIAL GUESS FOR P PASSED FROM CALLING ROUTINE

INPUT/OUTPUT:
P OPTIONAL INPUT INITIAL GUESS/OUTPUT CALCULATED
ABSOLUTE PRESSURE

ITMAX INPUT MAXIMUM NUMBER OF ITERATIONS PERMITTED.
THE RECOMMENDED VALUE IS 50.
OUTPUT ITMAX IS SET TO THE NUMBER OF ITERATIONS
ERFORMED

OOOOOOOOOOQOOOOOOOO0000000000000

"d

OUTPUT:
Y0 VECTOR MOLE FRACTIONS IN THE VAPOR PHASE
IERO VECTOR ERROR PARAMETERS .
IER(1)=1 AT LEAST ONE OF IER(2)-IER(11) WAS NOT ZERO
IER(2)=1 LIQUID ROOT PASSED FROM FUGI WAS NOT REAL
ON LAST ITERATION

OOOOOOO

125

C 1ER(3)=1 VAPOR ROOT PASSED FROM FUGI WAS NOT REAL
C ON LAST ITERATION ~
C IER(4)=4,5,6 TERMINAL ERROR RETURNED FROM FUGI
CALCULATION
C THE NUMBER TELLS WHICH COMPONENT OF FUGI'S
C ERROR VECTOR WAS SIGNIFICANT
C =4 NEGATIVE LOG CALCULATED
C =5 LOG OF FUGACITY COEFFICIENT CAUSES OVERFLOW
C =6 ITERATION ON COMPRESSIBILITY FACTOR DID NOT
CONVERGE
C IER(5)=1 CALCULATIONS DETERMINED VAPOR & LIQUID ROOTS
EQUAL
C (TRIVIAL SOLUTION)
C IER(6)=1 FAILED TO CONVERGE IN ITMAX LOOPS
C IER(7)=1 AN ITERATION WAS PERFORMED WITH NO
IMPROVEMENT
IN THE OBJECTIVE FUNCTION
IER(8)=1 AN ELEMENT OF TC, PC, OR x WAS SPECIFIED
INCORRECTLY.
IER(9)=1 THE T SPECIFIED WAS LESS THAN ZERO.
IER(10)=1 AN INITIAL GUESS FOR P WAS SPECIFIED BUT
IT WAS UNACCEPTABLE.
IER(11)=1 THE VALUE FOR NC WAS GREATER THAN 10.
IER(12)=1 THE VALUE OF ITMAX WAS LESS THAN 1.

NOTE: UNITS OF ALL TIE INPUTS SHOULD BE
CONSISTENT WITH UNITS OF RGAS. EXCEPT
FOR THIS, TI-E USER MAY CHOOSE HIS OWN UNITS.

REQD. ROUTINES:
PSTART, FUGI, ESTACT, SRICNR

SUBPROGRAM RESTRICTIONS:
AS WRITTEN, TI-E MAXIMUM NUMBER OF COMPONENTS
THAT CAN BE CONSIDERED IS TEN.

GENERAL COMMENTS:
THIS SUBROUTINE SOLVES THE OBJECTIVE FUNCTION GIVEN
IN TI-E LITERATURE REFERENCE BY A SECANT ITERATION
ON TI-E PRESSURE VARIABLE. THE SUBROUTINE CALLED
FOR FUGACITY CALCULATIONS ("FUGI") CONFORMS TO TIE
SPECIFICATIONS OF TI-E SOAVE EQUATION OF STATE GIVEN
IN PROCEDURE 8D1.1.

OOOOOOOOOOOOOOOOOOOOOOOOOOOOO

METHOD RELIABILITY:

126

TI-E AVERAGE ERRORS QUOTED BELOW ARE EXPECTED WHEN

E
C)

THE CORRELATIONS OF BINARY INTERACTION COEFFICIENT GIVEN
IN CHAP. 8 OF THE API TECHNICAL DATA BOOK.

SYSTEM TYPE AVERAGE PERCENT ERROR IN P
HYDROCARBON-HYDROCARBON 4.3
HYDROCARBON-HYDROGEN SULFIDE 4.8
HYDROCARBON-NITROGEN 14.0

HYDROCARBON-CARBON MONOXIDE 7.6
HYDROCARBON-CARBON DIOXIDE 7.4

REFERENCES:
ANDERSON, T.F.; PRAUSNITZ, J.M.; IND. ENG. CHEM. PROC.
DES. DEV., 199.14, (1930).

PROCEDURE 8D1.l OF TECHNICAL DATA BOOK.

******#**#*#*******t***#***#***#**t*******#t********¢**

OOOOOOOOOOOOOOOOOCO

IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION TC(*),Y(*),IER(*),XA1(10),XA2(10),ALPHA(NC)
DOUBLE PRECISION RHO,RHOV,RHOL,FUGLCL
COMMON/ETA/ETALETAVZLZV
COMMON/POSITION/ZLCL
COMMON/LOCAL/DELTAZ
COMMON FKEEP
C COMMON LOSFF,ZOSFF
DOUBLE PRECISION FUG,CHECKL,PLCL,CHECKV,LOSFF,ZOSFF
DIMENSION FUGCCALC(NC),IERF(6),FUGLCL(10),FUG(10),CHECKL(10),
lFUGL(10),FUGV(10),FUGCCALCL(NC),FUGCCALCV(NC),CHECKV(10),
XA(10)

C
C CHECK INPUTS FOR ERRORS.
C

DO 51=1,12

IER(I)=0

5 CONTINUE
IF(NC.GT.10)IER(1 1)=1
DO 81=8,12

8 IF(IER(11).NE.0)IER(1)=1
IF(IER(1).NE.0)GOTO 95

127

C IF(ISTART.EQ.0)THEN

101 D0 10 1=1,NC

10 FUGCCALC(I)=1

C ISTART=1

C ENDIF

C MAX=ITMAX

C BEGINNING OF MAIN ITERATION LOOP

C3 DO 40 I'I'ER=1,MAX
C ITMAX=ITER
BULK=0
C Y(1)=l.
C BULK IS USED To TELL SUBROUTINE ATTCALC TO USE EITHER THE
BULK .
C ATTRACTIVE TERM OR THE LCL ATTRACTIVE TERM.
C BULK=0 FOR YLCL
C BULK=1 FOR YBULK

C CALCULATE VAPOR COMPOSITION AND ERROR IN OBJECTIVE
FUNCTION.

SUMY=0
DO 20 1=1,NC
Y(I)=FUGLCL(I)/(FUGCCALC(I))/PLCL
SUMY=SUMY + Y(I)
C G=1.D0-SUMY
20 CONTINUE

1F(ITEREQ.1)THEN
FOR FIRST ITERATION, OLD VALUES ARE COMPUTED.

C

C

C POLD=PLCL
C GOLD=G
C PLCL=PLCL* .98D0
C GO TO 26
C ELSE
C FOR HIGHER ITERATIONS, A SECANT STEP IS TAKEN.
C IF(G.EQ.GOLD)THEN

C IER(7)=1

C GO TO 86

C END IF

128

CHNG=G*(PLCL-POLD)/(G-GOLD)
END IF

POLD=PLCL
GOLD=G

PNEw=POLD - CHNG
PLCL=PNEW
IF (PLCL.LT.0.0) THEN
CHNG=CHNG/2.0

CHNG=.3*POLD

WRITE(6,*)'PLCL IS -VE IN BUBPL'
GOTO 33
ENDIF

0000000800 00
u

C TEST FOR CONVERGENCE OF P AND Y. IF PRESSURE IS CONVERGED
C BUT Y IS NOT, SKIP RECALCULATION OF LIQUID FUGACITES AT
C TI-E START OF THE NEXT ITERATION.

C IF(DABS(SUMY-l).LE.1.D-9)TI-IEN
C GOTO 26
C END IF

C NOT CONVERGED. GET NEW VAPOR FUGACITES, CI-ECK FOR ERRORS.

26 DO 301=1,NC
Y(I)=Y(I)/SUMY

3o CONTINUE

.C START OF INSERT

C ALNFB=ALOG(FUGBULK(1))
C DO 222 1=1,NC

C222 ALNFL=DLOG(FUGBULK(I))
C FLCL=EXP(ALNFL)

C Calculate local a
C CALL ACALC(BETA,AB,ALCL,ZETA)

Using local parameters calculate VL and local density DENL.

A Newton-Raphson technique is used to converge the local fluid-
fluid pressure to the state where the local fluid-fluid fugacity

has the desired value. At each iteration in pressure, the local
volume (density) is determined. If zeta<1 the local density will

00000

129

zero because no fluidmolecule will fit in the slit

IF (ZETALE. 1) THEN

0000 0

RHOL(L) = o
GOTO 31
ENDIF
TESTOLD=1E10
LCOUNT=0
FSQD=0.
DO 25 J=1,400
C CALL VFCAL(ALCL,B,T,PLCL,VLCL,FLCALC)
SUMSQRs=o
LIQ=1

1 WRITE(60,*) 'PLCL = ',PLCL,' in BUBPL'
15 CALL FUGI(TC,RGAS,T,PLCL,Y,NC,LIQ,FUGCCALCL,
1 ZL,IERF,RHOL,BULK,XA1,FSQ,ALPHA,FSQD,LOSFF,ZOSFF)
IF(IERF(1).NE.0)THEN
DO 35 DE=4,6
IF(IERF(DE).EQ.1)IER(4)=DE
35 CONTINUE
PAUSE
END IF
DO 103 1=1,NC
103 FUGL(I)=Y(I)*PLCL*(‘FUGCCALCL(I))

SUMSQRSL=0
DO 37 1=1,NC
CHECKL(I)=(FUGLCL(I)-FUGL(I))**2
s7 CONTINUE
DO 88 1=1,NC
33 SUMSQRSL=CHECKL(I)+SUMSQRSL

LIQ=0
16 CALL FUGI(TC,RGAS,T,PLCL,Y,NC,LIQ,FUGCCALCV,
1 ZV,IERF,RHOV,BULK,XA2,FSQ,ALPHA,FSQD,LOSFF,ZOSFF)
IF(IERF(1).NE.0)THEN
DO 65 DE=4,6
IF(IERF(DE).EQ.1)IER(4)=DE
65 CONTINUE
PAUSE
END IF

DO 104 1=1,NC
104 FUGV(I)=Y(I)*PLCL*(FUGCCALCV(I))

130

94

93

SUMSQRSV=0

DO 94 1=1,NC
CI-IECKV(I)=(FUGLCL(I)-FUGV(I))* *2
CONTINUE

DO 93 1=1,NC
SUMSQRSV=CHECKV(I)+SUMSQRSV

DO 105 1=1,NC

FUG(I)=FUGV(I)

RHO=RHOV

FUGCCALC(I)=FUGCCALCV(I)

SUMSQRS=SUMSQRSV

XAUFXAZO)

z=zv

LIQ=0

IF(FUGL(I).LT.FUGV(I))THEN 1CR1TERIA IS NOW THE LOWEST

FUGACITY

91

105

C

FUG(I)=FUGL(I)

RHO=RHOL

FUGCCALC(I)=FUGCCALCL(I)

SUMSQRS=SUMSQRSL

XAUFXAIO)

Z=ZL

LIQ=1

ENDIF

CONTINUE

DO 70 1=1,NC
IF(FUG(I).LT.0) WRITE (13*) 'FLCL= '

, FUG(1),PLCL

FR=FUGLCL(I)/FUG(I)
TEST=ABS(FR-l.)

IF ITERATIONS ARE NOT CONVERGING, SUBDIVIDE TIE

INTERVAL BY 2

C

C
C

UP TO 3 TIMES.
IF(I.EQ.150)THEN
PRINT*,PLCL,RHO,FSQ
PRINT*,Z,LIQ,ALPHA(1),FUG(1),FUGLCL(1)
FSQD=0.3
GOTO 70
ENDIF
IF(J.EQ.250)THEN
FSQD=0.8
GOTO 70

131

ENDIF
11(1.eq.350)THEN
FSQD=1.

GOTO 70
ENDIF

IF (TESTOLDLTTEST .AND. LCOUNT.LT.6.)TIEN
LCOUNT=LCOUNT + 1

2 PLCL=PLCL-DELP
DELP=DELP* .3
GOTO 23
ENDIF
LCOUNT = 0
TESTOLD=TEST
IF (TEST.LT.0.001) GOTO 86
C USE A LIMITED NEWTON-RAPHSON TECHNIQUE BASED
ON RlenF = VdP.
C THE FACTOR OF 1.2 LIMITS THE INITIAL STEP SIZE.

IF (FR.LT.0) WRITE (*,*) 'FR= ', FR
ALNFR=DLOG(FR)
DELP=RGAS*T*ALNFR*RHO/1.2
IF(-DELP.GT.PLCL) DELP=-0.5*PLCL

23 PLCL = PLCL+DELP
C WRITE(*,*)'FUGLCL= ',FUGLCL(1), 'FUG= ’,FUG(1)
SUMY=0
DO 200 K=1,NC
Y(K)=FUGLCL(K)/(FUGCCALC(K))/PLCL
SUMY=SUMY + Y(K)
C G=1.D0-SUMY
200 CONTINUE

IF(ISTART.EQ.1)THEN
FOR FIRST ITERATION, OLD VALUES ARE COMPUTED.

POLD=PLCL
GOLD=G
PLCL=PLCL" .98D0
ISTART=ISTART+1
GOTO 9
ELSE

0 000000 00

FOR HIGHER ITERATIONS, A SECANT STEP IS TAKEN.

132

IF(G.EQ.GOLD)THEN
IER(7)=1
GO TO 86
END IF
CHNG=G*(PLCL-POLD)/(G-GOLD)
END IF

POLD=PLCL
GOLD=G
PNEW=POLD - CHNG
PLCL=PNEW
IF (PLCL.LT.0.0) THEN
CHNG=CHNG/2.0
CHNG=.3 *POLD
WRITE(6,*)'PLCL 1s -VE IN BUBPL'
GOTO 33
ENDIF
DO 300 K=1,NC
Y(K)=Y(K)/SUMY
300 CONTINUE

00000000800 000000
1.»

C DO 901 L=1,NC

C901 FR=FUGLCL(L)/FUG(L)
C TEST=ABS(FR-l.)

C IF(TESTOLD.LT.TEST)THEN

C CHNG=CHNG*.35

C PLCL=PLCL - CHNG

C GOTO 70

C ENDIF

C TESTOLD=TEST

C 1F (TEST.LT.0.001) GOTO 86

70 CONTINUE

25 CONTINUE
WRITE (6,*) 'DID NOT CONVERGE, RHO=',RHO,' PLCL=',PLCL
C30 CONTINUE
END OF INSERT

C 70 CONTINUE
C IF(SUMSQRS.LE.1.D-8)GOTO 86

C ELSE
C GOTO 40

133

C END IF
C IF(IERF(1).NE.0)THEN

C DO 35 DE=4,6

C IF(IERF(DE).EQ.1)IER(4)=DE
C35 CONTINUE

C
C END IF
C GOTO 86

C40 CONTINUE

C END OF MAIN ITERATION LOOP
C .
C ONLY WAY FOR PROGRAM TO REACH TIIIS NEXT STATEMENT IS FOR
C NUMBER OF ITERATIONS TO EXCEED ITMAX

C THEREFORE INDICATE ERROR

IER(6H
C PERFORM FINAL ERROR CHECKS.

86 CONTINUE 1WRITE(*,*)FUG(1),FUGLCL(1),PLCL
IF(IERF(2).EQ.1)IER(2)=1
IF(IERF(3).EQ.1)IER(3)=1
DO 90 DE=2,7
SAVEIT=PLCL
IF(IER(DE).NE.0)IER(1)-—-1
9o CONTINUE
95 RETURN
END

134

r:
2'.
O

 

000000 *000000000000 *0 *0 *0 *000000 *00000 *00000

&71 2 3 4 5 6 72

***************#***********************t***********#*#************

* SUBROUTINE CUBIC *

********##*******¥*#***********itt**t*ttt*******#*****#*#*##***ttt

* TIIIS SUBROUTINE FINDS TIE ROOTS OF A CUBIC EQUATION OF TIE

"' FORM X**3 + A2*X**2 + A1*X + A0 = 0 ANALYTICALLY. *

t*****#*¥****#****t****##*******¢it!*ttlii*t.**#*****¥$t**¢*¥*****

* VARIABLES *

**#*¥******$********iiiit*ttttl*ti**itttit*tt*#*¥*¢t*$**¥##*¢*t**¢

* A0 ---- TIE ZEROETH ORDER TERM OF TIE NORMALIZED CUBIC

* EQUATION *

"' A1 ---—- T'IE FIRST ORDER TERM OF TIE NORMALIZED CUBIC *

* EQUATION

* A2 ----- THE SECOND ORDER TERM OF T:IE NORMALIZED CUBIC *
"‘ EQUATION

* C1 ---- TIE COMPLEX ARGUMENT OF ROOT #1 OF TIE EQUATION

* C2 ---- TIE COMPLEX ARGUMENT OF ROOT #2 OF TIE EQUATION
* C3 ---- TIE COMPLEX ARGUMENT OF ROOT #3 OF THE EQUATION

* CCIECK -- TIE SAME AS "CIECK" BUT CONVERTED TO COMPLEX

* NUMBER FORMAT *

* CIECK —-- Q**3 + R**2, USED TO CHECK FOR TIE CASE OF TIE *
* SOLUTION AND IN FINDING TIE ROOTS OF TIE EQUATION, *
* DOUBLE PRECISION

* DAO ---- "A0" CONVERTED TO DOUBLE PRECSION "‘

"' DAl ---- "A1" CONVERTED TO DOUBLE PRECSION "‘

* DA2 --- "A2" CONVERTED TO DOUBLE PRECSION *

* ES] ---- AN INTERMEDIATE CALCULATION TO USED IN TIE *

* CALCULATION OF "S l "

* E82 ---- AN UNITERMEDIATE CALCULATION TO USED IN TIE "'

* CALCULATION OF "82"

* IFLAG --- A FLAG TO INDICATE TIE CASE OF TIE SOLUTION OF TIE

"' EQUATION: =1 ONE REAL + TWO COMPLEX ROOTS, *

"‘ =2 ALL REAL ROOTS, AT LEAST TWO TIE SAME *

* =3 THREE DISTINCT REAL ROOTS *

* P1 ---- AN INTERMEDIATE SUM USED IN TIE CALCULATION OF *
* "SSI"

"' P2 ---- AN INTERMEDIATE SUM USED IN TIE CALCULATION OF *

135

00000000000000000

* "582" *

 

* Q ----- AN INTERMEDIATE SUM USED IN CALCULATING "CHECK"
"‘ R —----- AN INTERMEDIATE SUM USED IN CALCULATING "CHECK"
* R1 ---- TIE REAL ARGUMENT OF ROOT #1 OF TIE EQUATION

* R2 ---- TIE REAL ARGUMENT OF ROOT #2 OF TIE EQUATION

* M ---- TIE REAL ARGUMENT OF ROOT #3 OF TIE EQUATION

" RECK --- TIE SAME AS "CIIECK", BUT SINGLE PRECISION REAL

* Sl ---- AN INTERMEDIATE VALUE USED TO FIND TIE ROOTS OF
* TIE EQUATION, COMPLEX NUMBER

* S2 -—-- AN INTERMEDIATE VALUE USED TO FIND TIE :ROOTS OF
* TIE EQUATION, COMPLEX NUMBER

* SSl ---- TIE SAME AS 51 BUT DOUBLE PRECSION REAL *

* SSZ ---- TIE SAME AS 82 BUT DOUBLE PRECSION REAL *

"' Zl --—- ROOT #1 OF THE EQUATION, COMPLEX NUMBER *

* 22 ---- ROOT #2 OF TIE EQUATION, COMPLEX NUMBER *

* Z3 ROOT #3 OF TIE EQUATION, COMPLEX NUMBER *

t
t
*

1:

*#***#**#**#******#*¢tit!*ttttttt*tfittfiittlttfifiitttttfitlttttfittiti

SUBROUTINE CUBIC(A2,A1,AO,R1,R2,R3,C1,C2,C3,IFLAG)
DOUBLE PRECISION A2,A1,AO,R1,R2,R3
DOUBLE PRECISION CHECK,DA0,DA1,DA2,P1,P2,Q,RSS1,552
COMPLEX Es 1,Esz,SI,Sz,ZI ,Z2,Z3,CCI-IECK
DAO = DBLE(A0)
DAl = DBLE(A1)
DA2 = DBLE(A2)
Q = DAl/3.DOO - DA2*DA2/9.DOO
R = (DA1*DA2 - 3.DOO"'DAO)/6.DOO - (DA2/3.DOO)"3
CHECK = Q**3 + R*R
IF (CHECK.GT.0.0) THEN
[FLAG = 1
P1 = R + DSQRT(CHECK)
IF(Pl.EQ.0.0)P1=1.D-7
P2 = R - DSQRT(CHECK)
IF(Pl .EQ.0.0)P1=1.D—7
IF (Pl.LT.0.0) THEN
881 = -DEXP((DLOG(-l.DOO*Pl))/3 .1300)
ELSE
881 = DEXP((DLOG(PI))/3.DOO)
ENDIF
IP(P2.EQ.0.0)P2=1.D-7
IF (Pz.LT.o.o) THEN
$32 = -DEXP((DLOG(-l.DOO*P2))/3.DOO)
ELSE
ssz = DEXP((DLOG(P2))/3.DOO)
ENDIF
R1 = 881 + $82 - DA2/3.DOO
R2 = —(SSl + 352) - DA2/3.DOO

136

a:
a:

R3 = R2
C1 = 0.0
C2 = (SQRT(3.))*(SSI - ssz)/2.Doo
C3 = -C2
ELSE IF (CHECK.LT.0.0) THEN
IFLAG = 3
RR = 1.*R
RECK = 1.*CIIECK
CCHECK = CMPLX(RECK,0.0)
E8] = CLOG(RR + CSQRT(CCHECK))/3.
E82 = CLOG(RR - CSQRT(CCHECK))/3.
S] = CEXP(ESl)
$2 = CEXPCESZ)
21 = ($1 + 52) - A2/3
22 = -(Sl + sz)/2 - A2/3 + (CMPLX(o.o,3".5))*(S1 - $2)/2
Z3 = -(81 + sz)/2 - A2/3 - (CMPLX(0.0,3**.5))"‘(Sl - sz)/2
R1 = REAL(Zl)
R2 = REAL(Z2)
R3 = REAL(Z3)
C1 = 0.0
C2=C1
C3 = C1
ELSE
C “IttittttttmflmMINIMUM”!IIIIIItIII*tlttttttttttttttttittittttttttttmflm
C * IF THE ROOTS OF THE EQUATION ARE VERY, VERY SMALL AND
VERY, *
C * VERY CLOSE TOGETHER THIS SUBROUTINE MAY ERRONEOUSLY
REPORT *
C * THAT THE EQUATION HAS ONLY ONE ROOT NEAR ZERO *
C thrills"!!!IttttttttttttttittttttMitttit!*tttttttttttttttttttttitanium
IFLAG = 2
IF (RLT.o.0) THEN
881 = -DEXP((DLOG(-l.DOO*R))/3.DOO)
ELSE IF (R.EQ.o.0) THEN
$81 = 0.0
ELSE
$51 = DEXP((DLOG(R))/3 .DOO)
ENDIF
$82 = 881
R1 = $81 + $32 - DA2/3.Doo
R2 = -(SSl + ssz)/2 - DA2/3.DOO
R3 = R2
C1 = 0.0
C2=C1
C3 = C2
ENDIF

137

RETURN

END
C &71

2

138

72

Fugggity Calculation

C CSFUGIFOR
C LATEST REVISION : 9/94 jre
C : 1/96 (switched to chempot, ADDED POLYETHYLENE jre)
c 7/96 PS, PPO, PEO, PIB (ram natarajan)
SUBROUTINE
FUGI(TC,RGAS,T,P,X,NC,LIQ,FUGC,Z,IER,RHO,BULK,XA,FSQ,
1ALPHA,FSQD,LOSFF,ZOSFF)
IMPLICIT DOUBLE PRECISION(A-H,K,O-Z)
PARAMETERCNMX=10)
DIMENSION TC(*),X(*),FUGC(NC),IER(6)
1
,YQVIJ(NMX,NMX),KVE(NMX),YQVI(NMX),Y(NMX,NMX),EOK(NI\D(,NMX)
1 ,CVI(N1VD(),CVIJ(NI\D(,N1\D(),QV(NMX,N1\D(),XA(NI\D()
COMMON/ESD/KCSTAR,DH,C,Q,VX,VLK,ND,EOKP
COWON/ETA/ETAL,ETAV,ZL,ZV
COMMON/DEPFUN/DUONKTDAONKT,DSONK,DHONKT
COMMON/POSITION/ZLCL
COWON/CONSTK/IUJWNIDQJNITIAL
COMMON FKEEP
C COMMON LOSFF,ZOSFF
DIMENSION

EOKmeo,KCSTAmeo,DH(Nwo,C(Nwo,Q(Nwo,VX(NI/DO

& ,ND(NI\D(),k1(nmx),VLK(N1\D(),ALPHA(NMX),RALPHA(NI\D()

DOUBLE PRECISION RHO,P,Z,ZOLD,B,ETA,LOSFF,ZOSFF

DOUBLE PRECISION FSQ,F,A2,AI,Ao,R1,R2,R3,FSQD

C ND IS THE DEGREE OF POLYMERIZATION OF EACH SPECIES
C EOKP IS THE DISPERSE ATTRACTION OVER BOLTz k FOR PURE SPECIES
C KCSTAR IS THE DIMENSIONLESS BONDING VOLUME FOR PURE
C DH 18 THE REDUCED BONDING ENERGY
C C,Q,Vx ARE THE PURE COMPONENT EOS PARAMETERS
C KU IS THE BINARY INTERACTION COEFFICIENT
C 2 IS PV/NOKT HERE
C IER= 1 -AT LEAST ONE ERROR .
2 - TOO MANY ALPHA ITERATIONS
3 -
4 - ERROR IN ALPHA CALCULATION, SQRT(ALPHA) OR ITERATIONS
5 -
6 - TOO MANY z ITERATIONS

000000

DATA K10,K2,ZM/1.7745,l.0617,9.5/
TRY=ZLCL
IF(INITIALEQ.0)THEN

C WRITE(*,*)'INPUT Id 1'

139

C READ(*,*)k11
DO 1 1=1,NC
Kl(I)=K10
C POLYMER DATA WAS REMOVED FROM THIS AREA!
I CONTINUE
INITIAL=1
END IF
DATA TBASE,TSLOPE/4oo,o/
DO 11 1=1,NC
11 FUGC(I)=1
IER(2)=o
IER(3)=O
YQVM=o.o
VM=o.o
CVM=0
KlYVM=0
DO 30 1=1,NC
IF(X(I).LT.O)PAUSE'ERROR 1N FUGI, x1<o'
DO 40 1=1,NC

EOK(I,J)=DSQRT(EOKP(I)*EOKP(J))*
& ( 1-KU(I,I)*(1+(T-TBASE)/I'BASE*TSLOPE))
CALL ATTCALC(EOK,Y,I,J,T,BULK,LOSFF,ZOSFF)
C Y(I,J) AND EOK(I,J) REPRESENT THE LOCAL TERMS CALCULATED IN
ATTCALC
QWLJ) = (Q(I)"‘VXU) + QU)*VX(1)) / 2-0
YQVU(1J)=QV(1J)*Y(U)
CVUOJ) = (C(I)*VX(J) + C(J)*VX(1))/ 2-0
YQVM=YQVM+YQVU(IJ)*X(I)*X(D
CVM = CVM + CVIJ(I,J)*X(I)*X(J)
4o CONTINUE
c kcstan(l)=0.075
KVE(I)=KCSTAR(I)*VX(I)*( DEXP(DH(I)/T*TC(I))-1 )
VM=VM+X(I)"'VX(I)
K1YVM=K1YVM+X(I)*KI(I)*Y(I,I)*VX(I)

30 CONTINUE
IF(KlYVM.eq.O)write(6,*)'klyvm in fugi=',k1yvm
C INITIATE SECANT ITERATION ON RHO

PORT=PIRGASIT
PVMORT=PORT*VM

C GUESS FOR RHO

140

ETA=PVMORT
IF(LIQ.EQ.I)ETA=.52
RHO=ETANM
IF(LIQ.EQ.0)RHO=PORT
WRITE(6,"‘)'Z IN FUGI :2

00000

C ALPSOL WILL USE TIIIS VALUE OF RHO TO CALCULATE A NEW VALUE
OF RIIO

C CALL ALPSOL(X,NC,KVE,ND,VM,
C & RHO,ZASSOC,XA,RALPH,FASSOC,IER)
C39 CONTINUE

IF(IER(4).EQ.1)GOTO 86
ZREP=4*CVM*RHO/(l-l.9DO*VM*RHO)
ZATT=-ZM*YQVM*RHO/(1+K1YVM*RHO)
Z=(1+ZREP+ZATT+ZAS SOC)

RHOLD=RHO
ERROLD=PORT-RHO*Z
RHO=RHO/1.05

IF (RHO.LT.O) WRITE(6,31)LIQ
NITER=0

C100 NITER=NITER+1

C IF(NITER.GT.2500)GO TO 86

00000 0000

C IER(2)=0

C 1ER(3)=O

C ETA=RHO*VM
C10 CALL

ALPSOL(X,NC,KVE,ND,VM,RHO,ZASSOC,XA,RALPH,FASSOC,IER)
IF(IER(4).EQ.1)GOTO 86
ZREP=4*CVM*RHO/(l-1.9DO*VM*RHO)

ZA'IT=-ZM*YQVM*RHO/(1+K1YVM*RHO)

Z=(1+ZREP+ZATT+ZAS SOC)

ERR=PORT-RHO*Z
CHNG=ERRI(ERR-ERROLD)“(RHO—RHOLD)
RHOLD=RHO

ERROLD=ERR

RHO=RHO-CHNG
IF(DABS(CHNG/RHO).GT.1.D-10)GO TO 100

0 000000 0000

WRITE(*,*)'ETA,Z',ETA,Z

141

C WRITE(*,*)'ZREP,ZATT',ZREP,ZATT

C IF (RHO.LT.0) WRITE(6,31)LIQ

C ITERATION ON RHO HAS CONCLUDED. GET DEPARTURES AND
FUGACITY COEFFS.

IF(FSQD.GT.0.01)THEN
FSQ=FSQD

GOTO 2

ENDIF

FSQ=0.03
IF(LIQ.EQ.1)FSQ=O.8
ZOLD=5.

DO 50 ITER=1,1000
B=VM*P/RGAS/T

A2=KlYVM*P/RGAS/T-l .-1.9"'B+FSQ

A1=-1.9*KIYVM/VM*B**2-K1YVM*BNM+1.9*B-4*CVM*BNM+
IZM*YQVM*P/RGASIT+FSQ*K1YVM*P/(RGAS*T)

AO=1.9*K1YVMNM*B"2-4“KlYVM’CVM/VM/VM‘B’Q-

ZM*1.9*YQVM/VM*B"2

81

CALL CUBIC(A2,A1,AO,R1,R2,R3,C1,C2,C3,IFLAG)

IF(IFLAGNE. 1) THEN
CALL ARRANGE(R1,R2,R3)
ZV=(Rl)

ZL=(R3)

IF(LIQ.EQ.1)THEN
RHO=PIZL/RGAS/T
ETA=VM*RHO
DO 81 1=1,NC
ALPHA(I)=KVE(I)/VM"‘ETA/(l-1.9*ETA)
1F(ALPHA(I).LT.0.)THEN
FSQ=FSQ*.3
GOTO 50
ENDIF
CONTINUE
CALL FCALC(X,NC,ALPHA,FSQ,F,XA,RALPHA,ND)
Z=ZL
ELSE

142

82

83

50
51

RHO=P/ZV/RGAS/T
ETA=VM*RHO
DO 821=1,NC
ALPHA(I)=KVE(I)/VM*ETA/(l-l .9*ETA)
IF(ALPHA(I).LT.0.)THEN
FSQ=FSQ/2.
GOTO 50
ENDIF
CONTINUE
CALL FCALC(X,NC,ALPHA,FSQ,F,XA,RALPHA,ND)
z=zv
ENDIF

ELSE
z=R1
RHO=P/Z/RGAS/T
ETA=VM*RHO
DO 83 1=1,NC
ALPHA(I)=KVE(I)/VM*ETA/(1-I .9*ETA)
IF(ALPHA(I).LT.0.)THEN
FSQ=FSQ/2.
GOTO 50
ENDIF
CONTINUE

CALL FCALC(X,NC,ALPHA,FSQ,F,XA,RALPHA,ND)
ENDIF
COMPARE=DABS(ZOLD/Z-l.)
IF(COMPARE.LE.0.001)GOTO 52
ZOLD=Z
CONTINUE
WRITE(*,*)'WARNING-NO. ITERATIONS EXCEEDED! (FUGD'
WRITE(*,*)'LIQ = ',LIQ
WRITE(*,*)'Z = ',z
WRITE(*,*)'ZV = ',ZV
WRITE(*,*)'ZL = ',ZL
WRITE(*,*)'COMPARE = ',COMPARE
WRITE(*,*)'P = ',P,'in FUGI'
WRITE(*,"‘)'T = ',T
WRITE(*,*)'ITER = ',ITER
IF(ETAGT0.52.ORETA.LT.0.)THEN
PR1NT*,'ETA= ',ETA
ETA=DABS(ETA/10.)
RHO=ETANM
z=.9
DO 19 L=1,NC
XA(L)=.9

143

19

77

80
9O

FUGC(L)=.9
CONTINUE
ENDIF
CONTINUE
WRITE(52,*)'Y=',Y

ZREP=4*CVM*RHO/(1-l.9DO*VM*RHO)
ZATT=—ZM*YQVM*RHO/(1+K1YVM*RHO)
Z=(1+ZREP+ZATT+ZASSOC)

ETA=RHO*VM
DO 77 1=1,NC
YQVI(I)=0.DO
CVI(I)=O.DO
CONTINUE

ASSOC=0
UNUMER=0
UDENOM=1
UATT=0
DO 90 1=1,NC
ASSOC=ASSOC+X(I)*ND(I)*( 2*DLOG(XA(I))+l-XA(I) )

UATT=UATT+X(I)*VX(I)*EOK(I,I)/T*(Y(I,I)+K2)*K1(I)
UFACTI=X(I)*ND(I)*RALPHA(I)*XA(I)*(2-XA(I))
UDENOM=UDENOM+X(I)*ND(I)*(RALPHA(I)*XA(I))**2
EBETHI=DH(I)/T*TC(I)
IF(EBETHI.GT.1.D-5)THEN

EBETHI=(EXP(EBETHI)—l)/EBETHI
ELSE

EBETHI=1
END IF
DO 80 1=1,NC

EBETHJ=DH(J)/T*TC(J)

IF(EBE'I‘IU.GT.1.D—5)THEN

EBETHJ=(EXP(EBETHJ)~1)/EBETHJ
ELSE
EBETHJ=1

END IF

QIJ=(EXP(DH(J)/T*TC(J))/EBETHJ+EXP(DH(I)/’l‘*TC(I))/EBETHI)/2

UNUMER=UNUMER+UFACTI*X(J)*ND(J)*RALPHA(J)*XA(J)*QIJ

YQVI(I)=YQVI(I)+YQVU(U)*X(J)

CVIOFCVKI) + CVUOJYXU)

CONTINUE

CONTINUE

IF (IER(Z).NE.O)WRITE(6,*)'ERROR IN LIQUID PHASE IN ALPSOL'
IF (IER(3).NE.O)WRITE(6,*)'ERROR IN VAPOR PHASE IN ALPSOL'

144

ETA=RHO*VM
IF (LIQ.EQ.1) THEN
ETAL=ETA
ZL=Z
ELSE
ETAV=ETA
zv=z
ENDIF

C QUEST=1-l.9DO*ETA
C IF (QUEST.LT.O)WRITE(6,*)'WARNING! (l-l .9*ETA) IS -VE IN FUGI'

DO 1201=1,NC
FREP=-4.DO/1.9DO*DLOG(1.0-l.9‘ETA)*CVM/VM
FUGREP=FREP"'( 2*CVI(I)/CVM-VX(I)NM ) + ZREP*VX(I)/VM

FATT=-ZM*YQVM/K1YVM*DLOG(1.DO+K1YVM*RHO)
FUGATT=ZATT*K1(I)*Y(I,I)*VX(I)/K1YVM+
+ FATT*( 2*YQVI(I)/YQVM-Kl(1)*Y(l,I)*VX(I)/K1YVM)

FUGASN=-1.9DO*RHO*VX(I)*FSQ/(1.DO-l.9DO*ETA)
FUGBON=2 *ND(I)*DLOG( XA(I) )+FUGASN

C IF (Z.LT.0.0)TIEN

C WRITE(6,*)'WARN1NG! z NEGATIVE IN FUGI!‘
C ENDIF

FUGC(I)=FUGREP+FUGATT+FUGBON-DLOG(Z)
FUGC(I)=DEXP(FUGC(I))

120 CONTINUE

UATT=-ZM*YQVM*RHO/(I+K1YVM*RHO)*UATT/KIYVM
UASSOC=-UNUMER/UDENOM

DUONKT=UATT+UASSOC
DAONKT=FREP+FATT+ASSOC-DLOG(Z)

DSONK =DUONKT-DAONKT

DHONKT=DUONKT+Z-1

RETURN
86 WRITE(6,686)

31 . FORMAT(lX,'LIQ=',1X,Il,2X,',’,'WARNING! RIIO (1) -VE IN FUGI')
33 FORMAT(lX,'LIQ=',1X,Il,2X,',','WARNING! RHO (2) -VE IN FUGI')

145

686 FORMAT( ERROR IN FUGI. ')

C IF(NITER.GE.2500)THEN

C WRITE(*,*)’TOO MANY z ITERATIONS'

C IER(6)=1

C END IF
IF(IER(4).EQ.1) WRITE(*,*)'ERROR IN ALPSOL'
IER(I)=1

RETURN
END

146

* *000000000000000000000000

Appendix C:
BINARY MIXTURE ADSORPTION PROGRAM

PROGRAM BINARYADSORPTION

DISCLAIMER

This code is distributed for educational purposes only.

The code is not to be sold.

There is no waITanty or guarantee of suitablitiy for any
application.

By using the code the user agrees to these terms.

Note: BULK PHASE MUST BE A GAS
Define Psi as a statement function
Use a 10—4 potential where PSI is the negative of the
intermolecular potential.
SIGGSl and 810682 are sigma fluid-wall for components 1 & 2, in Angstroms.
EPSGSl & 2 are the fluid-wall potential in K
EXCI & 2 are the excess in micro mol/m"2
AMTl & 2 are the amount adsorbed in mmng
Y1 & 2 are the compositions of the two components in the adsorbed phase
P1 & 2 are the fiIgacity coefficients for the two components
F 1 & 2 are the fiIgacities
p is the pressure in Mpa
TC, PC and W represent the critical temperature, pressure and omega's
ZETA = (slit width in A)/sigmaGG : entire slit width in terms
of the quantity of "diameters" of fluid particles
This program is to calculate the adsorption in a mixture for
a Slit

REAL K12b,K12,Vl,V2,W1,VV2,LWIDTH
CHARACTER RFILE*8O
INTEGER ILK,IFLAG,L,NROOT

Function statements for the Lennard-J ones potentials

PSIl (potential) for component A resulting from wall 1
PSIl(ETA)=4.0*3.1415926*O.382*SIGGS1*SIGGSI‘EPSGS1*(-O.2
l/ETA"10+0.5/ETA**4+O.5/(ETA+ALPHA)"4+O.5/(ETA+2.0*
I“ALPIIA)"""4+O.5/(ETA+3.0"‘ALPHA)"""4+O.5/(ETA+4.O"'ALPHA)
***4)

P812 (potential) for component B resulting from wall 1
PSIZ(ETA)=4.0*3.1415926*0382*SIGGSZ*SIGGSZ*EPSGSZ*(-0.2
l/ETA" 10+0. 5/ETA**4+O. 5/(ETA+ALPHA)*"‘4+0. 5/(ETA+2.0*
*ALPHA)‘*4+0.5/(ETA+3.0*ALPHA)**4+O.5/(ETA+4.0*ALPHA)
***4)

I47

C

C
C

000 00000000000

PSI3 (potential) for component A resulting from wall 2
PSI3(XI)=4.0*3.1415926*O.382*SIGGS1*SIGGSl*EPSGSl*(-O.2
l/XI"10+0.S/XI"4+O.5/(XI+ALPHA)**4+0.5/(XI+2.0*
*ALPHA)**4+O.S/(XI+3.0*ALPHA)**4+O.5/(XI+4.0*ALPHA)
*it4)

PSI4 (potential) for component B resulting from wall 2
PSI4(XI)=4.0*3.1415926*0382*SIGGSZ*SIGGSZ*EPSGSZ*(-0.2
1/XI"10+O.5/XI**4+O.5/(XI+ALPHA)*"‘4+O.5/(XI+2.0*
*ALPHA)"4+O.5/(XI+3.0*ALPHA)"4+O.S/(XI+4.0*ALPHA)
*tt4)

Total potential for component A
PSIA(ETA)=PSII(ETA)+PSIB(XI)
Total potential for component B
PSIB(ETA)=PSIZ(ETA)+PSI4(XI)

Fluid-fluid center-to-center distance averaged between components
"Note: according to Reid,Prausnitz, and Poling,

Methane SIGFF=3.758, Ethylene SIGFF=4. 163
SIGFF = 4

Solid-solid (carbon) center-tO-center distance
SIGWW = 3.4
SIGGS l=(SIGFF+SIGWW)/2
SIGGS2=SIGGSI

ALPHA=3.35/SIGGSI

Gas Constant
R = 8.314

Component critical properties
Ethylene

TC1=282.4

PC1=5.032

W1=0.085

ethane

TC2=3 05.4

PC2=4.880

W2=0.099

Mixing parameters (from experimental data)
COZ-TOLUENE

K12b=.1058

METHANE-ETHYLENE
K12b=.0209
METHANE-ETHANE

148

No

K12b=.0094

Unit conversion from MPa to Pa
PC2 = lOOOOOO’PC2
PCl = 1000000*PC1

User interface; HWIDTH is slit width in angstroms
PRINT“, 'ENTER T(K), P(MPA),EPSI, EPSZ, HWIDTH'
READ“, T,P,EPSGSl,EPSGSZ,HWIDTH
K12=K12b
LWIDTH=HWIDTH"‘ 1 .O- SIGWW
ZETA=LWIDTH* 1.0/SIGFF
Conversion to Pa
P=1000000‘P
Mole fi'action of component 1
PRINT“, 'ENTER Vl'
READ“, V1

Initial adsorption values set to 0
EXC1=0.0
EXC2=0.0
AMT1=0.0
AMT2=0.0

WRITE(*,*) 'Output will be submitted to test.txt and testZ.txt'
OPEN (UNIT=2,NAME='test.txt',TYPE=UNK1\IOWN')
OPEN(UNIT=5,NAME='test2.txt',TYPE='UNKNOWN')
OPEN(UNIT=12,NAME=‘test3 .txt',TYPE='UNKNOWN')
WRITE(‘,*)
WRITE(*,*) 'Enter any letter to proceed with calculations.’
READ(*,*) RFILE

Universal mole fraction variables: V1 and V2
Bulk mole fraction variables: VIIN and V21N
V2 = l-Vl
VIIN = V1
V21N = V2

Finite difference specifications for composition

DDVl: comp. 1 mole fraction change, DDV2: comp. 2
DDVI = 113-2
DDV2 = lE-2

Reduced temperatures

TRl = T/I‘Cl

TR2 = T/TC2

149

* SET DELTA P FOR PARTIAL DERIVATIVES; finite difference Spec.
C DDP=.01*P

C BULK CALCULATIONS
* CALCULATE COMPRESSIBILITY FACTOR (Z)
* CALCULATE A1,A2, TIEN A

 

 

Tabulated ESD parameters; b (cm"3/mol); EPSlff (K)
c,q dimensionless
See ESDparmS file
Component 1: Ethylene
Component 2: Ethane
Ethylene
EPSZfl‘=210.275
c2=l.305
q2=1.581
b2=15.013
Ethane
EPSfo=220.429
c2=l.3552
q2=1.6765
b2=l6.716
Methane
EPSlfi‘=178.082
c1=1.0382
q1=1.0728
b1=10.863
Toluene
EPSlfi=3327S2
c1=1.971
q1=2.849
b1=36.227
C Carbon Dioxide
EPSZfi’=178.269
c2=1.832
q2=2.585
b2=10.534

000000000000000000000

C F sq is association parameter not currently being applied
F sq=0.0
C Carbon surface area (m"""2/g)
Area=1300.
C Convert bl and b2 from cm"3/mole to m**3/mole
bl=bl/10"6
b2=b2/10**6

150

0000

C

C

Calculation of Y(bulk)
Ylb: bulk Y for component 1
Y2b: bulk Y for component 2
Y12b: intercomponent bulk Y
CALL YCALC(EPSIff,T,Y1b)
CALL YCALC(EPSfo,T,Y2b)
EPSIZ=((EPSlff‘EPSZfi)".5)*(l-K12b)
CALL YCALC(EPSIZ,T,Y12B)
EPS128=((EPSIII‘EPSZfl)".5)*(l-K12)
CALL YCALC(EPS12S,T,Y12BS)

Mixing of bulk parameters b (bb),c (cb), and q (qb)
CALL BCQMIX(V lIN,V21N,b1,b2,cl,c2,q1,q2,Y1b,Y2b,Y12b,bb,cb,
1 qub,Ybb)
WRITE(*,*)"Returned from BCQMIX"
WRITE(*,*)"bb = ",bb

Variables necessary for cubic calculation
capBb=(bb*P)/(R*T)
caquBb=qub*P/R/l‘
capYBb=Ybb*P/R/T

R1, R2, and R3: root values for Z (compressibility)
WRITE(*,*)"Calling SCUBIC"
CALL SCUBIC(capBb,cb,qu,caquBb,capYBb,Rl,R2,R3,IFLAG)
IF(IFLAGEQ. l) Z=R1
IF (IFLAG.EQ.2) THEN
Z=Rl
IF(R2.GT.R1) Z=R2
END IF
IF(IFLAG.EQ.3) THEN
CALL ARRANGE(R1,R2,R3)
Z=Rl
END IF

DENB is bulk molar density
DENB=P/Z/R/T
WRITE(*,*)"DENB = ",DENB
WRITE(*,*)"1.9N = ",DENB"'bb*1.9

Fugacity Coeflicient Calculation (Pl,P2)
Component 1
WRITE(*,*)"Calling FUG for P1 bulk"
CALL FUG(cb,c1,bb,b1,b2,Ylb,Y12b,DENB,VlIN,V21N,ql,q2,qub,Ybb,
1 Z,Pl)

151

C Component 2
WRITE(*,*)"Calling FUG for P2 bulk"
CALL FUG(cb,c2,bb,b2,bl,Y2b,Y12b,DENB,V211\I,VIII\I,q2,ql,qub,Ybb,
1 Z,P2)
C WRITE(*,*) 'A & B star ',AS,BS,'FUG',P1,P2

C Conversion to MPa from Pa
PBULK=PI1000000

WRITE(2,*) 'T =',T,' K ',‘P =', PBULK,’ MPa'

WRITE(2,*) 'BULK DEN =', DENB,‘ gmol/m"3',' Width = ',ZETA
WRITE(2,*)'Epsilon l & 2 =', EPSGS1,EPSG82

WRITE(2,*) 'inital mol. frac. Vl = ',V1,' V2 = ', V2

C Bulk Fugacity Calculation
F1BULK=V11N*P1*P
F2BULK=V21N*P2*P

PETA=P
Pold1=P
Pold2=P
Pold3=P
Pold4=P
WRITE(*,*)”BULK DONE"
C Bulk calculations done

C...

 

 

«Start iterating in ETA
Local slit calculations begin

ZETA: slit width in terms of a quantity of "fluid diameters"
DELETA: increment width in "fluid diameters" by which the
loop steps from the slit center to the edge of the region
able to hold complete particles
BETA: position in the slit relative to the edge of wall 1 in
terms of "fluid diameters";starts at center and works back
ETA: position in slit relative to wall 1 in terms of
"fluid-solid diameters"
XI: same as ETA except relative to wall 2

00000000000

INCREMENT=100
BETA=ZETA/2.+(ZETA/2.-0.48)/INCREMENT
DELETA=(ZETA/2.-0.48)/INCREMENT
L=0
WRITE (2,102)

152

C

WRITE(S, *)'DENB= ',DENB,‘ ','FlBULK=‘,F1BULK
WRITE (5, *) 'Lefi to Right: BETA,ETA,L, DEN,Pl,P2,F1,P,PSIA'
WRITE (5,*)

BEGINNING OF STEP LOOP ACROSS SLIT

5 BETA=BETA-DELETA

000

C

ETA=(BETA*2.*SIGFF+SIGwm/(SIGFF+SIGWW)
XI=(ZETA"‘SIGFF+SIGWW/2.-BETA*SIGFF)/SIGGSl
WRITE(12,*) 'L = ’,L

FACTOR = 1.0

SET DELTA P FOR PARTIAL DERIVATIVES (finite difference spec.)
DDP=.01’PETA

Initial value of Objective function for initial comparison
GOLD=1E5

Counter variables (H is for fugacity matching loop)
K=0
H = 0

Local calculation of attractive terms: Y12 and Y22
Reminder: Ylb and Y2b are bulk parameters
Yle: LOCAL intercomponent Y parameter

CALL YZCALC(BETAYlb,le,ZETA)
CALL YZCALC(BETA,Y2b,Y21,ZETA)
CALL YZCALC(BETA,Y12bs,Y122,ZETA)

Fluid-fluid fugacities for each component IN slit determined
from 10-4 potential

Fl = FlBULK*EXP(PSIA(ETA)/T)

F2 = F2BULK*EXP(PSIB(ETA)/T)

Pressure at a specific point inside the slit
(Redefined from bulk pressure P)

INITIAL PRESSURE GUESS FOR SLIT CENTER
P=F1/P1+F2/P2

BEGINNING OF ITERATIONS ON V1, V2, AND P TO MATCH LOCAL

FUGACITY

C

C7

EXPRESSIONS; 200 STEPS MAX.

IF ((Poldl-Pole) .GT. 3E5) THEN

153

DPDL1=(POld1 - Pold2)
DPDL2=(((POldl-Pold2)-(Pold2-Pold3))/(Pold2-Pold3))
DPDL3=DPDL2-(((Pold2-Pold3)-(Pold3-Pold4))/(Pold3-Pold4))
P = Poldl + FACTOR’DPDLI

ENDIF

IF (L .GT. 50) TIEN
WRITE(2,*) I GUESS = ',P

ENDIF

0 00000000

WRITE(*,*) 'F1,F2,P',F1,F2,P

0

'bz' and 'cz' are designated local variables for b,c,q parameters
10 CALL BCQMIX(V1,V2,bl,b2,cl,c2,ql,q2,le,Y22,Y12z,bz,cz,quz,sz)

IF (L .GT. 85) THEN
WRITE(12,‘) H = ',H
ENDIF

C Variables necessary for calculation of roots
casz=(bz*P)/(R*T)
caquBz=quz*P/R/T
capsz=sz*P/R/T

C Calculate roots and type of case
CALL SCUBIC(casz,cz,qu,caquBz,capYBz,Rl,R2,R3,IFLAG)

C Single root case
IF(IFLAG.EQ.1) TIEN
Z1=Rl
C DENg is local density
DENg=P/Zl/R/T
C P1 and P2 are local fugacity coefficients redefined fi'om bulk
CALL FUG(cz,cl,bz,b1,b2,le,Y122,DENg,Vl,V2,ql,q2,quz,sz,
1 Zl,Pl)
CALL FUG(cz,c2,bz,b2,bl,Y22,Y12z,DENg,V2,Vl,q2,q1,quz,sz,
l Zl,P2)
"' CALCULATE OBJECTIVE FUNCTIONS (to match fugacity expressions)
G] = V1*Pl*P/F1-l
G2 = V2*P2*P/FZ-l
GB = 1 - V1 - V2
G0=SQRT(G1*G1+ G2"'G2 + G3 *G3)
NROOT=1
GOTO 25
ENDIF

154

C

C

Double root case:
IF (IFLAGEQZ) THEN
Zl=R1
Z2=R2
IF(R2.GT.R1) THEN
Zl=R2
ZZ=R1
ENDIF
ENDIF

Triple root case: take root yielding smallest initial Objective
function
IF(IFLAG.EQ.3) TIEN
CALL ARRANGE(R1,R2,R3)
Zl=Rl
Z2=R3
ENDIF

MULTIPLE ROOT FUGACITY ANALYSIS
Fugacity calculation for root 1 (Z1)-gas; DENrgas density
DENg=P/Zl/R/T
PlA AND P2A ARE LOCAL FUG. COEFF. FOR GAS ROOT
CALL FUG(cz,cl,bz,bl,b2,Ylz,Y122,DENg,Vl,V2,ql,q2,quz,Yb7,
1 Zl,PlA)
CALL FUG(cz,c2,bz,b2,b1,Y21,Y122,DENg,V2,Vl,q2,q1,quz,sz,
1 Zl,P2A)

CALCULATE OBJECTIVE FUNCTIONS for gas phase

GlA = V1*PlA*P/Fl-l

G2A = V2‘P2A*P/F2-l

G3A = 1 - V1 - V2

GOA=SQRT(G1A*GlA+ G2A*G2A+ G3A*G3A)
WRITE(*,*)'G1A,G2A',G1A,G2A

UPPER ROOT TAKEN for 2 root case
G0=G0A
G1=GIA
G2=G2A
G3 =G3A
Z=Z1
P1=PlA
P2=P2A
NROOT=1

DEN=DENg

CRITERIA FOR DROPPING LOWER ROOT: WIEN (l-l .9‘bz‘DENl) < 0

155

C MEANING: REPULSIVE TERM BECOMES ATTRACTIVE DUE TO SIGN
CHANGE

C (ALSO, LOWER ROOT FUGACITY DOESN'T NEED TO BE

C CALCULATED FOR 2-ROOT REGION)

DEN1=P/ZZ/R/T

IF((l-l.9*bz*DENl) .LT. 0) THEN
C PRINT*, 'SKIPPING 22'
GOTO 25
ENDIF

C Fugacity calculation for low root (Z2)—liquid; DEN1=liq. density

C PlB AND P2B ARE FUG. COEFF. FOR LIQUID ROOT
CALL FUG(chl,bz,b1,b2,Ylz,Yl2z,DENl,V1,V2,ql,q2,quz,sz,
l 22,P1B)
CALL FUG(cz,c2,bz,b2,b1,Y22,Y122,DEN1,V2,V1,q2,q1,quz,sz,
1 22,P2B)

* CALCULATE OBJECTIVE FUNCTIONS for liquid phase
GIB = V1*P1B*P/Fl-1
G2B = V2*PZB*P/Fz-l
G3B = 1 - VI - V2
GOB=SQRT(G1B*GIB + G2B*G2B + G3B*GBB)

C WE TAKE TIE ROOT THAT YELDS AN INITIAL OBJECTIVE FUNCTION
C CLOSEST TO ZERO AS SHOWN BY NEXT STATEMENT; IT IS ASSUMED
C THAT THIS ROOT YELDS TIE 'CORRECT' FUGACITY REPRESENTATION
IN
C TIE SLIT (APPLES TO 3-ROOT REGION ONLY)
IF(GOA.GT.GOB) TIEN
DEN=DENI
G0=GOB
G1=G1B
G2=G2B
G3=G3B
Z=Z2
P1=P1B
P2=P2B
NROOT=2
ENDIF
25 CONTINUE
PRINT",'GO=',GOA,GOB,'NROOT=',NROOT,'IFG=',IFLAG

C Local Fugacity calc. wrt local composition and local pressure

156

F1CALC=V1*P1*P
F2CALC=V2*P2*P
VT=V1+V2

WRITE(*,"') FlCALC,Fl,P
WRITE(*,*) F2CALC,F2,H
PRINT“, 'ETA,Y1',ETA,V2

IF (IF K .LT. 1 .AND. GOLD.LT.G0) TIEN
This block circumvents Newton-Raphson by using fractional
portions of the changes calculated from the previous
Newton-Raphson iteration.
K=K+l
Vl=V1-DV1
V2=V2-DV2
P=P-DP
DVl = DV1/2
DV2 = DV2/2
DP = DP/2
GOTO 100
END IF
K=0
GOLD=G0
9999 FORMAT(4G15.5,1X,13)

0000000000000 000

IF (L .GT. 85) THEN
WRITE(12,*)'IFLAG = ',IFLAG,‘ NROOT = ',NROOT,’ z = ',z
WRITE(12,*)' G0 = ',G0,’ F1 = ',Fl,‘ F2 = ',F2
WRITE(12,*)'P = ',P,’ VI = ',Vl,' V2 = ',V2
WRITE(12,*)'G1 = ',Gl,’ G2 = ',GZ,’ G3 = ',G3
WRITE(12,*)
ENDIF

IF (P .GT. 1E8 .OR (GO-GOLD) .GT. 1) THEN
FACTOR=FACTOR+01
IF (FACTOR .GT. 5) THEN
STOP
ENDIF
GOTO 7
ENDIF

0000000

C WIEN CONVERGENCE CRITERIA ARE MET, THIS BLOCK

PROGRESSIVELY

C INTEGRATES ADSORPTION OVER THE WIDTH OF TIE SLIT
IF(G0.LT. lE-3) TIEN

C WRITE (*,*) 'Bubble Pressure Calculation'

157

C WRITE(*,101)
C WRITE (*,103) X1,X2,T

DEN=P/Z/8.3 14/1‘
PETA=P
C Component 1 local adsorption
EXC1=SIGGSI *(Vl *DEN-V11N*DENB)*DELETA* l .E-4+EXC1
AMT1=EXC1*Area/1000
C AMT 1=Vl *DEN” SIGGS] *DELETA*Area* l .E-7+AMT1
C Component 2 local adsorption
EXC2=SIGGS 1 *(V2 *DEN -V21N*DENB)*DELETA* 1 .E—4+EXC2
AMT2=EXC2*Area/1000
C AMT2=V2*DEN*SIGGSI*DELETA*Area*l.E-7+AMT2

WRITE (*,104) V1,DEN,P,BETA,H,L
WRITE (2,104) V1,DEN,P,BETA,H,L,AMT1,AMT2
WRITE (5,110) BETA,ETA,L,DEN,PI,P2,F1,P,PSIA(ETA)

101 FORMAT(IOX,'X1',10X,'X2',10X,'T-K')
102 FORMAT(4X,'Y1',7X,'DEN',9X,'PRESS',9X,'ETA',8X,'H, L')
104 FORMAT(Fs.5,2x,G12.5,2x,G12.5,2X,F5.3,2x14,2x,14)
110 FORMAT(Fs.1,1x,F5.1,1x,I3,1xF7.3,

1 Ix,F5.4,1x,F5.4,1X,G12.5,Ix,G12.5,1x,G12.5)

GOTO 200

ENDIF

C Incremental changes used in Newton-Raphson derivatives
W1 = V1 + DDVl
W2 = V2 + DDV2
PP = P + DDP

"' CALCULATE DERIVATIVES WRT P
C DENd: density in derivative variation
C PNl ,PN2: fug. coefi‘. in derivative variation

casz=(bz*PP)/(R*T)
caquBz=quz*PP/R/T
capYBz=sz*PP/R/T
CALL SCUBIC(casz,cz,qu,caquBz,capYBz,Rl,R2,R3,IFLAG)

C Single Root
IF(IFLAGEQ. l) TIEN
Z=R1
DENd=PP/Z/R/T
CALL FUG(cz,cl,bz,b1,b2,le,Y12z,DENd,Vl,V2,ql,q2,quz,sz,
1 Z,PNI)
CALL FUG(cz,c2,bz,b2,b1,Y22,Y122,DENd,V2,V1,q2,ql,quz,sz,

158

1 Z,PN2)
GOTO 30
ENDIF

C Double Root: largest root always taken
IF(IFLAG.EQ.2) THEN
BIG=R1
SMALL=R2
IF(Rl .LT.R2) TIEN
BIG=R2
SMALL=R1
ENDIF
R1=BIG
R3=SMALL
ENDIF

C Triple Root
IF(IFLAG.EQ.3) TIEN
CALL ARRANGE(R1,R2,R3)
ENDIF

IF(NROOT.EQ. 1) TIEN
C Ifupper (gaseous) root was taken in original loc. calculation,
C it will also be taken in this derivative calculation.
Z=R1
DENd=PP/Z/R/T
CALL FUG(cz,cl,bz,b1,b2,le,Y122,DENd,V1,V2,ql,q2,quz,sz,
1 Z,PNl)
CALL FUG(cz,c2,bz,b2,bl,Y22,Y122,DENd,V2,Vl,q2,ql,quz,sz,
1 Z,PN2)
ELSE
C Iflower (liquid) root was taken in original calculation,
C it will also be taken in this derivative calculation.
Z:
DENd=PP/Z/R/T
CALL FUG(cz,c1,bz,bl,b2,le,Y122.,DENd,Vl,V2,ql,q2,quz,sz,
1 Z,PNI)
CALL FUG(cz,c2,bz,b2,bl,Y22,Y122,DENd,V2,Vl,q2,ql,quz,sz,
l Z,PN2)
ENDIF

30 DIDP = (V1*PN1*PP-V1*P1*P)/DDP/FI
D2DP = (V2*PN2*PP-V2*P2"'P)/DDP/F2
D3DP = 0

* CALCULATE DERIVATIVES WRT V1

159

C
C

CALL BCQMIX(WI,V2,bl,b2,cl,c2,q1,q2,Ylz,Y22,Y122,bz,cz,quz,sz)

casz=(bz*P)/(R*T)
caquBz=quz*P/R/T
capsz=sz"'P/Rfl‘
CALL SCUBIC(casz,cz,qu,caquBz,capYBz,Rl,R2,R3,IFLAG)

Single root case
IF(IFLAG.EQ. l) TIIEN
Z=Rl
DENd=P/Z/R/T
CALL FUG(cz,cl,bz,b1,b2,le,Yl22,DENd,Wl,V2,ql,q2,quz,sz,
1 Z,PNl)
CALL FUG(cz,c2,bz,b2,bl,Y2z,Y12z,DENd,V2,Wl,q2,ql,quz,sz,
1 Z,PN2)
GOTO 40
ENDIF

Double root case:
IF(IFLAG.EQZ) TIEN
BIG=R1
SMALL=R2
IF (Rl .LT .R2) TIEN
BIG=R2
SMALL=R1
ENDIF
R1=BIG
R3=SMALL
ENDIF

Triple root case
IF(IFLAG.EQ.3) TIEN
CALL ARRANGE(RI,R2,R3)
ENDIF

IF(NROOT.EQ. 1) TIEN
Ifupper root was taken in original calculation, it will also
be taken in this derivative calculation.
Z=Rl
DENd=P/Z/R/T
CALL FUG(cz,c1,bz,bl,b2,Ylz,Y122,DENd,VV1,V2,ql,q2,quz,sz,
l Z,PNl)
CALL FUG(cz,c2,bz,b2,bl ,Y2z,Y12z,DENd,V2,W1,q2,ql,quz,sz,
1 Z,PN2)
ELSE

160

C Iflower root was taken in original calculation, it will also
C be taken in this derivative calculation.
Z=R3
DENd=P/Z/R/T
CALL FUG(cz,cl ,bz,b1 ,b2,Y1z,Y122,DENd,VVl ,V2,ql,q2,quz,sz,

l Z,PNl)
CALL FUG(cz,c2,bz,b2,b1,Y22,Yl2z,DENd,V2,VV1,q2,ql,quz,sz,

1 Z,PN2)

ENDIF

40 DlDl = (PNl*VV1*P-P1*V1*P)/DDVl/Fl
D2D1 = (PN2*V2*P-P2*V2*P)/DDV1/FZ
D3D1 = -1

"‘ CALCULATE DERIVATIVES WRT V2
CALL BCQMIX(V1,W2,b1,b2,cl,c2,ql,q2,le,Y2z,Y122,bz,cz,quz,sz)

casz=bz*P/(R*T)
caquBz=quz*P/R/T
capYBz=sz*P/R/I'

CALL SCUBIC(casz,cz,qu,caquBz,capYBz,Rl,R2,R3,IFLAG)

* Single root case

IF(IFLAG.EQ. 1) THEN

Z=Rl

DENd=P/Z/R/T

CALL FUG(cz,cl,bz,bl,b2,le,Y122,DENd,Vl,W2,ql,q2,quz,sz,
l Z,PNl)
CALL FUG(cz,c2,bz,b2,bl,Y22,Y122,DENd,W2,Vl,q2,ql,quz,sz,
1 Z,PN2)

GOTO 50
ENDIF

* Double root case: largest root always taken
IF(IFLAG.EQZ) THEN

BIG=R1

SMALL=

IF(R1.LT.R2) TIEN
BIG=R2
SMALL=R1

ENDIF

R1=BIG

R3=SMALL

ENDIF

161

C Triple root case
IF(IFLAG.EQ.3) THEN
CALL ARRANGE(R1,R2,R3)
ENDIF

IF(NROOT.EQ. l) TIEN
C Largest root is taken as in original calculation.
Z=Rl
DENd=P/Z/R/T
CALL FUG(cz,cl,bz,bl,b2,le,Y122,DENd,Vl,VV2,ql,q2,quz,sz,
l Z,PNl)
CALL FUG(cz,c2,bz,b2,bl,Y22,Y122,DENd,W2,Vl,q2,q1,quz,sz,
l Z,PN2)
ELSE
C Smallest root is taken as in original calculation.
Z=R3
DENd=P/Z/R/T
CALL FUG(cz,cl,bz,bl,b2,Y1z,YlZz,DENd,Vl,W2,ql,q2,quz,sz,
l Z,PNl)
CALL FUG(cz,c2,bz,b2,bl,Y22,Y12z,DENd,W2,Vl,q2,ql,quz,sz,
l Z,PN2)
ENDIF

50 mm = (PNI*V1*P-P1*V1*P)/DDV2/F1
mm = (PN2*W2*P-P2*V2*P)/DDV2/F2
mm = -1

* CALCULATE INCREMENTS
CALL INV(D1DP,D2DP,D3DP,D 1D1 ,D2D1,D3D1,D 1D2,D2D2,D3D2,DET)

DVI = -Gl"'D1D1 - G2*D1D2 - G3*D1DP
DV2 = -G1*D2D1 - G2*D2D2 - G3 *D2DP
DP = -G1"‘D3Dl - G2*D3D2 - G3 *D3DP
PRINT“, ETA,VI,DV1,P,DP
IF(ABS(DV1).GT.VI) THEN
PRINT*,'WARNING, LARGE DY',DV1
DVl=0.5*Vl"'DVl/ABS(DV1)
DV2=-DV1
C P=9.E7
C DP=0
END IF
100 VI = V1 + DVl
V2 = V2 + DV2
P = P + DP

162

IF (L .GT. 91) THEN

WRITE(12,*) 'IFLAG = ',IFLAG,‘ NROOT = ',NROOT,‘ G0 = ',60
WRITE(12,*) 'z =',z,'1=1 = ',F1,' F2 = ',F2

WRITE(12,*) 'P = ',P,‘ Vl = ',Vl,‘ V2 = ',V2

WRITE(12,*)'G1 = ',Gl,‘ G2 = ',GZ,’ G3 = ',G3

WRITE(12,*)

ENDIF

IF(VI.LT.0)V1=0
IF(V1.GT.I)V1=1.
IF(V2.LT.0)V2=0
IF(V2.GT.1)V2=1.

H=H+1

IF (H .LT. 200) GO T010
C END OF ITERATIVE LOOP FOR MATCHING FUGACITY

PRINT“, 'ITERATIONS EXCEEDED'
WRITE (12,*) 'ITERATIONS EXCEEDED'
C DDV1=lE-3
C DDV2=lE-3
C DDP=.01*P
C GOTO 10
200 CONTINUE
L=L+l
C Record previous 4 pressure values
Pold4=Pold3
Pold3=Pold2
Pold2=Pold1
Pold1=P
IF(L.LT.INCREMENT+1) GOTO 5
C END OF STEP LOOP ACROSS SLIT

WRITE(*,*) EXCESSI = ',EXC1,'EXCESSZ = ',EXC2
WRITE(*,*) 'AMTl = ’,AMTl,’ AMT2 = ',AMT2

WRITE(2,*) EXCESSI = ',EXC1,' EXCEssz = ',EXC2

WRITE(2,*)'AMT1 = ',AMT1,' AMT2 = ',AMT2

PRINT", ENTER 1 FOR SAME T, DIFFERENT x1'
PRINT“, ENTER 2 FOR NEW T AND NEW Xl'
PRINT", ENTER 0 To QUIT'

READ*, IDB

IF (IDB.EQ.1) GO To 3

I63

IF (IDB.EQ.2) GO TO 2
END

***********##8##******************#***#********#*tt*******#***¢**¥**t***

* SUBROUTINES BELOW 1!!!

* EPSX: fluid-fluid interaction parameter
SUBROUTINE YCALC(EPSX,T,YIj)
Y1j=exp(EPSX/T )-1.0617
RETURN
END

 

SUBROUTINE BCQMIX(V],V2,b1,b2,c1,c2,ql,q2,Y1,Y2,Y12,b,c,qu,Yb)
REAL Vl,V2
b=Vl"‘b1 + V2*b2
c=V1"‘cl + V2*c2
qY'b=(Vl "2*b1 *ql *Y1+Vl *V2*Y12*(bl *q2+b2*ql)+V2**2*b2"'q2*Y2)
Yb=V1*b1*Yl+V2*b2*Y2
RETURN
END

 

C
C This group of subroutines defines the coefiicients in the
C expression Z**3 + TM2*Z**2 + TM1"'Z + TMO = 0

SUBROUTINE TM2(B,qu,YB,T2)

T2 = l.7745*YB-1-1.9*B+qu

* WRITE (*,*) T2

RETURN

END

SUBROUTINE TM1(B,c,F sq,YB,qYB,T1)
T1 = -l.9*l.7745*B*YB-1.7745*YB+1.9*B-4*c*B+9.S*qYB+
l qu*l.7745*YB
* WRITE (*,*) T1
RETURN
END

SUBROUTINE TMO(B,c,YB,qYB,T0)
T0 = 1.9"“l.7745*B*YB-4*1.7745*c*B*YB-9.5*1.9*B*qYB
* WRITE(*,*) TO
RETURN
END

 

C ESD F ugacity coefficient calculation for component i
SUBROUTINE FUG(c,ci,b,bi,bj,Yi,Yij,rho,xi,xj,qi,qj,qu,Yb,Z,PHI)

164

C

REAL lnPHI,PHI,TERM1,TERM2,TERM3a,TERM3b,TERM3c

REAL TERM4a,TERM4b,TERM5
C WRITE(*,*)"IN FUG, CALCULATING TERM 1"
TERM] = -(4./l.9)*ci*log(1-l.9*b*rho)
C WRITE(*,*)" TERM 1 DONE"
TERM2 = 4*c*bi*rho/(1-(1.9*b*rho))
C WRITE(*,*)" TERM 2 DONE"

TERM3a = -(9.5/1.7745)*log(1+1.7745*Yb*rho)/(Yb)
TERMBb = xj*Yij*(bi*qj+bj*qi)+xi*Yi*2*bi*qi
TERMBc = -qu*Yi"'bi/(Yb)

C WRITE(*,*)" TERM 3 DONE"

TERM4a = -(9.5*qu/(Yb)) '
TERM4b = Yi*bi"‘rho/(l+l .7745*Yb*rho)

 

C WRITE(*,*)" TERM 4 DONE"
TERMS = -log(Z)
C WRITE(*,*)" TERM 5 DONE"
lnPHI =
TERM1+TERM2+TERM3 a*(TERM3b+TERM3 c)+TERM4a*TERM4b+TERMS
C WRITE(*,*)" SUMMATION DONE"
PHI = exp(lnPHI)
C WRITE(*,*)" PHI CALCULATED"
RETURN
END
C
SUBROUTINE SCUBIC(B,c,qu,qYB,YB,R1,R2,R3,IFLAG)
CALL TM2(B,qu,YB,T2)
CALL TMl(B,c,qu,YB,qYB,T1)
CALL TMO(B,c,YB,qYB,T0)

* WRITE(*,*)T2,T1,T0
CALL CUBIC(T2,T1,T0,R1,R2,R3,CI,C2,C3,1FLAG)
WRITE(*,*)IFLAG
WRITE(*,*)'R1 R2 R3'
WRITE(*,"') R1,R2,R3
RETURN
END

***

 

SUBROUTINE ARRANGE(R1,R2,R3)

165

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

IF(R2.GT.R1) THEN
TEMP=R1
R1=R2
=TEMP
ENDIF

IF(R3.GT.R2) THEN
TEMP=R2
R2=R3
R3=TEMP

ENDIF

20 CONTINUE

RETURN
END
C
C 'YB' is bulk Y parm.; 'Y' is local Y parm.
SUBROUTINE YZCALC(BETA,YB,Y,ZETA)

 

SA=ZETA-0. S-BETA

IF (BETA.LE.1.5) THEN
BETAL=BETA
SAL=SA
IF(BETALT.0.5) THEN
BETAL = 0.5
SAL = ZETA - 1.0
ENDIF
Y1=BETAL+05
Y2=I .-1./SAL"3
Y3=1./3."'Y2
Y=YB*(Y1+Y3)*3./8.
END IF

1F(BETA.GT.(ZETA-1.5)) THEN
Y1=SA+1.
Y2=1.-l./(BETA-0.5)**3
Y3=Y2/3.
Y=YB*(Y1+Y3)*3./8.

END IF

166

IF((BETAGT.1.5).AND.(BETA.LE.(ZETA-l.5))) THEN
Y1=7./3.-1./3./(BETA-0.5)**3
Y2=1.-1./SA**3
Y3=Y2/3.
Y=YB*(Y1+Y3)*3./8.
END IF

RETURN
END

 

SUBROUTINE
INV(DlDP,D2DP,D3DP,D1D1,D2D],D3D1,D1D2,D2D2,D3D2,DET)
C011 = D2D2‘D3DP-D3D2‘D2DP
C012 = -(DZD1*D3DP-D2DP*D3D1)
C013 = D2Dl "'D3D2-D3Dl *D2D2

C021 = -(D 1D2 *D3DP-D3D2’D1DP)
C022 = D1Dl‘D3DP-D3D1’D1DP
C023 = -(D1D1 *D3D2-D3D1 I"D1D2)

C031 = D1D2*D2DP-D2D2*D1DP
C032 = -(DlDl*D2DP-D2D1*D1DP)
C033 = D1D1*D2D2-D2D1*D1D2

DET = D1D1*COll + D1D2*C012 + D1DP*C013
WRITE(*,*) 'DET'
WRITE(*,*) DET

DlDl = COI 1/DET

DID2 = C021/DET

DIDP = CO31/DET

D2D1 = C012/DET
D2D2 = C022/DET
D2DP = CO32/DET

D3Dl = COl3/DET
D3D2 = C023/DET
D3DP = CO33/DET
RETURN
END

C &71 2 3 4 5 6 72

C tilttfittl*tttlttfiittit.$1flittttttttttt*tltlfiltfittfififitfifilfifl**#***##

C "' SUBROUTINE CUBIC *

167

00000000000 *000000000000 *0 *0 *0 *000000 *00000 *00

**#****************#**********#*******¢***#**l************t¥******

* TIIIS SUBROUTINE FINDS TIE ROOTS OF A CUBIC EQUATION OF TIE

* FORM X"3 + A2*X**2 + A1*X + A0 = 0 ANALYTICALLY. *

********#*******#***t*******#*tt**#******#****#****#**************

* VARIABLES *

**#*****#**##ttiittit*******titititti*****8*¥************#******¢*

* A0 ----- TIE ZEROETH ORDER TERM OF TIE NORMALIZED CUBIC

* EQUATION *

* A1 ---—- THE FIRST ORDER TERM OF THE NORMALIZED CUBIC *

* EQUATION

* A2 --—— THE SECOND ORDER TERM OF THE NORMALIZED CUBIC *
* EQUATION

* CI —-—- THE COMPLEX ARGUMENT OF ROOT #1 OF THE EQUATION

"‘ C2 —---- TIE COMPLEX ARGUMENT 0F ROOT #2 OF TIE EQUATION
* C3 ----- TIE COMPLEX ARGUMENT OF ROOT #3 OF TIE EQUATION

"‘ CCIECK -- TIE SAME AS "CIECK" BUT CONVERTED TO COMPLEX

* NUMBER FORMAT *

* CIECK -—- Q**3 + R**2, USED TO CIECK FOR TIE CASE OF TIE *
"‘ SOLUTION AND IN FINDING TIE ROOTS OF TIE EQUATION, *
"' DOUBLE PRECISION

* DAO ---- "A0" CONVERTED TO DOUBLE PRECSION *

* DAl ----- "A1" CONVERTED TO DOUBLE PRECSION "‘

* DA2 ---- "A2" CONVERTED T0 DOUBLE PRECSION *

* ESl ---- AN INTERMEDIATE CALCULATION TO USED IN TIE *

* CALCULATION OF "SI"

* E82 ---- AN INTERMEDIATE CALCULATION TO USED IN TIE *

* CALCULATION OF "S2"

* IFLAG --- A FLAG TO INDICATE TIE CASE OF TIE SOLUTION 0F TIE

* EQUATION: =1 ONE REAL + TWO COMPLEX ROOTS, *
* =2 ALL REAL ROOTS, AT LEAST TWO THE SAME *
* =3 THREEDISTINCT REAL ROOTS *

* P1 ----- AN INTERMEDIATE SUM USED IN TIE CALCULATION OF *
j «881"
* P2 ----- AN INTERMEDIATE SUM USED IN TIE CALCULATION OF *
* nsszn

* Q ----- AN INTERNEDIATE SUM USED IN CALCULATING "CIECK"
* R ---- AN INTERMEDIATE SUM USED IN CALCULATING "CIECK"
* R1 ----- TIE REAL ARGUMENT 0F ROOT #1 OF TIE EQUATION "'

* R2 ---- TIE REAL ARGUMENT OF ROOT #2 OF TIE EQUATION *

168

*
t

000000000000

" R3 ----- TIE REAL ARGUMENT OF ROOT #3 OF TIE EQUATION *

"‘ RECK ---- TIE SAME AS "CIECK", BUT SINGLE PRECISION REAL *

* Sl ---- AN INTERMEDIATE VALUE USED TO FIND TIE ROOTS OF
* TIE EQUATION, COMPLEX NUMBER "'

* $2 —---- AN INTERMEDIATE VALUE USED TO FIND TIE ROOTS OF
* TIE EQUATION, COMPLEX NUMBER *

* SS] -—-- TIE SAME AS 81 BUT DOUBLE PRECSION REAL *

"' SSZ ----- TIE SAME AS 82 BUT DOUBLE PRECSION REAL *

" Zl ----- ROOT #1 OF TIE EQUATION, COMPLEX NUMBER *
* 22 ----- ROOT #2 OF TIE EQUATION, COMPLEX NUMBER *
* Z3 ROOT #3 OF TIE EQUATION, COMPLEX NUMBER *

 

*ttt*#******titttittttfi*****#****t*ttfiifittttfi**#***#**********#***

SUBROUTINE CUBIC(A2,A1,A0,R1,R2,R3,C1,C2,C3,IFLAG)
DOUBLE PRECISION CHECK,DAO,DA1,DA2,P1,P2,Q,R,SSI,SSZ
COMPLEX ESl,E82,Sl,SZ,Zl,Z2,23,CCI-IECK
DAo = DBLE(AO)

DA1 = DBLE(A1)
DA2 = DBLE(A2)
Q = DA1/3.DOO - DA2*DA2/9.D00
R = (DA1*DA2 - 3.DOO*DAO)/6.D00 - (DA2/3.D00)**3
CHECK = Q**3 + R*R
IF (CHECK.GT.o.o) THEN
IFLAG = 1
P1 = R + DSQRT(CHECK)
P2 = R - DSQRT(CHECK)
IF (P1.LT.o.o) THEN
$81 = -DEXP((DLOG(-l.D00*P1))/3.DOO)
ELSE
881 = DEXP((DLOG(PI ))/3.DOO)
ENDIF
IF (P2.LT.0.0) THEN
$82 = -DEXP((DLOG(-1.DOO*P2))/3 .DOO)
ELSE
SS2 = DEXP((DLOG(Pz))/3.D00)
ENDIF
R1 = $81 + 532 - DA2/3.DOO
R2 = -(SSl + $32) - DA2/3.Doo
R3 = R2
C1 = 0.0
C2 = (SQRT(3.))*(SSI - ssz)/2.Doo
C3 = -C2
ELSE IF (CI-IECK.LT.0.0) THEN
IFLAG = 3
RR = 1 .*R
RECK = 1."'CHECK
CCHECK = CMPLX(RECK,0.0)

169

t

t

ESI = CLOO(RR + CSQRT(CCHECK))/3.
E82 = CLOG(RR - CSQRT(CCHECK))/3.
S] = CEXP(ESI)
$2 = CEXP(ES2)
21 = (81 + 32) - A2/3
22 = —(Sl + sz)/2 - A2/3 + (CMPLX(0.0,3".5))*(SI - mm
23 = -(Sl + S2)/2 - A2/3 - (CMPLX(0.0,3**.S))*(SI - sz)/2
R1 = REAL(ZI)
R2 = REAL(22)
R3 = REAL(23)
CI = 0.0
C2=C1
C3 = C1
ELSE
C **#****#*3***tfittliItttIl*********$¢*#****¥*I*¢***ititfifittttfitit
C * IF THE ROOTS OF THE EQUATION ARE VERY, VERY SMALL AND
VERY, *
C * VERY CLOSE TOGETHER, THIS SUBROUTINE MAY ERRONEOUSLY
REPORT *
C * THAT THE EQUATION HAS ONLY ONE ROOT NEAR ZERO *
C *1!!!itflit**#*********#*##*#*****tt*Iitttttltttfitfitttfitt***#******¥*
H1AG=2
IF (RLT.o.o) THEN
SSI = -DEXP((DLOG(-1.DOO*R))/3.DOO)
ELSE IF (REQ.o.o) THEN
881 = 0.0
ELSE
SSI = DEXP((DLOG(R))/3 .DOO)
ENDIF
ssz = 881
R1 = SS1 + SS2 - DA2/3.Doo
R2 = —(SSl + ssz)/2 - DA2/3.Doo
R3 = R2
C1 = 0.0
C2=C1
C3 = C2
ENDIF
RETURN
END
C &7 1 2 3 4 5 6 7 2

170

"IMINIMAL“TUTTI?