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