Tech Manual
User Manual:
Open the PDF directly: View PDF .
Page Count: 47
Download | ![]() |
Open PDF In Browser | View PDF |
DWSIM - Process Simulation, Modeling and Optimization Technical Manual Version 4.1, Revision 0 November 2016 License DWSIM is released under the GNU General Public License (GPL) version 3. Contact Information Author: Daniel Medeiros Website: http://dwsim.inforside.com.br / http://www.sourceforge.net/projects/ dwsim Contact: danielwag@gmail.com Contents 1 Introduction 2 2 Thermodynamic Properties 3 2.1 Phase Equilibria Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2.1.1 Fugacity Coefficient calculation models . . . . . . . . . . . . . . . . . . . 4 2.1.2 Chao-Seader and Grayson-Streed models . . . . . . . . . . . . . . . . . . 7 2.1.3 Calculation models for the liquid phase activity coefficient . . . . . . . . 7 2.1.4 Models for Aqueous Electrolyte Systems . . . . . . . . . . . . . . . . . . 12 2.2 Enthalpy, Entropy and Heat Capacities . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Speed of Sound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Joule-Thomson Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3 Transport Properties 18 3.1 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2 Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3 Surface Tension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4 Isothermal Compressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.5 Bulk Modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4 Thermal Properties 4.1 23 Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5 Aqueous Solution Properties 25 5.1 Mean salt activity coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.2 Osmotic coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.3 Freezing point depression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6 Specialized Models / Property Packages 25 6.1 IAPWS-IF97 Steam Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 6.2 IAPWS-08 Seawater . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 6.3 Black-Oil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.4 FPROPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.5 6.6 CoolProp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Sour Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 7 Reactions 30 7.1 Conversion Reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 7.2 Equilibrium Reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 7.2.1 7.3 Solution method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 Kinetic Reaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 8 Property Estimation Methods 8.1 32 Petroleum Fractions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 8.1.1 Molecular weight . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 8.1.2 Specific Gravity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Contents 8.2 Contents 8.1.3 Critical Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 8.1.4 Acentric Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 8.1.5 Vapor Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 8.1.6 Viscosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 Hypothetical Components . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 9 Other Properties 37 9.1 True Critical Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 9.2 Natural Gas Hydrates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 9.3 9.4 9.2.1 Modified van der Waals and Platteeuw (Parrish and Prausnitz) method . 38 9.2.2 Klauda and Sandler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 9.2.3 Chen and Guo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Petroleum Cold Flow Properties . . . . . . . . . . . . . . . . . . . . . . . . . . 40 9.3.1 Refraction Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 9.3.2 Flash Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 9.3.3 Pour Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 9.3.4 Freezing Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 9.3.5 Cloud Point . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 9.3.6 Cetane Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Chao-Seader Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 References DWSIM - Technical Manual 43 1 1 INTRODUCTION 1 Introduction The thermodynamic calculations are the basis of the simulations in DWSIM. It is important for a process simulator to cover a variety of systems, which can go from simple water handling processes to complex, more elaborated cases, such as simulations of processes in the petroleum/chemical industry. DWSIM is able to model phase equilibria between solids, vapor and up to two liquid phases where possible. External CAPE-OPEN Property Packages may have different equilibrium capabilities. The following sections describe the calculation methods used in DWSIM for the physical and chemical description of the elements of a simulation. DWSIM - Technical Manual 2 2 THERMODYNAMIC PROPERTIES 2 Thermodynamic Properties 2.1 Phase Equilibria Calculation In a mixture which finds itself in a vapor-liquid equilibria state (VLE), the component fugacities are the same in all phases, that is [1]: fiL = fiV (2.1) The fugacity of a component in a mixture depends on temperature, pressure and composition. in order to relate fiV with temperature, pressure and molar fraction, we define the fugacity coefficient, φi = fiV , yi P (2.2) which can be calculated from PVT data, commonly obtained from an equation of state. For a mixture of ideal gases, φi = 1. The fugacity of the i component in the liquid phase is related to the composition of that phase by the activity coefficient γi , which by itself is related to xi and standard-state fugacity fi0 by γi = fiL . xi fi0 (2.3) The standard state fugacity fi0 is the fugacity of the i-th component in the system temperature, i.e. mixture, and in an arbitrary pressure and composition. in DWSIM, the standard-state fugacity of each component is considered to be equal to pure liquid i at the system temperature and pressure. If an Equation of State is used to calculate equilibria, fugacity of the i-th component in the liquid phase is calculated by φi = fiL , xi P (2.4) with the fugacity coefficient φi calculated by the EOS, just like it is for the same component in the vapor phase. The fugacity coefficient of the i-th component either in the liquid or in the vapor phase is obtained from the same Equation of State through the following expressions RT ln φL i ˆ∞ " = ∂P ∂ni T,V,nj VL RT ln φVi ˆ∞ " = V V ∂P ∂ni T,V,nj RT − V RT − V # dV − RT ln Z L , (2.5) # dV − RT ln Z V , (2.6) where the compressibility factor Z is given by ZL = DWSIM - Technical Manual PV L RT (2.7) 3 2.1 Phase Equilibria Calculation 2 ZV = THERMODYNAMIC PROPERTIES PV V RT (2.8) 2.1.1 Fugacity Coefficient calculation models Peng-Robinson Equation of State The Peng-Robinson equation [2] is an cubic Equation of State (characteristic related to the exponent of the molar volume) which relates temperature, pressure and molar volume of a pure component or a mixture of components at equilibrium. The cubic equations are, in fact, the simplest equations capable of representing the behavior of liquid and vapor phases simultaneously. The Peng-Robinson EOS is written in the following form P = RT a(T ) − (V − b) V (V + b) + b(V − b) (2.9) where P pressure R ideal gas universal constant v molar volume b parameter related to hard-sphere volume a parameter related to intermolecular forces For pure substances, the a and b parameters are given by: a(T ) = [1 + (0.37464 + 1.54226ω − 0.26992ω 2 )(1 − Tr(1/2) )]2 0.45724(R2 Tc2 )/Pc b = 0.07780(RTc )/Pc (2.10) (2.11) where ω acentric factor Tc critical temperature Pc critical pressure Tr reduced temperature, T /T c For mixtures, equation 2.9 can be used, replacing a and b by mixture-representative values. a and b mixture values are normally given by the basic mixing rule, am = XX i xi xj q (ai aj )(1 − kij ) (2.12) x i bi (2.13) j bm = X i DWSIM - Technical Manual 4 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES where xi,j molar fraction of the i or j component in the phase (liquid or vapor) ai,j i or j component a constant bi,j i or j component b constant kij binary interaction parameter which characterizes the i-j pair The fugacity coefficient obtained with the Peng-Robinson EOS in given by The binary interaction parameters used by DWSIM are loaded from the databank and can be modified in the Property Package configuration window. P fi bi A ln = (Z − 1) − ln (Z − B) − √ xi P bm 2 2B xk aki bi − am bm k ln Z + 2, 414B Z − 0, 414B , (2.14) where Z in the phase compressibility factor (liquid or vapor) and can be obtained from the equation 2.9, Z 3 − (1 − B)Z 2 + (A − 3B 2 − 2B)Z − (AB − B 2 − 2B) = 0, (2.15) A= am P R2 T 2 (2.16) B= bm P RT (2.17) Z= PV RT (2.18) Soave-Redlich-Kwong Equation of State The Soave-Redlich-Kwong Equation [3] is also a cubic equation of state in volume, P = RT a(T ) − , (V − b) V (V + b) (2.19) The a and b parameters are given by: a(T ) = [1 + (0.48 + 1.574ω − 0.176ω 2 )(1 − Tr(1/2) )]2 0.42747(R2 Tc2 )/Pc b = 0.08664(RTc )/Pc (2.20) (2.21) The equations 2.12 and 2.13 are used to calculate mixture parameters. Fugacity is calculated by ln fi bi A = (Z − 1) − ln (Z − B) − xi P bm B P xk aki bi − am bm k ln Z +B Z (2.22) The phase compressibility factor Z is obtained from the equation 2.19, Z 3 − Z 2 + (A − B − B 2 )Z − AB = 0, DWSIM - Technical Manual (2.23) 5 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES A= am P R2 T 2 (2.24) B= bm P RT (2.25) Z= PV RT (2.26) The equations 2.15 and 2.23, in low temperature and pressure conditions, can provide three roots for Z. In this case, if liquid properties are being calculated, the smallest root is used. If the phase is vapor, the largest root is used. The remaining root has no physical meaning; at high temperatures and pressures (conditions above the pseudocritical point), the equations 2.15 and 2.23 provides only one real root. Peng-Robinson with Volume Translation Volume translation solves the main problem with two-constant EOS’s, poor liquid volumetric predictions. A simple correction term is applied to the EOS-calculated molar volume, v = v EOS − c, (2.27) where v =corrected molar volume, v EOS =EOS-calculated volume, and c =component-specific constant. The shift in volume is actually equivalent to adding a third constant to the EOS but is special because equilibrium conditions are unaltered. It is also shown that multicomponent VLE is unaltered by introducing the volume-shift term c as a mole-fraction average, EOS vL = vL − X xi ci (2.28) Volume translation can be applied to any two-constant cubic equation, thereby eliminating the volumetric defficiency suffered by all two-constant equations [4]. Peng-Robinson-Stryjek-Vera PRSV1 A modification to the attraction term in the Peng-Robinson equation of state published by Stryjek and Vera in 1986 (PRSV) significantly improved the model’s accuracy by introducing an adjustable pure component parameter and by modifying the polynomial fit of the acentric factor. The modification is: κ = κ0 + κ1 1 + Tr0.5 (0.7 − Tr ) (2.29) κ0 = 0.378893 + 1.4897153 ω − 0.17131848 ω 2 + 0.0196554 ω 3 (2.30) where κ1 is an adjustable pure component parameter. Stryjek and Vera published pure component parameters for many compounds of industrial interest in their original journal article. PRSV2 DWSIM - Technical Manual 6 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES A subsequent modification published in 1986 (PRSV2) [5] further improved the model’s accuracy by introducing two additional pure component parameters to the previous attraction term modification. The modification is: κ = κ0 + κ1 + κ2 (κ3 − Tr ) 1 − Tr0 .5 1 + Tr0.5 (0.7 − Tr ) (2.31) κ0 = 0.378893 + 1.4897153 ω − 0.17131848 ω 2 + 0.0196554 ω 3 (2.32) where κ1 , κ2 , and κ3 are adjustable pure component parameters. PRSV2 is particularly advantageous for VLE calculations. While PRSV1 does offer an advantage over the Peng-Robinson model for describing thermodynamic behavior, it is still not accurate enough, in general, for phase equilibrium calculations. The highly non-linear behavior of phase-equilibrium calculation methods tends to amplify what would otherwise be acceptably small errors. It is therefore recommended that PRSV2 be used for equilibrium calculations when applying these models to a design. However, once the equilibrium state has been determined, the phase specific thermodynamic values at equilibrium may be determined by one of several simpler models with a reasonable degree of accuracy. 2.1.2 Chao-Seader and Grayson-Streed models Chao-Seader ([6]) and Grayson-Streed ([7]) are older, semi-empirical models. The GraysonStreed correlation is an extension of the Chao-Seader method with special applicability to hydrogen. In DWSIM, only the equilibrium values produced by these correlations are used in the calculations. The Lee-Kesler method is used to determine the enthalpy and entropy of liquid and vapor phases. Chao Seader Use this method for heavy hydrocarbons, where the pressure is less than 10 342 kPa (1 500 psia) and the temperature is between the range -17.78 C and 260 C. Grayson Streed Recommended for simulating heavy hydrocarbon systems with a high hy- drogen content. 2.1.3 Calculation models for the liquid phase activity coefficient The activity coefficient γ is a factor used in thermodynamics to account for deviations from ideal behaviour in a mixture of chemical substances. In an ideal mixture, the interactions between each pair of chemical species are the same (or more formally, the enthalpy of mixing is zero) and, as a result, properties of the mixtures can be expressed directly in terms of simple concentrations or partial pressures of the substances present. Deviations from ideality are accommodated by modifying the concentration by an activity coefficient. . The activity coefficient is defined as γi = [ ∂(nGE /RT ) ]P,T,nj6=i ∂ni (2.33) where GE represents the excess Gibbs energy of the liquid solution, which is a measure of how far the solution is from ideal behavior. For an ideal solution, γi = 1. Expressions for GE /RT provide values for the activity coefficients. DWSIM - Technical Manual 7 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES The UNIQUAC equation considers g ≡ GE /RT formed by UNIQUAC and UNIFAC models two additive parts, one combinatorial term g C to take into account the size of the molecules, and one residual term g R , which take into account the interactions between molecules: g ≡ gC + gR (2.34) The g C function contains only pure species parameters, while the g R function incorporates two binary parameters for each pair of molecules. For a multicomponent system, gC = X xi ln φi /xi + 5 X i qi xi ln θi /φi (2.35) i and gR = − X X qi xi ln( θj τj i) i j (2.36) where X φi ≡ (xi ri )/( xj rj ) (2.37) j and θi ≡ (xi qi )/( X xj q j ) (2.38) j The i subscript indicates the species, and j is an index that represents all the species, i included. All sums are over all the species. Note that τij 6= τji . When i = j, τii = τjj = 1. In these equations, ri (a relative molecular volume) and qi (a relative molecular surface area) are pure species parameters. The influence of temperature in g enters by means of the τij parameters, which are temperature-dependent: τij = exp(uij − ujj )/RT (2.39) This way, the UNIQUAC parameters are values of (uij − ujj ). An expression for γi is found through the application of the following relation: ln γi = ∂(nGE /RT )/∂ni (P,T,nj6=i ) (2.40) The result is represented by the following equations: ln γi = ln γiC + ln γiR (2.41) ln γiC = 1 − Ji + ln Ji − 5qi (1 − Ji /Li + ln Ji /Li ) (2.42) ln γiR = qi (1 − ln si − X θj τij /sj ) (2.43) j where DWSIM - Technical Manual 8 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES X Ji = ri /( rj xj ) (2.44) j X L = qi /( q j xj ) (2.45) j si = X θl τli (2.46) l Again the i subscript identify the species, j and l are indexes which represent all the species, including i. all sums are over all the species, and τij = 1 for i = j. The parameters values (uij − ujj ) are found by regression of binary VLE/LLE data. The UNIFAC method for the estimation of activity coefficients depends on the concept of that a liquid mixture can be considered a solution of its own molecules. These structural units are called subgroups. The greatest advantage of this method is that a relatively small number of subgroups can be combined to form a very large number of molecules. The activity coefficients do not only depend on the subgroup properties, but also on the interactions between these groups. Similar subgroups are related to a main group, like “CH2”, “OH”, “ACH” etc.; the identification of the main groups are only descriptive. All the subgroups that belongs to the same main group are considered identical with respect to the interaction between groups. Consequently, the parameters which characterize the interactions between the groups are identified by pairs of the main groups. The UNIFAC method is based on the UNIQUAC equation, where the activity coefficients are given by the equation 2.40. When applied to a solution of groups, the equations 2.42 and 2.43 are written in the form: ln γiC = 1 − Ji + ln Ji − 5qi (1 − Ji /Li + ln Ji /Li ) ln γiR = qi (1 − X (θk βik /sk ) − eki lnβik /sk ) (2.47) (2.48) k The parameters Ji e Li are still given by eqs. 2.58 and ??. Furthermore, the following definitions apply: ri = X (i) (2.49) (i) (2.50) νk Rk k qi = X νk Qk k (i) eki = (νk Qk )/qi βik = X emk τmk (2.51) (2.52) m X X θk = ( xi qi eki )/( xj qj ) i DWSIM - Technical Manual (2.53) i 9 2.1 Phase Equilibria Calculation 2 sk = X THERMODYNAMIC PROPERTIES θm τmk (2.54) m si = X θl τli (2.55) l τmk = exp(−amk )/T (2.56) The i subscript identify the species, and j is an index that goes through all the species. The k subscript identify the subgroups, and m is an index that goes through all the subgroups. The (i) parameter νk is the number of the k subgroup in a molecule of the i species. The subgroup parameter values Rk and Qk and the interaction parameters −amk are obtained in the literature. Modified UNIFAC (Dortmund) model The UNIFAC model, despite being widely used in various applications, has some limitations which are, in some way, inherent to the model. Some of these limitations are: 1. UNIFAC is unable to distinguish between some types of isomers. 2. The γ − φ approach limits the use of UNIFAC for applications under the pressure range of 10-15 atm. 3. The temperature is limited within the range of approximately 275-425 K. 4. Non-condensable gases and supercritical components are not included. 5. Proximity effects are not taken into account. 6. The parameters of liquid-liquid equilibrium are different from those of vapor-liquid equilibrium. 7. Polymers are not included. 8. Electrolytes are not included. Some of these limitations can be overcome. The insensitivity of some types of isomers can be eliminated through a careful choice of the groups used to represent the molecules. The fact that the parameters for the liquid-liquid equilibrium are different from those for the vapor-liquid equilibrium seems not to have a theoretical solution at this time. One solution is to use both data from both equiibria to determine the parameters as a modified UNIFAC model. The limitations on the pressure and temperature can be overcome if the UNIFAC model is used with equations of state, which carry with them the dependencies of pressure and temperature. These limitations of the original UNIFAC model have led several authors to propose changes in both combinatorial and the residual parts. To modify the combinatorial part, the basis is the suggestion given by Kikic et al. (1980) in the sense that the Staverman-Guggenheim correction on the original term of Flory-Huggins is very small and can, in most cases, be neglected. As a result, this correction was empirically removed from the UNIFAC model. Among these modifications, the proposed by Gmehling and coworkers [Weidlich and Gmehling, 1986; Weidlich and Gmehling, 1987; Gmehling et al., 1993], known as the model UNIFAC-Dortmund, is one of the most promising. In this model, the combinatorial part of the original UNIFAC is replaced by: DWSIM - Technical Manual 10 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES ln γiC = 1 − Ji + ln Ji − 5qi (1 − Ji /Li + ln Ji /Li ) 3/4 Ji = ri /( X 3/4 rj xj ) (2.57) (2.58) j where the remaining quantities is defined the same way as in the original UNIFAC. Thus, the correction in-Staverman Guggenheim is empirically taken from the template. It is important to note that the in the UNIFAC-Dortmund model, the quantities Rk and Qk are no longer calculated on the volume and surface area of Van der Waals forces, as proposed by Bondi (1968), but are additional adjustable parameters of the model. The residual part is still given by the solution for groups, just as in the original UNIFAC, but now the parameters of group interaction are considered temperature dependent, according to: (0) (1) (2) τmk = exp(−amk + amk T + amk T 2 )/T (2.59) These parameters must be estimated from experimental phase equilibrium data. Gmehling et al. (1993) presented an array of parameters for 45 major groups, adjusted using data from the vapor-liquid equilibrium, excess enthalpies, activity coefficients at infinite dilution and liquidliquid equilibrium. enthalpy and entropy of liquid and vapor. Modified UNIFAC (NIST) model This model [8] is similar to the Modified UNIFAC (Dortmund), with new modified UNIFAC parameters reported for 89 main groups and 984 group–group interactions using critically evaluated phase equilibrium data including vapor–liquid equilibrium (VLE), liquid–liquid equilibrium (LLE), solid–liquid equilibrium (SLE), excess enthalpy (HE), infinite dilution activity coefficient (AINF) and excess heat capacity (CPE) data. A new algorithmic framework for quality assessment of phase equilibrium data was applied for qualifying the consistency of data and screening out possible erroneous data. Substantial improvement over previous versions of UNIFAC is observed due to inclusion of experimental data from recent publications and proper weighting based on a quality assessment procedure. The systems requiring further verification of phase equilibrium data were identified where insufficient number of experimental data points is available or where existing data are conflicting. NRTL model Wilson (1964) presented a model relating g E to the molar fraction, based mainly on molecular considerations, using the concept of local composition. Basically, the concept of local composition states that the composition of the system in the vicinity of a given molecule is not equal to the overall composition of the system, because of intermolecular forces. Wilson’s equation provides a good representation of the Gibbs’ excess free energy for a variety of mixtures, and is particularly useful in solutions of polar compounds or with a tendency to association in apolar solvents, where Van Laar’s equation or Margules’ one are not sufficient. Wilson’s equation has the advantage of being easily extended to multicomponent solutions but has two disadvantages: first, the less important, is that the equations are not applicable to systems where the logarithms of activity coefficients, when plotted as a function of x, show a maximum or a minimum. However, these systems are not common. The second, a little more DWSIM - Technical Manual 11 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES serious, is that the model of Wilson is not able to predict limited miscibility, that is, it is not useful for LLE calculations. Renon and Prausnitz [9] developed the NRTL equation (Non-Random, Two-Liquid) based on the concept of local composition but, unlike Wilson’s model, the NRTL model is applicable to systems of partial miscibility. The model equation is: n P ln γi = τji xj Gji j=1 n P + xk Gki n X j=1 k=1 xj Gij n P xk Gkj τij − k=1 n P τmj xm Gmj m=1 n P xk Gkj , (2.60) k=1 Gij = exp(−τij αij ), (2.61) τij = aij /RT, (2.62) where γi Activity coefficient of component i xi Molar fraction of component i aij Interaction parameter between i-j (aij 6= aji ) (cal/mol) T Temperature (K) αij non-randomness parameter for the i-j pair (αij = αji ) The significance of Gij is similar to Λij from Wilson’s equation, that is, they are charac- teristic energy parameters of the ij interaction. The parameter is related to the non-randomness of the mixture, i.e. that the components in the mixture are not randomly distributed but follow a pattern dictated by the local composition. When it is zero, the mixture is completely random, and the equation is reduced to the two-suffix Margules equation. For ideal or moderately ideal systems, the NRTL model does not offer much advantage over Van Laar and three-suffix Margules, but for strongly non-ideal systems, this equation can provide a good representation of experimental data, although good quality data is necessary to estimate the three required parameters. 2.1.4 Models for Aqueous Electrolyte Systems Revised LIQUAC model (LIQUAC*) In electrolyte systems, different properties, such as mean activity coefficients, osmotic coefficients, boiling point elevations, freezing point depressions and salt solubilities can be calculated using the new electrolyte models like LIQUAC and LIFAC. Common usage scenarios: Ù desalination processes Ù crystallization processes Ù waste water treatment DWSIM - Technical Manual 12 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES In the LIQUAC* model [10], the activity coefficient is calculated by three different terms: R ln γi = ln γ LR + ln γ M + ln γiSR i i (2.63) These three terms, the long range term (LR), the middle range term (MR) and the short range term (SR), consider the different kinds of interactions in electrolyte solutions. The long range term is taken into account by the Debye–Hckel theory as modified by Fowler and Guggenheim to consider different solvents and solvent mixtures. This term takes into account direct charge effects like attraction and repulsion between ions and the formation of a solvate shell in solution and is calculated differently for ions and solvents. The middle range term was developed from the semiempirical Pitzer model and takes into account the indirect charge effects such as interactions between dipoles–dipoles and dipoles–indirect dipoles. The short range term was developed from the corresponding local composition model and takes into account direct neighborhood effects of the compounds in solution. For the calculation of the short range term the part consists of a combinatorial (C) and a residual (R) part. While the combinatorial part takes into account the entropic interactions, i.e. the size and the form of the molecules the residual part considers the enthalpic interactions. Extended UNIQUAC [11] Sander et al. presented in 1986 an extension of the UNIQUAC model by adding a Debye-Hckel term allowing this Extended UNIQUAC model to be used for electrolyte solutions. The model has since been modified and it has proven itself applicable for calculations of vapor-liquid-liquid-solid equilibria and of thermal properties in aqueous solutions containing electrolytes and non-electrolytes. The model is shown in its current form here as it is presented by Thomsen (1997). The extended UNIQUAC model consists of three terms: a combinatorial or entropic term, a residual or enthalpic term and an electrostatic term ex ex Gex = Gex Combinatorial + GResidual + GExtended Debye−H ückel (2.64) The combinatorial and the residual terms are identical to the terms used in the traditional UNIQUAC equation. The electrostatic term corresponds to the extended Debye-Hckel law. The combinatorial, entropic term is independent of temperature and only depends on the relative sizes of the species: X Gex Combinatorial = xi ln RT i φi xi zX qi xi ln − 2 i φi θi (2.65) The two model parameters ri and qi are the volume and surface area parameters for component i. In the classical application of the UNIQUAC model, these parameters are calculated from the properties of non electrolyte molecules. In the Extended UNIQUAC application to multi component electrolyte solutions, this approach gave unsatisfactory results. The volume and surface area parameters were instead considered to be adjustable parameters. The values of these two parameters are determined by fitting to experimental data. Especially thermal property data such as heat of dilution and heat capacity data are efficient for determining the value of the surface area parameter q, because the UNIQUAC contribution to the excess enthalpy and excess heat capacity is proportional to the parameter q. The residual, enthalpic term is dependent on temperature through the parameter ψji : DWSIM - Technical Manual 13 2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES X X Gex Residual xi qi ln θj ψji =− RT i j (2.66) the parameter ψji is given by: uji − uii ψji = exp − T (2.67) uji and uii are interaction energy parameters. The interaction energy parameters are considered symmetrical and temperature dependent in this model uji = uij = u0ij + uTij (T − 298.15) (2.68) The values of the interaction energy parameters and are determined by fitting to experimental data. The combinatorial and the residual terms of the UNIQUAC excess Gibbs energy function are based on the rational, symmetrical activity coefficient convention. The Debye-Hckel electrostatic term however is expressed in terms of the rational, symmetrical convention for water, and the rational, unsymmetrical convention for ions. The electrostatic contributions to the water activity coefficients and the ionic activity coefficients are obtained by partial molar differentiation of the extended Debye-Hckel law excess Gibbs energy term. The term used for water is DH ln γw = 3 σ (x) = 3 x 2 Mw AI 3/2 σ bI 1/2 3 1 1+x− − 2 ln (1 + x) 1+x (2.69) (2.70) In this expression, b = 1.5 (kg/mol)½. The term used for ions is: ln γi∗DH = −Zi2 √ A I √ 1+b I (2.71) Based on table values of the density of pure water, and the relative permittivity of water, εr , the Debye-Hckel parameter A can be approximated in the temperature range 273.15 K < T < 383.15 K by: h i 2 A = 1.131 + 1.335E − 3 (T − 273.15) + 1.164E − 5 (T − 273.15) (2.72) The activity coefficient for water is calculated in the Extended UNIQUAC model by summation of the three terms: C R DH ln γw = ln γw + ln γw + ln γw (2.73) The activity coefficient for ion i is obtained as the rational, unsymmetrical activity coefficient according to the definition of rational unsymmetrical activity coefficients by adding the three contributions: ln γi∗ = ln DWSIM - Technical Manual γiC γiR + ln + ln γi∗DH γiC∞ γiR∞ (2.74) 14 2.2 Enthalpy, Entropy and Heat Capacities 2 THERMODYNAMIC PROPERTIES The rational, unsymmetrical activity coefficient for ions calculated with the Extended UNIQUAC model can be converted to a molal activity coefficient. This is relevant for comparison with experimental data. The temperature dependency of the activity coefficients in the Extended UNIQUAC model is built into the model equations as outlined above. The temperature dependency of the equilibrium constants used in the Extended UNIQUAC model is calculated from the temperature dependency of the Gibbs energies of formation of the species Parameters for water and for the following ions can be found in[11] H+, Na+, K+, NH4+, Cl-, SO42-, HSO4-, NO3-, OH-, CO32-, HCO3-, S2O82-. A significant advantage of the Extended UNIQUAC model compared to models like the Bromley model or the Pitzer model is that temperature dependence is built into the model. This enables the model to also describe thermodynamic properties that are temperature derivatives of the excess Gibbs function, such as heat of mixing and heat capacity. 2.2 Enthalpy, Entropy and Heat Capacities H id values are calculated from Peng-Robinson, Soave-Redlich-Kwong For the cubic equations of state, enthalpy, entropy the ideal gas heat capacity. For and heat capacities are calculated by the departure functions, which relates the phase properties mixtures, a molar average is in the conditions of the mixture with the same mixture property in the ideal gas state.This way, used. The value calculated by the following departure functions are defined [12], the EOS is for the phase, independently of the number of H − H id S − S id = X; =Y RT R components present in the mixture. (2.75) values for X and Y are calculated by the PR and SRK EOS, according to the table 1: Table 1: Enthalpy/Entropy calculation with an EOS H−H id RT PR Z −1− 1 21,5 bRT × ln SRK In DWSIM, Po = 1 atm. h S−S id R da a − T dT × i V +2,414b V +0,414b ln(Z − B) − ln PP0 − 21,5AbRT h i × ln VV +2,414b +0,414b da a − T dT × × ln 1 + Vb Z −1− 1 bRT A ln(Z − B) − ln PP0 − B × ln 1 + B Z T T da a dT da a dT × × Heat capacities are obtained directly from the EOS, by using the following thermodynamic relations: Cp − Cpid ˆV =T ∂2P ∂T 2 dV − ∞ Cp − Cv = −T DWSIM - Technical Manual T (∂P/∂T )2V −R (∂P/∂V )T ∂P 2 ∂T V ∂P ∂V T (2.76) (2.77) 15 2.2 Enthalpy, Entropy and Heat Capacities Lee-Kesler 2 THERMODYNAMIC PROPERTIES Enthalpies, entropies and heat capacities are calculated by the Lee-Kesler model [13] through the following equations: H − H id d2 b2 + 2b3 /Tr + 3b4 /Tr2 c2 − 3c3 /Tr2 + + 3E = Tr Z − 1 − − RTc Tr Vr 2Tr Vr2 5Tr Vr2 (2.78) S − S id + ln R (2.79) P P0 = ln Z − d1 b2 + b3 /Tr2 + 2b4 /Tr3 c1 − 2c3 /Tr3 + + 2E − 2 Vr 2Vr 5Vr5 Cv − Cvid 2 (b3 + 3b4 /Tr ) 3c3 = − 3 2 − 6E 2 R Tr Vr Tr Vr Cp − R Cpid = Cv − R Cvid − 1 − Tr ∂Pr ∂Tr ∂Pr ∂Vr (2.80) 2 Vr (2.81) Tr c4 γ γ E= β + 1 − β + 1 + exp − 2Tr3 γ Vr2 Vr2 Pr Vr D c4 B C Z= =1+ + 2+ 5+ 3 2 Tr Vr Vr Vr Tr Vr An iterative method is required to calculate Vr . The user γ β+ 2 Vr γ exp − 2 Vr (2.82) (2.83) should always watch the values generated by DWSIM in order to detect any issues in the B = b1 − b2 /Tr − b3 /Tr2 − b4 /Tr3 (2.84) C = c1 − c2 /Tr + c3 /Tr3 (2.85) D = d1 + d2 /Tr (2.86) compressibility factors generated by the Lee-Kesler model. Each property must be calculated based in two fluids apart from the main one, one simple and other for reference. For example, for the compressibility factor, Z = Z (0) + ω (r) (0) Z − Z , ω (r) (2.87) where the (0) superscript refers to the simple fluid while the (r) superscript refers to the reference fluid. This way, property calculation by the Lee-Kesler model should follow the sequence below (enthalpy calculation example): 1. Vr and Z (0) are calculated for the simple fluid at the fluid Tr and Pr . using the equation 2.78, and with the constants for the simple fluid, as shown in the table 3, (H − H 0 )/RTc (0) is calculated. This term is (H − H 0 )/RTc . in this calculation, Z in the equation 2.78 is Z (0) . 2. The step 1 is repeated, using the same Tr and Pr , but using the constants for the reference fluid as shown in table 3. With these values, the equation 2.78 allows the calculation of (r) (H − H 0 )/RTc . In this step, Z in the equation 2.78 is Z (r) . 3. Finally, one determines the residual enthalpy for the fluid of interest by DWSIM - Technical Manual 16 2.3 Speed of Sound 2 THERMODYNAMIC PROPERTIES (0) (H − H 0 )/RTc = (H − H 0 )/RTc + (r) (0) ω 0 0 (H − H )/RT − (H − H )/RT , c c ω (r) (2.88) where ω (r) = 0, 3978. Table 3: Constants for the Lee-Kesler model Constant Simple Fluid Reference Fluid b1 0.1181193 0.2026579 b2 0.265728 0.331511 b3 0.154790 0.027655 b4 0.030323 0.203488 c1 0.0236744 0.0313385 c2 0.0186984 0.0503618 c3 0.0 0.016901 c4 0.042724 0.041577 d1 × 104 0155488 0.48736 d2 × 104 0.623689 0.0740336 β 0.65392 1.226 γ 0.060167 0.03754 2.3 Speed of Sound The speed of sound in a given phase is calculated by the following equations: s c= K , ρ (2.89) where: c Speed of sound (m/s) K Bulk Modulus (Pa) ρ Phase Density (kg/m³) 2.4 Joule-Thomson Coefficient In thermodynamics, the Joule–Thomson effect (also known as the Joule–Kelvin effect, Kelvin–Joule effect, or Joule–Thomson expansion) describes the temperature change of a real gas or liquid when it is forced through a valve or porous plug while kept insulated so that no heat is exchanged with the environment. This procedure is called a throttling process or Joule–Thomson process. At room temperature, all gases except hydrogen, helium and neon cool upon expansion by the Joule–Thomson process. The rate of change of temperature with respect to pressure in a Joule–Thomson process is the Joule–Thomson coefficient. The Joule-Thomson coefficient for a given phase is calculated by the following definition: DWSIM - Technical Manual 17 3 µ= ∂T ∂P TRANSPORT PROPERTIES , (2.90) H The JT coefficient is calculated rigorously by the PR and SRK equations of state, while the Goldzberg correlation is used for all other models, 2 0.0048823Tpc 18/Tpr −1 µ= , Ppc Cp γ (2.91) for gases, and µ=− 1 , ρCp (2.92) for liquids. 3 Transport Properties 3.1 Density Liquid Phase Liquid phase density is calculated with the Rackett equation for non-EOS models when experimental data is not available [12], Vs = RTC [1+(1−Tr )2/7 ] Z , PC RA (3.1) where: Vs Saturated molar volume (m³/mol) Tc Critical temperature (K) Pc Critical pressure (Pa) Tr Reduced temperature ZRA Rackett constant of the component (or the mixture) R Ideal Gas constant (8,314 J/[mol.K]) For mixtures, the equation 3.1 becomes If T > Tcm , the Rackett method does not provide a X xi Tci [1+(1−Tr )2/7 ] Vs = R ZRA , Pc i value for Vs and, in this case, DWSIM uses the EOS-generated compressibility with (3.2) Tr = T /Tcm , and factor to calculate the density of the liquid phase. Tcm = XX φi φj Tcij , xi V c φi = P i , xi Vci 1/2 (3.4) 1/2 8 Vci Vcj , Tcij = 3 Tci Tcj 1/3 1/3 Vci + Vcj DWSIM - Technical Manual (3.3) (3.5) 18 3.1 Density 3 TRANSPORT PROPERTIES where: xi Molar fraction Vci Critical volume (m³/mol) If ZRA isn’t available, it is calculated from the component acentric factor, ZRA = 0.2956 − 0.08775ω, (3.6) If the component (or mixture) isn’t saturated, a correction is applied in order to account for the effect of pressure in the volume, β+P , V = Vs 1 − (0.0861488 + 0.0344483ω) ln β + Pvp (3.7) with β P 1/3 2/3 = −1 − 9.070217 (1 − Tr ) + 62.45326 (1 − Tr ) − 135.1102 (1 − Tr ) + 4/3 + exp 4.79594 + 0.250047ω + 1.14188ω 2 (1 − Tr ) , (3.8) where: V Compressed liquid volume (m³/mol) P Pressure (Pa) Pvp Vapor pressure / Bubble point pressure (Pa) Finally, density is calculated from the molar volume by the following relation: ρ= MM , 1000V (3.9) where: ρ Density (kg/m³) V Specific volume of the fluid (m³/mol) MM Liquid phase molecular volume (kg/kmol) For the Ideal Gas Property Vapor Phase Vapor phase density is calculated from the compressiblity factor generated by Package, the compressibility the EOS model, according with the following equation: factor is considered to be equal to 1. ρ= MM P , 1000ZRT (3.10) where: ρ Density (kg/m³) MM Molecular weight of the vapor phase (kg/kmol) DWSIM - Technical Manual 19 3.2 Viscosity 3 P Pressure (Pa) Z Vapor phase compressibility factor R Ideal Gas constant (8,314 J/[mol.K]) T Temperature (K) TRANSPORT PROPERTIES For ideal gases, the same equation is used, with Z = 1. Mixture If there are two phases at system temperature and pressure, the density of the mixture is calculated by the following expression: ρm = fl ρl + fv ρv , (3.11) where: ρm,l,v Density of the mixture / liquid phase / vapor phase (kg/m³) fl,v Volumetric fraction of the liquid phase / vapor phase (kg/kmol) 3.2 Viscosity Liquid Phase When experimental data is not available, liquid phase viscosity is calculated from ! ηL = exp X xi ln ηi , (3.12) i where ηi is the viscosity of each component in the phase, which depends on the temperature and is calculated from experimental data. Dependence of viscosity with the temperature is described in the equation η = exp A + B/T + C ln T + DT E , (3.13) where A, B, C, D and E are experimental coefficients (or generated by DWSIM in the case of pseudocomponents or hypotheticals). Vapor Phase Vapor phase viscosity is calculated in two steps. First, when experimental data is not available, the temperature dependence is given by the Lucas equation [12], ηξ = 0, .807Tr0,618 − 0.357 exp(−0.449Tr ) + 0.34 exp(−4.058Tr ) + 0.018 ξ = 0, 176 Tc M M 3 Pc4 (3.14) 1/6 , (3.15) where η Viscosity (µP ) Tc , Pc Component (or mixture) critical properties DWSIM - Technical Manual 20 3.3 Surface Tension 3 Tr Reduced temperature, T /Tc MM Molecular weight (kg/kmol) TRANSPORT PROPERTIES In the second step, the experimental or calculated viscosity with the Lucas method is corrected to take into account the effect of pressure, by the Jossi-Stiel-Thodos method [12], " (η − η0 ) Tc M M 3 Pc4 #1/4 1/6 +1 = 1.023 + 0.23364ρr + + 0.58533ρ2r − 0.40758ρ3r + 0.093324ρ4r , (3.16) where η, η0 Corrected viscosity / Lucas method calculated viscosity (µP ) Tc , Pc Component critical properties ρr Reduced density, ρ/ρc = V /Vc MM Molecular weight (kg/kmol) If the vapor phase contains more than a component, the viscosity is calculated by the same procedure, but with the required properties calculated by a molar average. 3.3 Surface Tension When experimental data is not available, the liquid phase surface tension is calculated by doing a molar average of the individual component tensions, which are calculated with the Brock-Bird equation [12], σ 2/3 1/3 Pc Tc 11/9 = (0.132αc − 0.279) (1 − Tr ) Tbr ln(Pc /1.01325) αc = 0.9076 1 + , 1 − Tbr (3.17) (3.18) where σ Surface tension (N/m) Tc Critical temperature (K) Pc Critical pressure (Pa) Tbr Reduced normal boiling point, Tb /Tc 3.4 Isothermal Compressibility Isothermal compressiblity of a given phase is calculated following the thermodynamic definition: β=− DWSIM - Technical Manual 1 ∂V V ∂P (3.19) 21 3.5 Bulk Modulus 3 TRANSPORT PROPERTIES The above expression is calculated rigorously by the PR and SRK equations of state. For the other models, a numerical derivative approximation is used. 3.5 Bulk Modulus The Bulk Modulus of a phase is defined as the inverse of the isothermal compressibility: K= DWSIM - Technical Manual 1 β (3.20) 22 4 THERMAL PROPERTIES 4 Thermal Properties 4.1 Thermal Conductivity Liquid Phase When experimental data is not available, the contribution of each component for the thermal conductivity of the liquid phase is calculated by the Latini method [12], λi A(1 − Tr )0.38 = 1/6 (4.1) Tr A A∗ Tb0.38 , M M β Tcγ = (4.2) where A∗ , α, β and γ depend on the nature of the liquid (Saturated Hydrocarbon, Aromatic, Water, etc). The liquid phase thermal conductivity is calculated from the individual values by the Li method [12], λL λij = = PP φi φj λij −1 −1 2(λ−1 i + λj ) xi V c φi = P i , xi Vci (4.3) (4.4) (4.5) where λL liquid phase thermal conductivity (W/[m.K]) Vapor Phase When experimental data is not available, vapor phase thermal conductivity is calculated by the Ely and Hanley method [12], λV = λ∗ + 1000η ∗ 3R 1.32 Cv − , MM 2 (4.6) where λV vapor phase thermal conductivity (W/[m.K]) Cv constant volume heat capacity (J/[mol.K]) λ∗ and η ∗ are defined by: λ∗ = λ0 H H= DWSIM - Technical Manual 16.04E − 3 M M/1000 1/2 (4.7) f 1/2 /h2/3 (4.8) 23 4.1 Thermal Conductivity 4 THERMAL PROPERTIES λ0 = 1944η0 (4.9) f= T0 θ 190.4 (4.10) h= Vc φ 99.2 (4.11) θ = 1 + (ω − 0.011)(0.56553 − 0.86276 ln T + − 0.69852/T + (4.12) φ = 1 + (ω − 0.011)(0.38650 − 1.1617 ln T + ) 0.288/Zc (4.13) If Tr 6 2, T + = Tr . If Tr > 2, T + = 2. h= η ∗ = η0 H η0 = 10−7 Vc φ 99.2 (4.14) M M/1000 16.04E − 3 (4.15) 9 X (n−4)/3 Cn T0 (4.16) n=1 T0 = T /f DWSIM - Technical Manual (4.17) 24 6 SPECIALIZED MODELS / PROPERTY PACKAGES 5 Aqueous Solution Properties 5.1 Mean salt activity coefficient The mean salt activity coefficient is calculated from the activity coefficients of the ions, ν+ ν− ∗,m ∗,m ∗,m ln γ± = ln γ+ ln γ− + ν ν (5.1) In this equation ν+ and v− are the stoichiometric coefficients of the cations and anions of the salt, while ν stands for the sum of these stoichiometric coefficients. With the mean salt activity coefficient the real behavior of a salt can be calculated and it can, e.g. be used for the calculation of electromotoric forces EMF. 5.2 Osmotic coefficient The osmotic coefficient represents the reality of the solvent in electrolyte systems. It is calculated by the logarithmic ratio of the activity and mole fraction of the solvent: Φ=− Ms ln a P i ion mion (5.2) 5.3 Freezing point depression The Schrder and van Laar equation is used: 4m hi ln ai = (1 − (Tm,i /T )) RTm,i (5.3) On the right hand side of the equation a constant factor is achieved, while on the left hand side the activity depends on temperature and composition. For a given composition the freezing point of the system can be calculated iteratively by varying the system temperature. The best way to do this is by starting at the freezing point of the pure solvent. This equation also allows calculating the freezing point of mixed solvent electrolyte systems. 6 Specialized Models / Property Packages 6.1 IAPWS-IF97 Steam Tables Water is used as cooling medium or heat transfer fluid and it plays an important role for air-condition. For conservation or for reaching desired properties, water must be removed from substances (drying). In other cases water must be added (humidification). Also, many chemical reactions take place in hydrous solutions. That’s why a good deal of work has been spent on the investigation and measurement of water properties over the years. Thermodynamic, transport and other properties of water are known better than of any other substance. Accurate data are especially needed for the design of equipment in steam power plants (boilers, turbines, condensers). In this field it’s also important that all parties involved, e.g., companies bidding for equipment in a new steam power plant, base their calculations on the same property data values because small differences may produce appreciable differences. DWSIM - Technical Manual 25 6.2 IAPWS-08 Seawater 6 SPECIALIZED MODELS / PROPERTY PACKAGES A standard for the thermodynamic properties of water over a wide range of temperature and pressure was developed in the 1960’s, the 1967 IFC Formulation for Industrial Use (IFC67). Since 1967 IFC-67 has been used for ”official” calculations such as performance guarantee calculations of power cycles. In 1997, IFC-67 has been replaced by a new formulation, the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam or IAPWS-IF97 for short. IAPWS-IF97 was developed in an international research project coordinated by the International Association for the Properties of Water and Steam (IAPWS). The formulation is described in a paper by W. Wagner et al., ”The IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam,” ASME J. Eng. Gas Turbines and Power, Vol. 122 (2000), pp. 150-182 and several steam table books, among others ASME Steam Tables and Properties of Water and Steam by W. Wagner, Springer 1998. The IAPWS-IF97 divides the thermodynamic surface into five regions: Ù Region 1 for the liquid state from low to high pressures, Ù Region 2 for the vapor and ideal gas state, Ù Region 3 for the thermodynamic state around the critical point, Ù Region 4 for the saturation curve (vapor-liquid equilibrium), Ù Region 5 for high temperatures above 1073.15 K (800 °C) and pressures up to 10 MPa (100 bar). For regions 1, 2, 3 and 5 the authors of IAPWS-IF97 have developed fundamental equations of very high accuracy. Regions 1, 2 and 5 are covered by fundamental equations for the Gibbs free energy g(T,p), region 3 by a fundamental equation for the Helmholtz free energy f(T,v). All thermodynamic properties can then be calculated from these fundamental equations by using the appropriate thermodynamic relations. For region 4 a saturation-pressure equation has been developed. In chemical engineering applications mainly regions 1, 2, 4, and to some extent also region 3 are of interest. The range of validity of these regions, the equations for calculating the thermodynamic properties, and references are summarized in Attachment 1. The equations of the high-temperature region 5 should be looked up in the references. For regions 1 and 2 the thermodynamic properties are given as a function of temperature and pressure, for region 3 as a function of temperature and density. For other independent variables an iterative calculation is usually required. So-called backward equations are provided in IAPWS-IF97 which allow direct calculation of properties as a function of some other sets of variables (see references). Accuracy of the equations and consistency along the region boundaries are more than sufficient for engineering applications. More information about the IAPWS-IF97 Steam Tables formulation can be found at http: //www.thermo.ruhr-uni-bochum.de/en/prof-w-wagner/software/iapws-if97.html?id= 172. 6.2 IAPWS-08 Seawater The IAPWS-08 Seawater Property Package is based on the Seawater-Ice-Air (SIA) library. The Seawater-Ice-Air (SIA) library contains the TEOS-10 subroutines for evaluating a wide range DWSIM - Technical Manual 26 6.3 Black-Oil 6 SPECIALIZED MODELS / PROPERTY PACKAGES of thermodynamic properties of pure water (using IAPWS-95), seawater (using IAPWS-08 for the saline part), ice Ih (using IAPWS-06) and for moist air (using Feistel et al. (2010a), IAPWS (2010)). TEOS-10 is based on a Gibbs function formulation from which all thermodynamic properties of seawater (density, enthalpy, entropy sound speed, etc.) can be derived in a thermodynamically consistent manner. TEOS-10 was adopted by the Intergovernmental Oceanographic Commission at its 25th Assembly in June 2009 to replace EOS-80 as the official description of seawater and ice properties in marine science. A significant change compared with past practice is that TEOS-10 uses Absolute Salinity SA (mass fraction of salt in seawater) as opposed to Practical Salinity SP (which is essentially a measure of the conductivity of seawater) to describe the salt content of seawater. Ocean salinities now have units of g/kg. Absolute Salinity (g/kg) is an SI unit of concentration. The thermodynamic properties of seawater, such as density and enthalpy, are now correctly expressed as functions of Absolute Salinity rather than being functions of the conductivity of seawater. Spatial variations of the composition of seawater mean that Absolute Salinity is not simply proportional to Practical Salinity; TEOS-10 contains procedures to correct for these effects. More information about the SIA library can be found at http://www.teos-10.org/ software.htm. 6.3 Black-Oil When fluids flow from a petroleum reservoir to the surface, pressure and temperature decrease. This affects the gas/liquid equilibrium and the properties of the gas and liquid phases. The black-oil model enables estimation of these, from a minimum of input data. The black-oil model employs 2 pseudo components: 1. Oil which is usually defined as the produced oil, at stock tank conditions. 2. Gas which then is defined as the produced gas at atmospheric standard conditions. The basic modeling assumption is that the gas may dissolve in the liquid hydrocarbon phase, but no oil will dissolve in the gaseous phase. This implies that the composition of the gaseous phase is assumed the same at all pressure and temperatures. The black-oil model assumption is reasonable for mixtures of heavy and light components, like many reservoir oils. The assumption gets worse for mixtures containing much of intermediate components (propane, butane), and is directly misleading for mixtures of light and intermediate components typically found in condensate reservoirs. In DWSIM, a set of models calculates properties for a black oil fluid so it can be used in a process simulation. Black-oil fluids are defined in DWSIM through a minimum set of properties: Ù Oil specific gravity (SGo) at standard conditions Ù Gas specific gravity (SGg) at standard conditions Ù Gas-to-oil ratio (GOR) at standard conditions Ù Basic Sediments and Water (%) DWSIM - Technical Manual 27 6.4 FPROPS 6 SPECIALIZED MODELS / PROPERTY PACKAGES Black oil fluids are defined and created through the Compound Creator tool. If multiple black-oil fluids are added to a simulation, a single fluid is calculated (based on averaged black-oil properties) and used to calculate stream equilibrium conditions and phase properties. The Black-Oil Property Package is a simplified package for quick process calculations involving the black-oil fluids described above. All properties required by the unit operations are calculated based on the set of four basic properties (SGo, SGg, GOR and BSW), so the results of the calculations cannot be considered precise in any way. They can exhibit errors of several orders of magnitude when compared to real-world data. For more accurate petroleum fluid simulations, use the petroleum characterization tools available in DWSIM together with an Equation of State model like Peng-Robinson or SoaveRedlich-Kwong. 6.4 FPROPS FPROPS is a free open-source C-based library for high-accuracy evaluation of thermodynamic properties for a number of pure substances. It makes use of published data for the Helmholtz fundamental equation for those substances. It has been developed by John Pye and others, can function as standalone code, but is also provided with external library code for ASCEND so that it can be used to access these accurate property correlations from within a MODEL. Currently FPROPS supports calculation of the properties of various substances. The properties that can be calculated are internal energy u, entropy s, pressure p, enthalpy h and Helmholtz energy a, as well as various partial derivatives of these with respect to temperature and density. FPROPS reproduces a limited subset of the functionality of commercial programs such as REFPROP, PROPATH, EES, FLUIDCAL, freesteam, SteamTab, and others, but is free open-source software, licensed under the GPL. More information about the FPROPS Property Package can be found in DWSIM’s wiki: http://dwsim.inforside.com.br/wiki/index.php?title=FPROPS_Property_Package 6.5 CoolProp CoolProp [14] is a C++ library that implements pure and pseudo-pure fluid equations of state and transport properties for 114 components. The CoolProp library currently provides thermophysical data for 114 pure and pseudo-pure working fluids. The literature sources for the thermodynamic and transport properties of each fluid are summarized in a table in the Supporting Information available in the above reference. For the CoolProp Property Package, DWSIM implements simple mixing rules based on mass fraction averages in order to calculate mixture enthalpy, entropy, heat capacities, density (and compressibility factor as a consequence). For equilibrium calculations, DWSIM requires values of fugacity coefficients at system’s temperature and pressure. In the CoolProp Property Package, the vapor and liquid phases are considered to be ideal. More information about CoolProp can be found at http://www.coolprop.org. 6.6 Sour Water The Sour Water Property Package is based on the SWEQ model described in the USEPA Report EPA-600/2-80-067: ”A new correlation of NH3, CO2, and H2S volatility data from DWSIM - Technical Manual 28 6.6 Sour Water 6 SPECIALIZED MODELS / PROPERTY PACKAGES aqueous sour water systems”, by Wilson, Grant M., available online at http://nepis.epa. gov/Exe/ZyPDF.cgi?Dockey=9101B309.PDF. In this model, chemical and physical equilibria of NH3, CO2, and H2S in sour water systems including the effects of release by caustic (NaOH) addition are considered. The original model is applicable for temperatures between 20 °C (68 °F) and 140 °C (285 °F), and pressures up to 50 psi. In DWSIM, use of the PR EOS to correct vapour phase non-idealities extends this range but, due to lack of experimental data, exact ranges cannot be specified. The Sour Water Property Package supports calculation of liquid phase chemical equilibria between the following compounds: Ù Water (H2O, ChemSep database) Ù Ammonia (NH3, ChemSep database) Ù Hydrogen sulfide (H2S, ChemSep database) Ù Carbon dioxide (CO2, ChemSep database) Ù Hydron (H+, Electrolytes database) Ù Bicarbonate (HCO3-, Electrolytes database) Ù Carbonate (CO3-2, Electrolytes database) Ù Ammonium (NH4+, Electrolytes database) Ù Carbamate (H2NCOO-, Electrolytes database) Ù Bisulfide (S-2, Electrolytes database) Ù Sulfide (HS-, Electrolytes database) Ù Hydroxide (OH-, Electrolytes database) Ù Sodium Hydroxide (NaOH, Electrolytes database) Ù Sodium (Na+, Electrolytes database) The following reactions in the liquid (aqueous) phase are taken into account by the SWEQ model: Ù CO2 ionization, CO2 + H2O <–> H+ + HCO3Ù Carbonate production, HCO3- <–> CO3-2 + H+ Ù Ammonia ionization, H+ + NH3 <–> NH4+ Ù Carbamate production, HCO3- + NH3 <–> H2NCOO- + H2O Ù H2S ionization, H2S <–> HS- + H+ Ù Sulfide production, HS- <–> S-2 + H+ Ù Water self-ionization, H2O <–> OH- + H+ Ù Sodium Hydroxide dissociation, NaOH <–> OH- + Na+ DWSIM - Technical Manual 29 7 REACTIONS 7 Reactions DWSIM includes support for chemical reactions through the Chemical Reactions Manager. Three types of reactions are available to the user: Conversion, where you must specify the conversion (%) of the limiting reagent as a function of temperature Equilibrium, where you must specify the equilibrium constant (K) as a function of temperature, a constant value or calculated from the Gibbs free energy of reaction (∆G/R). The orders of reaction of the components are obtained from the stoichiometric coefficients. Kinetic, where you should specify the frequency factor (A) and activation energy (E) for the direct reaction (optionally for the reverse reaction), including the orders of reaction (direct and inverse) of each component. For each chemical reaction is necessary to specify the stoichiometric coefficients of the compounds and a base compound, which must be a reactant. This base component is used as reference for calculating the heat of reaction. 7.1 Conversion Reaction In the conversion reaction it is assumed that the user has information regarding the conversion of one of the reactants as a function of temperature. By knowing the conversion and the stoichiometric coefficients, the quantities of the other components in the reaction can be calculated. Considering the following chemical reaction: aA + bB → cC, (7.1) where a, b and c are the stoichiometric coefficients of reactants and product, respectively. A is the limiting reactant and B is in excess. The amount of each component at the end of the reaction can then be calculated from the following stoichiometric relationships: NA = NB = NC = NA0 − NA0 XA b NB0 − NA0 XA a c NC0 + (NA0 XA ) a (7.2) (7.3) (7.4) where NA,B,C are the molar amounts of the components at the end of the reaction, NA0 ,B0 ,C0 are the molar amount of the components at the start of the reaction and XA is the conversion of the base-reactant A. 7.2 Equilibrium Reaction In the equilibrium reaction, the quantity of each component at the equilibrium is related to equilibrium constant by the following relationship: DWSIM - Technical Manual 30 7.3 Kinetic Reaction 7 K= n Y REACTIONS (qj )νj , (7.5) j=1 where K is the equilibrium constant, q is the basis of components (partial pressure in the vapor phase or activity in the liquid phase) ν is the stoichiometric coefficient of component j and n is the number of components in the reaction. The equilibrium constant can be obtained by three different means. One is to consider it a constant, another is considering it as a function of temperature, and finally calculate it automatically from the Gibbs free energy at the temperature of the reaction. The first two methods require user input. 7.2.1 Solution method For each reaction that is occurring in parallel in the system, we can define ξ as the reaction extent, so that the molar amount of each component in the equilibrium is obtained by the following relationship: Nj = Nj0 + X νij ξi , (7.6) i whereξi is the coordinate of the reaction i and νij is the stoichiometric coefficient of the j component at reaction i. Defining the molar fraction of the component i as xj = nj /nt , where nt is the total number of mols, including inerts, whe have the following expression for each reaction i: fi (ξ) = X ln(xi ) − ln(Ki ) = 0, (7.7) i where the system of equations F can be easily solved by Newton-Raphson’s method [15]. 7.3 Kinetic Reaction The kinetic reaction is defined by the parameters of the equation of Arrhenius (frequency factor and activation energy) for both the direct order and for the reverse order. Suppose we have the following kinetic reaction: aA + bB → cC + dD (7.8) The reaction rate for the A component can be defined as rA = k[A][B] − k 0 [C][D] (7.9) where DWSIM - Technical Manual 31 8 PROPERTY ESTIMATION METHODS k = A exp (−E/RT ) (7.10) k0 = A0 exp (−E 0 /RT ) (7.11) The kinetic reactions are used in Plug-Flow Reactors (PFRs) and in Continuous-Stirred Tank Reactors (CSTRs). In them, the relationship between molar concentration and the rate of reaction is given by ˆV FA = FA0 + rA dV, (7.12) where FA is the molar flow of the A component and V is the reactor volume. 8 Property Estimation Methods 8.1 Petroleum Fractions 8.1.1 Molecular weight Riazi and Al Sahhaf method [16] MM = 3/2 1 (6.97996 − ln(1080 − Tb ) , 0.01964 (8.1) where MM Molecular weight (kg/kmol) Tb Boiling point at 1 atm (K) If the specific gravity (SG) is available, the molecular weight is calculated by MM = 42.965[exp(2.097 × 10−4 Tb − 7.78712SG + +2.08476 × 10−3 Tb SG)]Tb1.26007 SG4.98308 (8.2) Winn [17] M M = 0.00005805P EM e2.3776 /d150.9371 , (8.3) where P EM e Mean Boiling Point (K) d15 Specific Gravity @ 60 °F DWSIM - Technical Manual 32 8.1 Petroleum Fractions 8 PROPERTY ESTIMATION METHODS Riazi[17] MM 42.965 exp(0.0002097P EM e − 7.78d15 + 0.00208476 × P EM e × d15) × = ×P EM e1.26007 d154.98308 (8.4) Lee-Kesler[17] t1 t2 = −12272.6 + 9486.4d15 + (8.3741 − 5.9917d15)P EM e = (1 − 0.77084d15 − 0.02058d15 ) × ×(0.7465 − 222.466/P EM e) × 107 /P EM e t3 = (8.6) 2 (1 − 0.80882d15 − 0.02226d15 ) × ×(0.3228 − 17.335/P EM e) × 1012 /P EM e3 MM (8.5) 2 = t1 + t2 + t3 (8.7) (8.8) Farah MM = exp(6.8117 + 1.3372A − 3.6283B) MM = exp(4.0397 + 0.1362A − 0.3406B − 0.9988d15 + 0.0039P EM e), (8.9) (8.10) where A, B Walther-ASTM equation parameters for viscosity calculation 8.1.2 Specific Gravity Riazi e Al Sahhaf [16] SG = 1.07 − exp(3.56073 − 2.93886M M 0.1 ), (8.11) where SG Specific Gravity MM Molecular weight (kg/kmol) 8.1.3 Critical Properties Lee-Kesler [16] Tc = 189.8 + 450.6SG + (0.4244 + 0.1174SG)Tb + (0.1441 − 1.0069SG)105 /Tb ln Pc = (8.12) 5.689 − 0.0566/SG − (0.43639 + 4.1216/SG + 0.21343/SG2 ) × ×10−3 Tb + (0.47579 + 1.182/SG + 0.15302/SG2 ) × 10−6 × Tb2 − −(2.4505 + 9.9099/SG2 ) × 10−10 × Tb3 , DWSIM - Technical Manual (8.13) 33 8.1 Petroleum Fractions 8 PROPERTY ESTIMATION METHODS where Tb NBP (K) Tc Critical temperature (K) Pc Critical pressure (bar) Farah Tc = 731.968 + 291.952A − 704.998B (8.14) Tc = 104.0061 + 38.75A − 41.6097B + 0.7831P EM e (8.15) Tc = 196.793 + 90.205A − 221.051B + 309.534d15 + 0.524P EM e (8.16) Pc = exp(20.0056 − 9.8758 ln(A) + 12.2326 ln(B)) (8.17) Pc = exp(11.2037 − 0.5484A + 1.9242B + 510.1272/P EM e) (8.18) Pc = exp(28.7605 + 0.7158 ln(A) − 0.2796 ln(B) + 2.3129 ln(d15) − 2.4027 ln(P EM (8.19) e)) Riazi-Daubert[17] Tc = 9.5233 exp(−0.0009314P EM e − 0.544442d15 + 0.00064791 × P EM e × d15) × ×P EM e0.81067 d150.53691 Pc = (8.20) 31958000000 exp(−0.008505P EM e − 4.8014d15 + 0.005749 × P EM e × d15) × ×P EM e−0.4844 d154.0846 (8.21) Riazi[17] Tc 35.9413 exp(−0.00069P EM e − 1.4442d15 + 0.000491 × P EM e × d15) × = ×P EM e0.7293 d151.2771 (8.22) 8.1.4 Acentric Factor Lee-Kesler method [16] ω= Pc 6 − ln 1.10325 − 5.92714 + 6.09648/Tbr + 1.28862 ln Tbr − 0.169347Tbr 6 15.2518 − 15.6875/Tbr − 13.472 ln Tbr + 0.43577Tbr (8.23) Korsten[17] ω = 0.5899 × ((P EM V /Tc )1.3 )/(1 − (P EM V /Tc )1.3 ) × log(Pc /101325) − 1 DWSIM - Technical Manual (8.24) 34 8.1 Petroleum Fractions 8 PROPERTY ESTIMATION METHODS 8.1.5 Vapor Pressure Lee-Kesler method[16] ln Prpv = 6 5.92714 − 6.09648/Tbr − 1.28862 ln Tbr + 0.169347Tbr + (8.25) 6 +ω(15.2518 − 15.6875/Tbr − 13.4721 ln Tbr + 0.43577Tbr ), where Prpv Reduced vapor pressure, P pv /Pc Tbr Reduced NBP, Tb /Tc ω Acentric factor 8.1.6 Viscosity Letsou-Stiel [12] η = ξ0 = ξ1 = ξ = ξ0 + ξ1 ξ 2.648 − 3.725Tr + 1.309Tr2 7.425 − 13.39Tr + 5.933Tr2 1/6 Tc 176 M M 3 Pc4 (8.26) (8.27) (8.28) (8.29) where η Viscosity (Pa.s) Pc Critical pressure (bar) Tr Reduced temperature, T /Tc MM Molecular weight (kg/kmol) Abbott[17] t1 log v100 = 4.39371 − 1.94733Kw + 0.12769Kw2 + 0.00032629AP I 2 − 0.0118246KwAP I + +(0.171617Kw2 + 10.9943AP I + 0.0950663AP I 2 − 0.869218KwAP I t1 , = AP I + 50.3642 − 4.78231Kw t2 log v210 = = (8.30) (8.31) −0.463634 − 0.166532AP I + 0.000513447AP I 2 − 0.00848995AP IKw + +(0.080325Kw + 1.24899AP I + 0.19768AP I 2 t2 , AP I + 26.786 − 2.6296Kw (8.32) (8.33) where DWSIM - Technical Manual 35 8.2 Hypothetical Components v100 Viscosity at 100 °F (cSt) v210 Viscosity at 210 °F (cSt) Kw Watson characterization factor AP I Oil API degree 8 PROPERTY ESTIMATION METHODS 8.2 Hypothetical Components The majority of properties of the hypothetical components is calculated, when necessary, using the group contribution methods, with the UNIFAC structure of the hypo as the basis of calculation. The table 4 lists the properties and their calculation methods. Table 4: Hypo calculation methods. Property Symbol Method Critical temperature Tc Joback [12] Critical pressure Pc Joback [12] Critical volume Vc Joback [12] Normal boiling point Tb Joback [12] Vapor pressure P pv Lee-Kesler (Eq. 8.25) Acentric factor ω Lee-Kesler (Eq. 8.23) Vaporization enthalpy ∆Hvap Vetere [12] Ideal gas heat capacity Cpgi ∆Hf298 Harrison-Seaton [18] Ideal gas enthalpy of formation DWSIM - Technical Manual Marrero-Gani [19] 36 9 OTHER PROPERTIES 9 Other Properties 9.1 True Critical Point The Gibbs criteria for the true critical point of a mixture of n components may be expressed of various forms, but the most convenient when using a pressure explicit cubic equation of state is A11 A12 A21 .. . A22 An1 ... ... A11 A12 ... A21 .. . A22 An−1,1 ... ... An−1,n ∂L ∂n1 ... ... ∂L ∂nn L= M= ... A1n =0 (9.1) Ann A1n = 0, (9.2) where A12 = ∂2A ∂n1 ∂n2 (9.3) T,V All the A terms in the equations 9.1 and 9.2 are the second derivatives of the total Helmholtz energy A with respect to mols and constant T and V . The determinants expressed by 9.1 and 9.2 are simultaneously solved for the critical volume and temperature. The critical pressure is then found by using the original EOS. DWSIM utilizes the method described by Heidemann and Khalil [20] for the true critical point calculation using the Peng-Robinson and Soave-Redlich-Kwong equations of state. 9.2 Natural Gas Hydrates The models for natural gas hydrates equilibrium calculations are mostly based in statistical thermodynamics to predict in which temperature and pressure conditions there will be formation or dissociation of hydrates. In these conditions, fwi = fwH , (9.4) that is, the fugacity of water in hydrate is the same as in the water in any other phase present at equilibria. The difference in the models present in DWSIM is mainly in the way that water fugacity in the hydrate phase is calculated. In the modified van der Waals and Platteeuw model, the isofugacity criteria is used indirectly through chemical potentials, which must also be equal in the equilibria: µiw = µH w (9.5) remembering that DWSIM - Technical Manual 37 9.2 Natural Gas Hydrates 9 OTHER PROPERTIES fi = xi P exp((µi − µgi i )/RT ). (9.6) 9.2.1 Modified van der Waals and Platteeuw (Parrish and Prausnitz) method The classic model for determination of equilibrium pressures and temperatures was developed by van der Waals and Platteeuw. This model was later extended by Parrish and Prausnitz [21] to take into account multiple ”guests” in the hydrate structures. The condition of equilibrium used in the vdwP model is the equality of the chemical potential of water in the hydrate phase and in the other phases, which can be liquid, solid or both. Chemical potential of water in the hydrate phase In the modified var der Waals method, the chemical potential of water in the hydrate phase is calculated by: β µH w = µw + RT X νm ln(1 − m X θmj ), (9.7) j where µβw is the chemical potential of water in the empty hydrate lattice (something like an ”ideal” chemical potential) and νm is the number of m cavities by water molecule in the lattice. The fraction of cavities m-type cavities occupied by the gaseous component l is θml = (Cml fi )/((1 + X Cmj fj )), (9.8) j where Cmj is the Langmuir constant and f i is the fugacity of the gaseous component l. The Langmuir constant takes into account the interactions between the gas and the molecules of water in the cavities. Using the Lennard-Jones-Devonshire cell theory, van der Waals e Platteeuw showed that the Langmuir constant can be given by ˆ ∞ exp[(−w(r))/kT ]r2 dr, C(T ) = 4π/kT (9.9) 0 where T is the absolute temperature, k is the Boltzmann constant and w(r) is the spherically symmetric potential which is a function of the cell radius, the coordination number and the nature of the gas-water interaction. In this method, the Kihara potential with a spherical core is used, w(r) = 2ze[σ 12 /(R11 r)(δ 10 + a/Rδ 11 ) − σ 6 /(R5 r)(δ 4 + a/Rδ 5 )], (9.10) δ N = [(1 − r/R − a/R)−N ) − (1 + r/R − a/R)−N )]/N, (9.11) where N is equal to 4, 5, 10 or 11; z and R are, respectively, the coordination number and the cavity cell radius. Supported hydrate formers CH4, C2H6, C3H8, iC4H10, H2S, N2 and CO2. Model applicability range DWSIM - Technical Manual Temperature: 211 to 303 K; Pressure: 1 to 600 atm. 38 9.2 Natural Gas Hydrates 9 OTHER PROPERTIES 9.2.2 Klauda and Sandler The model proposed by Klauda and Sandler [22] uses spherically symmetric Kihara potentials determined from viscosity data and the second virial coefficient, in opposition to the traditional models which adjust these parameters to experimental hydrate data. In general, this method predicts hydrate formation data more precisely than the other models. Fugacity of water in the hydrate phase fwH = exp(Aβg ln T + (Bgβ )/T + 2, 7789 + Dgβ T ) × exp Vwβ [P − exp(Aβg ln T + (Bgβ )/T + 2, 7789 + Dgβ T )]/RT × X X X exp[ νm ln(1 − (Cml fl )/(1 + Cmj fj ))] m (9.12) j The A, B and D constants are specific for each hydrate former and represent the vapor pressure of the component in the empty hydrate lattice. Vwβ represents the basic hydrate molar volume (without the presence of guests) and the Langmuir constant (C ) is calculated by the following equation: ˆ R−a exp[(−w(r))/kT ]r2 dr C(T ) = 4π/kT (9.13) 0 In the Klauda e Sandler method the spherically symmetric Kihara potential is also used, w(r) = 2ze[σ 12 /(R11 r)(δ 10 + a/Rδ 11 ) − σ 6 /(R5 r)(δ 4 + a/Rδ 5 )] (9.14) δ N = [(1 − r/R − a/R)−N ) − (1 + r/R − a/R)−N )]/N (9.15) with a modifications in the potential to include the effects of the second and third cell layers, w(r) = w(r)[1] ) + w(r)[2] ) + w(r)[3] ). (9.16) Supported hydrate formers CH4, C2H6, C3H8, iC4H10, H2S, N2 and CO2. Model applicability range Temperature: 150 to 320 K; Pressure: 1 to 7000 atm 9.2.3 Chen and Guo Chen and Guo [23] developed a model based in a formation mechanism based in two steps, the first being a quasi-chemical reaction to form the ”basic hydrate” and the second as being a small gas absorption process in the linking cages of the basic hydrate. The results showed that this model is capable of predict hydrate formation conditions for pure gases and mixtures. Fugacity of water in the hydrate phase In the Chen and Guo model, a different approximation is used for the equilibrium condition. Here the equilibrium is verified by means of an isofugacity criteria of the hydrate formers in the hydrate and vapor phase. The fugacity of the component in the vapor phase is calculated by: DWSIM - Technical Manual 39 9.3 Petroleum Cold Flow Properties 9 OTHER PROPERTIES fiH = fi0 (1 − θi )α , (9.17) where α depends on the structure of the hydrate formed, and fi0 = f 0 (P )f 0 (T )f 0 (xw γw ), (9.18) f 0 (P ) = exp(βP/T ), (9.19) 0 0 0 f 0 (T ) = A exp(B /(T − C )), (9.20) f 0 (xw γw ) = (xw γw )(−1/λ2 ) , (9.21) where β and λ2 depend on the structure of the hydrate formed and A’, B’ and C’ depends on the hydrate former. xw and γm are, respectively, the water molar fraction and activity coefficient in the liquid phase. In the Chen and Guo model, the Langmuir constants are calculated with an Antoine-type equation with parameters obtained from experimental data, for a limited range of temperature: C(T ) = X exp(Y /(T − Z)) (9.22) Supported hydrate formers CH4, C2H6, C3H8, iC4H10, H2S, N2, CO2 and nC4H10. Model applicability range Temperature: 259 to 304 K, Pressure: 1 to 700 atm. 9.3 Petroleum Cold Flow Properties 9.3.1 Refraction Index API Procedure 2B5.1 I = 0.02266 exp(0.0003905 × (1.8M eABP ) + 2.468SG − 0.0005704(1.8M eABP ) × SG) × ×(1.8M eABP )0.0572 SG−0.72 (9.23) r= 1 + 2I 1−I 1/2 (9.24) where r Refraction Index SG Specific Gravity M eABP Mean Averaged Boiling Point (K) DWSIM - Technical Manual 40 9.3 Petroleum Cold Flow Properties 9 OTHER PROPERTIES 9.3.2 Flash Point API Procedure 2B7.1 P F = {[0.69 × ((t10AST M − 273.15) × 9/5 + 32) − 118.2] − 32} × 5/9 + 273.15 (9.25) where PF Flash Point (K) t10AST M ASTM D93 10% vaporized temperature (K) 9.3.3 Pour Point API Procedure 2B8.1 P F L = [753 + 136(1 − exp(−0.15v100 )) − 572SG + 0.0512v100 + 0.139(1.8M eABP )] /1.8 (9.26) where PFL Pour Point (K) v100 Viscosity @ 100 °F (cSt) 9.3.4 Freezing Point API Procedure 2B11.1 P C = −2390.42 + 1826SG + 122.49K − 0.135 × 1.8 × M eABP (9.27) where PC Freezing Point (K) K API characterization factor (API) 9.3.5 Cloud Point API Procedure 2B12.1 h i 0.315 −0.133SG) P N = 10(−7.41+5.49 log(1.8M eABP )−0.712×(1.8M eABP ) /1.8 (9.28) where PN Cloud Point (K) 9.3.6 Cetane Index API Procedure 2B13.1 IC = 415.26 − 7.673AP I + 0.186 × (1.8M eABP − 458.67) + 3.503AP I × × log(1.8M eABP − 458.67) − 193.816 × log(1.8M eABP − 458.67) DWSIM - Technical Manual (9.29) 41 9.4 Chao-Seader Parameters 9 OTHER PROPERTIES where IC Cetane Index AP I API degree of the oil 9.4 Chao-Seader Parameters The Chao-Seader parameters needed by the CS/GS models are the Modified Acentric Factor, Solubility Parameter and Liquid Molar Volume. When absent, the Modified Acentric Factor is taken as the normal acentric factor, either read from the databases or calculated by using the methods described before in this document. The Solubility Parameter is given by δ= ∆Hv − RT VL 1/2 (9.30) where ∆Hv Molar Heat of Vaporization VL Liquid Molar Volume at 20 °C DWSIM - Technical Manual 42 References References References [1] J. Smith, Intro to Chemical Engineering Thermodynamics. New York: McGraw-Hill Companies, 1996. [2] D.-Y. Peng and D. B. Robinson, “A new two-constant equation of state,” Industrial and Engineering Chemistry Fundamentals, vol. 15, no. 1, pp. 59–64, 1976. [3] G. Soave, “Equilibrium constants from a modified redlich-kwong equation of state,” Chemical Engineering Science, vol. 27, no. 6, pp. 1197–1203, June 1972. [4] C. H. Whitson and M. R. Brule, Phase Behavior (SPE Monograph Series Vol. 20), S. of Petroleum Engineers, Ed. Society of Petroleum Engineers, 2000. [5] R. Stryjek and J. H. Vera, “Prsv2: A cubic equation of state for accurate vaporliquid equilibria calculations,” The Canadian Journal of Chemical Engineering, vol. 64, pp. 820– 826, 1986. [6] K. C. Chao and J. D. Seader, “A general correlation of vapor-liquid equilibria in hydrocarbon mixtures,” AICHE Journal, vol. 7, pp. 599–605, 1961. [7] H. G. Grayson and C. W. Streed, “Vapor-liquid equilibria for high temperature, high pressure hydrogen-hydrocarbon systems,” N/A, vol. VII, 1963. [8] J. W. Kang, V. Diky, and M. Frenkel, “New modified {UNIFAC} parameters using critically evaluated phase equilibrium data,” Fluid Phase Equilibria, vol. 388, no. 0, pp. 128 – 141, 2015. [Online]. Available: http://www.sciencedirect.com/science/article/pii/ S0378381214007353 [9] H. Renon and J. M. Prausnitz, “Local compositions in thermodynamic excess functions for liquid mixtures,” AICHE Journal, vol. 14, pp. 135–144, 1968. [10] A. Mohs and J. Gmehling, “A revised {LIQUAC} and {LIFAC} model (liquac*/lifac*) for the prediction of properties of electrolyte containing solutions,” Fluid Phase Equilibria, vol. 337, no. 0, pp. 311 – 322, 2013. [Online]. Available: http: //www.sciencedirect.com/science/article/pii/S0378381212004426 [11] K. Thomsen, “Aqueous electrolytes: model parameters and process simulation,” Ph.D. dissertation, Department of Chemical Engineering, Technical University of Denmark, 1997. [12] R. Reid, The Properties of Gases and Liquids. New York: McGraw-Hill, 1987. [13] B. I. Lee and M. G. Kesler, “A generalized thermodynamic correlation based on threeparameter corresponding states,” AIChE Journal, vol. 21, no. 3, pp. 510–527, 1975. [14] I. H. Bell, J. Wronski, S. Quoilin, and V. Lemort, “Pure and pseudo-pure fluid thermophysical property evaluation and the open-source thermophysical property library coolprop,” Industrial & Engineering Chemistry Research, vol. 53, no. 6, pp. 2498–2508, 2014. [Online]. Available: http://pubs.acs.org/doi/abs/10.1021/ie4033999 [15] M. Michelsen and J. Mollerup, Thermodynamic Models: Fundamentals and Computational Aspects. Denmark: Tie-Line Publications, 2007. DWSIM - Technical Manual 43 References References [16] M. R. Riazi, H. A. Al-Adwani, and A. Bishara, “The impact of characterization methods on properties of reservoir fluids and crude oils: options and restrictions,” Journal of Petroleum Science and Engineering, vol. 42, no. 2-4, pp. 195–207, April 2004. [17] M.-R. Riazi, Characterization and properties of petroleum fractions, M.-R. Riazi, Ed. ASTM, 2005. [18] B. K. Harrison and W. H. Seaton, “Solution to missing group problem for estimation of ideal gas heat capacities,” Industrial and Engineering Chemistry Research, vol. 27, no. 8, pp. 1536–1540, 1988. [19] J. Marrero and R. Gani, “Group-contribution based estimation of pure component properties,” Fluid Phase Equilibria, vol. 183-184, pp. 183–208, July 2001. [20] R. A. Heidemann and A. M. Khalil, “The calculation of critical points,” AIChE Journal, vol. 26, no. 5, pp. 769–779, 1980. [21] W. R. Parrish and J. M. Prausnitz, “Dissociation pressures of gas hydrates formed by gas mixtures,” Industrial and Engineering Chemistry Process Design and Development, vol. 11, no. 1, pp. 26–35, 1972. [22] J. Klauda and S. Sandler, “A fugacity model for gas hydrate phase equilibria,” Industrial and Engineering Chemistry Research, vol. 39, no. 9, pp. 3377–3386, 2000. [23] G.-J. Chen and T.-M. Guo, “A new approach to gas hydrate modelling,” Chemical Engineering Journal, vol. 71, no. 2, pp. 145–151, December 1998, url. DWSIM - Technical Manual 44
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : No Page Count : 47 Page Mode : UseOutlines Author : Title : Subject : Creator : LaTeX with hyperref package Producer : pdfTeX-1.40.17 Create Date : 2016:12:01 19:29:07-04:00 Modify Date : 2016:12:01 19:29:07-04:00 Trapped : False PTEX Fullbanner : This is MiKTeX-pdfTeX 2.9.6050 (1.40.17)EXIF Metadata provided by EXIF.tools