THESiS / ICHiGAN STATE UNIVERSITY UBRAR mini mum\\\\\\\\\\g\:;\\\um“win 3 1293 016 026 This is to certify that the thesis entitled Adsorption of Water on Carbon: A Study Using the ESD Equation of State presented by Cassandra A. Smith has been accepted towards fulfillment of the requirements for M.S. degree in Chemicai Engr (with Major professor Date ’Zj'q/97 0-7639 MS U i: an Affirmative Action/Equal Opportunity Institution LIBnAm? Michigan State University PLACE IN RETURN Box to remove this checkout from your record. TO AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE "091514 ADSORPTON OF WATER ON CARBON: A STUDY USING THE ESD EQUATION OF STATE By Cassandra A. Smith A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE Department of Chemical Engineering 1997 ABSTRACT ADSORPTON OF WATER ON CARBON: A STUDY USING THE ESD EQUATION OF STATE By Cassandra A. Smith Several models have been developed to predict the adsorption of water vapor on BPL carbon. The Elliott, Suresh, and Donohue (ESD) equation of state is used along with the Simplified Local Density (SLD) assumptions to predict the thermodynamic properties of adsorbing water on carbon surfaces of varying geometries. These models attempt to predict the trends observed experimentally for the adsorption of water at temperatures ranging from 298 K to 398 K. The inaccuracy of the calculated isotherms dictated the evolution of the model from a basic slit pore to the inclusion of a square well potential, a slit size distribution, and an active site on the surface of the slit pore. The slit model predicts the behavior of adsorbing ethylene very well and each model can predict the threshold pressure of adsorbing water, but none are able to predict the magnitude of the water isotherms that we see experimentally. The difficulties associated with this problem lead us to believe that there is a physical effect of hydrogen bonding that we have not yet identified that is not captured in the models presented here and is not as important when studying the behavior of non-hydrogen bonding fluids. To my family, my strength. And especially Emily, Auntie’s favorite little girl. iii ACKNOWLEDGMENTS I would like to thank my adviser, Dr. Carl T. Lira, for all of his support and patience. I do appreciate every moment of the time you spent helping me. I would also like to thank Caroline Roush and Dr. Christian Lastoskie for their efforts and input for the comparisons made to Monte Carlo simulations. They provided all of the simulation results contained in this thesis as well as contributed to the writing of chapter 4. Thank you Candy, Faith, and Julie, my best source of information. Your help, patience and smiles mean a lot. Thanks to all of my Friends. You’ve helped me manage through the rough times of graduate school and for that you are recognized faithfully in my prayers. iv TABLE OF CONTENTS LIST OF TABLES- - - - -- ...... - - - - _ ..... VI LIST OF FIGURES-.-" - - ......... - ............... . VII NOMENCLATURE - _ ................... - - _ --....IX CHAPTER 1 : INTRODUCTION - - ..1 CHAPTER 2 : ADSORPTION MODEL DEVELOPMENT ...................................... 7 General ............................................................................................................................................ 7 Slits ................................................................................................................................................ 13 Square Well .................................................................................................................................... l4 Slit Distribution 22 Active Site ...................................................................................................................................... 27 CHAPTER 3 : SUMMARY AND CONCLUSIONS ................................................. 38 Summary ......................................................................................................................................... 38 Conclusions .................................................................................................................................... 38 CHAPTER 4 : COMPARISON OF THE SLD MODEL TO MC SIMULATIONS 40 MC Model ...................................................................................................................................... 40 SLD Model ..................................................................................................................................... 41 Discussion ...................................................................................................................................... 44 Conclusions .................................................................................................................................... 53 CHAPTER 5 : CY ANIDE REMOVAL FROM WASTE WATER .......................... 57 Background .................................................................................................................................... 57 Treatment Options .......................................................................................................................... 58 Apparatus ....................................................................................................................................... 59 Procedures ..................................................................................................................................... 59 Results ............................................................................................................................................ 61 APPENDIX A .................................................................................................................................... 63 APPENDIX B .................................................................................................................................... 67 APPENDIX C .................................................................................................................................... 69 APPENDIX D .................................................................................................................................... 85 LIST OF REFERENCES .......................................................................................... 102 LIST OF TABLES TABLE 1: EXPERIMENTAL AND CALCULATED VALUES OF SATURATION PRESSURES OF PURE WATER. ................................................................................................................... 14 TABLE 2: EQUILIBRIUM BULK PRESSURES CALCULATED BY SLD FOR REDUCED DENSITIES. THE VALUE OF SATURATION DENSITY IS CALCULATED AT 211.22 K. .......................... 44 TABLE 3: CYANIDE REMOVAL RESULTS FROM AIR STRIPPER. ......................................... 62 TABLE 4: PARAMETERS USED FOR THE DETERMINATION OF WATER DIAMETER. ............... 64 TABLE 5: FLUID-SOLID INTERACTIONS AND POLARIZABIIJI‘IES. ........................................ 64 LIST OF FIGURES FIGURE 1-1 2 IUPAC CLASSIFICATION OF ADSORPTION ISOTI-IERMS. ................................... 2 FIGURE 1-22 EXPERIMENT AL DATA OF RUDISILL ET AL. (1992) AT TEMPERATURES RANGING FROM 298 K TO 398 K. .............................................................................................. 6 FIGURE 2-1: LAYOUT OF SLIT PORE MODEL. ..................................................................... 9 FIGURE 2-2: ACTIVE SITE MODEL ILLUSTRATION. ........................................................... 12 FIGURE 2-3: WATER ISOTHERMS IN SLIT MODEL WITH DISPERSION FORCES ONLY ............. 15 FIGURE 24: ILLUSTRATION OF POSITION OF SQUARE WELL POTENTIAL. ........................... 16 FIGURE 2-5: WATER ADSORPTION IN A SLIT WITH SQUARE WELL POTENTIAL ONLY. ......... 18 FIGURE 2-6: WATER ADSORPTTON IN A SINGLE 3er OF WIDTH 13.4 A. ............................. 19 FIGURE 2-7: EFFECTS OF ADDING A SQUARE WELL TO DISPERSION IN SLIT MODEL ............ 20 FIGURE 2-8: RELATIVE MAGNITUDE OF SQUARE WELL POTENTIAL. .................................. 21 FIGURE 2-9: WATER ADSORPTION USING A PORE SIZE DISTRIBUTION AT 373 K. ................ 23 FIGURE 2-10: INDIVIDUAL Ser ADSORPTTON IN THE PORE SIZE DISTRIBUTION MODEL. ..... 24 FIGURE 2-1 1: DENSITY PROFILE FOR 13 .4 A SLIT WIDTH. ................................................ 25 FIGURE 2-12: DENSITY PROFILE FOR 1 1.591 A SLIT WIDTH. ............................................ 26 FIGURE 2-13: RELATIONSHIP BETWEEN RADIAL AND RECTANGULAR POSITIONS. ............. 29 FIGURE 2-142 INDIVIDUAL CONTRIBUTIONS TO WATER ADSORPTION AT 373 K FOR ACTIVE SITE MODEL. ............................................................................................................ 30 FIGURE 2-15: DENSITY PROFILES AT e=0°, 45°, AND 90°. .............................................. 32 FIGURE 2-16: RELATIVE MAGNITUDE OF FLUID-SOLID AND FLUID-SITE POTENTIALS ......... 3 5 FIGURE 2-172 ISOTHERMS CALCULATED FOR ETHYLENE ON CARBON USING THE SLIT MODEL WITH Eps/K=95 KANDH=13.4 A. ............................................................................. 36 FIGURE 2-18: ETHYLENE ADSORPTION USING THE 5er MODEL AT 21 1.22 K WITH VARIOUS VALUES OF FLUID-SOLID INTERACTION PARAMETER. ................................................. 3 7 FIGURE 4-1: COMPARISON OF PENG-ROBINSON SLD METHOD To EXPERIMENTAL DATA OF REICHETAL. (1980)WITHEps/K=125 KANDASLTTWIDTHOF13.4 A .................... 42 FIGURE 4-2: COMPARISON OF ENERGY PROFILES AT L=9.1943 AND T=21 1.22 K. ........... 46 FIGURE 4-3: COMPARISON OF ENERGY PROFILES AT L=9.1943 AND T=339.42 K. ........... 47 FIGURE 4-4: COMPARISON OF ENERGY PROFILES AT L=4.5 AND T=21 1.22 K. ................. 49 FIGURE 4-5: COMPARISON OF ENERGY PROFILES AT L=4. 5 AND T=339.42 K. ................. 50 FIGURE 4-6: COMPARISON OF ENERGY PROFILES AT L=2.37 AND T=21 1.22 K. ............... 51 FIGURE 4-7: COMPARISON OF ENERGY PROFILES AT L=2.37 AND T=339.42 K. ............... 52 FIGURE 4-8: COMPARISON OF ENERGY PROFILES AT L=1.0 AND T=21 1.22 K. ................. 54 FIGURE 4-9: COMPARISON OF ENERGY PROFILES AT L=1.0 AND T=339.42 K. ................. 55 FIGURE 5-1: AIR STRIPPER APPARATUS. ........................................................................ 60 FIGURE A-l: DETERWATION OF WATER MOLECULAR DIAMETER. ................................. 65 FIGURE A-2: DETERMINATION OF Eps/K FOR WATER/CARBON SYSTEM. ............................ 66 viii NOMENCLATURE a - Peng- Robinson attractive energy parameter b - Peng-Robinson size parameter Eta - distance from first wall in slit model/fluid-solid diameter f - fugacity ' H - distance from carbon center to carbon center HCN - hydrogen cyanide k - Boltzmann’s constant L - distance from carbon surface / diameter of fluid molecule N - number of molecules P - pressure ppb - parts per billion ppm - parts per million r - radial position R - gas constant SQWL - square well T - temperature U - internal energy Xi - distance from second wall in slit/fluid-solid diameter Y - ESD attractive energy parameter 2 - direction of finite length 2 - compressibility factor ZR - radial distance from active site/fluid-solid diameter Greek or - spacing between carbon planes 8 - interaction parameter p — density a - molecular diameter I“ - excess ‘1’ - fluid-solid potential Subscripts atoms - property dependent on number of atoms attr - attractive term b - bulk property fa - fluid-active site property if - fluid-fluid property fs - fluid-solid property ix LCL - local property mean - arithmetic mean r - reduced property rep - repulsive term 55 - solid-solid property 2 - position dependent property Superscripts ig - ideal gas sat - saturation property Chapter1 : INTRODUCTION There is a need to better understand the behavior of water adsorption onto carbon surfaces. Water vapor is present as a contaminant in many systems where adsorption is used as a separation and purification technique as well as the regeneration of porous carbon. The shape of the Type V water isotherms suggest that the water molecules are more attracted to each other than they are to the porous carbon surface (Talu and Meunier, 1996). Water also takes on the Type III shape when adsorbed onto non-porous surfaces (Carrott, 1992). Figure 1-1 illustrates the five types of adsorption behavior according to the IUPAC classifications. Water tends not to adsorb very strongly on carbon surfaces at low pressures, but a sharp rise in the adsorbed amount is noticed in the region immediately below the saturation pressure (Rudisill, Hacskaylo, and LeVan, 1992). This phenomenon seems to be due to the fact that water does not wet the carbon surface when there is a lack of functional groups on the surface of the adsorbent. There is loss of hydrogen bonding that would otherwise occur in bulk water as a result of the lack of polar groups, and the water/carbon interaction is significantly smaller than the water/water interaction. Models such as those presented by Gubbins and Ulberg (1995) and Dubinin (1980) assume that the total adsorption of water onto these surfaces is very weakly dependent on the dispersion forces which can, therefore, be neglected as a first approximation. Most of the surface excess is due to the clustering of water molecules around such surface impurities as 5-0 "U (D '53 n 0 IE U) ‘5 < E a ..J o \ E < 14 151 P / Psa“ Figure 1-1: IUPAC classification of adsorption isotherms. 3 oxygen groups as well as other water molecules. The density of these surface groups also plays a role in the adsorption of most species (Segarra and Glandt, 1994). There have been few theories developed to help understand the mechanisms involved in water vapor adsorption. While the opinions on the mechanism may vary, most agree that there is a difference in the way that water vapor is attracted to carbon surfaces that is not seen with other common non-hydrogen bonding adsorbates. Gubbins and Ulberg (1995) propose that the bulk of the adsorption is due to functional groups in the carbon. Here, the authors model slits in graphitic carbon using grand canonical Monte Carlo simulations. They use a water diameter slightly smaller than that of carbon (3.15 AB .4 A). The graphite walls consist of 3 layers of carbon and no functional groups, even though this is the proposed mechanism of adsorption. In further work (Maddox et al., 1995), the active sites are assumed to be acetic acid (carboxyl) groups placed on the surfaces of a 20 A slit (center to center) with surfaces 30x30 A. The number of active sites used were 0, 2, 8 and 32 which were evenly spaced throughout the carbon surfaces. This means that a slit pore with 32 sites on its surface yields 28.125 A/ site. There is no mention of any comparisons made to the trends seen experimentally to indicate which site density may be more realistic. Dubinin (1980) reports similar adsorption mechanisms for water on activated carbons. He concludes that the dispersion forces that are so important in the adsorption of other fluids can be neglected in water adsorption. The assumption made here is that the active sites on the surfaces are unspecified oxygen complexes and that they are available to the water molecules to form hydrogen bonds. There are 1.43 mmol of these active sites 4 per gram of oxidized active carbon from sucrose as determined graphically by the theory of volume filling of micropores (TVFM) equation azaoch/(l—ch) (1) where a is the adsorption value for a relative pressure h=P/P“‘, a0 is the number of active sites, and c is a ratio of equilibrium constants kz/kl. Dubinin also states that as water clusters around the active sites, the fluid molecules become secondary sites of adsorption due to the formation of hydrogen bonds. Talu and Meunier (1996) propose a similar theory of water adsorption around oxygen containing sites, but suggest that there are three regimes to the Type V isotherm. There is the low loading region where water is adsorbed onto the active sites, intermediate loading where hydrogen bonds are largely responsible for adsorption, and the high loading region where the pore begins to fill and a plateau is observed. The development of this model is an attempt to characterize all Type V isotherms using water to check the results. The theory does a good job of predicting isotherms for different carbons but uses five adjustable parameters to generate data for each type of carbon. Work done by Huggahalli and Fair (1995) indicates that the shape of the water isotherms is due primarily to capillary condensation. Here, they take work done by Doong and Yang (1987) on volume filling theory from the original Dubinin-Radushkevich equation and decouple and linearize for use in their own model. This model is then compared to the experimental data of Rudisill et a1. (1992) for temperature ranges from 298 K to 373 K with some success, although this model is more dependable at higher temperatures. 5 This thesis focuses on the adsorption of gas phase water in various slits. The isotherms predicted by each of the methods are compared to experimental data by Rudisill et a1. (1992) seen in Figure 1-2 at temperatures ranging from 298 K to 398 K on BPL carbon manufactured by Calgon Carbon Corporation (Pittsburgh, PA). Water is the focus of this work because of its ability to form a network of hydrogen bonds with other water molecules which makes it difficult to model. The Elliott, Suresh, and Donohue equation of state is used in this work because of its ability to predict the thermodynamic properties of non-spherical and hydrogen bonding fluids such as water. This EOS introduces accuracy for linear chains of bonding fluids. Its advantage over more simple equations is that it incorporates a shape factor into the repulsive term of the equation. This allows for better agreement with molecular simulations for hard spheres and chains (Suresh and Elliott, 1992). The Simplified Local Density approximation is used in conjunction with the ESD EOS. This approximation assumes that the local density determines the interactions between fluid molecules and that density can be calculated analytically using any equation of state. om TITlITIITITIlr 5 5 0 5 0 11 =25... ac Ensues .9355. ac l7 c w 0 TlTerlili-ll- z m .5. m 5 0 3.2.5. .9— nU F .— Ensoem 11 5.2.5. .9. 0 Figure 1-2: Experimental data of Rudisill et a1. (1992) at temperatures ranging from 298 K to 398 K. Chapter 2 : ADSORPTION MODEL DEVELOPMENT Four major attempts have been made at modeling the adsorption behavior of water vapor. Each of the models developed included many of the same theories with specific ideas incorporated to make it unique. The models included a simple slit pore geometry, a slit with a square well potential, a slit distribution and an active site on otherwise smooth walls. M The Elliott, Suresh and Donohue equation of state was the backbone of each model. As stated previously, it accurately predicts the thermodynamic properties of hydrogen bonding fluids such as water. It takes the form Z=l+Z'°"+Z'"’, where "P: 4c” (2) 1-l9n Zdttr = qunY (3) 1+klnY Here, c is a shape factor for the repulsive term, q is a shape factor for the attractive temI, n is the reduced number density ( q=bp), Y is the attractive energy parameter, and z", and k, are constants for the equation of state. The fluid diameter as of 3.4 A is calculated using the van der Waals (vdW) volume of water and the Lennard-J ones diameter. Appendix A contains an explanation of the calculations and the figures used to illustrate the relationships between the vdW volume, the Lennard-J ones molecular diameters, and the ESD size parameter b. 8 A 10-4 interaction potential which makes use of the fluid-solid interaction potential at. was used in each of the models. It incorporates five layers in the form of: ‘ I 0.2 _ 0.5 _ 0.5 . Elam £104 Eta + a 4 \P(Z) : 47rpatomso'kgj3I O 5 ( O 5) O S I (4) f (Eta + 2a)“ _ (Eta + 3a)“ ' (Eta +4a)‘ I where ais the plane Spacing between the solid particles (3.35A/ 01,), 0]} is the average of the fluid and solid molecular diameters [(3, = ( O'fl‘ + 0..) / 2], Eta is the ratio of the distance from the carbon center to of, and paw, represents the number of atoms per square Angstrom (0.382atoms/A2). The potential in relation to the second wall can be calculated by replacing Eta with Xi, which is the distance from the second wall divided by the fluid- solid diameter (see Figure 2-1). For the active site model, the water-site potential was calculated using the Lennard-Jones 12-6 potential. ‘I’(r) = 4.03,,[-Z—;;+ER}—n-) (5) The radial distance is represented by r and ZR is the distance from the center of the active site reduced by the fluid-solid diameter. The interaction parameters Efi/k for water/carbon and af/k (k is Boltzmann’s constant) for water/active site are both fitted the same way for each potential to yield calculated isotherms that mimic the behavior seen experimentally. The SLD approach can be used with any equation of state. The thermodynamic properties of the fluid of interest can be estimated with Equations 6-8 below derived from superimposing the potentials according to SLD (Rangarajan, Lira and Subramanian, 1995). The chemical potential of the fluid is calculated assuming that the fluid-fluid fugacity can be estimated using the fluid-solid potential ‘I’. Figure 2-1: Layout of slit pore model. 10 tub =tufi(ztr)+,ufs(ztr) (6) #170) = p.(n+RT1n[f’Tm-] <7) -‘I’T(z,r)+SQWL fflr(z,r) = f. eXP[ kT ] (8) In these equations [1,, and fl, are the bulk chemical potential and fugacity, f0 is 1 bar, ,uo is the Chemical potential at 1 bar, and ya is the fluid-solid contribution to the chemical potential. Equation 8 and the rearrangement of Equation 6 allow the fluid chemical potential to be calculated from Equation 7. Note that ‘I’T(z,r) is ‘I’(z)+‘P(r) and that ‘I’(r) is only used in the active site model, as explained in the Active Site section below. In all other models, it has a value of zero and the fugacity and chemical potential are functions of 2 (position) only. In some cases, a square well potential SQWL is added to the slit model to increase the surface energy. This model is detailed in the Square Well section. For the ESD equation, the density and local pressure are obtained when the calculated fiigacity of Equation 9 equals the fugacity of Equation 8. f..t.(z,r) = y.P,.a¢.- ' (9) where the local pressure is represented by PLCL, y,- is the vapor mole fraction and is equal to 1 for a pure component calculation, and ¢,~ is the calculated fugacity coeflicient. Another constant throughout the models is the excluded volume of the slit. In previous SLD work of Rangarajan, Subramanian and Lira (1995) the region from the wall to Z/O'fi =O.5 has been excluded, assuming that no molecular centers could be contained in this area. This work assumes that the cutoff should be the point where the local fiigacity l 1 ft; is one tenth of a percent of the bulk fiigacity fb. This yields a difi‘erent calculation for each slit size. From Equation 8, one can calculate ‘I’(z,r)/k as approximately 2500 K when fg(z,r)/fb=0.001 at 373 K (assuming no square well). 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. For a slit of 10 A fi'om carbon surface to surface, the difference is slight. Using the 12-6 solute- solvent potential the cutoff is Z/Ufl‘ =O.416 at 373 K. The method of integration used for each model was the modified Simpson’s rule. In this method, the difl‘erence between the local and bulk densities is integrated in the correct geometric form over the distance to yield excess, I“. In the case of adsorption on homogeneous walls, this integration over the entire slit is divided by two to represent the adsorption in a slit for the two boundary walls. It is also divided by two in the active site model because we want to represent a particle on a wall where only halfis exposed to the fluid (see Figure 2-2) between the walls. The attractive term Y in the ESD equation (Elliott, Suresh, Donohue, 1990) is replaced with the local attractive term YLCL(z,r) which is dependent on the position from the wall or active site (Rangarajan, Lira, and Subramanian, 1995 and Subramanian, Pyada, and Lira, 1995). This requires a difference in the method of calculation for each model type. In slits, the local attraction is predicted based on the bulk value of Y and the position in rectangular coordinates. The calculations for the active site method were a preliminary test to see if the model could accurately predict water vapor adsorption. These calculations are explained in the Active Site section and equations are included in Appendix B. 13 These changes allow different concepts to be included in the same general theory. The major difl‘erences will now be discussed for the four models of adsorption along with the results obtained for each. Appendix C contains a copy of the main programs used for each of the methods along with the subroutines that are unique to the specific geometry. Appendix D contains all of the common subroutines. Slits This program was developed to use the ESD EOS to calculate the adsorption of fluids in slits. This model focused on physical adsorption of a gas phase fluid in the pores of BPL activated carbon. The adsorbent has a surface area of 1050-1150 m2 per gram of carbon as specified by the manufacturer (Granular Carbon data sheet #27-119, June 1986). Previous work by Reich et a1. (1980) specify a surface area of 988 mz/g as determined by the nitrogen BET surface area method. In this thesis, a surface area of 1100 mz/g was used as an average value from the manufacturer. The mean pore size is estimated as 13.4 Angstroms from carbon center to carbon center based on the fit to experimental data of ethylene (Reich et al., 1980) in Chapter 4. The geometry of the surface is assumed to be slit-like with smooth, homogeneous walls. The code is a modification of work originally done by Elliott et al. used to predict the phase behavior of single components and mixtures for BSD. During modification, several techniques for calculating thermodynamic properties were changed. The Newton- Raphson technique replaced the secant iteration for predicting local pressure because of its inefficiency when performing the diflicult calculations required in the modified program. Also, the iteration on density was replaced by analytically solving the cubic equation for the three roots. By setting the criteria for selection to be the root with the Fri IL, lb Ll FIFK Vi UK “a l4 smallest fugacity value, the possibility of obtaining the wrong root was eliminated. The densities of liquid water calculated by this equation of state are consistently lower than the experimental values. The bulk density calculated by the ESD is 0.72495 g/cm3 at 373 K and 0.102115 MPa. This is 24.4% lower than the experimental value of 0.95838 g/cm3 at 373 K and 0.101325 MPa. The results obtained from this model were not as accurate as we had hoped initially. Figure 2-3 shows the calculated isotherms for a temperature of 373 K with H=13.4 A and ea/k =35 K, 50 K, and 100 K. Increasing the value of Salk moves the threshold pressure to lower values of P/P“t and clearly shows that this model does not accurately predict the shape nor magnitude of water adsorption. Given in Table l are the saturation pressures for each temperature as predicted by the ESD equation of state along with the experimental values. Table 1: Experimental and calculated values of saturation pressures of pure water. Temperature P‘" (MPa) P'” (MPa) (K) calculated experimental 298 0.00369 0.0031690 323 0.0133 0.012344 348 0.03981 0.038563 373 0.1021 0.10132 398 0.23123 0.23201 Sguaflefl The slit program was modified to include a square well potential up to one fluid molecular diameter out from the surface of the carbon wall as seen in Figure 2-4. This was done in an attempt to represent the localized active sites in the first layer. The 24 T 1 T=373 K ‘ H=13.4 A 1 19 4 a 14 4 3 E g . ah I 9 L sf,/k=100 K x “" I 4 T a,/k=50 K .. , ,, J l w- J ‘ x ‘ f ‘ f' stlk=35 K / "-m-—H‘—A-I-I-I-I----I--- .‘11 A A; 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 PIPsat Figure 2-3: Water isotherms in slit model with dispersion forces only. aloe-COIIOOCCIIIIIOIIcan...one...an...oooouoounoonun-IOIOIOIOOIIIOIOIIOOII-OCIOIIOIQOIO coo-eguo-aonaouc oucllonuoaoonocuoo cocoooooooo.no-o-ocoooootouan-OIIIOI [It'll-IO...DUO....ICOI.OOIOIODICODCOIOODOOIIO 2 c0 = anon-000...... Figure 2-4: Illustration of position of square well potential. V = Water 17 expectation was to increase the height of the isotherm and possibly correct the shape to better resemble that of the experimental data. The use of square well potential in the model tends to create the same effect as active sites in that it gives the water molecules more of an attraction to the region near the surface of the adsorbent and allows them to pack more closely to one another within one molecular diameter of the carbon surface. Toward the center of the slit, the molecules are more dispersed and less dense due to a lack of functional groups and lack of other water molecules to bond with. Without the addition of the square well potential, the water excess predicted in a slit was insignificant for the same dispersion forces. The addition of this potential predicted approximately the same amount of adsorption as the slit model. Illustrated in Figure 2-5 are the effect that the square well alone has on the predicted behavior. The adsorption is linear and very small. Compare this to the results observed from the combination of the square well with the dispersion forces in Figure 2-6. The sum of the two is greater than the adsorption of either potential singularly. The effect of the square well potential, in conjunction with the dispersion forces, is to the move the threshold pressure to higher values with increasing square well potential SQWL as well as increase the total amount adsorbed on the surface. This trend can be seen in Figure 2-7 with constant physical fluid-solid interaction value of 35 K. The layer of water adsorbed is contained within the region where the square well potential is active. This indicates that the adsorption is due almost completely to the addition of this potential. The magnitude of SQWL is on the order of the fluid- solid potential, ‘I’(z) (see Equation 8), and can be seen in Figure 2-8 for Sfilk=35 K, SQWL=275 K, and H=13.4 Angstrom. This value of SQWL is chosen because of the 18 T=373 K H=13.4 A SQWL only 0.03 T 0.025 «A 002 - - - SQWL=150K ‘ — —sow1.=275 K sowr=400 K A m = O 2 0.015 - v 3 / L-t / / / I 0.01 + 0.005» .,.«*" 0 "‘ 7‘ II 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 PIP‘“ Figure 2-5: Water adsorption in a slit with square well potential only. 19 9 .. l 8 T 7 T 6 T A T=373 K g 5 - H=13.4 A g sfs/k=35 K E SQWL=275 K r 4 + Q Li 3 “r- 2 I l 1 + l 0 4 t t I - t t T T . 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 HP“ Figure 2-6: Water adsorption in a single slit of width 13 .4 A. 10 . I"’x (mmollg) 34L 20 T=373 K H=13.4 A ayk=35 K - - - SQWL=150K 1 — -SQWL=275 K SQWL=400 K J —--T fin]. I I 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 HP“ “nan-I'll - I 4 l I T T 1 Figure 2-7: Effects of adding a square well to dispersion in slit model. 21 2000 4 1500 4» =373 K ‘ H=13.4 A ’ aB/k=35 K SQWL=275 K 1000 4» + 2 El 0 G 500 «~ (I) A 0 N v 9* 0 0.: -500 . -1000 i 2,6“ Figure 2-8: Relative magnitude of square well potential. +fluid—solid Potential +SQWL Potential 22 placement of the threshold pressure near 0.5 as seen experimentally. The ath of 35 K is based on its relationship to the polarizability as detailed in Appendix A. Slit Distribution The effect of including a pore size distribution for the modeling is simply to slope the region above P/P“‘=0.8 and create less of a plateau as the one observed in Figure 2-6 for a slit. The individual pores fill at difl‘erent rates with the smaller pores filling last at higher pressures. Figure 2-9 shows a single isotherm for the combination of 5 slit sizes at 9.78, 11.591, 13.4, 15.209, and 17.018 Angstrom. These individual slit contributions are distributed as 10, 20 and 40 percent of the adsorption with the smallest percentage going to the smallest and largest pore size, 40 % to the median slit, and 20% to the remaining slits. Figure 2-10 illustrates the individual slit isotherms. With each of the larger slit sizes, the pores exhibit a plateau region at higher pressures. This behavior is not observed for a slit size of 9.782 A or 11.591 A. The initial thought was that the square well overlapped for these slits and the pore was able to retain more fluid due to the increased energy from opposing walls. However, it was observed that the 11.591 A slit has a space in the center of the slit of approximately 1.4 A that did not rely on the energy from the square well. It is also observed that the density profiles of 13.4 A (Figure 2-11) and 11.591 A (Figure 2-12) slits with the same parameters looks very Similar, but the larger slit has a thicker adsorbed layer. Since the 11.591 A slit has not yet filled the region where the square well is turned on, this would indicate that the pore still has room to accept more fluid at a pressure of 0.1 MPa. This also supports the idea that more energy is required to fill the smaller slits due to the limited space available in them. 23 a -. 7 - 6 A 5 J. - T=373 K g Slit Distribution 2 afs/k=35 K E 4 T SQWL=275 K 5' i— 3 3 2 . 1 3 0 r t t t 4.1 1 *fi 1 t t 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PIP‘“ Figure 2-9: Water adsorption using a pore size distribution at 373 K. 24 9 _. 3 I 1 =373 K I 8 H=13.4 A ,» . t 7 a,/k=35 K g ‘- SQWL=275 K - 6 -<>- A 51 = 5 'Ji- 8 E. x 4 A» 0 i... 3 — - H=9.782A ‘ -~ - H=11.591 A r r '- H=13.4A W ”H=15209A 2 + H=17.018 A 1 I 0 W ‘7 I I I I - - —i - - _J; - I T 4I 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 PIP‘“ Figure 2-10: Individual slit adsorption in the pore size distribution model. 25 0.04 ~~ 0.035 T 0.03 4i I T=373 K «5" 0°25 * H=13.4 A 0 sit/1635 K 2 SQWL=275 K g P=0.1MPa v 0.02 1 '3 c C) O 0.015 . 0.01 A» 0.005 A» 0 J——-‘ + _A 0 0.5 1 1.5 2,6" Figure 2-11: Density profile for 13 .4 A slit width. 26 0.04 «r 0.035 «A l 0.03 + T=373 K at" 0°25 " H=11.591 A E e,,/k=35 K = SQWL=275 K g P=O.1MPa v 0.02 r E in r: o O 0.015 A - 0.01 I» 0.005 A 0 -————d I H 0 0.5 1 1 5 2,0“ Figure 2-12: Density profile for 11.591 A slit width. 2.5 27 In each of the three larger slits, the bulk of the fluid adsorbed is within the square well region. Active Site An active site model was developed to act as functional groups on the surface of a carbon slit. This method integrates the density of the water vapor radially around the active site between the walls of the slit and considers the number of similar sites that may exist in a sample of activated carbon. It is a modification of work done previously by Subramanian et a1. (1995) for clustering in supercritical fluids. Figure 2-2 illustrates how such a model would look. Such a complicated geometry warrants complicated integrations to calculate the attractive term YLCL(z,r) which depends on both the radial and rectangular positions. Before attempting this double integration, it was first necessary to determine if this model would correctly duplicate the trends seen experimentally. One way to estimate the outcome of the isotherms is to do separate integrations of the attractive term over the slit and the spherical active site and develop a relationship between the two terms. The exact equations for both of these geometries are presented in Appendix B from work done by Rangarajan, Lira, and Subramanian (1995) and Subramanian, Pyada, and Lira (1995). The calculation of YLCL(r) for three dimensional clustering takes the form _th;(r):1+c (10) b where C is a correction factor to account for the volume of the active site itself. For the estimation of the integration over I and z, the rectangular equation can simply be substituted into Equation 10 to yield YLCL in the form 28 YLCL(r):YLCL(z)+C (11) Y. I; As C approaches zero (radius of active site-)0), Equation 10 approaches a bulk fluid and Equation 11 approaches the value seen in a slit. This estimation at difl‘erent values of 0 (angle between the active site and the molecule of interest as seen in Figure 2-13) should provide some clue as to the shape of the isotherm. The calculation of r in relation to z is based on the sine of the angle where sin(0)=z/r. To test the form of the isotherms, the density profile at a fixed angle was assumed to hold at all angles for the purposes of integration. Therefore, the isotherm for an angle with the thickest adsorbed layer will over estimate the true isotherm if thinner adsorbed layers exist at other angles. At angles of 0°, 10°, 45°, 70°, 80° and 90°, the trends become clear for isotherms generated at 373 K with Sfi/k =40 K, a;./k=550 K and site density SD of 2.5 mmol/g activated carbon. From Figure 2-14 we can predict that the integration over 0=0-)1t would yield a shape that is not similar to the experimental data of Rudisill et al. at 373 K based on the individual adsorption isotherms. The isotherm predicted at 0=0° is an over-estimation based on the density at the wall being continuous throughout the slit. This is not the case. The adsorption decreases toward the center of the slit and away from the active site. This idea is validated by the isotherm generated at 0=10° which shows how quickly the adsorption drops off at such a small angle away for the surface. The isotherms are almost identical between 0=20° and 65°, which constitutes approximately half of the 180° range since the same effect is expected for 0 between 115° 29 Figure 2-13: Relationship between radial and rectangular positions. 30 35 373 K I ayk=40 K ayk=550 K SD=2.5 mmol/g 30 .. g .. - ' 25 + a 20 A " " ‘ 9=0 E ..... ”0:10 E ——0=45 g - - -- 9:70 x ——0=80 1515“ ----0=90 - Rudisill et al. (1992) 10 i 5 .. l’ ’ - / f 0 1 ‘ ”‘4“ L” I ““7““ i 1 t t t o 0.1 0.2 0.3 0.4 0.5 0.5 0.7 0.8 0.9 1 PIPsat Figure 2-14: Individual contributions to water adsorption at 373 K for active site model. and The fig: ads: 110 caltt incre. 011 lhi where shows from c POIaIiz Physic: the cor is fit to along II 3:15.. (1992) ( lhé‘ SUI‘fa density 01 “mined 31 and 160°. The isotherms for 0<20° would comprise less than 25% of the angle range. Therefore, the contribution to the isotherm from the largest curves at 0° and 180° in Figure 2-14 would add little to the actual calculations. The shape and magnitude of the adsorption curve would primarily resemble that at 45°. Increasing the width of the slit demonstrates an interesting result. The value of YLa/Yb approaches unity in the center of the slit for both the radial and rectangular calculations. Therefore, as the slit size increases, the density in the center of the slit increases as well. The density is driven by the value of the attractive term which depends on the position in the slit. The fluid is just as likely to condense in the center of the slit where YLCI/Yb-9 1 as it is to condense near the wall and the active site. Figure 2-15 shows the density profiles for this model with 0=0°, 45°, and 90° for a slit width of 24.0 A from carbon center to center. As noted previously, the water/carbon parameter Sfilk is assumed based on the polarizabilities. Appendix A contains the relationship between the values of efi/k and this physical property. Values of 35-40 K are reasonable for water. They vary slightly to yield the correct threshold pressure seen experimentally. The water-active site parameter ea/k is fit to give the correct threshold pressure and to adjust the magnitude of adsorption along with the site density. A value of site density SD=2.5 mmng converts to a square of 8.55 A on each side. Compare this value to the arbitrary values reported by Gubbins et al. (1992) of 28. 125 A per site as well as the 1.43 mmng reported by Dubinin (1995). Using the surface area of BPL carbon (1100 mz/g), Dubinin’s site density yields 11.3 A/Site. The density of 8.55 A/site will likely cause overlap of active sites in a slit of 13.4 A (as estimated from data of Reich et al.), yet the magnitude of the isotherm is well below T=373 K H=24. A =40 K s,./k=550 K P=0.1MPa a, 32 l - -0=0 ---0=45 -U'J“ IIIIIIIIIIIIIIIIIIIIIIIIII 7...! llllllllllll I IIIIIIIIIIIIIIIIII .- l ./ q d [+ 4 1‘ d + “I + iii—l 5 4 5 3 5 2 5 1 M 0 4 0 3 O 2 0 1 0 0 0 0. 0 0 0 0. 0 0. 0 0 0 0 0 Ana—8:25 33:00 5.5 3.5 ['1fo 0.5 Figure 2-15: Density profiles at 0=0°, 45°, and 90°. 33 experimental adsorption. This indicates that the model is unrealistic in its predictions. The site density would have to exceed values seen in actual carbons in order to obtain the magnitudes needed, but the shape would still be inaccurate. An explanation as to why the active site model may not have worked well enough to provide the trends expected is that the addition of the Lennard-J ones potential ‘I’(r) in conjunction with ‘I’(z) act similarly to the results of including a square well potential illustrated in Figure 2-8. The difference between the two is that the potential in rectangular coordinates for the active site model has taken on the magnitude and approximate shape as the square well potential in the previous model. Figure 2-16 shows what the two potentials together look like for 0=45°. The two figures are very similar in both shape and magnitude. This means, essentially, that the only difl‘erence between the square well model and the active site model is the radial integration over density to yield adsorbed excess. It’s also important to emphasize that the slit model works very well for the adsorption of ethylene on carbon as seen in Figure 2-17. There is no need to manipulate the model by adding energy in the form of a square well potential. The fluid-solid interaction parameter is 95 K and the slit width is 13.4 A. The ESD equation of state does a good job of predicting this non-hydrogen bonding fluid with the SLD assumptions. This indicates that there is a physical feature that is missing fiom our model that is obviously important when attempting to model a hydrogen bonding fluid such as water. The isotherms calculated for ethylene adsorption in slits supports this conclusion. With various values of the interaction parameter from 95 K to 25 K in Figure 2-18 at 34 211.22 K a gradual rise is seen when Sfi/k is decreased to 25 K. There is no rapid increase as seen with the same model for water isotherms. This, again, indicates that there is a physical phenomenon that we have not accounted for in the hydrogen bonding isotherms. 35 2500 ~~ T=373 K H=13.4 A chlk=40 K 2000 __ Sfa/k=550 K SD=2.5 mmoI/g 1500 + 5 1°00 " —fluid-solid Potential (2) 1': +fluid—active site Potential (r) 9 74‘ V 500 ~~ 9' 0 4 A 1 1 1 1 3 1.4 15 1 1 9 \ . -500 1» -1000 ‘5 ZIOff Figure 2-16: Relative magnitude of fluid-solid and fluid-site potentials. f“. x X X 4"" u”'~—. p a." a“ . I 5' .l I I 1"" (mmollg) I Reich et al. at 301.4 K x Reich et al. at 260.2 K A Reich et al. at 212.7 K I - - ~ Slit calculation at 301.4 K T; t ‘ W *Slit calculation at 260.2 K i Slit calculation at 212.7 K 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 P (MPa) Figure 2-17: Isotherrns calculated for ethylene on carbon using the slit model with Sfilk=95 K and H=13.4 A. 37 1‘” (mmollg) 0 0.1 0.2 0.3 0.4 0. 5 0.6 P (MPa) Figure 2-18: Ethylene adsorption using the slit model at 211.22 K with various values of fluid-solid interaction parameter. Chapter 3 : SUMMARY AND CONCLUSIONS Summary Modeling the adsorption of water on carbon is an industrially important problem. This work has applied the ESD equation to model the water-water interaction to several modeling techniques developed to represent the adsorption surface. The basic slit pore model had very little success in predicting the behavior of adsorbing water. The addition of a square well potential to the slit model was unsuccessful in following the trends seen experimentally as well. The additional energy of the square well was not enough to overcome the unlikely geometry of the model. Incorporating slit pores of differing widths served to eliminate the plateau region seen in the previous model above the reduced pressure of 0.8, but still did not increase the magnitude or improve the shape of the isotherms. The active site model predicts the behavior of water vapor by assuming that the fluid adsorbs in mounds around active centers on the carbon surfaces. Despite the more suitable geometry of the surface, this model did not yield any better agreement with experimental data than the previous attempts. Conclusions The active site model, which many suggest is the correct mechanism for water adsorption, is incapable of simulating the behavior of water vapor on carbon surfaces in the current form. Even with three adjustable parameters, the model does not correct the 38 39 Shape or magnitude of the previous models and is actually very similar to the square well model. Also, the site density is physically unrealistic and cannot be seen in actual carbons. In essence, we have taken on the challenge of a very dificult problem. We believe that the models presented in this thesis lack an important physical feature which prevent them from capturing the characteristics of water adsorption. We have been unable to identify the missing factor, but there are several possibilities that may be explored in future work. The first option includes exploring the non-local efl‘ects that may be important for water adsorption because of its unique ability to form hydrogen bonds, but negligible for non-hydrogen bonding fluids. Also, adsorbing water may hydrogen bond difl‘erently than what we see in homogeneous fluids, which may mean that it is necessary to represent the bonds at the surface explicitly and allow the equation of state, which we know to represent the bonding in homogeneous fluids well, to continue to predict the behavior away from the walls. Chapter 4 : COMPARISON OF THE SLD MODEL TO MC SIMULATIONS In an effort to evaluate the accuracy of the Simplified Local Density (SLD) assumption, a study was done to make comparisons with Monte Carlo (MC) simulations using the adsorption of ethylene molecules on carbon as the system of interest. The SLD approach was introduced as a means of modeling gas adsorption onto carbon surfaces by Lira et al. (1995). It, along with the MC method, assumes a slit-like structure for the carbon pore. These slits are assumed to be infinite in the x and y directions with a finite 2 length as seen in Figure 2-1. It is necessary to go into the background for both of the methods in order to clearly explain the basis for comparison. Both models use a 12-6 fluid-fluid or fluid-solid interaction potential as modeled after the Lennard-J ones pairwise potential. ~13. (r) = 4% [(3,1) ' (3%)] Here, one is the diameter of a fluid molecule, r is the separation distance and at, is a fitted parameter for the bulk ethylene well depth where x is either fluid for MC or solid for SLD. MC Model The MC simulation was based on a canonical ensemble. The initial configuration of fluid particles corresponds to molecules placed in a fee lattice. The number of molecules, N, is directly proportional to the density of the system p. Smaller values of N were 40 41 needed for lower values of L, due to the smaller system sizes. Simulations for systems of 108 to 864 ethylene molecules were performed for ethylene adsorption in carbon slits with widths ranging from 4.22 A to 38.8 A from carbon surface to surface. Monte Carlo configurational sampling was performed until it was determined that the original lattice had been relaxed and equilibrium conditions established. Additional sampling was then carried out to obtain ensemble-averaged density and molar energy profiles. SLD Model The SLD approach assumes that fluid-fluid interaction is most greatly determined by the local fluid density. This assumption makes it possible to analytically solve for local density. The Peng-Robinson equation of state is used with the SLD assumption to predict the thermodynamic properties of the fluid. The model has been refined slightly from the above citations. In the previous publications for slits using this SLD model, the fluid between the wall surface and Z/Ofi = 1 was ignored. In this study, the cutoff distance is selected to be the point where the local fugacity of the fluid was approximately one tenth of a percent of the value of the bulk firgacity. This cutoff point has a different value for each slit size. For the SLD model, the parameters used were chosen based on experimental data for ethylene adsorption onto activated carbon (Reich et al., 1980). The fluid-solid interaction parameter Sfi/k and the slit size were varied in order to obtain the trends seen experimentally at temperatures of 212.7 K, 260.2 K, and 301.4 K. A value for Sfi/k of 125 K and a slit width of 13.4 A from carbon center to center were selected to fit the data reasonably and can be seen in Figure 4-1. The Sfi/k is not required for the SLD model. 42 A A A I g 5- , . x E * " E x 4 4- G t_. hi 3 r I Reich et al. at 212.7 K A Reich et al. at 260.2 K 2 x Reich et al. at 301.4 K { PR-SLD at 212.7 1 1* -— -PR-SLD at 260.2 K Jr PR-SLD at 301.4 K 1 0 l ‘ A a i i o 50 100 150 200 250 P (psia) Figure 4-1: Comparison of Peng-Robinson SLD method to experimental data of Reich et a1. (1980) with wk =125 K and a slit width of 13.4 A. 43 For the MC calculations, a value of ea/k=224.7 K was selected. The MC method will not match the critical temperature and pressure predicted by the SLD method due to the mean-field assumption of the SLD model. Further, the SLD and MC methods predict different state densities for non-ideal gas conditions, making direct comparisons impossible. Therefore, to conduct meaningful comparisons between the two models, a common basis between the two methods was established in which the energy profiles were computed at the same reduced densities pan/p”, where p‘“ is the liquid saturation density at the subcritical temperature of 21 1.22 K. The mean density is calculated on the basis of L (free volume). Comparisons at two temperatures are provided, one subcritical and one supercritical. The comparison of energy values was obtained in terms of the dimensionless energy per molecule for each case. For the SLD approach, the internal energy departure firnction was used in the form of (U -LI‘)/RT (U— ___U_R“)_ H52)” dp b where Z=1+Z'°"+Z'"’— -1+——p—a-&- . It rs interesting to note that the 1— -bp RT 1+ 2bp- -b p repulsive energy does not contribute to the departure function for the Pang-Robinson equation of state. To discuss the general trends observed, several slit widths and reduced densities were chosen for the two temperatures. Table 2 illustrates the values of the reduced densities used for the comparison and the pressures that they correspond to using the Peng-Robinson equation of state. In this table, the saturation density p“It at the subcritical temperature (211.22 K) is used to reduce the densities at both the subcritical and supercritical temperatures. In order to model realistic pressure values, a limit of 50 MPa bulk pressure was chosen for each temperature and slit width, as calculated with the SLD model. Table 2: Equilibrium bulk pressures calculated by SLD for reduced densities. The value 44 of saturation density is calculated at 211.22 K. T=211.22 K T=339.42 K Pmcm/put P (MP3) P (MPa) L=1.0 0.2 1.526E-6 5.859E-3 0.25 4.688E-3 1.289 0.27 0.3516 50 0.278 50 L=2.37 0.1 1.146E-5 7.7S4E-3 0.3 2.148E-3 0.2202 0.539 0.125 50 0.583 50 L=4.5 0.1 2.836E-4 7.076E-2 0.3 0.1265 2.708 0.579 0.3256 50 0.674 50 L=9.194 0.1 2.312E-2 0.7920 0.3 0.4647 7.556 0.579 0.6353 50 0.718 50 Discussion Each pair of the following figures may be used to compare the reduced energy profiles calculated by the two methods. The MC method is well-known to produce layers in the density profiles. We found that both MC attractive and repulsive configurational 45 energies exhibit oscillations due to these layers. However, the effects of attraction and repulsion largely cancel each other and the total energy profiles are relatively structureless. L=9. 1943 For a slit width corresponding to L=9. 1943 (this corresponds to (L + 03/65) = 10, where G..=3 .4 A and Ga=4.22 A ) two distinctive behaviors were noted for each of the two temperatures as illustrated in Figure 4-2 and Figure 4-3. It should be noted that the energy profiles pertaining to the MC method have been averaged at each position to result in the symmetrical profiles shown. At the lower temperature of 21 1.22 K and reduced density of 0. 1, a thicker surface “energy” layer of fluid and a flatter energy profile are observed for the MC method with both models showing a lack of fluid in the center of the slit. At a reduced density of 0.3, the difl‘erence in the thickness of contact layers decreased and both models predict low fluid density in the slits center. Again, the MC profile is flatter. At p,=0.579, the SLD method does not predict fluid in the center of the slit, whereas the MC method predicts a more uniform density profile. At the highest reduced density, corresponding to a value of 0.71 8, similar behavior is observed for both methods. This behavior corresponds to fluid being present in the center of the pore, but the fluid is much more structured in the MC method. The energy near the wall is found to be more density (pressure) dependent for the MC method. At the supercritical temperature of 339.42 K, the same trends in behavior are observed for all reduced densities. However, a slightly thicker contact layer is observed in the MC method for the lowest reduced density of 0. 1. At the highest reduced density, a uniform energy profile is observed for the SLD approach, whereas the MC simulation 0 . :v.-:VA~ -2 I. ar"’“"”~n... 0" " -4 L I “‘1‘ g I [xv] kd"\ 0 -6 I...” I .m— C I \-.1 La! 1 m "n‘ ‘0', \' \ I c .811 \‘ ..... I ...... ‘a..-‘CQ‘.-.. 1] _p___0.1 “" '-‘v I’ .\ .1..§ -' ‘3“. ' - A. .p- l". .r- .‘\ I ‘u'.’ II. “9:03 “‘0 x.’ “\p‘ ‘" “" '" "’ ‘c"‘\./ * " ~ p=0.579 ""'" " p=0.718 -12 4 i i i A. + 4% a i 0 5 1O 15 20 25 30 35 4O 45 Distance from carbon center (Angstroms) a.) MC energy profiles. Reduced density is represented by p. o - I i r "‘ ' " "I -05 -l_ 1 l I : j i . . l -1 at I i , : : I -1.5 w >. i z : l E, -2 I I a I 345, I . . , "J ' l : 2 -3 - :5, . . l —p=0.1 -3 5 A \‘2‘ I r z ‘ f.’ m “hp-=03 -.-I .....: .._ "-9.055 -4 °-._ ._.-"'I "" ' p=0.718 4.5 A" . 4. 0 5 10 15 20 25 30 35 40 45 Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-2: Comparison of energy profiles at L=9.1943 and T=21 1.22 K. Ji 1"”'~"~c l r “I -2 I “ ’ x .3 f \ I >‘ i 1’ “ 94 f \ i 3 I ‘ .5 J m R. I . ...... . \ ,, -6 i 1" "n’ “ ' fl... km“ | I .‘. ’0..- .-..‘ _pfl.1 -7T ‘\ .’ ‘I‘ a. f “ ‘vu’ 0a.. ’0 - -M.3 '8+ " ’0 O ‘0 I _ X . " ' " p=0.579 -9 I I ' I I I I I '1 fiI 4] 0 5 1O 15 20 25 30 35 40 45 Distance from carbon center (Angstroms) a.) MC energy profiles. Reduced density is represented by p. ”m-.. ' ”0.” ~. ',cr ‘5‘ ~ " .'.. —pi).l ’. 'g' '“ "*p=0.3 -_..-‘ " " " p=0.579 15 20 25 30 35 4O 45 Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-3: Comparison of energy profiles at L=9.1943 and T=339.42 K. 48 shows a slight increase in the dimensionless energy profile. This may be due, again, to slightly different filling pressures. The MC energy is more density dependent at the wall. L=4.5 For a slit width corresponding to L=4.5, the profiles at the subcritical temperature show similar behavioral patterns as those observed in the larger slit at all densities. This behavior is illustrated in Figure 4-4. At the supercritical temperature, Figure 4-5 shows similar trends to those observed in the larger slit. One other note is that contact layers in MC method are hard to distinguish due to the more uniform density profiles at this point, whereas the SLD method gives consistent patterning of these layers with large mean density decreases in the center. L=2.37 For a slit width corresponding to L=2.37, the greatest discrepancy between the two methods is observed at the reduced densities of 0.1 and 0.3 as seen in Figure 4-6. There is a layer of fluid observed at the walls with an absence in the center of the slit and an accompanying lack of energy for the SLD method. This is not the case for the MC method which indicates that there is density in the center of the slit even though the energy profile contradicts this observation. The behavior patterns between the two theories for the higher reduced densities of 0.539, and 0.583 predict similar energy profile and very dissimilar density profiles. For the supercritical temperature, Figure 4-7 shows similar trends as Figure 4-6. 49 0; -1.L adv- 2W r, “s. ' I I \ .3 ‘ ~5le \\"\ 5“ “""""' m», .. -5 I , C I l.l.l '6 iv .' ....... n. . ”J _p=0.l -7, \v “~‘." ~‘ . 10‘ ~" In. “2.:79 \ l I a a -8 -‘ I. .5 r . ‘ ' . - -- .674 -gi g" \'_-l \u- 9‘0 107 % i fl) 0 5 10 15 20 25 Distance from carbon center (Angstroms) a.) MC energy profiles. Reduced density is represented by p. 0 to.” r.' "-\ -0.5 i . II I , I -1 —I- l . I I a ' i -1 5 l . I > ' | i 5’ '2 i ' ' ' i :45 ~ " ' I.” \‘zr.o) : : \ .a' —p=0.l '3 ’- .\. 5" . | {,7 .- "p=0.3 . ~-,_ . . - - - .579 -3'5 ._ ~“""" L'”:’/. '" ' 3.674 -4 T s, “I...” -4.5 I + I I I H 0 5 10 15 20 25 Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-4: Comparison of energy profiles at L=4.5 and T=211.22 K. .‘ZI‘L’,( I-l '._-—r—-"‘a \ L f I -2 I I \ 9 l \ ' 3-3 , w} \A g 5 Lm-erwvr ’ “A's/Ar-w-i 1:,4 ’ - I W , : —p=o.1 -5 s ‘ .: _... -p=03 6 ' ' “g" ---p=o.579 \ ‘l \. '“v -7f .- .'+ g 4 0 5 1O 15 20 25 Distance from carbon center (Angstroms) a.) MC energy profile. Reduced density is represented by p. Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-5: Comparison of energy profiles at L=4.5 and T=339.42 K. r! 51 .1 -o- -2 -. > 9-3 3 m .4 db —p=().l -5 .. """" “p=0.3 " " " p=0.539 F .5 .t "" ' p=0.583 ‘7 T 0 2 4 6 8 1O 12 14 Distance from carbon center (Angstroms) a.) MC energy profile. Reduced density is represented by p. n 0» __‘ I F” ~. 1 I -o.5 4 I ‘ 5 . I I I ~ 9 -1 -I d g i ll! 0 i r C "I i i, In -15 I ti g I —p=0.l 5‘ U I \ U y - ~p=o.3 2 i I \ l! - - - p=0.539 \\ I _ . -2.5 I I I a I I I o 2 4 6 8 1o 12 14 Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-6: Comparison of energy profiles at L=2.37 and T=211.22 K. 52 0 2 -1q. 1 5 2 i i ‘1‘ > :’K\ 2'” I 9-34- . K l-P‘ r o . \”\~/“ . 1:.4 , . m ." "’. _p=0.1 5? §“ I" mmp=0.3 .64. ~."'o-v'*"" ""p=0.539 -7 ‘L I I I . . I O 2 4 6 8 10 12 14 Distance from carbon center (Angstroms) a.) MC energy profile. Reduced density is represented by p. ——w&1 -~fim --wm$9 I N r F 4 0 2 4 6 8 1O 12 14 Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-7: Comparison of energy profiles at L=2.37 and T=339.42 K. 53 For a slit width corresponding to L=l, the behavior is illustrated in Figure 4-8 and Figure 4-9. At both the lower and supercritical temperatures, quite similar behavior is observed for the designated densities. This is the opposite of what was expected at the outset of the comparison project when we expected the smallest pores might show the greatest differences between models It should be noted that the reduced densities at this slit width cover a much smaller range than at the other slit sizes. The pressure ranged from a value of 1.526E-6 MPa to 50 MPa as noted in Table 2 above. In both methods, similar energy profiles are observed. Conclusions Several trends are common for the slit widths considered. The SLD energy profile is “thinner” and more pronounced near the wall than the MC profile for lower densities. At high densities, the profiles are extremely similar. The SLD energy profile doesn’t have as much density dependence at the wall. The Peng-Robinson equation does not have a contribution to repulsive energy, so it is not possible to directly compare repulsive and attractive energies of the two models. This work has not directly uncovered the reason for the smaller density dependence of the SLD energy profile at the wall. Comparisons of the MC and SLD density profiles are difiicult due to the layered structure of the MC results. The energy profile via SLD is directly related to the local density, whereas in MC, nonlocal effects arise. Ifthe efi‘ect is due to nonlocal effects, then it is a shortcoming of the SLD method, and the deficiencies are more severe for larger pores. The other possibility is that the lack of repulsive fluid-fluid interaction in the Peng-Robinson is the primary shortcoming. Currently, the repulsive contribution to the Peng-Robinson equation is 54 r 9y a: O '— — p=0.2 — p=0.25 -— p=0.27 ------ p=0.278 #L i L l l I T f T 1 2 3 4 5 6 7 8 Distance from carbon center (Angstroms) a.) MC energy profile. Reduced density is represented by p. -o.1 ~ l -o.2 .- >to.4 -- a: -o.5 «» r: ll] -0.6 1’- -o.7 «» -O.9 «~ b.) Figure 4-8 - - p='0.2 — p=o.25 — p=0.27 ------ p=0.278 J , 3 i 5 e % é Distance from carbon center (Angstroms) SLD energy profile. Reduced density is represented by p. : Comparison of energy profiles at L=1.0 and T=211.22 K. 55 — -p=o.2 _p==0.25 — p=0.27 1 2 3 4 5 6 7 8 Distance from carbon center (Angstroms) a.) MC energy profile. Reduced density is represented by p. Energy I I I f: , S S o ~L~———+——-+———~+———+——+——4——+——1 .C'> u -0.5 "' "" p=O.2 -0.6 — p=0.25 — .270 -O.7 p=0 '08 i i fi 41 i I I 0 1 2 3 4 5 6 7 8 Distance from carbon center (Angstroms) b.) SLD energy profile. Reduced density is represented by p. Figure 4-9: Comparison of energy profiles at L=l.0 and T=339.42 K. 56 modeled using the bulk fluid mean field term which is known to be deficient even in bulk fluids. An increased fluid-fluid repulsion near the wall would result in a flatter density profile at a fixed mean density, improving agreement between the models. However, rather than modifying the Peng-Robinson equation, it would be preferable to adapt the SLD to a difl‘erent equation of state that has a repulsive contribution to energy, and to develop a greater understanding of the repulsive energies near the wall. There is a significant cancellation of attractive and repulsive contributions in the MC energy profile, which is quite structureless relative to the density profile. This cancellation evidently extends to the local entropy, since the local chemical potential is invariant with respect to position, and is given by a combination of energy and entropy. Chapter 5 : CYANIDE REMOVAL FROM WASTE WATER Background Work was started on a Phase I Small Business Innovation Research (SBIR) proposal granted by the United States Department of Agriculture (USDA) to investigate the feasibility of recovering cherry flavoring from the fi'uit’s pits. This grant was a source of funding for six months and an opportunity to work with Natura, Inc, one of the businesses associated with Michigan Biotechnology International (MBI). During Phase 1, several objectives were set to be accomplished in this 6 month period. They included the following: 6. 7. Screening adsorbents to be used for recovery. Generating breakthrough curves for the adsorbents and optimizing flow rates and bed volumes. Modeling the regeneration of the adsorbent. Determining cyanide treatment and disposal methods along with experimental validation. Economic and engineering analysis of the above. Preliminary design of apparatus. Preparation of a final report. It was necessary to develop a method to remove hydrogen cyanide from the waste water of the proposed plant. Processing of the cherry pits requires a one hour hydrolysis in order to get the desired products, benzaldehyde and benzyl alcohol, into solution. Since 57 58 cyanide naturally occurs in the hydrolysis of cherry pits, extraction of the desired products leaves the cyanide in an aqueous solution. The proposed waste output of this process is 622 gallons per hour with approximately 10 mg/L of cyanide. The proposed area of the plant is near Traverse City, Michigan, where eighty percent of the cherries produced in the United States are grown. Disposal limits for this area call for HCN concentrations of 30 pg per liter of solution or less. Therefore, this contaminant must be reduced in the waste stream, which can then be used to spray irrigate. Treatment Options Several conventional methods were researched in order to determine the best method of removal. Activated carbons are typically used to reduce metal cyanide concentrations. This method requires a pH between 6 and 8.5 and the addition of copper chloride. There is also a problem with organics competing with the cyanide for active sites, which requires more carbon. Since there is a large amount of organic matter in the aqueous solution, which has a pH of 5.7, activated carbons would not suit our needs and the copper chloride content may need to be reduced before disposal. Cyanide is also not complexed with any heavy metal in this process. Alkaline chlorination and ozonation require basic solutions and would be very costly to implement due to large volumes of waste water. While ozone is a strong oxidizing agent, it is non-selective and will oxidize any organic available. In order to implement such large changes in the pH for any of these methods, the use of chemicals that may be as strictly regulated as cyanide is required which would lead to more waste disposal problems. This clearly indicates that a safe alternative is necessary. The best method for the removal of HCN for this flavor recovery waste was 59 determined to be an air stripper. This device would allow cyanide concentrations to be reduced to within disposal limits without the use of other hazardous materials. The stripped cyanide could also be recovered using a solution of sodium hydroxide which would yield a smaller amount of more concentrated waste. Apparatus A miniature column with 15 bubble cap trays was used in a hood to perform the experiments. The total height of the column was 3 ft with a diameter of 2 inches. A 1000 ml round bottom flask with an air inlet line was used as the liquid reservoir. Figure 5-1 illustrates the entire setup. The cap of the column contained an inlet line for the feed as well as an outlet for the air. The air outlet led to a gas scrub bottle containing the IN sodium hydroxide for cyanide recovery. A peristaltic pump was used to control the liquid flow rate into the column and a flow meter controlled the air flow from the house source. The column was wrapped with heating tape and attached to a variable autotransforrner to maintain the temperature when necessary. Procedures Two conditions were chosen for the feed, room temperature and 50 ° C. This data would illustrate any temperature sensitivity of the cyanide removal. A simple procedure was followed for each experiment. The hydrolyzate used to recover the cherry flavoring from the pits was prepared (~ 1 liter) and filtered to remove all solid particles. During preparation, the solution is heated to 50 ° C and allowed to cool while the initial cyanide concentration was determined for the room temperature runs. Otherwise, the solution continued to be heated in a water bath at 50 ° C. The steps involved in measuring the HCN concentration include diluting 5 ml of the sample in 500 ml of reverse osmosis (R0) water. A 25 ml sample is required for 60 Sample Line EU .jlzll=fl '.__ll.__ll Peristaltic Pump —-— Gas Scrub Bottle r—N Air Column supply h—J Figure 5-1: Air Stripper Apparatus. 61 the colorimetric test kit manufactured by Hach (Model CYN-3). A series of reagents are added to the sample and, after 20 minutes, the color of the sample is matched to an untreated sample of the same solution. The actual concentration is then read from the color disc provided with the test kit. The column is primed with the feed solution prior to introducing air to the system. Once the feed has reached every stage, the flow rate is adjusted to 2.54 ml/min and the air turned on to 1.54 L/min. The column is allowed to operate for approximately 1.5 hours, making sure there are no flooding problems every 20-30 minutes. The first sample is then taken from the bottom stage (approximately 7 ml), and the cyanide test is repeated. The dilution rate is adjusted in order to obtain sample concentrations from 0-0.3 mg/L, the limits of the colorimetric disc. Since the minimum sample size for this particular test is 25 ml, it is always necessary to add approximately 18 ml R0 water. This sample is mixed thoroughly with a magnetic stirrer. The cyanide test is repeated every 30 minutes to determine concentration with time and the time required to reach steady state. For room temperature distillations, the time required to reach steady state was greater than 3 hours versus approximately 3 hours for experiments performed at 50 ° C. The same procedure is repeated for the higher temperature runs with the start time being recorded after the feed and column liquid reach a temperature of 50 " C. The feed is kept heated in a water bath and the column with heating tape. Results The results indicate that this method of separation is feasible for the process. Table 3 gives the concentration of cyanide versus time for each of the experiments performed. It is clear that the stripper works better with higher feed temperatures. The 62 final test at 50 ° C shows that the concentration can be reduced to slightly higher than 50 ppb which is a significant reduction, but not quite low enough for direct disposal. The error associated with such colorimetric tests can be great due to the visual determination of a color match. More accurate tests would give a better determination of the final concentration remaining in the waste stream as well as more sophisticated distillation equipment and tighter temperature control. However, these results are sumcient for the feasibility stage of the work. Table 3: Cyanide removal results from Air Stripper. Experiment 1 (room T) 2 (50 C) 3 (room T) 4 (50 C) 5 (room T) 6 (50 C) Feed (ppm) 12 9 9 6 10 5-10 time(min) 95 86 95 120 120 130 HCN (ppm) 3 0.77 1.43 0 0.42 0.037 time(min) 1 15 125 152 165 HCN (ppm) 1.19 0.85 0.29 0.051 time(min) 188 201 HCN (ppm) 0.134 0.054 APPENDICES APPENDIX A Fluid Diameters The determination of the diameter of a water molecule as used in this thesis is based on the van der Waals volume and the Lennard-J ones molecular diameter. ' The vdW volumes for several non hydrogen bonding fluids are plotted against the L-J diameters to illustrate the linear relationship between the two. Table 4 contains these values as well as the ESD size parameter b. The L-J diameter of water does not have the same linear relationship as the rest of the molecules as seen in Figure A-l. This molecular diameter falls short of the line created by the other fluids. The molecular diameter used throughout this thesis is calculated assuming that water falls on the same line as the non hydrogen bonding fluids. Assuming that the vdW volume is accurate at 12.37 cm3/mol, the cube of the new molecular diameter is simply read from the plot as 3.4 Angstroms. The relationship between the ESD parameter b and the Lennard-Jones diameter is similar to the trends seen in Figure A-l. Polarizabilitr'es The fluid-solid interaction parameter eg/k for the water/carbon system is estimated in a similar way to the molecular diameter. The relationship between Sfi/k as determined by Subramanian (1995) and Tan et al., 1997 and polarizabilities is used to estimate the 63 64 interaction parameter for water based on other fluids. The interactions are determined from the fit to experimental data in work done by the previously mentioned authors. Table 5 contains the fluids used and the values of their parameters. Note that the value of 35 K for the interaction parameter for water is included in the table as an estimate used in this work. Figure A-2 illustrates the relationship between these properties. Values of 35 K and 40 K are both used in this thesis. Table 4: Parameters used for the determination of water diameter. 1D Component L-J o (A) 0301?) ESD b vdW volume (cm’lmol) (cm3/mol) 1 METHANE 3.758 53.07259551 10.863 17.05 2 ETHANE 4.443 87.70592631 16.716 27.34 3 PROPANE 5.118 134.060503 22.921 37.57 201 ETHYLENE 4.163 72.14715875 15.013 23.88 202 PROPYLENE 4.678 102.3718738 20.89 34.08 909 CARBON DIOXIDE 3.941 61.20956662 10.534 19.7 914 ARGON 3.542 44.43 709609 7.998 919 NEON 2.82 22.425768 4.346 920 KRYPTON 3 .655 48.82723638 9.896 1921 WATER 2.641 18.42066072 9.411 12.37 Table 5 : Fluid-solid interactions and polarizabilities. Component Name a *10025 (61113) Sa/k (K) H20 WATER 15.9 35 CH4 METHANE 26 107 C02 CARBON DIOXIDE 26.5 80 C2114 ETHYLENE 42.6 140 C3113 PROPANE 62.9 1 10 65 401 35¢ 30+ y = 0.3333): - 0.6899 25 i R2 = 0.988 20~ 15‘ vdw volume (cm 3) 10W 54 o T T T T I 1 1 0 20 4O 60 80 1 00 1 20 1 40 03 (A3) Figure A-l: Determination of water molecular diameter. 66 140 w 0 120 * 100 w 804» o Stalk (K) 60 e 40 «- 20 » 0 + i i i + t 0 10 20 30 40 50 POLARIZABILITY (1025 cm3) Figure A-2: Determination of 8ka for water/carbon system. 60 -L 70 APPENDIX B Calculating Y with geometric variations Rectangular Calculations: L/03= (H-cfi/oa z =distance from carbon surface For U05 >3 Z/O'g (0.5 OR Z/Ufl' >(L/O'g -0.5) YLCL=Yb*3./8.*(0.5+5./6.-1./(3.*(L/0'g -l.0)"'*3)) (B. 1) 0.1 S 2/0’55 1.5 YLCL= Yb*3./8.*(Z/C55 +5./6.-1./(3.*(L/O’g -Z/Og -0.5)**3)) (B.2) 1.532103 S(L/O'g -l.5) YLCL= Yb*3./8.*(8./3.-1./(3.*(Z/O'g- -0.5)**3)-1./(3."‘(L/Cg 4105 -0.5)**3)) (3.3) (L/O’fl' -1.5) S 2105 SW03 '05) YLCL= Yb*3./8.*(L/Og -Z/0'3' +5./6.-l ./(3."‘(z/O'g -0.5)**3)) (3.4) 2 S L/O’g <3 2105 <0.5 0R Z/O’g <(L/ca -0.5) YLCL= Yb*3./8.*(0.5+5./6.-1./(3.*(L/03-1.0)**3)) (B.5) 0.552103 S(L/O'g -1.5) YLCL= Yb*3./8.*(ZJGg +5./6.-1./(3.*(L/0’g ~z/O’g -0.5)**3)) (3.6) (Lloa -l.5) S 2163 $1.5 YLCL= Yb*3./8.*(L/O'g '1.) (3.7) 1552/65 S(L/O’g «0.5) YLCL= Yb*3./8.*(L/O'5 -Z/0'g +5./6.-1./(3.*(Z/O'fl' -O.5)**3)) (B.8) 1.55L/Cg S2 YLCL= Yb*3./8.*(L/Og -1.) (3.9) L/ca 51.5 YLCL= Yb*3./16. (3.10) 67 68 Spherical Calculations: Z =EI3*O’{, Rl=( of} £52+ZZ)/(ZZ) R/off $1.5 A0=8./3.+8./3. A1=2.*(Z-Rl)/Ug A2= 0.3/(a 0.2+ 22-2R1*Z)) A3=( 053)/(Z( 0.2+ 2320.2» A4=2653/(3(Z+ om A5=-8./3. R/O’ff > 1. 5 A0= 16/3 A 1 =-2/3 0.3/(24:93 A2=2/3 Gas/(2+ O'f,)3 A3= 0.3/(2% 0,2 412-220..» A4= (ma/(2% of} +zz+2Zofi)) A5=0 Relationship between 2 and r (3.11) (3.12) (3.13) (3.14) (3.15) (3.16) (3.17) (3.18) (3.19) (3.20) (3.21) (3.22) Equations B. l-B. 10 are used to calculate the value of Ym/Yb in rectangular coordinates. This value is then substituted into Equation 3.23 to find the value in radial coordinates. The resulting equations take the form of : YLCL(z,r) = 3/16);[Y£;(—Z) A0+(A1+ A2 + A3+ A4+ A5):l b (3.23) APPENDIX C Computer Codes: Main Program and Varying Subroutines SL113 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C THIS PROGRAM CALCULATES EXCESS IN A SLIT FOR MIXTURES USING C THE ESD EQUATION OF STATE. A SQUARE UELL POTENTIAL HAS BEEN ADDED. C C 2 IS FROM CARBON SURFACE. THEREFORE: ZOSFF=0 AT THE SURFACE OF UALL. C ZLCL IS DISTANCE FROM CARBON CENTER IN FLUID-SOLID DIAMETERS. C ZLCL IS EQUIVALENT TO ETA IN PREVIOUSLY URITTEN PROGRAMS. IMPLICIT DOUBLE PRECISION(A-H1KaO-Z) CHARACTER*1 ANSUER CHARACTER*E TYPE CHARACTERKED NAME PARAMETER (NMX=ID) DIMENSION TC(NMX)1PC(NMX)1ACEN(NMX)aID(NMX) + aZFEED(NMX)1S(NMX)1NAME(NMX) COMMON/ESD/KCSTAR(NMX)aDH(NMX)1C(NMX)aQ(NMX)1VX(NMX)1 8 VLK(NMX)1ND(NMX)1EOKP(NMX) COMMON/ETA/ETALaETAV1ZL1ZV COMMON/EQN/TYPEaAPPX COMMON/CONSTK/KIJ(NMXaNMX)aINITIAL COMMON/FEED/ZFEED COMMON/KVALUES/S COMMON/NAME/NAME COMMON/LOCAL/DELTAZ COMMON FKEEP DIMENSION LCLDENS(EUD)aFUGBULK(NMX)1PSI(NMX)1FUGLCL(NMX)1 lFUG(NMX)aFUGCCALC(NMX)aXA(NMX)aSIGMA(NMX)1FKEEP(NMX) COMMON/POSITION/ZLCL DOUBLE PRECISION LCLDENS:EXCEaLOSFFaRHODaFKEEP DOUBLE PRECISION FUGaRHOaPBaSIGFFsFUGLCLalOSFF COMMON/INTERACTIONS/SIGFFa SIGFS:SIGSS DIMENSION X(NMX)1IER(IE)1Y(NMX)aPSIL(NMX)1PSIB(NMX)1RHOD(EUU) OPEN(UNIT=521FILE='SLTOUTPUT.DAT') OPEN(UNIT=53aFILE='SLTPROFILE-DAT') OPEN(UNIT=SH«FILE='IMPTDATAI.DAT') OPEN(UNIT=SS~FILE='IMPTDATAE.DAT') OPEN(UNIT=551FILE='IMPTDATAB.DAT') URITE(53121) 69 El EB 29 1100 ID 15 30 SOL 505 501 70 F0R"AT(6X1'RHO'1ISX1'RHOB'1LSX1'PB'wllX1'pLCL' URITE(53a*)' ' URITE(59122) FORMAT(6X:'PLCL'1LSXa'ZLCL'allXa'ZOSFF'1LLX1'PSI') URITE(591*>' ' URITE(SL~EH) FORMAT(BXa'EPSFS'alsxa'XA'alsxa'H'ulSXa'ZOSFF'alSXa'PSI') URITE(551*)' ' URITE(SE:*)' ZOSFF'a' FUGLCL'a' RHO' URITE(ba*)'ENTER NUMBER OF COMPONENTS (0 TO QUIT) ' READ(Sa*)NC INITIAL=D IF (NC.EQ.D)STOP URITE(E~*)'ENTER ID NOS- OF THE COMPONENTS OR ZEROS FOR LIST ' READ(51¥)(ID(I)1I=14NC) CALL GETCRIT(NCaIDaNAMEaTCaPCaACEN) IF(ID(l)-EQ-U)GO TO 1100 CALL GETESD(NC:IDaEOKPaKCSTARaDHaC1Q1VX1ND) URITE(*1*)'THE Kij MATRIX IS ESTIMATED AS' URITE(*1'(BX13IB)')(ID(I)1I=I1NC) DO 15 I=21NC DO l0 J=laI-l CALL ESTACT(IDaIaJ1EOKP1VXaKIJ) CONTINUE URITE(L~LDI)IaID(I)1(KIJ(IaJ)ad=laI-l) CONTINUE CONTINUE URITE(E:*)'ENTER TYPE OF PHASE EQUILIBRIUM CALCULATION' URITE(B~*)'1. BP FOR BUBBLE POINT PRESSURE' URITE(E:*)'E. QT TO QUIT' URITE(ba*)' ' READ(SaSDS)TYPE FORMAT