.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
Sï¬iuéu es O'QC? as Mdgcxr 1mm of]
AC‘kH/OU‘PA Caf‘éah ama Ohfrrs F/avor
ï¬ccovery From C146â€)! P
presentedby
Aaf‘on D So UK l 6
has been accepted towards fulï¬llment
of the requirements for
Ma$1€r o {NSC Pncedegreein Chem C 11 Eng nei’r 03
(Made
V
Major professor
Date ï¬ll/‘7?
0-7639 MS U is an Aï¬maliw 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 fulï¬llment 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 ï¬rst thrust considers the application of
adsorption, as well as some other separations processes, to the recovery of natural cherry
flavoring ï¬'om cherry pits. Some benchmark work has been established relating to the
ï¬ltration, 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 Simpliï¬ed
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 proï¬les 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 ï¬t, the model can only
demonstrate qualitative properties at this point in time.
I dedicate this ï¬nished 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 ï¬nished product from my education. I would also like to thank Dr. Charles
Petty for giving me a second chance to complete the transport phenomena ï¬nal 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 Proï¬le 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 Proï¬le 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 Deï¬ne 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 ï¬rst 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 ï¬nite length
2 - compressibility factor
Greek
or - spacing between carbon planes
8 - interaction parameter
1] - reduced density
p - density
¢ - fugacity coefï¬cient
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 deï¬ned by the equation of state
ff - fluid-fluid property
fs - fluid—solid property
xi
Chapter 1:
INTRODUCTION AND BACKGROUND
Adsorption is deï¬ned 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 signiï¬cance of this work lies in the fact that natural cherry flavoring
(benzaldehyde) is much more valuable than artiï¬cial flavoring. Even though the natural
and artiï¬cial 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 beneï¬cial 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 Simpliï¬ed Local Density model
ï¬rst developed in the dissertation of Bharath Rangarajan. The SLD model was ï¬rst 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 ï¬ts of ethane-methane and
ethylene-methane. However, a ï¬t 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 modiï¬ed to
deal with the different particle sizes present in mixtures. Secondly, it must be modiï¬ed
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 artiï¬cial flavors. While being
identical in taste quality, artiï¬cial 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 puriï¬cation 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 ï¬ltration 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
conï¬rming 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 difï¬cult 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 speciï¬c 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 conï¬rm 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 ï¬rst
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
identiï¬es 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) identiï¬es 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 efï¬cient ï¬ltration a challenging task—especially for large quantities of
hydrolyzate. A couple of stages were necessary in removing particulates.
The ï¬rst 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 ï¬lter (bottom to top): 4 grams of glass wool, 32 grams of diatomaceous
earth (normally used in swimming pool ï¬lters), 130 grams of sand, and 50 grams of
cherry pits for the purpose of spreading liquid flow over the entire diameter of the ï¬lter.
Figure 2-1 shows a schematic of this ï¬lter. A pump pushed the hydrolyzate through the
ï¬lter from a feed tank. While this ï¬lter effectively removed all noticeable particulates
from the hydrolyzate, the ï¬lter 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 ï¬ltration. While the diatomaceous earth increased the
volumetric buildup of the ï¬lter cake, it made the ï¬lter cake sufï¬ciently 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
ï¬lled 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 identiï¬ed 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 signiï¬cant 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 superï¬cial 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 superï¬cial
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, superï¬cial 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 reï¬lled 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, superï¬cial 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, superï¬cial 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, superï¬cial 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—Cï¬n) 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, superï¬cial 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, superï¬cial
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 efï¬ciency 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 ï¬rst step is to ï¬nd 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 ï¬rst 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 signiï¬cant 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 inefï¬ciencies 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 Masï¬gg) % 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 ï¬rst 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 Proï¬le 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 ï¬rst 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 ï¬ns 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 efï¬ciency 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 signiï¬cantly 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 proï¬le 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 proï¬le of the ethyl
acetate layer, and Table 2-19 gives a GC proï¬le 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 Proï¬le 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 ï¬rst 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 proï¬le 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 ï¬rst 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
signiï¬cant 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 modiï¬ed 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 efï¬cient methods for obtaining good
approximations. The best methods incorporate as few adjustable parameters as possible.
The Langmuir, Toth, and Freundlich models are easy to ï¬t, but require temperature-
dependent parameters. The simpliï¬ed 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 simpliï¬ed 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 ï¬t both the
Henry’s law region and the adsorbent capacity. Further, when ï¬tting data, it is difï¬cult
to maintain good ï¬ts 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 ï¬ve 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'ï¬â€˜ + 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 ï¬rst 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)+flï¬(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 ï¬igacity 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 proï¬le 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 modiï¬ed 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 ï¬t 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 ï¬tted: eg/k (fluid-solid interaction potential in Kelvin) and
H (slit width in angstroms). The parameters were ï¬t to optimize the simultaneous
representation of all data in a given ï¬gure rather than optimization of individual
isotherms.
Figure 3-1: Model of a slit-shaped pore showing the variables used to deï¬ne 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 ï¬tted with
the ESD version of the simpliï¬ed 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 ï¬t over 167 K. Other ï¬ts 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 ï¬ts 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 ï¬tting low pressure data.
Better ï¬ts seem to be generated when the pressure ranges are wider. Nitrogen shows
inaccurate predictions of temperature dependence. The ï¬t of propane is weak in that, at
the highest temperature, the knee region is overpredicted. While butane, propylene, and
propane deviate signiï¬cantly from the sphericity assumption of the model, it is not
known why the nitrogen data are not ï¬t 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 speciï¬ed, 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 eï¬lk=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 eï¬/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 ï¬tted 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 ï¬ne-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 ï¬tted 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 ï¬tted 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
ï¬tted. Figure 4-2 shows carbon dioxide adsorption on Columbia Grade L carbon. The
isotherms are more straight because the data set does not go signiï¬cantly beyond the
Henry’s Law region. The optimized parameters give a good ï¬t to this data, but these
68
same parameters might not give a good ï¬t 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 ï¬ts 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 signiï¬cantly 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
eï¬/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 ï¬rgacity 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 ï¬rgacity coefï¬cient. 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
coeï¬'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=ï¬l—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 deï¬ned 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 ï¬t 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 deï¬nitions of the local position, Eta (in fluid-solid diameters).
Also, separate exclusion regions for each component would have to be deï¬ned 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, ï¬rture
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 deï¬ne the 10-4 Lennard-
Jones potentials relative to the ï¬rst 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) = 4npmoiflsï¬j 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 deï¬ne 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 deï¬ned 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 coefï¬cients 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 deï¬ne 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 diï¬â€™erent mixtures have been ï¬tted using
this theory. While the results are signiï¬cantly 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 ï¬t: 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 ï¬ts 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 ï¬t 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 ï¬ts 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 ï¬t 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 ï¬t at 212.7 K on the same BPL carbon.
Once again, methane is under-predicted, but the model ï¬ts 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 ï¬tting these adjustable parameters to a data
set. Before the advent of this OPTIMIZATION program, the adjustable parameters were
ï¬tted manually by trial and error. All of the ï¬tted 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 ï¬ts 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 ï¬rst. This reasoning will be discussed in the next section.
A data ï¬le 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 ï¬le 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 ï¬le.
RATIONALE / METHOD
The program reads in the pressures from the input ï¬le and runs them through a
modiï¬ed 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 ï¬t is.
The optimization is done one parameter at a time. The user chooses which one to
optimize ï¬rst. 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 ï¬nely 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 ï¬t. After H is successï¬illy ï¬tted, this new H value is held constant
while 8ka is being adjusted to minimize the mean square error of the data. Sï¬/k is ï¬rst
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 ï¬tted parameter values no longer
change.
97
OUTPUT
Once the optimal parameters are obtained, the user is given the option to plot a
ï¬tted curve to the data. For the ï¬tted 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 ï¬tted data is then exported into an output ï¬le
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 Deï¬nitions 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 ï¬le 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 ï¬le will be read from; the data sets of interest
! should be pasted into this ï¬le from different ï¬les
! Read in experimental data
! Each data set is assigned a temperature
! Pressure (MPa) is in the ï¬rst 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 proï¬les 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 speciï¬ed 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 modiï¬ed version of the Slit program is called to calculate
! surface excess at the pressures speciï¬ed by the above data
! ï¬le-> 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 ï¬rst?"
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 ï¬rll
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 modiï¬ed 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
modiï¬ed 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 ï¬tted data points will now be plotted
! Side-by-side on the same output ï¬le.
1
!
***************************************************#**************
1500 FORMAT(A,F6.2)
1510 WRITE(*,*)
WRITE(*,1500) " EPSOVERK = ",EPSOVERK
WRITE(*,1500) " H = ",H
WRITE(*,*)
1600 Response = 1
WRITE(*,*) "Generate complete ï¬tted 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 ï¬rst column is pressure in MPa."
WRITE(32,*) "The second column is experimental adsorption (mmol/g)."
WRITE(32,*) "The third column is ï¬tted 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 ï¬rst column is pressure in MPa."
WRITE(33,*) "The second column is experimental adsorption (mmol/g)."
WRITE(33,*) "The third column is ï¬tted 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 ï¬t another data set (Y=1/N=2)?"
READ(*,*) Response
IF (Response == 1) THEN
GOTO l
ENDIF
! CLOSE(60)
5000 END
109
Modiï¬ed 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 ï¬t 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*tï¬ttï¬ittlttï¬ï¬itttttï¬tlttttï¬ttiti
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
Deï¬ne 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 ï¬Igacity coefï¬cients for the two components
F 1 & 2 are the ï¬Igacities
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 ï¬'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; ï¬nite 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 ï¬le
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
EPSlï¬â€˜=178.082
c1=1.0382
q1=1.0728
b1=10.863
Toluene
EPSlï¬=3327S2
c1=1.971
q1=2.849
b1=36.227
C Carbon Dioxide
EPSZï¬â€™=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‘EPSZï¬)".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, *) 'Leï¬ 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 (ï¬nite 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 speciï¬c point inside the slit
(Redeï¬ned 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 coefï¬cients redeï¬ned ï¬'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. coeï¬â€˜. 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 deï¬nes the coeï¬icients 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 coefï¬cient 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 tilttï¬ttl*tttlttï¬ittit.$1flittttttttttt*tltlï¬ltï¬ttï¬ï¬ï¬tï¬ï¬lï¬ï¬‚**#***##
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*#******titttittttï¬*****#****t*ttï¬iï¬ttttï¬**#***#**********#***
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***tï¬ttliItttIl*********$¢*#****Â¥*I*¢***ititï¬ï¬ttttï¬tit
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*Iitttttltttï¬tï¬tttï¬tt***#******Â¥*
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?