MODELING AND SIMULATION OF SOLID-LIQUID EQUILIBRIUM BY PERTURBED-CHAIN STATISTICAL ASSOCIATING FLUID THEORY THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF TECHNOLOGY In CHEMICAL ENGINEERING By SUNIL KUMAR MAITY UNDER THE GUIDANCE OF
Prof. Saibal Ganguly & Prof. Sirshendu De
DEPARTMENT OF CHEMICAL ENGINEERING INDIAN INSTITUTE OF TECHNOLOGY, KHARAGPUR INDIA January 2003
Indian Institute of Technology Kharagpur
CERTIFICATE
India
________________________________________________________________________
This is to certify that SUNIL KUMAR MAITY, final year M.Tech student in the Department of Chemical Engineering, Indian Institute of Technology, Kharagpur has carried out his project work under our guidance during the academic session 2002-03, and herewith submitted the thesis entitled “MODELING AND SIMULATION OF SOLID-LIQUID EQUILIBRIUM BY PERTURBED-CHAIN STATISTICAL ASSOCIATING FLUID THEORY”, in partial fulfillment of the requirements for the award of the degree of “MASTER OF TECHNOLOGY” in “CHEMICAL ENGINEERING”.
Prof. S. Ganguly
Prof. S. De
DEPARTMENT OF CHEMICAL ENGINEERING INDIAN INSTITUTE OF TECHNOLOGY KHARAGPUR, INDIA
ACKNOWLEDGEMENTS
I express my deep sense of gratitude and indebtedness to Dr. S. GANGULY and Dr. S. DE, Associate Professor, Department of Chemical Engineering for their invaluable guidance, inspiration, and encouragement given to me at all stages of my M.Tech thesis work. I am also indebted to Mr. K. Gayen, Mr. M. K. Purakait, and Mr. D. P. Chakraborty, for their encouragement and help given to me carrying out the study for the dissertation work in their company. I also like to thank Mr. B. De and Mr. S. Bayen for their constant support to me at all stages of my work. Finally, I thank one and all who directly or indirectly rendered their help for the successful completion of my project work.
I.I.T., Kharagpur 13/01/2003
(SUNIL KUMAR MAITY)
CONTENTS Page Chapter-1
INTRODUCTION
1
Chapter-2
LITERATURE REVIEW
4
Chapter-3 3.1
MODELING OF SOLID-LIQUID EQUILIBRIUM Theory of Solid-Liquid Equilibrium
9 10
3.2
Solid-Liquid Equilibrium Phase Diagram
13
Chapter-4 4.1
PERTURBED-CHAIN STATISTICAL ASSOCIATING FLUID THEORY EQUATION OF STATE PC-SAFT Equation of State
14 15
4.2
Summary of Equation for Calculating Thermo-Physical Properties Using Perturbed Chain SAFT Equation of State
25
Chapter-5
REGRESSION ANALYSIS OF SOLUBILITY DATA
34
5.1
Solid-Liquid Equilibrium of n-Alkanes
35
5.2
Solid-Liquid Equilibrium of Aromatic Compounds
40
5.3
Effect of Pressure on Solid-Liquid Equilibrium
42
5.4
Effect of Solvent on Solid-Liquid Equilibrium
44
5.5
Effect of Molecular Weight and Melting Temperature on Solid-Liquid Equilibrium Discussion
44
46
6.1
SENSITIVITY STUDY SYSTEM Results of Sensitivity Study
6.2
Discussion
5.6 Chapter-6
FOR
POLYETHYLENE
45
47 50
Chapter-7
RESULTS OF SOLUBILITY OF POLYETHYLENE
51
7.1
Experimental Determination of Solubility
53
7.2
Conclusions and Future Scope of Work
55
References
57
Nomenclature
58
Appendix A
Derivation of the Pure-Solute F02l/F02s
61
Appendix B
Program of Solid-Liquid Equilibrium Calculation
64
List of Tables Sl no.
Name of Tables
Page
4.1.1
Model constants for the integrals I1( ,m) and I2( ,m) of square-well chains used in 4.1.20 and 4.1.21.
21
5.1
PC-SAFT Parameters of Organic Solutes and Solvents
35
5.1.1 5.1.2 5.1.3
5.1.4 5 .1.5 5.1.6
Experimental SLE Data for System n-Dodecane and nHeptane Experimental SLE Data for System n-Hexadecane and nHeptane Experimental SLE Data for System n-Octadecane and nHeptane Solid-Solid Transition and Melting Properties of Aromatic Compounds Solid-Solid Transition and Melting Properties of Normal Alkanes Experimental SLE Data for System n-Dotriacontane and nHeptane
36 37 37 38 38 39
5.2.1
Experimental SLE Data for System Biphenyl and Benzene
40
5.2.2
Experimental SLE Data for System ε-Caprolactone and Toluene
41
5.3.1
Experimental Data for System n-Octacosane and Decane
42
5.3.2 5.3.3
6.1
Correlation of Molar Volume (Cm3/Mol) and Temperature for n-Octacosane Experimental Data for System n-Octacosane, P-Xylene, and n-Decane PC-SAFT Parameters of Polyethylene
43 43
47
7.1
Experimental SLE Data for System Polyethylene and mXylene
52
7.1.1
Properties of Polyethylene
53
7.1.2
Experimental SLE Data for Grade1 Polyethylene in Xylene Experimental SLE Data for Grade2 Polyethylene in Xylene
7.1.3 7.1.4
PC-SAFT Parameters of Xylenes
53 53 55
List of Figures Sl no. 2.1 2.2
List of figures Schematic Diagram of PVT Cell Apparatus The Schematic Diagram of the High-Pressure Optical Vessel with a Video Microscope
3.2.1
Solid Liquid Equilibrium Phase Diagram
5.1.1
SLE for System n-Dodecane + n-Heptane
5.1.2 5.1.3 5.1.4 5.2.1 5.2.2
SLE for System n-Hexadecane and n-Heptane SLE for System n-Octadecane and n-Heptane SLE for System n-Dotriacontane and n-Heptane SLE for System Biphenyl and Benzene
5.3.1 5.3.2 5.4.1 5.5.1
6.1.1 6.1.2 6.1.3 6.1.4 6.1.5
SLE for System ε-Caprolactone fnd Toluene Effect of Pressure on Binary SLE of System n-Octacosane and n-Decane for Different Composition. Effect of Pressure on SLE of System n-Decane+P-Xylene + n- Octacosane Effect of Solvent on SLE Effect of Molecular Weight on SLE in n-Heptane Effect of Pressure on Solubility of Polyethylene in mXylene Effect of Crystallizability Fraction on Solubility of Polyethylene in M-Xylene Effect of Melting Point on Solubility of Polyethylene in mXylene Effect of Solvent on Solubility of Polyethylene in mXylene Effect of Kij on Solubility of Polyethylene in m-Xylene
Page 6 7
13 36 37 38 39 40 41 42 43 44 44
47 48 48 49 49
7.1 7.1.1 7.1.2
Solubility of Polyethylene in m-Xylene at One Bar Pressure and Prediction By PC-SAFT Model. Solubility of Polyethylene (PE30398.9) in Xylene Solubility of Polyethylene (PE32599.8) in Xylene
52 54 55
CHAPTER 1
INTRODUCTION
Chapter 1: Introduction The study of solid–liquid equilibrium, SLE is of great technical interest for developing and designing separation processes, such as crystallization and fractionation. Crystallization processes are used for separation of mixtures. Knowledge of solid-liquid equilibrium behavior is also important for pipeline design where undesirable crystallization can cause safety problem. In oil production, solubility of normal-alkanes as well as other materials such as aromatics and naphthalene are important. But most of the studies deal with solubility mostly at atmospheric pressure. In the production of crude oil, pressure is generally elevated and precipitation can be a troublesome problem. So the study of solubility at elevated pressure has huge industrial importance. Linear polyethylene is composed of a distribution of n-alkanes of different molecular weight. Lower molecular weight oligomer fractions have reasonably high solubility in various solvents. The study of the solubility of low molecular weight alkanes in various solvents and in polymer themselves will lead to their partition coefficients when a polymer is in contact with any of the solvents. Partition coefficients are useful in predicting the maximum or equilibrium levels of migration of these low molecular weight components from the polymer into contacting solvents when estimating the migration from food packaging materials into food. Polyethylene coming out of the reactor is separated from the solvents in a flush drum. Then cooling and crystallization is used to purify it. So crystallization on the surface of heat exchanger and flush drum and clogging of pipeline due to crystallization are the typical industrial problem. In the polyethylene production, reactor is operated at very high pressure. So the study of solubility at high pressure is industrially very important. The solubility data for low molecular weight n-alkanes and aromatic compounds is enormous in this regard. Using these data parameters of different equation was determined and hence solubility was predicted. But solubility data for polymers such as polyethylene is scarce even in atmospheric pressure. Again considerable experimental effort is generally required to study the high-pressure phase equilibrium for polymer system. The equations that are used for predicting solubility of low molecular weight aromatic compounds are not always applicable to highly non-ideal mixtures such as high molecular weight chain like polymers.
Modeling and Simulation of Solid-Liquid Equilibrium
2
Chapter 1: Introduction
In this work, PC-SAFT equation of state is used to model solid- liquid equilibrium since it has wide applicability starting from low molecular weight organic compounds to high molecular weight polymer system, highly non-ideal system to associating compounds. PC-SAFT equation of state is used for predicting thermo physical properties and liquidliquid, vapor-liquid equilibrium since last few years. This equation of state requires three pure component parameters: segment no (m), segment diameter (σ), and energy parameter (ε/K), and it has one adjustable solvent-solute binary interaction parameter. The simplest case of SLE is that of a pure crystalline totally crystallizable solute and liquid Solvent, where the solute has a finite solubility in the solvent, but the solvent solubility in the solid is zero (X2S=0). In case of polyethylene, it is not a totally crystalline solute. Crystallinity fraction of polyethylene varies from 0.4 to 0.6. This model is applicable to both totally crystalline to partially crystalline solutes. Here a model has been developed based on PC-SAFT equation of state, which is applicable to homopolymer system. This model requires melting point, crystallinity fraction, no of repeating unit data of polymer. In this work this model is initially tested with literature solubility data of low molecular weight n-alkanes and aromatic compound both at atmospheric and elevated pressure for different solvent systems. Then sensitivity study is done for polyethylene system to understand the effects of different parameters on solubility. Lastly solubility study is done for different grades of polyethylene in xylene at atmospheric pressure. This data is used to determine the model parameters such as crystallinity fraction and adjustable solvent-solute binary interaction parameter (Kij).
Modeling and Simulation of Solid-Liquid Equilibrium
3
CHAPTER 2
LITERATURE REVIEW
Chapter 2: Literature Review This chapter deals with the solubility study by previous workers and their findings, experimental set up, and modeling. These are discussed one by one in subsequent paragraphs. • Cheng Pan, Maciej Radosj (3) developed a solid-liquid equilibrium (SLE) model based on copolymer SAFT (Statistical Associating Fluid Theory). Copolymer SAFT was derived from well known homopolymer version of SAFT. This equation of state is applicable to heterosegmented chains; chains composed of segments varying in size, energy, and hence connected with different kinds of bonds. Copolymer-SAFT is used to calculate fugacity coefficients of solutes in the liquid mixture. This equation of state is applicable to totally crystalline to partially crystalline solutes. Initially this model was regressed and tested on solubility data for naphthalene, n-alkanes, and polyethylene. Then the model was used in sensitivity study to understand the effects of crystallizability, melting point, molecular weight, and pressure on SLE of polyethylene in supercritical and sub critical propane. The results of their sensitivity analysis are as follows: With increase in pressure, solubility decreases at fixed temperature. With increase in crystallinity fraction (C), solubility decreases at fixed temperature and pressure. With increase in melting point, solubility decreases at fixed temperature and pressure. • Shu-Sing Chang, John R. Maurey, and Walter J. Pummer (7) determined the solubility and phase equilibrium of two n-alkanes namely n-Octadecane (C18) and ndotriacontane (C32) in different solvent systems: n-heptane, ethanol, ethanol/water mixture, tributyrin, trioctanoin, and mixed triglycerides. Solubility was determined by visual observing the dissolution temperature of a mixture of solvent and solute of known composition. Magnetic stirrer enhanced the mixing. For lower solubility where visual method become impractical, 14C labeled tracers were used. With a sensitive liquid scintillation counter operating at 20-30 cpm background, it was possible to detect the presence of 10-10 g of the labeled alkanes in aliquots taken from the solution. Differential Scanning Calorimeter (DSC) measured heat of fusion and melting temperature of the two n-alkanes, required for modeling of SLE.
Modeling and Simulation of Solid-Liquid Equilibrium
5
Chapter 2: Literature Review • Hyo-Guk Lee, Frank R. Groves, and Joanne M. Walcott (6) measured the effect of pressure on binary (SLE) for system: n-decane + n-octacosane (C28) and ternary SLE for systems: n-decane + p-xylene + n-octacosane (C28), and n-decane + p-xylene + phenanthrene mixtures. Their measurements correspond to 10 mole% solid content and pressure up to 200 bars.
Fig 2.1: Schematic diagram PVT cell apparatus. A, Ruska pump; B, pressure gauge; C, PVT cell; D, air bath; E, cathetometer; F, mercury reservoir; G, CO2 reservoir; H, Flash separator; I, wet test meter; TC, temperature controller; TI, temperature indicator. Pressure effect was measured by a Ruska Pressure-Volume-Temperature (PVT) cell, in which equilibrium condition was observed visually through the sight glass of the cell as shown in fig 2.1. The results of their study showed that solubility of n-octacosane in n-decane is decreased by around 40% as pressure is increased from atmospheric to 200 bars. The results were correlated by calculating activity coefficient of solute in liquid mixture from Flory-Huggins plus regular solution equation including a pressure correction term.
Modeling and Simulation of Solid-Liquid Equilibrium
6
Chapter 2: Literature Review •
Roland Witting, Dana Constantineseu, and Jurgen Gmchling
(5)
measured the
solubility of ε-caprolactone by visual technique in the following solvent systems: benzene, toluene, cyclohexane, 1-propanol, methanol, water, and 2-pentanone. For the description of the activity coefficient the NRTL model was used. Using the experimental solubility data the interaction parameters of the NRTL model were evaluated. Also activity coefficient of the solutes were predicted using the group contribution method modified UNIFAC (Dortmund) using the available group interaction parameter of the “acyclic esters” group (COO). • Hyo-Guk Lee, Philip A. Schenewerk, and Joanne Walcott, and Frank R. Groves (12) Jrs studied the pressure effect on solubility for n-octacosane in a mixture of n-decane and carbon dioxide. Their measurement was based on PVT cell apparatus as described earlier. Using perturbed hard sphere chain equation of state with empirical mixing rule then they correlated the solubility data. • Y. Tanaka and M. Kawakami (11) measured high-pressure solid-liquid binary phase equilibrium for the following four systems: benzene +n-tetradecane, benzene +nhexadecane, cyclohexane + n-tetradecane, cyclohexane + n-hexadecane. They measured saturation condition using high pressure optical vessel with the aid of a video microscope as shown below in Fig 2.2
Fig 2.2: The schematic diagram of the high-pressure optical vessel with a video microscope. A, high-pressure optical vessel; B, CCD camera probe; C, color monitor; D, thermocouple; E, thermostat; F, pressure transducer; G, pressure and temperature indicator; H, sample inlet; I, sample outlet; J, oil pump; K, pressure intensifier; L, exchange of oil path.
Modeling and Simulation of Solid-Liquid Equilibrium
7
Chapter 2: Literature Review
• In the late sixties, E. McLaughlin and H.A. Zainal (8,9) studied the solubility of a large no aromatic compounds such as biphenyl, o-, m-, p-terphenyl, naphthalene, anthracene etc in benzene and carbon tetrachloride solvent. • Joachim Gross and Gabriele Sadowski (4) developed Perturbed-Chain Statistical Associating Fluid Theory (PC-SAFT) equation of state based on modified square well potential. It requires three pure components parameters: segment no., segment diameter, energy parameter. Also it uses one adjustable binary interaction parameter. It is widely applicable to non-spherical chain like molecules like polymer. This equation of state has capability of predicting pure and mixture density; excess enthalpy, entropy, free energy of mixtures; vapor-liquid, liquid-liquid and solid-liquid equilibrium; and vapor pressure.
Modeling and Simulation of Solid-Liquid Equilibrium
8
CHAPTER 3
MODELING OF SOLID-LIQUID EQUILIBRIUM
Chapter 3: Modeling of Solid-Liquid Equilibrium 3.1 Theory of Solid-Liquid Equilibrium At equilibrium, the crystalline-solute fugacity in the liquid phase is equal to that in the solid phase: fL2=fS2
(3.1.1)
Where 2 means solute. Further, the solute fugacities in both liquid and solid phases are: fL2=
L L 2x 2P
fS2=fS02
(3.1.2) (3.1.3)
Where, as usual, the solid phase is assumed to be pure crystalline solute. Substituting Eq. 3.1.2 and Eq. 3.1.3 into Eq. 3.1.1, we have L L S 2x 2P=f 02
(3.1.4)
Let us divide both sides of Eq. 3.1.4 by f02L, which is the fugacity of pure sub cooled liquid solute at constant T and P: L L L S L 2x 2P/f 02=f 02/f 02
(3.1.5)
Since f02L is the fugacity of pure liquid, x2L = 1, so f02L can be expressed as fL02= 0P Where
0
(3.1.6)
is the fugacity coefficient of pure sub-cooled liquid solute at constant T and P.
Next, let us substitute Eq. 3.1.6 into Eq. 3.1.5 and take the natural logarithm of both sides of Eq. 3.1.5: (3.1.7)
The fugacity ratio of pure solute on the right hand side of Eq. 3.1.7, derived from a thermodynamic cycle given in Appendices “A” is as follows:
Modeling and Simulation of Solid-Liquid Equilibrium
10
Chapter 3: Modeling of Solid-Liquid Equilibrium
(3.1.8)
Where Psat is the solute saturated-vapor pressure at its melting temperature. volume difference of liquid and solid solute defined as v=vL-vS.
v is the
The first three terms on the right-hand side of Eq. 3.1.8 are not of equal importance, as suggested by Prausnitz et al (13); the first term is dominant. The other two terms have opposite signs, and hence, have a tendency approximately to cancel each other. These two heat-capacity terms, therefore, are neglected. The last term in Eq. 3.1.8, accounts for the pressure effect. The pressure effects tend to be negligible at low pressures. At higher pressures, however, and in general for compressible solutions, the pressure effects can be significant. So, after neglecting the heat capacity effect, Eq. 3.1.8 becomes
(3.1.9) Substituting Eq. 3.1.9 into Eq. 3.1.7, we get our working equation: (3.1.10) For low-pressure systems, the second term in the right side of Eq. 3.1.10 vanishes. Eq. 3.1.10 is applicable to crystalline solutes, that is solutes with 100% crystallizability, such as pure normal alkanes and aromatic hydrocarbons. We extend this approach to crystallizable polymers, that is, macromolecular solutes with partial crystallizability that are usually referred to as semi-crystalline polymers; the amorphous polymers will have 0% crystallizability. We quantify crystallizability in terms of a crystallinity fraction c, where 0 c 1, assuming that the polymer contains c crystalline fraction and (1-c) amorphous fraction. Following Harismiadis and Tassios (14), who assume that the log of the ratio of fugacities is proportional to c, we express the effective log of the ratio of fugacities as follows:
Modeling and Simulation of Solid-Liquid Equilibrium
11
Chapter 3: Modeling of Solid-Liquid Equilibrium
(3.1.11)
Where u is the number of ethyl units in the backbone. The fugacity ratio of the crystal unit for the sub-cooled liquid and solid can be estimated as for the crystalline molecules, using Eq. 3.1.9, except the enthalpy of melting, Hm, is exchanged for Hu, and Psat is set equal to zero because Psat is very low for polymers:
(3.1.12)
Where Hu is the enthalpy of melting per mole of crystal unit. For the ethyl unit, Hu=8.22 kJ/mol, as reported by Van Krevelen (15). The polymer-volume change, v, is determined from the densities of an amorphous polymer, a, and a crystalline polymer, c: v=1/ a-1/ c. For polyethylene, a=0.853 g/cm3 and c=1.004 g/cm3. Combining Eq. 3.1.7, Eq. 3.1.11 and Eq. 3.1.12, we get (3.1.13) Where subscript p stands for polymer. The fugacity coefficients of polymer in solution, p
L
and of pure-liquid polymer,
p
0
are calculated by PC-SAFT equation of state.
Solid crystalline normal alkanes such as n-dotriacontane exhibits a solid-solid (ss) phase transition a few degrees below its melting point. Such phase transition is typical of nalkanes larger than n-C20; from a higher-temperature phase of a hexagonal geometry to more stable crystalline structures at lower temperatures. Furthermore, the two solid phases are in a state of thermodynamic equilibrium. Similar to the above approach, we include the effect of the ss phase transition as follows:
(3.1.14)
Modeling and Simulation of Solid-Liquid Equilibrium
12
Chapter 3: Modeling of Solid-Liquid Equilibrium Where Tss is the ss transition temperature and
Hss is the enthalpy of the ss transition.
3.2 Solid-Liquid Equilibrium Phase Diagram SLE data in this work are presented in the form of temperature-solubility, T–x, phase diagrams. The solubility is defined as the equilibrium mole fraction (for small molecules) or weight fraction (for macromolecules) of the solid solute in solution.
C L T A LS
Wt or mole fraction (x) Fig 3.2.1: solid liquid equilibrium phase diagram Continuous SL curve (AC) divides T-X phase plain into two regions, a liquid region (L) at higher temperature and a solid-liquid (SL) region at lower temperature. In SL region two phases coexists, the crystalline solid and the liquid solution.
Modeling and Simulation of Solid-Liquid Equilibrium
13
CHAPTER 4
PERTURBED-CHAIN STATISTICAL ASSOCIATING FLUID THEORY EQUATION OF STATE
Chapter 4: PC-SAFT Equation of State
4.1 PC-SAFT EQUATION OF STATE Introduction For the correlation and prediction of phase equilibrium in macromolecular systems, the equations of state for chain molecules have been successfully used for more than two decades. In many recent investigations, non-spherical molecules are conceived to be chains comprised of freely jointed spherical segments. Several routes have been established to obtain descriptions for those chain fluids. One particularly successful equation of state concept for chain molecules is based on Wertheim's theory of associating fluids. Applying Wertheim's first-order perturbation theory (TPT1), Chapman et al. derived an equation of state for chain mixtures, known as the statistical associating fluid theory (SAFT). Initially the chain structure was not accounted for in the dispersion term of the SAFT equation, since a hard-sphere reference was used within the chain term; the dispersion contribution of each segment in a chain was assumed to be equal to a nonbonded spherical molecule of the same diameter. Numerous investigators have subsequently examined the use of a square-well reference and a Lennard¯Jones reference fluid in the chain-term, leading to equations of state for square-well chains and Lennard¯Jones chains, respectively. These expressions are lengthy, and thus many of the most commonly applied engineering equations of state still utilize square-well dispersion terms, which do not account for the connectivity of the segments. PC-SAFT equation uses the same chain term and association term as the earlier SAFT equations. Because a hard chain fluid serves as a reference for perturbation theory, rather than the spherical molecules as in the SAFT modifications, the proposed model is referred to as perturbed-chain SAFT (PC-SAFT). This model is applicable to real chain molecules of any length, from spheres to polymers. Molecular Model: Modified Square Well Potential In the proposed equation molecules are conceived to be chains composed of spherical segments. Pair potential for the segments of a chain is given by modified square well potential as suggested by Chen and Kreglewski. ∝ r < (σ -s1) U(r) =
3ε
(σ -s1) ≤ r < σ
-ε
σ ≤ r < λσ
0.0
(4.1.1)
r ≥ λσ
Modeling and Simulation of Solid-Liquid Equilibrium
15
Chapter 4: PC-SAFT Equation of State Where U(r) is pair potential, r is radial distance between two segments, σ is temp independent segment diameter, ε is the depth of potential well, λ is reduced well width. As suggested by Chen and Kreglewski a ratio of s1 /σ = 0.12 is assumed. Any specific interactions, like hydrogen bonding or dipole-dipole forces have been neglected. Contributions to the Helmholtz free energy due to such interactions may be implemented separately. According to this model, nonassociating molecules are characterized by three pure component parameters: the temp independent segment diameter, σ; the depth of potential, ε; and the no of segments per chain, m. PC-SAFT Equation of state The properties of the square-well chain fluid are calculated applying a perturbation theory; where the structure of the fluid system is assumed dominate by the repulsive interactions. According to the perturbation theories, the interaction of molecules can be divided into a repulsive part and a contribution due to the attractive part of the potential. To calculate the repulsive contribution, a reference fluid in which no attractions are present is defined. Attractive interactions are treated as a perturbation to the reference system. In the frame work of Barker and Henderson’s perturbation theory a temperature dependent segment diameter, d(T) can be used to describe the soft repulsion of the molecules, where
d(T) = ∫ 0σ 1 – exp -U(r) /κT
dr
(4.1.2)
For the potential function used in Eq. 4.1.1, integration leads to the temperature dependent hard segment diameter, di (T) of component i, according to
ε ⎡ ⎛ d i = σ i ⎢1 − 0.12 exp⎜ − 3 i kT ⎝ ⎣
⎞⎤ ⎟⎥ ⎠⎦
(4.1.3)
The complete equation of state is given as an ideal gas contribution (id), a hard chain contribution (hc), and a perturbation contribution, which accounts for the attractive interactions (disp).
Modeling and Simulation of Solid-Liquid Equilibrium
16
Chapter 4: PC-SAFT Equation of State Z = Zid + Zhc + Zdisp
(4.1.4)
Where Z is the compressibility factor, with Z=PV/(RT) and Zid =1.0, P is the pressure, V is the molar volume, R denotes the gas constant. Hard Chain Reference Equation of State Based on Wertheim’s first-order thermodynamic perturbation theory Chapman et al. developed an equation of state, which for homonuclear hard-sphere chains is given by
(4.1.5) (4.1.6) Where xi is the mole fraction of chains of component i, mi is the number of segments in a chain of component i, is the total number density of molecules, giihs is the radial pair distribution function for segments of component i in the hard sphere system, and superscript hs indicates quantities of the hard-sphere system. Expressions of Boublik and Mansoori et al. are used in Eq. 4.1.3 for mixtures of the hard-sphere reference system:
(4.1.7) (4.1.8) Where
(4.1.9) And
i
is the segment diameter of component i.
Note that 3 is equal to the packing fraction , i.e., represents a reduced segment density.
3
= . The packing fraction
Perturbation Theory for Pure Chain Molecules The basic idea of the perturbation approach is that the Helmholtz free energy of a system can be expressed as an expansion in inverse temperature around the free energy of a Modeling and Simulation of Solid-Liquid Equilibrium
17
Chapter 4: PC-SAFT Equation of State reference system. The perturbation expansion is fast convergent and can be truncated after the second term, so that the perturbation contribution to the Helmholtz free energy of the system is given by
(4.1.10) Where A1 and A2 are the first-and second-order perturbation terms. The perturbation theory of Barker and Henderson was derived for spherical molecules. If square-well chains are to be treated within this theory, all intermolecular segmentsegment interactions between two molecules have to be considered. The appropriate equations become as
(4.1.11)
(4.1.12) With
(4.1.13) Where gαβhc (m; x,
) is the site-site radial distribution function of chains, which
represents the radial distribution function for a segment
of one chain and a segment
of another chain separated by the radial distance χαβ= x. In Eq. 4.1.11 and Eq. 4.1.12 homonuclear chains are assumed, where any two segments
and
on different chains
interact with the same depth of the pair potential, αβ = , and same well width, αβ = . Note, that in Eq. 4.1.11 and Eq. 4.1.12 already and N were introduced on a permolecule basis, i.e. the relations s= ·m and Ns=Nm were applied, where s is the segment density ( s = Ns/ V), and Ns is the number of segments in the system. The perturbation expressions originally required segment-based quantities ( s and Ns); the conversion to molecular quantities ( and N) leads to a factor m on the right hand side of
Modeling and Simulation of Solid-Liquid Equilibrium
18
Chapter 4: PC-SAFT Equation of State Eq. 4.1.12. The superscript "hc" in Eq. 4.1.11, Eq. 4.1.12, and Eq. 4.1.13 indicates quantities of the hard-sphere chain reference fluid. Besides an expression for the compressibility factor Zhc of the reference hard-sphere chain fluid, the perturbation theory of Barker and Henderson requires the site-site radial distribution function gαβhc (m; x, ) of the reference fluid. Chiew has derived equations for gαβhc (m; x, ) of chains from integral equation theory by applying the Percus¯Yevick closure and obtained an approximation for the average intermolecular radial distribution function gαβhc (m; x, ), given by
(4.1.14) The position of segments and within the appropriate chains has considerable influence on the site¯site radial distribution function. It is important for example, whether refers to a terminal segment (segment located at the end of a chain) or to a non-terminal segment. By means of the averaging (Eq. 4.1.14), the segments of chain molecules are non-distinguishable. They are characterized by an average intermolecular radial distribution function, which is also given on a segment¯segment basis. In the formulation of the perturbation terms given above in Eq. 4.1.11 and Eq. 4.1.12, all segment¯segment interactions have to be considered individually. However, it is fair to introduce an averaging into the perturbation theory analogous to that proposed by Chiew, i.e., Eq 4.1.14. For a pure fluid of homonuclear chains, the perturbation terms become
(4.1.15)
(4.1.16) Where for pure chain fluids the compressibility term can be obtained from Eq. 4.1.5 in the form
Modeling and Simulation of Solid-Liquid Equilibrium
19
Chapter 4: PC-SAFT Equation of State
(4.1.17) A Simple Mathematical Representation of the Perturbation Terms Since the integrations over reduced radius x in Eq. 4.1.15 and Eq. 4.1.16 have to be performed numerically, it is desirable to find a simple mathematical representation for those integrals. Let us therefore introduce the abbreviations
(4.1.18)
(4.1.19) Although the average radial distribution function depends upon radius, density and segment number, the integrations over radius yields expressions I1( , m) and I2( , m) which are functions of density and segment number only. Gulati and Hall have taken advantage of that fact by representing I1( , m) as a simple power series in density for the case of a spherical square-well fluid (m=1) and of a square-well dimer fluid (m=2). Those authors have obtained values of the radial pair distribution function of monomers and dimers from molecular dynamics simulations and upon integration obtained I1( ,m=1) and I1( , m=2)-values. They have subsequently fit coefficients of a power series to both sets of I1-values over a range of densities (0.025 0.475) in order to obtain simple expressions for I1( , m=1) and I1( , m=2). Hino and Prausnitz have also used a power series in density to substitute I1( ,m) for the case of a spherical square-well fluid (m=1). They have obtained the power series coefficients by fitting an analytic expression for the integral in Eq. 4.1.18, which was derived by Chang and Sandler. In the present study, it is aimed at a simple function which can represent I1( , m) in Eq. 4.1.19 for the entire range of segment numbers. As in the previous investigations, a power series in density is assumed, where the appropriate coefficients are now functions of the segment number. The integrals I1( , m) and I2( , m) can accurately be represented by a power series in density of sixth order:
Modeling and Simulation of Solid-Liquid Equilibrium
20
Chapter 4: PC-SAFT Equation of State
(4.1.20) Where ai(m) and bi(m) are coefficients of the power series in density. From Eq. 4.1.19, it becomes apparent that bi=(i+1)ai. However, this simple relation does not hold if more realistic pair potentials would be adopted. Let us now concern with the dependence of the power series coefficients ai(m) and bi(m) on the segment number. Only the function ai(m) will be considered in the following paragraph; the function bi(m) can be treated analogously. First coefficients are regressed ( ai(m)) of the power series i6ai· i to the integral in Eq. 4.1.18 for different segment numbers (m=1, 1.5, 2, 3, 4, 5, 6, 7, 8, 10, 100, 1000)). It was found, that the dependence of each of the power series coefficients on segment number can accurately be described with a relation proposed by Liu and Hu.
(4.1.21) Table 4.1.1. Model constants for the integrals I1( , m) and I2( , m) of square-well chains used in Eq. 4.1.20 and Eq. 4.1.21. i 0 1 2 3 4 5 6
a0i 0.91056314452 0.63612814495 2.68613478914 -26.5473624915 97.7592087835 -159.591540866 91.2977740839
a1i -0.30840169183 0.18605311592 -2.50300472587 21.4197936297 -65.2558853304 83.3186804809 -33.7469229297
a2i -0.09061483510 0.45278428064 0.59627007280 -1.72418291312 -4.13021125312 13.7766318697 -8.67284703680
i 0 1 2 3 4 5 6
b0i 0.72409469413 2.23827918609 -4.00258494846 -21.0035768149 26.8556413627 206.551338407 -355.602356122
b1i -0.57554980753 0.69950955214 3.89256733895 -17.2154716478 192.672264465 -161.826461649 -165.207693456
b2i 0.09768831158 -0.25575749816 -9.15585615297 20.6420759744 -38.8044300521 93.6267740770 -29.6669055852
With Eq. (4.1.21), the model-constants a0i, a1i, and a2i (with i=0,...6) are introduced. The constants a0i can easily be obtained by setting a0i=ai(m=1).The model-constants a1i and a2i were obtained by fitting them to a matrix of I1(m, ) values, where ranges of =0,...0.6
Modeling and Simulation of Solid-Liquid Equilibrium
21
Chapter 4: PC-SAFT Equation of State and m=1,...1000 were chosen. All I1(m, ) values were calculated from Eq. 4.1.18 using the average radial distribution function proposed by Chiew. The values of the modelconstants a0i, a1i, and a2i as well as b0i, b1i, and b2i are given in Table 4.1.1. Using Eq. 4.1.20 and Eq. 4.1.21, the perturbation terms of first and second order can be rewritten to the simple form
(4.1.22) (4.1.23) The three pure-component parameters required by the equation of state are those which entirely characterize square-well chain molecules: the segment number, m; the segment diameter, ; and the depth of the pair potential, /k. Extension to the Mixtures A rigorous application of Barker and Henderson's perturbation theory to mixtures (within the above described formalism) requires expressions for the average radial pair distribution function gijhc(mk , k; xij, ) of mixtures. These must be given for any pair of molecules i and j in the system of molecules with segment numbers mk and segment diameter k of all k components. O'Lenick and Chiew have derived a set of equations for the radial pair distribution function of mixtures; however, those equations are not given in analytical form. Therefore, Van der Waals one fluid mixing rules are adopted here to extend the perturbation terms to mixtures
(4.1.24)
(4.1.25) Conventional combining rules are employed to determine the parameters between a pair of unlike segments
Modeling and Simulation of Solid-Liquid Equilibrium
22
Chapter 4: PC-SAFT Equation of State
(4.1.26) (4.1.27) The one-fluid mixing concept of the compressibility term in Eq. 4.1.25 were applied, i.e., similarly to Eq. 4.1.17 it is
(4.1.28) The approximation within the Van der Waals one-fluid mixing rule, that the radial pair distribution function can be averaged for reduced radii, is a widespread approach for mixtures of spherical molecules. Since the radial pair distribution function for chain molecules in the above formalism is given on a per-segment basis, the one-fluid mixing rule is also applicable to chain molecules. The terms I1( , ) and I2( , ) in Eq. 4.1.24 and Eq. 4.1.25 are then evaluated for the mean segment number of the mixture, which is given by Eq. 4.1.6. The equation of state can be written in terms of the compressibility factor applying the relation
(4.1.29) The compressibility factor is calculated as Z = Zhc + Zpert, where the perturbation contribution is given by Zpert=Z1+Z2
(4.1.30)
and the perturbation terms of first-and second-order are given by
(4.1.31) With
Modeling and Simulation of Solid-Liquid Equilibrium
23
Chapter 4: PC-SAFT Equation of State
(4.1.32) and
(4.1.33) Where
(4.1.34) and where C1 and C2 are abbreviations defined as
(4.1.35) (4.1.36)
Modeling and Simulation of Solid-Liquid Equilibrium
24
Chapter 4: PC-SAFT Equation of State
4.2 SUMMARY OF EQUATIONS FOR CALCULATING THERMOPHYSICAL PROPERTIES USING PERTURBED-CHAIN SAFT EQUATION OF STATE This section provides a summary of equations for calculating thermo-physical properties using the perturbed-chain SAFT equation of state. The Helmboltz free energy A res is the starting point in this paragraph, as all other properties can be obtained as derivatives of A res . In the following, a tilde (~ ) will be used for reduced quantities, and caret symbols (^ ) will indicate molar quantities. The reduced Helmholtz free energy, for example, is given by ~ res
a
=
Ares NkT
(4.2.1)
At the same time, one can write in terms of the molar quantity ^ res
a
~ res
a = RT
(4.2.2)
Helmholtz Free Energy The residual Helmholtz free energy consists of the hard-chain reference contribution and dispersion. ~ res
a
~ hc
~ disp
=a +a
(4.2.3)
Hard-Chain Reference Contribution ~ hc
a
~ ~ hs
= m a − ∑ x1 (mi − 1) ln g iihs (σ ii )
(4.2.4)
i
~
Where m is the mean segment number in the mixture _
m = ∑ xi mi
(4.2.5)
i
The Helmholtz free energy of the hard-sphere fluid is given on a per-segment basis
Modeling and Simulation of Solid-Liquid Equilibrium
25
Chapter 4: PC-SAFT Equation of State ~ hs
a
=
⎤ ⎛ ζ 23 ⎞ ζ 23 A hs 1 ⎡ 3ζ 1ζ 2 ⎜⎜ 2 − ζ 0 ⎟⎟ ln (1 − ζ 3 )⎥ = + + ⎢ 2 N s kT ζ 0 ⎢⎣ (1 − ζ 3 ) ζ 3 (1 − ζ 3 ) ⎥⎦ ⎝ζ 3 ⎠
(4.2.6)
and the radial distribution function of the hard-sphere fluid is
g
hs ij
⎛ did j 1 = +⎜ (1 − ζ 3 ) ⎜⎝ d i + d j
⎛ did j ⎞ 3ζ 2 ⎟ +⎜ 2 ⎜d +d ⎟ (1 − ζ ) j 3 ⎝ i ⎠
2
⎞ 2ζ 22 ⎟ ⎟ (1 − ζ )3 3 ⎠
(4.2.7)
With ζ n defined as
π
n ∈ {0 ,1 , 2 , 3 } 6 i The temperature-dependent segment diameter d i of component i is given by
ζn =
ρ ∑ xi mi d in
ε ⎡ ⎛ d i = σ i ⎢1 − 0.12 exp⎜ − 3 i kT ⎝ ⎣
⎞⎤ ⎟⎥ ⎠⎦
(4.2.8)
(4.2.9)
Dispersion Contribution The dispersion contribution to the Helmholtz free energy is given by ~ disp
a
( )
− ⎛ −⎞ = −2πρl1 ⎜η , m ⎟m 2 ∈ σ 3 − πρ m C1l 2 η , m m 2 ∈2 σ 3 ⎠ ⎝
(4.2.10)
Where an abbreviation C1 is introduced for the compressibility expression, which is defined as
⎛ ∂Z hc C1 = ⎜⎜1 + Z hc + p ∂p ⎝
⎞ ⎟⎟ ⎠
−1
− 8η − 2η 2 − ⎛ 20η − 27η 2 + 12η 3 − 2η 4 ⎞ ⎛ ⎞ ⎜ ⎟⎟ = ⎜1 + m + ⎜1 − m ⎟ 4 2 ( ) [ ( )( ) ] − − − 1 η 1 η 2 η ⎝ ⎠ ⎝ ⎠
−1
(4.2.11)
Another set of abbreviations
⎛ ∈ij m 2 ∈ σ 3 = ∑ ∑ x i x j m i m j ⎜⎜ i j ⎝ kT
⎞ 3 ⎟⎟σ ij ⎠
Modeling and Simulation of Solid-Liquid Equilibrium
(4.2.12)
26
Chapter 4: PC-SAFT Equation of State ⎛ ∈ij m ∈ σ = ∑ ∑ x i x j m i m j ⎜⎜ i j ⎝ kT 2
2
3
2
⎞ 3 ⎟⎟ σ ij ⎠
(4.2.13)
Conventional combining rules are employed to determine the parameters for a pair of unlike segments.
σ ij =
1 (σ i + σ j ) 2
(4.2.14)
∈ij = ∈i ∈ j (1 − k ij )
(4.2.15)
The integrals of the perturbation theory are substituted by simple power series in density − − 6 I 1 ⎛⎜η , m ⎞⎟ = ∑ ⎛⎜ i 1 m ⎞⎟η i ⎝ ⎠ i =0 ⎝ ⎠
(4.2.16)
− − 6 I 2 ⎛⎜η , m ⎞⎟ = ∑ bi ⎛⎜ m ⎞⎟η i ⎝ ⎠ i =0 ⎝ ⎠
(4.2.17)
Where the coefficients a i and bi depend on the chain length according to −
−
−
− m− 1 m− 1 m− 2 a i ⎛⎜ m ⎞⎟ = a 0i + − a li + − a 2i − ⎝ ⎠ m m m −
−
(4.2.18)
−
− m− 1 m− 1 m− 2 bi ⎛⎜ m ⎞⎟ = b0i + − b1i + − b2 i − ⎝ ⎠ m m m
(4.2.19)
The universal model constants for a 0i , a1i , a 2i , b0i , b1i , and b2i are available in literature.
Density The density at a given system pressure psys must be determined iteratively by adjusting the reduced density η until psys = psys . A suitable staring value for a liquid phase is
[
( )]
η = 0.5; for a vapor phase, η = 10 −10 . Values of η > 0.7405 = .τ / 3 2 are higher then the closest packing of segments and have no physical relevance. The number of molecules p is calculated from η through
ρ=
−1 6 ⎛ η ⎜ ∑ xi mi d i3 ⎞⎟ ⎠ π ⎝i
Modeling and Simulation of Solid-Liquid Equilibrium
(4.2.20)
27
Chapter 4: PC-SAFT Equation of State
The quantities ζ n given in Eq 4.2.8 can now be calculated. For a converged value of η , ^
we obtain the molar density ρ , in units of kmol/m 3 , from 3
0 ρ ⎛⎜ 10 Α ⎞⎟ ⎛ −3 kmol ⎞ 10 ρ= ⎜10 ⎟ N ΑV ⎜⎜ m ⎟⎟ ⎝ mol ⎠ ⎝ ⎠ ^
(4.2.21)
0
Where p is, according to Eq 4.2.20, given in units of Α 3 and N ΑV = 6.022 × 10 23 mol1 denotes Avogadro,s number.
Pressure Equations for the compressibility factor will be derived using the thermodynamic relation
⎛ ~ res ⎜∂a Z = 1 + η⎜ ⎜ ∂η ⎝
⎞ ⎟ ⎟ ⎟ ⎠ T , x1
(4.2.22)
The pressure can be calculated in units of Ρa = Ν / m 2 by applying the relation 0 ⎛ ⎞ ⎜ 10 Α ⎟ Ρ = ZkTρ 10 ⎜⎜ m ⎟⎟ ⎝ ⎠ From Eqs (4.2.22) and (4.2.3), it is Z = 1 + Z hc + Z disp
3
(4.2.23)
(4.2.24)
Hard-Chain Reference Contribution The residual hard-chain contribution to the compressibility factor is given by Z hc = m Z hs − ∑ xi (mi − 1)(g iihs ) ρ −
−1
i
∂g iihs ∂ρ
(4.2.25)
Where Z hs is the residual contribution of the hard-sphere fluid, given by
Z
hs
ζ3 3ζ 23 − ζ 3ζ 23 3ζ 1ζ 2 = + + (1 − ζ 3 ) ζ 0 (1 − ζ 3 )2 ζ 0 (1 − ζ 3 )3
Modeling and Simulation of Solid-Liquid Equilibrium
(4.2.26)
28
Chapter 4: PC-SAFT Equation of State
ρ
∂g ijhs ∂p
=
⎛ d i d j ⎞⎛ 3ζ 2 ζ3 6ζ 2ζ 3 ⎞ ⎟⎜ ⎟+ + +⎜ 2 2 (1 − ζ 3 ) ⎜⎝ d i + d j ⎟⎠⎜⎝ (1 − ζ 3 ) (1 − ζ 3 )3 ⎟⎠
(4.2.27)
2
⎛ d i d j ⎞ ⎛ 4ζ 22 6ζ 22ζ 3 ⎞ ⎟ ⎜ ⎜ ⎟ + ⎜ d + d ⎟ ⎜ (1 − ζ )3 (1 − ζ )4 ⎟ i j 3 3 ⎠ ⎠ ⎝ ⎝ hs and g ij was given in Eq. 4.2.7.
Dispersion Contribution The dispersion contribution to the compressibility factor can be written as Z disp = −2πρ
− ⎡ ⎤ ∂ (ηI 1 ) 2 ∂ (ηI 2 ) + C 2ηI 2 ⎥ m 2 ∈2 σ 3 m ∈ σ 3 − πρ m ⎢C1 ∂η ∂η ⎣ ⎦
(4.2.28)
Where
∂ (ηI 1 ) 6 ⎛ − ⎞ = ∑ a1 ⎜ m ⎟( j + 1)η j j =0 ∂η ⎝ ⎠
(4.2.29)
∂ (ηI 2 ) 6 ⎛ − ⎞ = ∑ b1 ⎜ m ⎟( j + 1)η j i =0 ⎝ ∂η ⎠
(4.2.30)
and where C 2 is an abbreviation defined as
C2 =
3 2 − ⎛ − − 4η 2 + 20η + 8 ⎛ ∂C1 ⎞ 2η + 12η − 48η + 40 ⎞⎟ m 1 = −C12 ⎜⎜ m + − ⎜ ⎟ ∂η ⎝ ⎠ (1 − η )5 [(1 − η )(2 − η )]3 ⎟⎠ ⎝
(4.2.31)
Fugacity Coefficient The fugacity coefficient ϕ k (T , P ) is related to the residual chemical potential according to lnϕ k =
µ kres (T ,υ ) kT
− ln Z
(4.2.32)
The chemical potential can be obtained from
⎛ − res ⎜∂a µ (T , v ) = a + (Z − 1) + ⎜ kT ⎜ ∂x k ⎝ res k
− res
⎡ ⎛ − res ⎞ N ⎢ ⎜∂a ⎟ − ∑ ⎢x j ⎜ ⎟ j =1 ⎟ ⎢ ⎜⎝ ∂x j ⎠ T ,u , xijk ⎣
Modeling and Simulation of Solid-Liquid Equilibrium
⎤ ⎞ ⎥ ⎟ ⎥ ⎟ ⎟ ⎠ T ,u , xijk ⎥⎦
(4.2.33)
29
Chapter 4: PC-SAFT Equation of State
Where derivatives with respect to mole fraction are calculated regardless of the summation relation ∑ j x j = 1 . For convenience, one can define abbreviations for derivatives of Eq. 4.2.8 with respect to mole fraction. ⎛ ∂ζ ⎞ ζ n , xk = ⎜⎜ n ⎟⎟ ⎝ ∂x k ⎠ T . ρ , x
=
π 6
ρmk (d k )n
n ∈ {0,1,2,3}
(4.2.34)
j/k
Hard-Chain Reference Contribution ⎛ ~ hc ⎜∂a ⎜ ∂x ⎜ k ⎝
⎞ ⎛ ~ hs hs ~ ⎟ ⎜∂a = mk a + m⎜ ⎟ ⎟ ⎜ ∂x k ⎠T ,ρ , x j / k ⎝
⎞ hs ⎟ hs −1 ⎛ ∂g ii ∑ x1 (mi − 1)(g ii ) ⎜⎜ ⎟ i ⎟ ⎝ ∂xk ⎠T ,ρ , x j / k
⎞ ⎟⎟ ⎠T ,ρ ,x j / k
(4.2.35)
With
⎛ ~ hs ⎜∂a ⎜ ∂x ⎜ k ⎝
⎞ ζ 0, xk ~ hs 1 ⎟ a + = − ⎟ ζ0 ζ0 ⎟ ⎠ T , p, x j / k
⎡ 3(ζ 1, xk ζ 2 + ζ 1ζ 2, xk ) 3ζ 1ζ 2ζ 3, xk 3ζ 22ζ 2, xk ζ 23ζ 3, xk (3ζ 3 − 1) ⎤ + + + +⎥ ⎢ (1 − ζ 3 ) (1 − ζ 3 )2 ζ 3 (1 − ζ 3 )2 ζ 32 (1 − ζ 3 )3 ⎥ ⎢ ⎥ ⎢⎛ 3ζ 2ζ ζ − 2ζ 3ζ 3 ⎞ ⎛ ⎞ ζ 3, xk ζ 2 3, xk 2 ⎥ ⎢⎜ 2 2, xk 3 − ζ 0, xk ⎟ ln (1 − ζ 3 ) + ⎜⎜ ζ 0 − 2 ⎟⎟ 3 ⎟ ⎥ ⎢⎜⎝ ( ) 1 ζ − ζ3 ⎠ ζ3 3 ⎝ ⎠ ⎦ ⎣
⎛ ∂g ijhs ⎜ ⎜ ∂x k ⎝
⎞ ⎛ did j ζ 3, xk ⎟ = +⎜ 2 ⎜d +d ⎟ j ⎝ i ⎠ T , p , x j / k (1 − ζ 3 )
⎛ did j ⎜ ⎜d +d j ⎝ i
⎞ ⎟ ⎟ ⎠
2
⎛ 4ζ 2ζ 2, xk 6ζ ζ 3, xk ⎜ + ⎜ (1 − ζ )3 (1 − ζ )4 3 3 ⎝ 2 2
⎞⎛ 3ζ 2, xk 6ζ 2ζ 3, xk ⎞ ⎟⎜ ⎟+ + 2 ⎜ ⎟ (1 − ζ ) (1 − ζ 3 )3 ⎟⎠ 3 ⎠⎝
⎞ ⎟ ⎟ ⎠
Modeling and Simulation of Solid-Liquid Equilibrium
(4.2.36)
(4.2.37)
30
Chapter 4: PC-SAFT Equation of State
Dispersion Contribution
⎛ ~ hs ⎜∂a ⎜ ∂x ⎜ k ⎝
⎞ ⎟ = −2πp I 1. xk m 2 ∈ σ 3 + I 1 m 2 ∈ σ 3 ⎟ ⎟ ⎠ T , p, x j / k
[
⎧
(
) ]− xk
(4.2.38)
(
πp ⎨ ⎡⎢m k C1 I 2 + m C1, xk I 2 + m C1 I 2, xk ⎤⎥ m 2 ∈2 σ 3 + m C1 I 2 m 2 ∈2 σ 3 −
−
−
⎩⎣
⎦
)
xk
⎫ ⎬ ⎭
With
(m
2
(m
2
⎛ ∈kj ∈ σ 3 )xk = 2m k ∑ x j m j ⎜⎜ j ⎝ kT ∈ σ 2
3
)
xk
⎞ 3 ⎟⎟σ kj ⎠
⎛ ∈kj = 2m k ∑ x j m j ⎜⎜ j ⎝ kT
(4.2.39)
2
⎞ 3 ⎟⎟ σ kj ⎠
⎧ 8η − 2η 2 20η − 27η 2 + 12η 3 − 2η ⎫ C1, xk = C 2ζ 3, xk − C12 ⎨m k m − ⎬ k (1 − η )4 [(1 − η )(2 − η )]2 ⎩ ⎭
(4.2.40) (4.2.41)
− 6 ⎡ ⎤ I 1, xk = ∑ ⎢a i ⎛⎜ m ⎞⎟iζ 3, xkη i −1 + a i , xη i ⎥ i =0 ⎣ ⎝ ⎠ ⎦
(4.2.42)
− 6 ⎡ ⎤ I 2, xk = ∑ ⎢bi ⎛⎜ m ⎞⎟iζ 3, xkη i −1 + bi , xη i ⎥ i =0 ⎣ ⎝ ⎠ ⎦
(4.2.43)
−
−
a i , xk
⎞ 4 mk ⎛ = 2 a1i + 2 ⎜ 3 − − a 2i ⎟ − − ⎜ ⎟ m ⎠ m m ⎝
bi , xk
⎞ 4 mk ⎛ = 2 a1i + 2 ⎜ 3 − − b2i ⎟ − − ⎜ ⎟ m ⎠ m m ⎝
mk
−
(4.2.44)
−
mk
(4.2.45)
Enthalpy and Entropy ^ res
The molar enthalpy h is obtained from a derivative of the Helmholtz free energy with respect to temperature, according to
Modeling and Simulation of Solid-Liquid Equilibrium
31
Chapter 4: PC-SAFT Equation of State ^ res ⎛ ~ res ⎞ ⎜ ∂a ⎟ h = −T ⎜ (4.2.46) ⎟ + (Z − 1) RT ⎜ ∂T ⎟ ⎝ ⎠ p , xi Unlike the enthalpy of an ideal gas, witch is a function of temperature only; the entropy of a real gas is a function of both temperature and pressure (or density). Hence, the residual entropy in the variables P and T is different from the residual entropy for the specified conditions u and T . It is ^ res
s
^ res
(P, T ) = s (v, T ) + ln(Z ) R
(4.2.47)
R
~ res
All of the equations for a entropy can be written as ^ res
s
⎡⎛
are given in the variables v and T , so that the residual
~ res ⎤ ⎞ ⎟ a ⎥ + + ln (Z ) ⎢⎜ ∂T ⎟ T ⎥ ⎟ ⎜ ⎢⎝ ⎥ ⎠ p , xi ⎣ ⎦ ~ res
(P, T ) = −T ⎢⎜ ∂ a R
^ res
The residual molar Gibbs free energy g ^ res
^ res
^ res
g h S = − RT RT
(4.2.48)
(P, T ) is defined as
(P, T ) R
(4.2.49)
or simply as ^ res
~ res g = a + (Z − 1) − ln (Z ) RT
(4.2.50)
~ res
The temperature derivative of a in Eqs 4.2.46 and 4.2.48 is again the sum of two contributions. ⎛ ~ disp ⎞ ⎛ ~ hc ⎞ ⎛ ~ res ⎞ ⎜∂a ⎟ ⎜∂a ⎟ ⎜∂a ⎟ +⎜ (4.2.51) =⎜ ⎟ ⎟ ⎜ ∂T ⎟ ⎜ ∂T ⎟ ⎜ ∂T ⎟ ⎟ ⎜ ⎠ p , xi ⎠ p , xi ⎝ ⎠ p , xi ⎝ ⎝ With abbreviations for two temperature derivatives
∂d i ∈ ⎞⎤ ⎛ ∈ ⎞⎡ ⎛ (4.2.52) = σ i ⎜ 3 i 2 ⎟ ⎢ − 0.12 exp⎜ − 3 i ⎟⎥ kT ⎠⎦ ∂T ⎝ kT ⎠ ⎣ ⎝ π ∂ζ ζ n ,T = n = ρ ∑ xi mi nd i ,T (d i )n −1 n ∈ {1,2,3} (4.2.53) ∂T 6 i The hard-chain contribution and the dispersion contribution can conveniently be written. d i ,T =
Modeling and Simulation of Solid-Liquid Equilibrium
32
Chapter 4: PC-SAFT Equation of State
Hard-Chain Reference Contribution ⎛ ~ hs ⎞ ⎛ ~ hc ⎞ hs −⎜∂a ⎟ ⎜∂a ⎟ hs −1 ⎛ ∂g ii ⎜ ( ) ( ) − − 1 m x m g = ∑ i i ii ⎜ ∂T ⎟ ⎜ ∂T ⎟ ⎜ ∂T i ⎟ ⎜ ⎟ ⎜ ⎝ ⎠ p , xi ⎝ ⎠ p , xi ⎝ ⎛ ⎜∂a ⎜ ∂T ⎜ ⎝
~ hs
⎞ ⎟ ⎟ ⎟ ⎠ p , xi
⎞ ⎟⎟ ⎠ p , xi
(4.2.54)
⎡ 3(ζ 1,T ζ 2 + ζ 1ζ 2,T ) 3ζ 1ζ 2ζ 3,T 3ζ 22ζ 2,T ζ 23ζ 3,T (3ζ 3 − 1)⎤ + + + ⎥ ⎢ (1 − ζ 3 ) (1 − ζ 3 )2 ζ 3 (1 − ζ 3 )2 ζ 32 (1 − ζ 3 )3 ⎥ 1 ⎢ = ⎥ (4.2.55) 3 ζ 0 ⎢⎛ 3ζ 22ζ 2,T ζ 3 − 2ζ 23ζ 3,T ⎞ ζ ⎞ ⎛ ζ ⎥ ⎢⎜ ⎟ ln(1 − ζ 3 ) + ⎜ζ 0 − 2 ⎟ 3,T 3 2 ⎟ ⎜ ⎟ ⎥ ⎢⎜⎝ ζ3 ζ 3 ⎠ (1 − ζ 3 ) ⎝ ⎠ ⎦ ⎣
Equation 4.2.54 requires only the i − i pairs in temperature derivative of the radial pair distribution function g ijhs . For simplicity, one can restrict oneself to the i − i pairs in Eq 4.2.7 by equating ⎛ dd ⎞ 1 d i = ⎜⎜ i i ⎟⎟ 2 ⎝ di + di ⎠ The temperature derivative of the radial pair distribution function g ijhs is then
ζ 3,T 6ζ 2ζ 3,T ⎞ ⎛ 1 ∂g iihs ⎛ 1 ⎞⎛ 3ζ 2,T ⎛1 ⎞ 3ζ 2 ⎟ + ⎜ d i d i ,T ⎞⎟ + ⎜ d i ⎟⎜⎜ + = + ⎜ d i ,T ⎟ 3 ⎟ 2 2 ( ) ∂T (1 − ζ 3 ) ⎝ 2 ⎠ (1 − ζ 3 ) ⎝ 2 ⎠⎝ 1 − ζ 3 (1 − ζ 3 ) ⎠ ⎝ 2 ⎠ 2ζ
2 2
6ζ 22ζ 3,T ⎛ 1 ⎞ ⎛⎜ 4ζ 2ζ 2,T + + ⎜ di ⎟ 3 (1 − ζ 3 )4 ⎝ 2 ⎠ ⎜⎝ (1 − ζ 3 ) 2
(1 − ζ 3 )3
⎞ ⎟ ⎟ ⎠
Dispersion Contribution ⎛ ~ disp ⎞ ⎜∂a ⎟ I ⎞ ⎛ ∂I = −2πρ ⎜ 1 − 1 ⎟ m 2 ∈ σ 3 − ⎜ ∂T ⎟ ⎝ ∂T T ⎠ ⎜ ⎟ ⎝ ⎠ p , xi −
(4.2.56)
(4.2.57)
(4.2.58)
I ⎤ ∂I ⎡ ∂C1 I 2 + C1 2 − 2C1 2 ⎥ m 2 ∈2 σ 3 T⎦ ∂T ⎣ ∂T
πp m ⎢ With
∂I 1 6 ⎛ − ⎞ = ∑ a i ⎜ m ⎟iζ 3,T η i −1 ∂T i = 0 ⎝ ⎠ − 6 ∂I 2 = ∑ bi ⎛⎜ m ⎞⎟iζ 3,T η i −1 ∂T i = 0 ⎝ ⎠ ∂C1 = ζ 3,T C 2 ∂T
Modeling and Simulation of Solid-Liquid Equilibrium
(4.2.59) (4.2.60) (4.2.61)
33
CHAPTER 5
REGRESSION ANALYSIS OF SOLUBILITY DATA
Chapter 5: Regression Analysis of Solubility Data Plenty of solubility data are available for low molecular weight pure crystalline n-alkanes and aromatic compounds in different solvent system. These experimental data are collected from different literature to show the suitability of the developed PC-SAFT model for low molecular weight system. PC-SAFT parameters of these solvents and solutes are shown in Table 5.1. Table 5.1 PC-SAFT Parameters of Organic Solutes and Solvents Hydrocarbon Segment no Segment Diameter 0
(m)
Energy Parameter
(σ, A)
(ε/κ, K)
3.8959 3.9592 3.9668 4.0217 4.0217 3.8151 3.9902
249.21 254.70 256.20 252.0 252.0 327.42 255.05
3.7729 3.8049 3.8384 3.6478 3.7169
231.20 238.40 243.87 287.35 285.69
Solutes n-dodecane(C12) n-hexadecane(C16) n-octadecane(C18) n-octacosane(C28) n-dotriacontane(C32) biphenyl ε-caprolactone
5.3060 6.6485 7.3271 10.3622 11.835 3.8877 2.8149
n-pentane n-heptane n-decane benzene toluene
2.6896 3.4831 4.6627 2.4653 2.8149
Solvents
5.1: SOLID-LIQUID EQUILIBRIUM OF n-ALKANES Illustrated in Fig 5.1.1 and Fig 5.1.2 are solubility of n-dodecane and n-hexadecane in nheptane solvent respectively, on the basis of experimental data taken from Hoerr, Harwood (3) as presented in tabular form in Table 5.1.1 & 5.1.2. These data fits very well with PC-SAFT model with Kij = 0.0. Fig 5.1.3 represents the solubility of n-octadecane in n-heptane solvent. Corresponding experimental data are taken from S. Chang, J. R. Maurey, and W. J. Pummer (7). These data fit very well with the model with only small Kij (= 0.0003).
Modeling and Simulation of Solid-Liquid Equilibrium
35
Chapter 5: Regression Analysis of Solubility Data Table 5.1.1 Experimental SLE Data (3) for System n-Dodecane (2) and n-Hexane (1) Wt Fraction (W2)
T, K
Wt Fraction (W2)
T, K
0.05 0.1 0.15 0.2 0.3 0.4
215.5 223 228 233 240 245.5
0.5 0.6 0.7 0.8 0.9 1.0
249 252.5 256 259 261.5 263.6
260
Temperature(K)
250
240
230
220 (3)
Experimental PC-SAFT(Kij=0.0) 210 0.0
0.2
0.4
0.6
0.8
1.0
Weight Fraction of n-Dodecane
Fig 5.1.1:SLE for system n-dodecane (C12) + n-hexane at 1 bar
Modeling and Simulation of Solid-Liquid Equilibrium
36
Chapter 5: Regression Analysis of Solubility Data Table 5.1.2 Experimental SLE Data (3) for System n-Hexadecane (2) and n-Hexane (1) Wt Fraction (W2)
T, K
Wt Fraction (W2)
T, K
0.025 0.05 0.1 0.15 0.2 0.3 0.4
236 245 252 256 260 265.8 270
0.5 0.6 0.7 0.8 0.9 1.0
274.7 278.4 282 284.6 287.6 291.2
290
Temperature(K)
280 270 260 250 240 (3)
Experimental PC-SAFT(K ij=0.0)
230 0.0
0.2
0.4
0.6
0.8
1.0
W eight Fraction of n-Hexadecane
Fig 5.1.2: SLE for system n-hexadecane (C16) and n-hexane at 1 bar Table 5.1.3 Experimental SLE Data (7) for System n-Octadecane (2) and n-Heptane (1) Wt Fraction (W2)
T, K
0.213 0.37
273 281
Wt. Fraction (W2) 0.38 0.395
Modeling and Simulation of Solid-Liquid Equilibrium
T, K 281.4 282
37
Chapter 5: Regression Analysis of Solubility Data
300
Temperature(K)
290
280
270
260 (7)
E xp e rim e n ta l P C -S A F T (K ij= 0 .0 0 0 3 )
250 0 .0
0 .2
0 .4
0 .6
0 .8
1 .0
W e ig h t F ra c tio n o f n -O c ta d e c a n e
Fig 5.1.3: SLE for system n-octadecane (C18) and n-heptane at 1 bar For n-alkanes larger than C20 such as n-octacosane and n-dotriacontane exhibits solidsolid phase transition a few degree below it’s melting point as shown in Table 5.1.5. These compounds form a higher temperature hexagonal geometry to a more stable crystalline structure at lower temperature. Effect of phase transition is included by model equation 3.1.14. Table 5.1.4 Solid-Solid Transition and Melting Properties of Aromatic Compounds Hydrocarbon
Mol. Weight
Tm, K
∆Hm, J/Mol
biphenyl ε-caprolactone
154.211 96
342.1 272.13
18732 13821
Table 5 .1.5 Solid-Solid Transition and Melting Properties of Normal Alkanes Hydrocarbon Tss, K ∆Hss, J/mol Tm, K n-dodecane (C12) n-hexadecane (C16) n-octadecane (C18) n-octacosane (C28) n-dotriacontane (C32)
331.2 338.9
35447 42700
Modeling and Simulation of Solid-Liquid Equilibrium
263.6 291.2 301.1 334.4 342.1
∆Hm, J/mol 36977 53563 59400 64658 76000
38
Chapter 5: Regression Analysis of Solubility Data
Experimental solubility data for system n-Dotriacontane in n-Heptane is taken from S. Chang, J. R. Maurey, and W. J. Pummer (7) as listed in Table 5.1.6. Comparison of these data with theoretical model is shown in Fig 5.1.4. Table 5.1.6 Experimental SLE Data (7) for System n-Dotriacontane (2) and n-Heptane (1) Wt Fraction (W2)
T, K
0.0491 0.0604 0.0976 0.201
302.9 304.5 308.1 314.0
Wt Fraction (W2)
T, K
0.332 0.499 0.666 0.802
319.4 324.7 329.3 333.9
340
Temperature(K)
330
320
310
300
290
280 0.0
(7)
Experimental PC-SAFT(Kij=0.0004) 0.2
0.4
0.6
0.8
1.0
Weight Fraction of n-Dotriacontane
Fig 5.1.4: SLE for system n-dotriacontane (C32) and n-heptane at 1 bar
Modeling and Simulation of Solid-Liquid Equilibrium
39
Chapter 5: Regression Analysis of Solubility Data 5.2: SOLID-LIQUID EQUILIBRIUM OF AROMATIC COMPOUNDS E. McLaughlin and H.A. Zainal (8) studied the solubility of biphenyl in benzene for higher mole fraction range as shown in Table 5.2.1. PC-SAFT model fits well with these data as shown in Fig 5.2.1. Table 5.2.1 Experimental SLE Data (8) for System Biphenyl (2) and Benzene (1) Mole Fraction (X2)
T, K
Mole Fraction (X2)
T, K
0.5118 0.6478
310 320.6
0.8195 0.8916
332.2 336.2
340 330 320
Temperature(K)
310 300 290 280 270 260 (8)
250
Experimental PC-SAFT(Kij=0.0)
240 0.0
0.2
0.4
0.6
Mole Fraction of Biphenyl
0.8
1.0
Fig 5.2.1: SLE for system biphenyl and benzene at 1 bar Illustrated in Table 5.2.2 are the saturation solubility data for system ε-caprolactone in toluene as taken from R. Witting, D. Constantineseu, and J. Gmchling (5). Corresponding figure and comparison with model predicted results are shown in Fig 5.2.2.
Modeling and Simulation of Solid-Liquid Equilibrium
40
Chapter 5: Regression Analysis of Solubility Data
Table 5.2.2 Experimental SLE Data (5) for System ε-Caprolactone (2) and Toluene (1) Mole Fraction (X2)
T, K
Mole Fraction (X2)
T, K
0.1009 0.1486 0.1942 0.2472 0.2981 0.3501 0.3975 0.4518 0.5039
210.83 218.64 223.74 229.07 233.61 237.57 240.92 244.43 247.63
0.5468 0.5896 0.6384 0.6878 0.7914 0.8474 0.8959 0.9369 1.0
250.02 252.41 254.98 257.42 262.44 265.08 267.4 269.26 272.18
280 270 260
Temperature(K)
250 240 230 220 210 200 190 (5)
Experimental PC-SAFT(K ij=0.0)
180 170 0.0
0.2
0.4
0.6
0.8
1.0
Mole Fraction of Caprolactone
Fig 5.2.2: SLE for system ε-caprolactone and toluene at 1 bar
Modeling and Simulation of Solid-Liquid Equilibrium
41
Chapter 5: Regression Analysis of Solubility Data 5.3: EFFECT OF PRESSURE ON SOLID-LIQUID EQUILIBRIUM H. G. Lee, F. R. Groves, and J. M. Walcott (6) measured the saturation condition for noctacosane + n-decane and n-octacosane + p-xylene + n-decane for approximately 10 mole% solid content and pressure up to 200 bar. Their experimental data are shown in Table 5.3.1 & 5.3.3. These data are plotted along with the model prediction for both binary and ternary system as shown in Fig 5.3.1 & 5.3.2. Data fits well with model with negligible Kij. Table 5.3.1 Effect of Pressure on Binary SLE: Experimental Data (6) for System n-Octacosane (2) and Decane (1) (X2 represents Mole Fraction of n-Octacosane) X2 =0.06013 X2=0.08198 X2=0.1074 T, K
P, Bar
T, K
P, Bar
T, K
P, Bar
307.2 308 308.9 309.7 310.4 311.6 312.5
1 45 77 110 143 189 224
310 310.3 311 312.2 313.4 314.4 314.7
1 18 46 97 155 195 212
312.5 312.7 313.5 314.4 315.3 316.3 317.4
2 11 45 89 126 169 216
250
---op en sym b ol-P C -S AF T(K ij= 0.0006) --closed sym b ol--E xp erim en tal
(6)
Pressure(Bar)
200
150
100
50
X 2= 0.06013 X 2= 0.08198 X 2= 0.1074
0 308
310
312
314
316
318
Tem p eratu re(K )
Fig 5.3.1: Effect of pressure on binary SLE of system n-octacosane (C28) and n-decane for different composition.
Modeling and Simulation of Solid-Liquid Equilibrium
42
Chapter 5: Regression Analysis of Solubility Data Correlations of solid and liquid molar volume with temperature for n-octacosane are listed in Table 5.3.2 Table 5.3.2 Correlation of Molar Volume (Cm3/Mol) and Temperature (K) for n-Octacosane Liquid molar volume, VL = 0.42203T + 365.588 Solid molar volume, VS = 0.11828T + 381.623 Table 5.3.3 Effect of Pressure on Ternary SLE: Experimental Data (6) for System n-Octacosane (3), P-Xylene (2) and n-Decane (1) (X1/X2 =2.0, X3 =0.09803) T, K P, Bar T, K P, Bar 310.5 1 313.3 126 311 23 314.6 194 311.6 45 315.3 218 312.3 83
250
200
Pressure(Bar)
150
100
50 (6)
Experimental PC-SAFT(K ij=-0.0004)
0 310
311
312
313
314
315
316
Temperature(K)
Fig 5.3.2: Effect of pressure on SLE of system n-decane (1)+p-xylene (2) + n-octacosane(3) (x1/x2=2.0,x3=0.09803)
Modeling and Simulation of Solid-Liquid Equilibrium
43
Chapter 5: Regression Analysis of Solubility Data 5.4: EFFECT OF SOLVENT ON SOLID-LIQUID EQUILIBRIUM 350 340
Temperature(K)
330 320 310 300 290
n -h ep tan e+ n -octacosan e(C 28) n -p en tan e+ n -octacosan e(C 28)
280 270 0.0
0.2
0.4
0.6
0.8
1.0
Mole F raction of n -O ctacosan e
Fig 5.4.1: Effect of solvent on SLE 5.5: EFFECT OF MOLECULAR WEIGHT AND MELTING TEMPERATURE ON SOLID-LIQUID EQUILIBRIUM 340
Temperature(K)
320 300 280 260 240 220
n -d od ecan e(C 12),(K ij= 0.000) n -h exad ecan e(C 16),(K ij= 0.000) n -d otriacon tan e(C 32),(K ij= 0.0004)
200 0.0
0.2
0.4
0.6
0.8
1.0
W eig h t F raction of S olu te
Fig 5.5.1: Effect of molecular weight on SLE in n-hexane Modeling and Simulation of Solid-Liquid Equilibrium
44
Chapter 5: Regression Analysis of Solubility Data 5.6 DISCUSSION Here some important conclusions are drawn regarding the results of low molecular weight organic compounds. •
As shown in Fig 5.1.1 to 5.1.3, PC-SAFT equation of state gives good agreement with the experimental results for n-alkanes. The requirement of adjustable binary interaction parameter (Kij) increases with increasing chain length for same solvent system.
•
For aromatic compounds, model prediction matches very well with literature experimental data even with Kij=0.0 as shown in Fig 5.2.1 and Fig 5.2.2.
•
Illustrated in Fig 5.3.1 and Fig 5.3.2 are the effects of pressure on solubility for binary and ternary systems. These figures show that with increase in pressure solubility decreases. So for saturated system if we increase the pressure at fixed temperature some solid will crystallize. This necessitates the study of solubility at elevated pressure.
•
Effect of solvent on solubility is small as shown in Fig 5.4.1. However with increasing the boiling point of solvent solubility decreases for same solute.
•
With increasing molecular weight of solute, melting point increases and consequently its solubility decreases for same solvent.
Modeling and Simulation of Solid-Liquid Equilibrium
45
CHAPTER 6
SENSITIVITY STUDY FOR POLYETHYLENE SYSTEM
Chapter 6: Sensitivity Study for polyethylene System Initially PC-SAFT model is regressed and tested on the solubility data for n-alkanes, and aromatic compounds as discussed in chapter 5. Also in the same chapter effect of pressure, solvent, melting point, and molecular weight are discussed rigorously for low molecular weight organic compounds. In this chapter, sensitivity study is performed using the same model as used for low molecular weight system to understand the effects of crystallizability, melting temperature, solvent, adjustable binary interaction parameter (kij) and pressure on solid– liquid equilibrium of polyethylene in m-xylene as shown in the following figures. The PC-SAFT parameters of polyethylene used for this study is listed below in Table 6.1 Table 6.1 PC-SAFT Parameters of Polyethylene Polyethylene
Segment No (m/M, mol/g)
LDPE HDPE
0.0263 0.0263
Segment Diameter (σ, 0A) 4.0217 4.0217
Energy Parameter (ε/κ, K) 249.5 252.0
6.1 RESULTS OF SENSITIVITY STUDY 420 410
Temparature(K)
400 390 380 370 360 350
1.0 bar 100 bar --------Kij=0.0,Tm=415K, C=0.4
340 330 0.0
0.2
0.4
0.6
0.8
W eight Fraction of Polyethylene(PE120K)
Fig 6.1.1: Effect of pressure on solubility of polyethylene in m-xylene
Modeling and Simulation of Solid-Liquid Equilibrium
47
Chapter 6: Sensitivity Study for polyethylene System
420 400
Temperature(K)
380 360 340
C=0.2 C=0.4 C=0.6 ---------Kij=0.0,Tm=415K
320 300 0.0
0.2
0.4
0.6
0.8
1.0
W t. Fraction of Polyethylen e(PE120K)
Fig 6.1.2: Effect of crystallizability fraction (C) on solubility of polyethylene in m-xylene at 1 bar
440 430 420
Tempetature(K)
410 400 390 380 370 360
Tm=415K Tm=440K -----1 bar,Kij=0.0,C=0.4
350 340 0.0
0.2
0.4
0.6
0.8
1.0
Wt. Fraction of Polyethylene(PE120K)
Fig 6.1.3: Effect of melting point (Tm) on solubility of polyethylene in m-xylene Modeling and Simulation of Solid-Liquid Equilibrium
48
Chapter 6: Sensitivity Study for polyethylene System
420 410 400
Temperature(K)
390 380 370 360 350
m -xy le n e c y c loh e xan e ----1 b ar,K ij= 0.0,Tm = 415K ,C = 0.4
340 330 0.0
0 .2
0.4
0 .6
0.8
1.0
W t. F rac tio n o f P oly e th y le n e(P E 1 20K )
Fig 6.1.4: Effect of solvent on solubility of polyethylene in m-xylene 420 410 400
Temperature(K)
390 380 370 360
kij=0.012 kij=0.01 kij=-0.0 kij=-0.008 ---1bar,Tm=415K,C=0.4
350 340 330 0.0
0.2
0.4
0.6
0.8
1.0
Wt. Fraction of Polyethylene(PE120K)
Fig 6.1.5: Effect of Kij on solubility of polyethylene in m-xylene
Modeling and Simulation of Solid-Liquid Equilibrium
49
Chapter 6: Sensitivity Study for polyethylene System 6.2 DISCUSSION Here some important conclusions are drawn regarding the results of sensitivity study of polyethylene system. In all figures solubility represent the saturation solubility. Here the study is based on molecular weight of polyethylene of 120000. But the phase diagram can also be drawn for other molecular weights as well. Here intentionally the effect of molecular weight is not shown because with change in molecular weight other property of polyethylene will also change such as melting point, binary interaction parameter, enthalpy of fusion etc. 1. With increase in pressure, solubility decreases at fixed temperature as shown in Fig 6.1. But the effect of pressure on solubility is not much for only a small change in pressure. But enormous change in pressure will decrease solubility of polyethylene to a great extent sufficient to cause industrial problem like clogging of pipelines due to crystallization of polyethylene. 2. With increase in crystallizability fraction (C), solubility decreases at fixed temperature and pressure as shown in Fig 6.2. Here influence of crystallizability fraction on solubility is much. 3. With increase in melting point, solubility decreases at fixed temperature, pressure and other fixed properties as shown in Fig 6.3. 4. Effect of solvent on solubility of polyethylene is not much; rather effect of solvent can be neglected when all other properties are same for both solvents as shown in Fig 6.4.But binary interaction parameter will be different for different solvents even for same solute. Solubility of polyethylene will be different for different solvents if this effect is accounted. 5. With increase in binary interaction parameter, Kij (adjustable parameter of PCSAFT EOS), solubility of polymer decreases at fixed temperature and pressure as shown in Fig 6.5. Here Kij influences much on solubility. 6. In all the figures, solute is totally soluble in the solvent at the melting point of solute. Extrapolation of the figures to the melting point gives the same result. 7. Sensitivity study for polyethylene system gives similar trend as it is seen for low molecular weight organic compounds.
Modeling and Simulation of Solid-Liquid Equilibrium
50
CHAPTER 7
RESULTS OF SOLUBILITY OF POLYETHYLENE
Chapter 7: Results of Solubility of Polyethylene Experimental solubility data for polyethylene system is limited. In the year 1946, Richards (3) studied the solubility of polyethylene (of molecular weight 17000 and melting point 387.5K) in m-xylene at atmospheric pressure. Their experimental data is tabulated in Table 7.1.
Table 7.1 Experimental SLE Data (3) for System Polyethylene (2) and m-Xylene (1) Wt Fraction (W2)
T, K
0.012 0.03 0.1 0.32
345 349 354 363.75
Wt Fraction (W2)
T, K
0.53 0.7 1
371.75 377 387.5
Experimental data fits very well with PC-SAFT model with C = 0.8 as shown in Fig 7.1
390
Temperature(K)
380
370
360
350 (3)
Experimental PC-SAFT (Kij=0.0055) 340 0.0
0.2
0.4
0.6
0.8
1.0
Wt. Fraction of Polyethylene
Fig 7.1: Solubility of polyethylene in m-xylene at 1 atmospheric pressure and prediction by PC-SAFT model
Modeling and Simulation of Solid-Liquid Equilibrium
52
Chapter 7: Results of Solubility of Polyethylene
7.1 EXPERIMENTAL DETERMINATION OF SOLUBILITY Solubility of polyethylene in xylene is determined for two commercial grade samples at atmospheric pressure. SLE was measured by visual technique. Since there is no standard apparatus for measurement of solubility, two opening conical flux fitted with a condenser serves the purpose. Temperature is measured using resistance thermometer inserted through side opening of the conical flux. Measured amount of solvent and solute was taken in a conical flux and heated slowly with occasional manual shaking. The temperature at which all solutes disappear is the saturation solubility temperature. The properties of two samples are tabulated below in Table 7.1.1 Table 7.1.1 Properties of Polyethylene Sample
Melting Point, K
Grade1 Grade2
405 410
Molecular Weight, g/mol 30398.9 32599.8
Experimental results of the study are shown below in Table 7.1.2 & Table 7.1.3 Table 7.1.2 Experimental SLE Data for Grade1 Polyethylene (2) in Xylene (1) Wt Fraction (W2)
T, K
0.019 0.04623 0.0674 0.106
360.3 362.1 365.6 366.7
Wt Fraction (W2) 0.1276 0.173 0.245
T, K 367.6 369.2 371.5
Table 7.1.3 Experimental SLE Data for Grade2 Polyethylene (2) in Xylene (1) Wt Fraction (W2)
T, K
0.0198 0.046 0.074
377.9 380 81.1
Wt Fraction (W2) 0.209 0.3028 0.4032
Modeling and Simulation of Solid-Liquid Equilibrium
T, K 386.2 388.8 392.2
53
Chapter 7: Results of Solubility of Polyethylene 0.101
82.5
0.504
394.7
Model prediction and comparison with experimental results are shown in subsequent figures. Crystallizability fractions (C) of the polyethylene samples are not known. So crystallizability fractions are determined from the model. It is found that for C = 0.8, model prediction fits very well with the experimental results for both grades of samples. Solubility measurement is done with xylene. But proportion of individual xylenes (o-, m-, and p-xylenes) is not known. So PC-SAFT parameters of solvent is not known though parameters of individual xylenes are available in literature. For simplicity parameters of m-xylene are used for modeling of solubility of polyethylene, since the parameters are close to one another. PC-SAFT parameters for polyethylene used for modeling are listed in Table 6.1. Parameters of xylenes are listed in Table 7.1.4.
Modeling and Simulation of Solid-Liquid Equilibrium
54
Chapter 7: Results of Solubility of Polyethylene
385
Temperature (K)
380
375
370
365
Experimental PC-SAFT C=0.8, Kij=0.002,Tm=405K
360
355 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
Wt. Fraction of Polyethylene
Fig 7.1.1: Solubility of polyethylene (PE30398.9) in xylene at atmospheric pressure Table 7.1.4 PC-SAFT Parameters of Xylenes Xylenes
Segment No (m)
m-xylene o-xylene p-xylene
3.1861 3.1362 3.1723
Segment Diameter (σ, 0A) 3.7563 3.7600 3.7781
Modeling and Simulation of Solid-Liquid Equilibrium
Energy Parameter (ε/κ, K) 291.05 283.77 288.13
55
Chapter 7: Results of Solubility of Polyethylene
410
Temperature (K)
400
390
380
370 0.0
Experimental PC-SAFT (Kij=0.0135) C = 0.8, Tm= 410 K 0.2
0.4
0.6
Wt. Fraction of Polyethylene
Fig 7.1.2: Solubility of polyethylene (PE32599.8) in xylene at atmospheric pressure
7.2 CONCLUSIONS AND FUTURE SCOPE OF WORK •
From the previous figures it is clear that PC-SAFT model can be used to predict the solubility of polyethylene with negligible error if crystallizability fraction, melting point, and molecular weight of polyethylene are known and binary interaction parameter is adjusted accordingly to fit the experimental data of polyethylene-solvent system.
Modeling and Simulation of Solid-Liquid Equilibrium
56
Chapter 7: Results of Solubility of Polyethylene •
Here crystallizability fraction is determined from the model. But if the values of the same are known or measured then from the model it is possible to determine correct value of binary interaction parameter. From this it is possible to correlate binary interaction parameter with molecular weight of polyethylene for same solvent system.
•
In this work the study is only limited to single solvent. This study can be extended to different solvents as well. In order to overcome modeling difficulty and to get parameters of the solvent, laboratory grade pure solvent is to be used for the experimental study.
•
The study can be extended to large variety of laboratory grade polymers of different molecular weight. Study can also be done with narrow molecular weight distribution polyethylene in order to correlate the properties of polyethylene with molecular weight.
•
Also the solubility study can be extended to high pressure as well in order to get idea about the pressure effect on solubility. Solubility of polymer under pressure is important industrially because most of the polymer is manufactured under high pressure.
•
Melting temperature and its relation with pressure is important for the study of solubility at high pressure. This can be measured by Differential Scanning Calorimeter (DSC).
•
No. average and weight average molecular weight of polymer and its molecular weight distribution is important in order to extend this model to that based on pseudocomponent approach.
Modeling and Simulation of Solid-Liquid Equilibrium
57
REFERENCES
1. Chemical Engineering Thermodynamics By Sandlar. 2. Molecular Thermodynamics of Fluid Phase Eqiulibria, 2nd edition by J. M. Prausnitz. 3. Cheng Pan, Maciej Radosz, Fluid Phase Equilibria 155, (1999) 57-73 4. Joachim Gross and Gabriele Sadowski, Ind. Engg. Chem. Res., 2001, 40, 1244-1260. 5. Roland Witting, Dana Constantineseu, and Jurgen Gmchling, J. Chem. Engg. Data, 2001, 46, 1490-1493. 6. Hyo-Guk Lee, Frank R. Groves, and Joanne M. Walcott, J. Chem. Engg. Data, 1993, 38, 257-259. 7. Shu-Sing Chang, John R. Maurey, and Walter J. Pummer, J. Chem. Engg. Data, 1983, 28, 187-189. 8. E. McLlaughlin and H.A. Zainal, J. Chem. Soc. 1959, 863 9. E. McLlaughlin and H.A. Zainal, J. Chem. Soc. 1960, 2485 10. U. Domanska, F. R. Groves, Jr., and E. McLlaughlin, J. Chem. Engg. Data, 1993, 38, 88-94 11. Y. Tanaka and M. Kawakami, Fluid Phase Equilibria, 125 (1996) 103-114 12. Hyo-Guk Lee, Philip A. Schenewerk, and Joanne Walcott, and Frank R. Groves Jr, Fluid Phase Equilibria, 128 (1997) 229-240. 13.J.M. Prausnitz, R.N. Lichtenthaler, E.G. Azevedo, Molecular Thermodynamics of Fluid Phase Equilibria, 2nd edn., Prentice-Hall, 1986. 14. V.I. Harismiadis, D.P. Tassios, Ind. Eng. Chem. Res. 35 1996 4667–4681. 15. D.W. Van Krevelen, Properties of Polymers, Elsevier, 1990, p. 120. 16. C.W. Hoerr, H.J. Harwood, J. Org. Chem. 16 1951 799
Nomenclature
NOMENCLATURE
C f
crystallizability fugacity
k
Boltzmann’s constant, 1.38 ×10 -23 J/K
N P sat P PE R
Avogadro’s number, 6.023 ×10-23 molecules/mol saturated pressure of solute at melting temperature pressure, bar polyethylene gas constant
SAFT
statistical association fluid theory
PC-SAFT SLE SL SS T u X, x Z
perturbed-Chain statistical association fluid theory solid–liquid equilibrium solid–liquid solid-solid temperature, K number of ethyl units in backbone mole fraction compressibility factor
ρ ρ a, ρ c φ φ0 ∆CP
density, mol/ml density of completely amorphous and completely crystalline, g / cm 3 fugacity coefficient fugacity coefficient of pure sub cooled liquid solute at given T and P CpL - CpS the heat capacity change of solute, J /K
∆G ∆H ∆ HU
enthalpy of melting per crystal unit, for ethyl unit, ∆HU =8.22 kJ / mol
∆S
molar entropy change
∆v
vL-vS volume change, cm3/mol
A
Helmholtz free energy, J
A1
Helmholtz free energy of first-order perturbation term, J
molar Gibbs free energy change molar enthalpy change
Modeling and Simulation of Solid-Liquid Equilibrium
58
Nomenclature A2
Helmholtz free energy of second-order perturbation term, J
a 01 , a 02 , a 03
model constants; defined in Eq. 4.2.18
a j (m )
functions defined by Eqs 4.2.18 and 4.2.19
b01 , b02 , b03
model constants defined in Eq. 4.2.19
d
temperature-dependent segment diameter, A
g hc
average radial distribution function of hard-chain fluid
g hc
site-site radial distribution function of hard-chain fluid
I1 , I 2
abbreviations defined by Eqs. 4.2.16, 4.2.17
k
Boltzmann constant, J/K
k ij
binary interaction parameter
Κ
Κ factor, Κ i = y i / x i
m
number of segments per chain
0
−
m M N P
mean segment number in the system, defined in Eq. 4.2.5 molar mass, g/mol total number of molecules pressure, Pa
R r
gas constant, J mol -1 Κ −1 radial distance between two segments, A
s1
constant defining the pair potential, defined in Eq.4.1.1,
T
temperature, Κ
u (r )
υ
0
pair potential function, J, defined in Eq.4.1.1
W x
molar volume, m 3 / mol weight fraction reduced radial distance between two segments
xi
mole fraction of component i
Z
compressibility factor
Greek Letters
∈
Depth of pair potential, J
η
packing fraction, η = ζ n , defined in Eqs.4.2.20 & 4.2.8
λ
reduced well width of square-well potential
Modeling and Simulation of Solid-Liquid Equilibrium
59
Nomenclature 0 3
p
total number density of molecules, 1 / A
σ
segment diameter, A
ζn
abbreviation (n = 0,...,3) defined by Eq. 4.2.8, A
0
0 n3
Superscripts calc crit disp exp hc hs id sat L ref res S i
calculated property critical property contribution due to dispersive attraction experimental property residual contribution of hard-chain system residual contribution of hard-sphere system ideal gas contribution property at saturation condition liquid phase reference residual solid phase component I
Subscripts 2 O2 i m p ss
solute pure solute component i melting condition polymer solid–solid phase transition
Modeling and Simulation of Solid-Liquid Equilibrium
60
Appendix A: Derivation of the Pure-Solute F02l/F02s
APPENDIX A: DERIVATION OF THE PURE-SOLUTE F02L/F02S The fugacity ratio of pure solute is derived from a thermodynamic cycle. Similar to what Prausnitz et al did for the temperature effect, we generate a loop for both, temperature and pressure effects. The detailed derivations are shown below.
The molar Gibbs energy change for solute from state a to state f is related to the fugacities of solid and sub cooled liquid.
(A-1) It is also related to the corresponding enthalpy and entropy changes by
G(a
f)= H(a
f)-T S(a
f)
(A-2)
For the enthalpy change from state a to state f, we have
H(a f)= H(a b)+ H(b H(d e)+ H(e f)
c)+ H(c
Modeling and Simulation of Solid-Liquid Equilibrium
d)+
(A-3)
61
Appendix A: Derivation of the Pure-Solute F02l/F02s
Using the Maxwell relations,
(A-4) Eq. A-3 can be rewritten in terms of the heat capacity, CP, and the enthalpy of melting, Hm, as follows:
(A-5)
For the entropy change from state a to state f, we have
S(a f)= S(a b)+ S(b S(d e)+ S(e f)
c)+ S(c
d)+
(A-6)
Using the Maxwell relations again,
(A-7) We rewrite the entropy change given by Eq. A-6 as
(A-8)
Modeling and Simulation of Solid-Liquid Equilibrium
62
Appendix A: Derivation of the Pure-Solute F02l/F02s
Substituting Eq. A-5 and Eq. A-8 into Eq. A-2, we obtain
(A-9)
Where CP=CPL-CPS, and v=vL-vS. Since pressure has little effect on taken out of the integration directly.
v,
v can be
Substituting Eq. A-9 into Eq. A-1, we get the fugacity ratio given in Eq. 3.1.8.
Modeling and Simulation of Solid-Liquid Equilibrium
63
Appendix B: Program of Solid-Liquid Equilibrium Calculation
APPENDIX B: PROGRAM OF SOLID-LIQUID EQUILIBRIUM CALCULATION IMPLICIT NONE INTEGER nc, nph PARAMETER (nc=20,nph=2) INTEGER ncomp, nphas DOUBLE PRECISION t,p,densys(nph),x_sys(nc) DOUBLE PRECISION parame(nc,5),kij(nc,nc) DOUBLE PRECISION d_sta c-----local variables------------------------------------------------INTEGER i,j,ph DOUBLE PRECISION phi(nc),phio ,tmax,z,phi2,r,tm,ts,hm,hs,alnx,x2 DOUBLE PRECISION xtrans(nc),dstart,parsys(nc,25),fugcoe(nc),den DOUBLE PRECISION paramel(nc,5),v,c,u,z1,mp,ms c--------------------------------------------------------------------PARAMETER (r=8.314) c--------------------------------------------------------------------OPEN (85,FILE='pe_et.inp') OPEN (86,FILE='output.txt') READ(85,*) Read(85,*)p READ(85,*) Read(85,*)t READ(85,*) Read(85,*)tmax READ(85,*) Read(85,*)tm READ(85,*) Read(85,*)ts READ(85,*) Read(85,*)hm READ(85,*) Read(85,*)hs READ(85,*) Read(85,*)c READ(85,*) Read(85,*)u READ(85,*) Read(85,*)z READ(85,*) Read(85,*)v READ(85,*) Read(85,*)mp READ(85,*) Read(85,*)ms
Modeling and Simulation of Solid-Liquid Equilibrium
64
Appendix B: Program of Solid-Liquid Equilibrium Calculation READ(85,*) do i=1,2 read(85,*)paramel(i,1),paramel(i,2),paramel(i,3) enddo close(85) p=p*1.0d5 10 ncomp=1 parame(1,1)= paramel(1,1) parame(1,2)= paramel(1,2) parame(1,3)= paramel(1,3) kij(1,1)=0.0 x_sys(ncomp)=1.0 d_sta=0.4 dstart = d_sta DO 11 i = 1,ncomp xtrans(i) = x_sys(i) parsys(i,1) = parame(i,1) parsys(i,2) = parame(i,2) parsys(i,3) = parame(i,3) 11 CONTINUE CALL PHIEOS (fugcoe,xtrans,t,p,parsys,kij,ncomp,dstart,den) phio=fugcoe(1)
12
ncomp=2 do i=1,ncomp do j=1,ncomp if(i.eq.j)then kij(i,j)=0.0 else kij(i,j)=0.00 endif enddo enddo do i=1,ncomp parsys(i,1) = paramel(i,1) parsys(i,2) = paramel(i,2) parsys(i,3) = paramel(i,3) enddo xtrans(1) = z xtrans(2) =1-z CALL PHIEOS (fugcoe,xtrans,t,p,parsys,kij,ncomp,dstart,den) phi2=fugcoe(1) alnx= -((hm/(r*tm))*((tm/t)-1.0)+(hs/(r*ts))*((ts/t)-1.0) 1 +v*p*1.0d-5/(82.03*t))*c*u -(phi2-phio) x2=exp(alnx) if(dabs(x2-z).lt.1.d-10)then
Modeling and Simulation of Solid-Liquid Equilibrium
65
Appendix B: Program of Solid-Liquid Equilibrium Calculation c
conversion of mole fraction to weight fraction z1=z*mp/(z*mp+(1-z)*ms) write(86,*)z1,t else z=z-(z-x2)*0.008 go to 12 endif if(t.lt.tmax)then t=t+2.0 go to 10 else endif STOP END c********************************************************************* SUBROUTINE PHIEOS (phi,x,t,p,parame,kij,ncomp,densta,dense) IMPLICIT NONE c-----variables used in the parameter list of subroutine-------------INTEGER nc PARAMETER (nc=20) INTEGER ncomp DOUBLE PRECISION phi(nc) DOUBLE PRECISION kij(nc,nc) DOUBLE PRECISION pges,pgesdz,gij(nc, nc) DOUBLE PRECISION fres DOUBLE PRECISION x(nc),t,p,parame(nc,25),mseg(nc) DOUBLE PRECISION densta,dense,dap_dx(nc,7),dbp_dx(nc,7) DOUBLE PRECISION order1,order2,apar(7),bpar(7) DOUBLE PRECISION A,B,C,D,AA,BB(nc),CC(nc),DD(nc),EFF(nc), & QQ(nc,nc),PIJK(nc,nc,nc) DOUBLE PRECISION PI, RGAS, NA, KBOL, TAU DOUBLE PRECISION uij(nc,nc),d00ij(nc,nc),d0(nc) c-----local variables------------------------------------------------INTEGER i, k, m DOUBLE PRECISION zms, rho, m_mean, term1, term2 DOUBLE PRECISION mhs(nc), mdisp(nc), mcha(nc), mpart(nc), & myres(nc), myresq, lnphi(nc) DOUBLE PRECISION dgijdx(nc, nc, nc),ddendx(nc) DOUBLE PRECISION zres, zges DOUBLE PRECISION I1, I2, I1_dx, I2_dx, & ord1dx, ord2dx, ddrdx, & drdpkt, te1,te2,te3, dte1dx,dte2dx,dte3dx c------obtain parameters and density independent expressions---------CALL PERTPAR (kij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK,
Modeling and Simulation of Solid-Liquid Equilibrium
66
Appendix B: Program of Solid-Liquid Equilibrium Calculation 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) c------density iteration---------------------------------------------CALL DENSITR (pges,pgesdz,gij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) c-----residual Helmholtz free energy---------------------------------CALL F_EOS (fres,gij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) zms = 1.d0 - dense rho = 6.d0 * dense / (PI*D) m_mean = A c-----compressibility factor z = p/(kT*rho)--------------------------zges = (p * 1.d-30)/(KBOL*t*rho) zres = zges - 1.d0 c-----calcul. the derivatives of f to mole fraction x ( d(f)/d(x) )--DO 1 k = 1,ncomp c-------d(f)/d(x) : hard sphere contribution-------------------------term1 = DLOG(zms)*(-1.d0+3.d0*CC(k)**2.d0 -2.d0*CC(k)**3.d0) & + (3.d0*CC(k)**2.d0 - 3.d0*CC(k)**3.d0) & * dense / zms**2.d0 term2 = (3.d0*CC(k)+3.d0*EFF(k)-CC(k)**3.d0 + AA*DD(k))/zms & + dense*3.d0 * BB(k) / zms**2.d0 & + 2.d0 * CC(k)**3.d0 / zms**3.d0 mhs(k) = mseg(k)* (term1 + term2 * dense) c-------d(f)/d(x) : chain term---------------------------------------c-------d(dense)/d(x)=ddendx(p,k) and d(gij)/d(x)-------------------ddendx(k) = dense/D*mseg(k)*d0(k)**3.d0 DO 10 i = 1, ncomp dgijdx(i,i,k) = ((1.d0+ 3.d0*PIJK(i,i,k)) & /(zms**2.d0) + 2.d0*dense*QQ(i,i)* & (3.d0+2.d0*PIJK(i,i,k))/zms**3.d0+ & dense**2.d0*6.d0*QQ(i,i)**2.d0/ & zms**4.d0)*ddendx(k) 10 CONTINUE mcha(k) = 0.d0 DO 12 i = 1, ncomp mcha(k) = mcha(k) + x(i) * (1.d0-mseg(i)) & * (1.d0/gij(i,i)) * dgijdx(i,i,k) 12 CONTINUE mcha(k) = mcha(k)+( 1.d0-mseg(k))*DLOG(gij(k,k)) c-------d(f)/d(x) : dispersion contribution--------------------------I1 = 0.d0
Modeling and Simulation of Solid-Liquid Equilibrium
67
Appendix B: Program of Solid-Liquid Equilibrium Calculation I2 = 0.d0 I1_dx = 0.d0 I2_dx = 0.d0 DO 14 m = 0,6 I1 = I1 + apar(m+1)*dense**DBLE(m) I2 = I2 + bpar(m+1)*dense**DBLE(m) I1_dx = I1_dx + apar(m+1)*DBLE(m)*dense**DBLE(m-1)*ddendx(k) & + dap_dx(k,m+1)*dense**DBLE(m) I2_dx = I2_dx + bpar(m+1)*DBLE(m)*dense**DBLE(m-1)*ddendx(k) & + dbp_dx(k,m+1)*dense**DBLE(m) 14 CONTINUE ord1dx = 0.d0 ord2dx = 0.d0 DO 16 i = 1,ncomp ord1dx= ord1dx & + 2.d0*mseg(k) *x(i)*mseg(i)*d00ij(i,k)**3.d0 *uij(i,k)/t ord2dx= ord2dx & + 2.d0*mseg(k) *x(i)*mseg(i)*d00ij(i,k)**3.d0*(uij(i,k)/t)**2.d0 16 CONTINUE te1 = zms**4.d0 *(2.d0-dense)**2.d0 te2 = zms**2.d0 * ( zms**2.d0*(2.d0-dense)**2.d0 & +dense*(-2.d0*dense**3.d0+12.d0*dense**2.d0 & -27.d0*dense+20.d0) ) te3 = dense* ( (2.d0-dense)**2.d0*(8.d0-2.d0*dense) & -zms**2.d0*(-2.d0*dense**3.d0 & +12.d0*dense**2.d0-27.d0*dense+20.d0) ) drdpkt = te1 / (te2 + m_mean*te3) dte1dx = ( - 4.d0*zms**3.d0 *(2.d0-dense)**2.d0 & - 2.d0*zms**4.d0 *(2.d0-dense) )*ddendx(k) dte2dx = dte1dx+ ( (3.d0*dense**2.d0 - 2.d0*2.d0*dense + 1.d0) & *(-2.d0*dense**3.d0+12.d0*dense**2.d0 & -27.d0*dense+20.d0) & + dense*zms**2.d0*(-6.d0*dense**2.d0 & +24.d0*dense -27.d0) )*ddendx(k) dte3dx = ddendx(k)* ( te3/dense & + dense*( -2.d0*(2.d0-dense)*(8.d0-2.d0*dense) & +(2.d0-dense)**2.d0*(-2.d0) & + 2.d0*zms*(-2.d0*dense**3.d0 & +12.d0*dense**2.d0-27.d0*dense+20.d0) & -zms**2.d0*(-6.d0*dense**2.d0 & +24.d0*dense -27.d0) ) ) ddrdx = ( dte1dx*(te2+m_mean*te3) & -te1*(dte2dx +(m_mean*dte3dx+mseg(k)*te3)) ) & / (te2 + m_mean*te3)**2.d0
Modeling and Simulation of Solid-Liquid Equilibrium
68
Appendix B: Program of Solid-Liquid Equilibrium Calculation mdisp(k)= -2.d0*PI*rho*(order1*I1_dx+ord1dx*I1) & - PI*rho*drdpkt*m_mean*(order2*I2_dx+ord2dx*I2) & - PI*rho*(drdpkt*mseg(k)+ddrdx*m_mean)*order2*I2 c-----d(f)/d(x) : summation of all contributions---------------------mpart(k) = mhs(k) + mcha(k) + mdisp(k) 1 CONTINUE myresq = 0.d0 DO i = 1, ncomp myresq = myresq - x(i)* mpart(i) END DO myresq = myresq + fres + zres DO k = 1, ncomp myres(k) = myresq + mpart(k) lnphi(k) = myres(k) - DLOG(zges) phi(k) = lnphi(k) END DO RETURN END c******************************************************************** SUBROUTINE P_EOS (pges,pgesdz,gij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) IMPLICIT NONE c-----variables used in the parameter list of subroutine-------------INTEGER nc PARAMETER (nc=20) INTEGER ncomp DOUBLE PRECISION pges,pgesdz,gij(nc, nc) DOUBLE PRECISION x(nc),t,p,parame(nc,25),mseg(nc) DOUBLE PRECISION densta,dense,dap_dx(nc,7),dbp_dx(nc,7) DOUBLE PRECISION order1,order2,apar(7),bpar(7) DOUBLE PRECISION A,B,C,D,AA,BB(nc),CC(nc),DD(nc),EFF(nc), & QQ(nc,nc),PIJK(nc,nc,nc) DOUBLE PRECISION PI, RGAS, NA, KBOL, TAU DOUBLE PRECISION uij(nc,nc),d00ij(nc,nc),d0(nc) c-----local variables------------------------------------------------INTEGER i, j, m DOUBLE PRECISION rho,zms,m_mean,ddendv DOUBLE PRECISION term1,term2,term3 DOUBLE PRECISION phs,pdisp,pcha,pideal DOUBLE PRECISION dgijdz(nc,nc),dgijd2(nc,nc),zmsdz, & dvzdz,te1dz,te2dz, & te3dz,pdspdz,fdspdz,fdspd2,pchadz,fchdz,fchd2, & piddz,phsdz DOUBLE PRECISION I2, edI1dz, edI2dz, edI1d2, edI2d2,
Modeling and Simulation of Solid-Liquid Equilibrium
69
Appendix B: Program of Solid-Liquid Equilibrium Calculation & z00, z00_dz, z00_d2, z00_d3, & drdpkt, drdpdz, drdpd2 c------------------------------------------------------------------rho = 6.d0 * dense / (PI*D) zms = 1.d0 - dense m_mean = A C-----gij , the derivative dgijdz=d(gij)/d(dense) -----------------C-----and dgijd2 = dd(gij)/d(dense)**2 ----------------------------DO 2 i = 1, ncomp j=i gij(i,j) = 1.d0/zms + 3.d0*QQ(i,j)*dense/zms**2.d0 & +2.d0*(QQ(i,j)*dense)**2.d0/zms**3.d0 dgijdz(i,j) = (1.d0+3.d0*QQ(i,j)) /zms**2.d0 & +dense*(6.d0*QQ(i,j)+4.d0*QQ(i,j)**2.d0)/zms**3.d0 & +dense**2.d0*(6.d0*QQ(i,j)**2.d0) /zms**4.d0 dgijd2(i,j) = 2.d0*qq(i,j)*(( 1.d0/qq(i,j) & + 6.d0+2.d0*qq(i,j))/zms**3.d0+ (9.d0*dense & + 12.d0*qq(i,j)*dense)/zms**4.d0 +(12.d0*qq(i,j) & * dense**2.D0)/ zms**5.d0) 2 CONTINUE c-----derivations of dense to volume (ddendv = -d(dense)/d(volume)--ddendv = dense**2.d0 *6.d0/ PI /D c-----p : ideal gas contribution------------------------------------pideal = rho c-----p : hard sphere contribution----------------------------------term1 = (A-C**3.d0/D**2.d0) / zms term2 = (3.d0*B*C/D-C**3.d0/D**2.d0)/zms**2.d0 term3 = (2.d0*C**3.d0/D**2.d0)/zms**3.d0 phs = ddendv*(term1 + term2 + term3) c-----p : chain term------------------------------------------------fchdz = 0.d0 DO 4 i= 1, ncomp fchdz= fchdz + (x(i)*(1.d0-mseg(i))) & * (1.d0/gij(i,i)) * dgijdz(i,i) 4 CONTINUE pcha = ddendv * fchdz c------p : dispersion contribution----------------------------------c------edI1dz is equal to d(dense*I1)/d(dense)----------------------c------edI2dz is equal to d(dense*I2)/d(dense)----------------------I2 = 0.d0 edI1dz = 0.d0 edI2dz = 0.d0 DO 6 m=0,6 I2 = I2 + bpar(m+1)*dense**DBLE(m) edI1dz = edI1dz + apar(m+1)*DBLE(m+1)*dense**DBLE(m) edI2dz = edI2dz + bpar(m+1)*DBLE(m+1)*dense**DBLE(m)
Modeling and Simulation of Solid-Liquid Equilibrium
70
Appendix B: Program of Solid-Liquid Equilibrium Calculation 6
CONTINUE
z00 = 1.d0 + m_mean*(4.d0*dense-2.d0*dense**2.d0)/ZMS**3.d0 & + (1.d0 - m_mean)*(5.d0*dense-2.d0*dense**2.d0) & /(ZMS*(2.d0-dense)) z00_dz = m_mean*(-2.d0*dense**2.d0+4.d0*dense+4.d0)/ZMS**4.d0 & + (1.d0 - m_mean)*(dense**2.d0-8.d0*dense+10.d0) & /(ZMS*(2.d0-dense))**2.d0 drdpkt = 1.d0/(Z00 + dense*Z00_dz) c--z00_d2 = m_mean*(-4.d0*dense**2.d0+8.d0*dense+20.d0)/ZMS**5.d0 & + (1.d0 - m_mean) & *(-2.d0*dense**3.d0+24.d0*dense**2.d0-60.d0*dense+44.d0) & /(ZMS*(2.d0-dense))**3.d0 drdpdz = (-2.d0*z00_dz -z00_d2*dense)*drdpkt**2.d0 c be aware: rho=6.d0*dense/(PI*D) fdspdz = -2.d0*PI *6.d0/(PI*D) * edI1dz * order1 & -( PI *6.d0*dense/(PI*D) *drdpdz* m_mean *I2 *order2 & +PI *6.d0/(PI*D) *drdpkt* m_mean *edI2dz *order2 ) pdisp = ddendv* fdspdz c-----p summation, p is obtained in unit [Pa] -----------------------pges = pideal+phs+ pdisp +pcha pges = pges * (KBOL*t) / 1.d-30 c-----calcul. the derivatives of P to dense-----[d(p)/d(dense)]------c-----the derivatives are usefull for the density-iteration ---------c-----using a gradient (Newton) algorithm----------------------------c-----abbreviations--------------------------------------------------zmsdz = -1.d0 dvzdz = 2.d0*dense*6.d0/pi/D c-----d(p)/d(dense) : ideal gas contribution-------------------------piddz = 6.d0/(pi*d) c-----d(p)/d(dense) : hard sphere contribution-----------------------te1dz = (a-c**3.d0/d**2.d0)/(-zms**2.d0)* zmsdz te2dz = (3.d0*b*c/d - c**3.d0/d**2.d0) * 2.d0/(-zms**3.d0)*zmsdz te3dz = (2.d0 * c**3.d0 /d**2.d0) * (-3.d0)/zms**4.d0 * zmsdz phsdz = dvzdz*(term1+term2+term3) + ddendv*(te1dz+te2dz+te3dz) c-----d(p)/d(dense) : chain term-------------------------------------fchd2 = 0.d0 DO 8 i= 1,ncomp fchd2 = fchd2 + ( x(i)*(1.d0-mseg(i))) & *(dgijd2(i,i)/gij(i,i)-dgijdz(i,i)**2.d0/ gij(i,i)**2.d0) 8 CONTINUE pchadz = ddendv*fchd2 + dvzdz*fchdz c-----d(p)/d(dense) : dispersion contribution------------------------edI1d2 = 0.d0
Modeling and Simulation of Solid-Liquid Equilibrium
71
Appendix B: Program of Solid-Liquid Equilibrium Calculation edI2d2 = 0.d0 DO 10 m=0,6 edI1d2 = edI1d2 + apar(m+1)*DBLE(m+1)*DBLE(m)*dense**DBLE(m-1) edI2d2 = edI2d2 + bpar(m+1)*DBLE(m+1)*DBLE(m)*dense**DBLE(m-1) 10 CONTINUE z00_d3 = m_mean*12.d0*(-dense**2.d0+2.d0*dense+9.d0)/ZMS**6.d0 & + (1.d0-m_mean)* 6.d0 & *(dense**4.d0-16.d0*dense**3.d0+60.d0*dense**2.d0 & -88.d0*dense+46.d0) & /(ZMS*(2.d0-dense))**4.d0 drdpd2 = (-3.d0*z00_d2 -z00_d3*dense)*drdpkt**2.d0 & + 2.d0/drdpkt*(drdpdz)**2.d0 fdspd2 = - 2.d0*PI *6.d0/(PI*D) * edI1d2 * order1 & - PI*6.d0/(PI*D)*order2*m_mean* ( drdpd2 *I2*dense & +2.d0 *drdpdz *edI2dz & + drdpkt *edI2d2 ) pdspdz = dvzdz*fdspdz + ddendv*fdspd2 c-----total d(p)/d(dense) with unit [Pa] ----------------------------pgesdz = (piddz+phsdz+pdspdz+pchadz) * (kbol*T)/1.d-30 RETURN END c******************************************************************** SUBROUTINE DENSITR (pges,pgesdz,gij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) IMPLICIT NONE c-----variables used in the parameter list of subroutine-------------INTEGER nc PARAMETER (nc=20) INTEGER ncomp DOUBLE PRECISION pges,pgesdz,gij(nc, nc) DOUBLE PRECISION x(nc),t,p,parame(nc,25),mseg(nc) DOUBLE PRECISION densta,dense,dap_dx(nc,7),dbp_dx(nc,7) DOUBLE PRECISION order1,order2,apar(7),bpar(7) DOUBLE PRECISION A,B,C,D,AA,BB(nc),CC(nc),DD(nc),EFF(nc), & QQ(nc,nc),PIJK(nc,nc,nc) DOUBLE PRECISION PI, RGAS, NA, KBOL, TAU DOUBLE PRECISION uij(nc,nc),d00ij(nc,nc),d0(nc) c-----local variables------------------------------------------------INTEGER i,start,max_i DOUBLE PRECISION x1,y1,dydx,acc_i c--------------------------------------------------------------------acc_i = 1.d-10 max_i = 50
Modeling and Simulation of Solid-Liquid Equilibrium
72
Appendix B: Program of Solid-Liquid Equilibrium Calculation i =0 x1 = densta 1 CONTINUE i=i+1 dense = x1 CALL P_EOS (pges,pgesdz,gij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) y1 = (pges / p ) - 1.d0 dydx = pgesdz/p x1 = x1 - y1/ dydx IF (x1.GT.0.9d0) x1 = 0.6d0 IF (x1.LE.0.d0) x1 = 1.d-10 start = 1 IF (DABS(y1).LT.acc_i) start = 0 IF (i.GT.max_i) THEN start = 0 write (*,*) 'density iteration failed' c stop ENDIF IF (start.EQ.1) GOTO 1 dense = x1 RETURN END c*********************************************************************** SUBROUTINE F_EOS (fres,gij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) IMPLICIT NONE c-----variables used in the parameter list of subroutine-------------INTEGER nc PARAMETER (nc=20) INTEGER ncomp DOUBLE PRECISION fres,gij(nc,nc) DOUBLE PRECISION x(nc),t,p,parame(nc,25),mseg(nc) DOUBLE PRECISION densta,dense,dap_dx(nc,7),dbp_dx(nc,7) DOUBLE PRECISION order1,order2,apar(7),bpar(7) DOUBLE PRECISION A,B,C,D,AA,BB(nc),CC(nc),DD(nc),EFF(nc), & QQ(nc,nc),PIJK(nc,nc,nc) DOUBLE PRECISION PI, RGAS, NA, KBOL, TAU DOUBLE PRECISION uij(nc,nc),d00ij(nc,nc),d0(nc) c-----local variables------------------------------------------------INTEGER i,m
Modeling and Simulation of Solid-Liquid Equilibrium
73
Appendix B: Program of Solid-Liquid Equilibrium Calculation DOUBLE PRECISION zms,rho,m_mean,term1,term2 DOUBLE PRECISION I1,I2,z00,z00_dz,drdpkt DOUBLE PRECISION fhs,fdisp,fcha c-----abbreviations--------------------------------------------------zms = 1.d0 -dense rho = 6.d0*dense/(PI*D) m_mean = A c-----Helmh. free energy : hard sphere contribution------------------term1 = -(A-C**3.d0/D**2.d0)*DLOG(zms) term2 =(C**3.d0/D**2.d0)/(zms**2.d0)+(3.d0*B*C/D*dense & -(C**3.d0/D**2.d0))/zms fhs = term1 + term2 c-----Helmh. free energy : chain term--------------------------------fcha = 0.d0 DO i = 1,ncomp fcha = fcha + x(i) *(1.d0- mseg(i)) *DLOG(gij(i,i)) END DO c-----Helmh. free energy : dispersion contribution-------------------I1 = 0.d0 I2 = 0.d0 DO m=0,6 I1 = I1 + apar(m+1)*dense**DBLE(m) I2 = I2 + bpar(m+1)*dense**DBLE(m) END DO Z00 = 1.d0 + m_mean*(4.d0*dense-2.d0*dense**2.d0)/(zms**3.d0) & + (1.d0 - m_mean)*(5.d0*dense-2.d0*dense**2.d0) & /(zms*(2.d0-dense)) z00_dz = m_mean*(-2.d0*dense**2.d0+4.d0*dense+4.d0)/zms**4.d0 & + (1.d0 - m_mean)*(dense**2.d0-8.d0*dense+10.d0) & /(zms*(2.d0-dense))**2.d0 drdpkt = 1.d0/(Z00 + dense*z00_dz) fdisp = -2.d0*PI*rho*I1*order1 - PI*rho*drdpkt*m_mean*I2*order2 c-----resid. Helmholtz free energy-----------------------------------fres = fhs + fcha + fdisp RETURN END c********************************************************************* SUBROUTINE PERTPAR (kij, 1 ncomp,x,t,p,parame,mseg,densta,dense,dap_dx,dbp_dx, 2 order1,order2,apar,bpar,A,B,C,D,AA,BB,CC,DD,EFF,QQ,PIJK, 3 PI,RGAS,NA,KBOL,TAU,d0,uij,d00ij) IMPLICIT NONE c-----variables used in the parameter list of subroutine-------------INTEGER nc PARAMETER (nc=20)
Modeling and Simulation of Solid-Liquid Equilibrium
74
Appendix B: Program of Solid-Liquid Equilibrium Calculation INTEGER ncomp DOUBLE PRECISION kij(nc,nc) DOUBLE PRECISION x(nc),t,p,parame(nc,25),mseg(nc) DOUBLE PRECISION densta,dense,dap_dx(nc,7),dbp_dx(nc,7) DOUBLE PRECISION order1,order2,apar(7),bpar(7) DOUBLE PRECISION A,B,C,D,AA,BB(nc),CC(nc),DD(nc),EFF(nc), & QQ(nc,nc),PIJK(nc,nc,nc) DOUBLE PRECISION PI, RGAS, NA, KBOL, TAU DOUBLE PRECISION uij(nc,nc),d00ij(nc,nc),d0(nc) c-----local variables------------------------------------------------INTEGER i,j,k,m DOUBLE PRECISION m_mean DOUBLE PRECISION ap(7,3),bp(7,3) DOUBLE PRECISION d00(nc),u0k(nc) c-----constants------------------------------------------------------PI = 3.14159265359d0 RGAS = 8.31441d0 NA = 6.022045d23 KBOL = RGAS/NA TAU = PI/3.d0/DSQRT(2.d0) c-----dispersion term constants--------------------------------------ap(1,1)= 0.91056314451539d0 ap(1,2)= -0.30840169182720d0 ap(1,3)= -0.09061483509767d0 ap(2,1)= 0.63612814494991d0 ap(2,2)= 0.18605311591713d0 ap(2,3)= 0.45278428063920d0 ap(3,1)= 2.68613478913903d0 ap(3,2)= -2.50300472586548d0 ap(3,3)= 0.59627007280101d0 ap(4,1)= -26.5473624914884d0 ap(4,2)= 21.4197936296668d0 ap(4,3)= -1.72418291311787d0 ap(5,1)= 97.7592087835073d0 ap(5,2)= -65.2558853303492d0 ap(5,3)= -4.13021125311661d0 ap(6,1)= -159.591540865600d0 ap(6,2)= 83.3186804808856d0 ap(6,3)= 13.7766318697211d0 ap(7,1)= 91.2977740839123d0 ap(7,2)= -33.7469229297323d0 ap(7,3)= -8.67284703679646d0 bp(1,1)= 0.72409469413165d0 bp(1,2)= -0.57554980753450d0 bp(1,3)= 0.09768831158356d0
Modeling and Simulation of Solid-Liquid Equilibrium
75
Appendix B: Program of Solid-Liquid Equilibrium Calculation bp(2,1)= 1.11913959304690d0 *2.d0 bp(2,2)= 0.34975477607218d0 *2.d0 bp(2,3)= -0.12787874908050d0 *2.d0 bp(3,1)= -1.33419498282114d0 *3.d0 bp(3,2)= 1.29752244631769d0 *3.d0 bp(3,3)= -3.05195205099107d0 *3.d0 bp(4,1)= -5.25089420371162d0 *4.d0 bp(4,2)= -4.30386791194303d0 *4.d0 bp(4,3)= 5.16051899359931d0 *4.d0 bp(5,1)= 5.37112827253230d0 *5.d0 bp(5,2)= 38.5344528930499d0 *5.d0 bp(5,3)= -7.76088601041257d0 *5.d0 bp(6,1)= 34.4252230677698d0 *6.d0 bp(6,2)= -26.9710769414608d0 *6.d0 bp(6,3)= 15.6044623461691d0 *6.d0 bp(7,1)= -50.8003365888685d0 *7.d0 bp(7,2)= -23.6010990650801d0 *7.d0 bp(7,3)= -4.23812936930675d0 *7.d0 c-----pure component parameters--------------------------------------DO 1 i = 1,ncomp mseg(i) = parame(i,1) d00(i) = parame(i,2) u0k(i) = parame(i,3) parame(i,5) = 0.12d0 d0(i)= d00(i)*(1.d0-parame(i,5)*DEXP(-3.d0*parame(i,3)/t)) 1 CONTINUE c-----combination rules----------------------------------------------DO 2 i = 1, ncomp DO 21 j = 1, ncomp d00ij(i,j)=0.5d0*( d00(i) + d00(j) ) uij(i,j)= (1.d0-kij(i,j))*(u0k(i)*u0k(j))**0.5d0 21 CONTINUE 2 CONTINUE c-----abbreviations--------------------------------------------------A = 0.d0 B = 0.d0 C = 0.d0 D = 0.d0 DO i = 1,ncomp A = A + x(i) * mseg(i) B = B + x(i) * mseg(i) * d0(i) C = C + x(i) * mseg(i) * d0(i)**2.d0 D = D + x(i) * mseg(i) * d0(i)**3.d0 END DO AA = A DO 4 k = 1, ncomp
Modeling and Simulation of Solid-Liquid Equilibrium
76
Appendix B: Program of Solid-Liquid Equilibrium Calculation BB(k) = B*C / D**2.d0 * d0(k)**3.d0 CC(k) = C/D * d0(k) DD(k) = d0(k)**3.d0 / D EFF(k) = B/D * d0(k)**2.d0 DO 41 i = 1, ncomp QQ(k,i)= C/D * (d0(k)*d0(i)) / (d0(k) + d0(i)) DO 411 j=1, ncomp PIJK(i,j,k) = (d0(i)*d0(j))/(d0(i)+d0(j)) / d0(k) 411 CONTINUE 41 CONTINUE 4 CONTINUE c-----dispersion term parameters for chain molecules-----------------m_mean = A DO m=1,7 apar(m) = ap(m,1) + (1.d0-1.d0/m_mean)*ap(m,2) & + (1.d0-1.d0/m_mean)*(1.d0-2.d0/m_mean)*ap(m,3) bpar(m) = bp(m,1) + (1.d0-1.d0/m_mean)*bp(m,2) & + (1.d0-1.d0/m_mean)*(1.d0-2.d0/m_mean)*bp(m,3) END DO c-----derivatives of apar, bpar to mole fraction ( d(apar)/d(x) )----DO k=1,ncomp DO m=1,7 dap_dx(k,m) = mseg(k)/m_mean**2.d0*ap(m,2) & +(3.d0*mseg(k)/m_mean**2.d0 & - 4.d0*mseg(k)/m_mean**3.d0)*ap(m,3) dbp_dx(k,m) = mseg(k)/m_mean**2.d0*bp(m,2) & +(3.d0*mseg(k)/m_mean**2.d0 & - 4.d0*mseg(k)/m_mean**3.d0)*bp(m,3) END DO END DO c-----van der Waals mixing rules for perturbation terms--------------order1 = 0.d0 order2 = 0.d0 DO i = 1,ncomp DO j = 1,ncomp order1 = order1 + x(i)*x(j)* mseg(i)*mseg(j) & *d00ij(i,j)**3.d0 * uij(i,j)/t order2 = order2 + x(i)*x(j)* mseg(i)*mseg(j) & *d00ij(i,j)**3.d0 * (uij(i,j)/t)**2.d0 END DO END DO RETURN END
Modeling and Simulation of Solid-Liquid Equilibrium
77
Appendix C: Saturation data Binary Systems (11): Table 1:
Table 2:
Table 3:
Table 4:
Appendix D: Saturation data for CO2/n-decane/n-octacosane (12).
