Tech Manual

User Manual:

Open the PDF directly: View PDF PDF.
Page Count: 47

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 SpeedofSound ................................... 17
2.4 Joule-Thomson Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3 Transport Properties 18
3.1 Density........................................ 18
3.2 Viscosity....................................... 20
3.3 SurfaceTension ................................... 21
3.4 Isothermal Compressibility . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.5 BulkModulus .................................... 22
4 Thermal Properties 23
4.1 ThermalConductivity................................ 23
5 Aqueous Solution Properties 25
5.1 Mean salt activity coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
5.2 Osmoticcoefficient ................................. 25
5.3 Freezing point depression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6 Specialized Models / Property Packages 25
6.1 IAPWS-IF97 Steam Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
6.2 IAPWS-08Seawater................................. 26
6.3 Black-Oil....................................... 27
6.4 FPROPS....................................... 28
6.5 CoolProp....................................... 28
6.6 SourWater ..................................... 28
7 Reactions 30
7.1 ConversionReaction................................. 30
7.2 EquilibriumReaction ................................ 30
7.2.1 Solutionmethod............................... 31
7.3 KineticReaction................................... 31
8 Property Estimation Methods 32
8.1 PetroleumFractions................................. 32
8.1.1 Molecularweight .............................. 32
8.1.2 SpecificGravity ............................... 33
Contents Contents
8.1.3 CriticalProperties.............................. 33
8.1.4 AcentricFactor ............................... 34
8.1.5 VaporPressure ............................... 35
8.1.6 Viscosity................................... 35
8.2 HypotheticalComponents.............................. 36
9 Other Properties 37
9.1 TrueCriticalPoint.................................. 37
9.2 NaturalGasHydrates ................................ 37
9.2.1 Modified van der Waals and Platteeuw (Parrish and Prausnitz) method . 38
9.2.2 KlaudaandSandler............................. 39
9.2.3 ChenandGuo................................ 39
9.3 Petroleum Cold Flow Properties . . . . . . . . . . . . . . . . . . . . . . . . . . 40
9.3.1 RefractionIndex............................... 40
9.3.2 FlashPoint ................................. 41
9.3.3 PourPoint.................................. 41
9.3.4 FreezingPoint................................ 41
9.3.5 CloudPoint ................................. 41
9.3.6 CetaneIndex................................. 41
9.4 Chao-Seader Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
References 43
DWSIM - Technical Manual 1
1 INTRODUCTION
1 Introduction
The thermodynamic calculations are the basis of the simulations in DWSIM. It is impor-
tant 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 capa-
bilities.
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]:
fL
i=fV
i(2.1)
The fugacity of a component in a mixture depends on temperature, pressure and composi-
tion. in order to relate fV
iwith temperature, pressure and molar fraction, we define the fugacity
coefficient,
φi=fV
i
yiP,(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 icomponent in the liquid phase is related to the composition of that
phase by the activity coefficient γi, which by itself is related to xiand standard-state fugacity
f0
iby
γi=fL
i
xif0
i
.(2.3)
The standard state fugacity f0
iis the fugacity of the i-th component in the system temper-
ature, 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 iat 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=fL
i
xiP,(2.4)
with the fugacity coefficient φicalculated 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=
∞
ˆ
VL"∂P
∂niT,V,nj−RT
V#dV −RT ln ZL,(2.5)
RT ln φV
i=
∞
ˆ
VV"∂P
∂niT,V,nj−RT
V#dV −RT ln ZV,(2.6)
where the compressibility factor Zis given by
ZL=P V L
RT (2.7)
DWSIM - Technical Manual 3
2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES
ZV=P V 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
(V−b)−a(T)
V(V+b) + b(V−b)(2.9)
where
Ppressure
Rideal gas universal constant
vmolar volume
bparameter related to hard-sphere volume
aparameter related to intermolecular forces
For pure substances, the aand bparameters are given by:
a(T) = [1 + (0.37464 + 1.54226ω−0.26992ω2)(1 −T(1/2)
r)]20.45724(R2T2
c)/Pc(2.10)
b= 0.07780(RTc)/Pc(2.11)
where
ωacentric factor
Tccritical temperature
Pccritical pressure
Trreduced temperature, T/T c
For mixtures, equation 2.9 can be used, replacing aand bby mixture-representative values.
aand bmixture values are normally given by the basic mixing rule,
am=X
iX
j
xixjq(aiaj)(1 −kij )(2.12)
bm=X
i
xibi(2.13)
DWSIM - Technical Manual 4
2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES
where
xi,j molar fraction of the ior jcomponent in the phase (liquid or vapor)
ai,j ior jcomponent aconstant
bi,j ior jcomponent bconstant
kij binary interaction parameter which characterizes the i-j pair
The binary interaction
parameters used by DWSIM are
loaded from the databank and
can be modified in the Property
Package configuration window.
The fugacity coefficient obtained with the Peng-Robinson EOS in given by
ln fi
xiP=bi
bm
(Z−1) −ln (Z−B)−A
2√2BPkxkaki
am−bi
bmln Z+ 2,414B
Z−0,414B,(2.14)
where Zin the phase compressibility factor (liquid or vapor) and can be obtained from the
equation 2.9,
Z3−(1 −B)Z2+ (A−3B2−2B)Z−(AB −B2−2B) = 0,(2.15)
A=amP
R2T2(2.16)
B=bmP
RT (2.17)
Z=P V
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
(V−b)−a(T)
V(V+b),(2.19)
The aand bparameters are given by:
a(T) = [1 + (0.48 + 1.574ω−0.176ω2)(1 −T(1/2)
r)]20.42747(R2T2
c)/Pc(2.20)
b= 0.08664(RTc)/Pc(2.21)
The equations 2.12 and 2.13 are used to calculate mixture parameters. Fugacity is calculated
by
ln fi
xiP=bi
bm
(Z−1) −ln (Z−B)−A
BPkxkaki
am−bi
bmln Z+B
Z(2.22)
The phase compressibility factor Zis obtained from the equation 2.19,
Z3−Z2+ (A−B−B2)Z−AB = 0,(2.23)
DWSIM - Technical Manual 5
2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES
A=amP
R2T2(2.24)
B=bmP
RT (2.25)
Z=P V
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=vEOS −c, (2.27)
where v=corrected molar volume, vEOS =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
cas a mole-fraction average,
vL=vEOS
L−Xxici(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+κ11 + T0.5
r(0.7−Tr)(2.29)
κ0= 0.378893 + 1.4897153 ω−0.17131848 ω2+ 0.0196554 ω3(2.30)
where κ1is 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−T0
r.51 + T0.5
r(0.7−Tr)(2.31)
κ0= 0.378893 + 1.4897153 ω−0.17131848 ω2+ 0.0196554 ω3(2.32)
where κ1,κ2, and κ3are 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 Grayson-
Streed correlation is an extension of the Chao-Seader method with special applicability to hy-
drogen. 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 )
∂ni
]P,T,nj6=i(2.33)
where GErepresents 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
UNIQUAC and UNIFAC models The UNIQUAC equation considers g≡GE/RT formed by
two additive parts, one combinatorial term gCto take into account the size of the molecules,
and one residual term gR, which take into account the interactions between molecules:
g≡gC+gR(2.34)
The gCfunction contains only pure species parameters, while the gRfunction incorporates
two binary parameters for each pair of molecules. For a multicomponent system,
gC=X
i
xiln φi/xi+ 5 X
i
qixiln θi/φi(2.35)
and
gR=−X
i
qixiln(X
j
θjτji)(2.36)
where
φi≡(xiri)/(X
j
xjrj)(2.37)
and
θi≡(xiqi)/(X
j
xjqj)(2.38)
The isubscript indicates the species, and jis 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 genters 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 γiis 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 γC
i+ ln γR
i(2.41)
ln γC
i= 1 −Ji+ ln Ji−5qi(1 −Ji/Li+ ln Ji/Li)(2.42)
ln γR
i=qi(1 −ln si−X
j
θjτij /sj)(2.43)
where
DWSIM - Technical Manual 8
2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES
Ji=ri/(X
j
rjxj)(2.44)
L=qi/(X
j
qjxj)(2.45)
si=X
l
θlτli (2.46)
Again the isubscript identify the species, jand lare 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 γC
i= 1 −Ji+ ln Ji−5qi(1 −Ji/Li+ ln Ji/Li)(2.47)
ln γR
i=qi(1 −X
k
(θkβik/sk)−ekilnβik/sk)(2.48)
The parameters JieLiare still given by eqs. 2.58 and ??. Furthermore, the following
definitions apply:
ri=X
k
ν(i)
kRk(2.49)
qi=X
k
ν(i)
kQk(2.50)
eki = (ν(i)
kQk)/qi(2.51)
βik =X
m
emkτmk (2.52)
θk= (X
i
xiqieki)/(X
i
xjqj)(2.53)
DWSIM - Technical Manual 9
2.1 Phase Equilibria Calculation 2 THERMODYNAMIC PROPERTIES
sk=X
m
θmτmk (2.54)
si=X
l
θlτli (2.55)
τmk = exp(−amk)/T (2.56)
The isubscript identify the species, and jis an index that goes through all the species. The k
subscript identify the subgroups, and mis an index that goes through all the subgroups. The
parameter ν(i)
kis the number of the ksubgroup in a molecule of the ispecies. The subgroup
parameter values Rkand Qkand 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 equilib-
rium.
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 γC
i= 1 −Ji+ ln Ji−5qi(1 −Ji/Li+ ln Ji/Li)(2.57)
Ji=r3/4
i/(X
j
r3/4
jxj)(2.58)
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 Rkand Qkare 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:
τmk = exp(−a(0)
mk +a(1)
mkT+a(2)
mkT2)/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 liquid-
liquid equilibrium. enthalpy and entropy of liquid and vapor.
Modified UNIFAC (NIST) model This model [8] is similar to the Modified UNIFAC (Dort-
mund), 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), infi-
nite 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 consis-
tency 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 gEto 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:
ln γi=
n
P
j=1
τjixjGji
n
P
k=1
xkGki
+
n
X
j=1
xjGij
n
P
k=1
xkGkj
τij −
n
P
m=1
τmj xmGmj
n
P
k=1
xkGkj
,(2.60)
Gij =exp(−τij αij ),(2.61)
τij =aij /RT, (2.62)
where
γiActivity coefficient of component i
xiMolar fraction of component i
aij Interaction parameter between i-j (aij 6=aji)(cal/mol)
TTemperature (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 depres-
sions 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:
ln γi= ln γLR
i+ ln γMR
i+ ln γSR
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
Gex =Gex
Combinatorial +Gex
Residual +Gex
Extended Debye−H¨uckel (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:
Gex
Combinatorial
RT =X
i
xiln φi
xi−z
2X
i
qixiln φi
θi(2.65)
The two model parameters ri and qi are the volume and surface area parameters for com-
ponent 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 sur-
face 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
Gex
Residual
RT =−X
i
xiqiln
X
j
θjψji
(2.66)
the parameter ψji is given by:
ψji = exp −uji −uii
T(2.67)
ujiand uii are interaction energy parameters. The interaction energy parameters are considered
symmetrical and temperature dependent in this model
uji =uij =u0
ij +uT
ij (T−298.15) (2.68)
The values of the interaction energy parameters and are determined by fitting to experi-
mental 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 co-
efficients are obtained by partial molar differentiation of the extended Debye-Hckel law excess
Gibbs energy term. The term used for water is
ln γDH
w=2
3MwAI3/2σbI1/2(2.69)
σ(x) = 3
x31 + x−1
1 + x−2 ln (1 + x)(2.70)
In this expression, b = 1.5 (kg/mol) . The term used for ions is:
ln γ∗DH
i=−Z2
i
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 Acan be approximated in the temperature range 273.15 K <T
<383.15 K by:
A=h1.131 + 1.335E−3 (T−273.15) + 1.164E−5 (T−273.15)2i(2.72)
The activity coefficient for water is calculated in the Extended UNIQUAC model by sum-
mation of the three terms:
ln γw= ln γC
w+ ln γR
w+ ln γDH
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 γC
i
γC∞
i
+ ln γR
i
γR∞
i
+ ln γ∗DH
i(2.74)
DWSIM - Technical Manual 14
2.2 Enthalpy, Entropy and Heat Capacities 2 THERMODYNAMIC PROPERTIES
The rational, unsymmetrical activity coefficient for ions calculated with the Extended UNI-
QUAC 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
Hid values are calculated from
the ideal gas heat capacity. For
mixtures, a molar average is
used. The value calculated by
the EOS is for the phase,
independently of the number of
components present in the
mixture.
Peng-Robinson, Soave-Redlich-Kwong For the cubic equations of state, enthalpy, entropy
and heat capacities are calculated by the departure functions, which relates the phase properties
in the conditions of the mixture with the same mixture property in the ideal gas state.This way,
the following departure functions are defined [12],
H−Hid
RT =X;S−Sid
R=Y(2.75)
values for Xand Yare calculated by the PR and SRK EOS, according to the table 1:
Table 1: Enthalpy/Entropy calculation with an EOS
H−Hid
RT
S−Sid
R
PR Z−1−1
21,5bRT a−Tda
dT ×ln(Z−B)−ln P
P0−A
21,5bRT T
a
da
dT ×
×ln hV+2,414b
V+0,414bi×ln hV+2,414b
V+0,414bi
SRK Z−1−1
bRT a−Tda
dT ×ln(Z−B)−ln P
P0−A
BT
a
da
dT ×
×ln 1 + b
V×ln 1 + B
Z
In DWSIM, Po= 1 atm. Heat capacities are obtained directly from the EOS, by using the following thermodynamic
relations:
Cp−Cid
p=T
V
ˆ
∞∂2P
∂T 2dV −T(∂P/∂T )2
V
(∂P/∂V )T−R(2.76)
Cp−Cv=−T∂P
∂T 2
V
∂P
∂V T
(2.77)
DWSIM - Technical Manual 15
2.2 Enthalpy, Entropy and Heat Capacities 2 THERMODYNAMIC PROPERTIES
Lee-Kesler Enthalpies, entropies and heat capacities are calculated by the Lee-Kesler model
[13] through the following equations:
H−Hid
RTc
=TrZ−1−b2+ 2b3/Tr+ 3b4/T 2
r
TrVr−c2−3c3/T 2
r
2TrV2
r
+d2
5TrV2
r
+ 3E(2.78)
S−Sid
R+ ln P
P0= ln Z−b2+b3/T 2
r+ 2b4/T 3
r
Vr−c1−2c3/T 3
r
2V2
r
+d1
5V5
r
+ 2E(2.79)
Cv−Cid
v
R=2 (b3+ 3b4/Tr)
T2
rVr−3c3
T3
rV2
r−6E(2.80)
Cp−Cid
p
R=Cv−Cid
v
R−1−Tr∂Pr
∂Tr2
Vr
∂Pr
∂VrTr
(2.81)
E=c4
2T3
rγβ+ 1 −β+1+ γ
V2
rexp −γ
V2
r (2.82)
An iterative method is required
to calculate Vr. The user
should always watch the values
generated by DWSIM in order
to detect any issues in the
compressibility factors
generated by the Lee-Kesler
model.
Z=PrVr
Tr
= 1 + B
Vr
+C
V2
r
+D
V5
r
+c4
T3
rV2
rβ+γ
V2
rexp −γ
V2
r(2.83)
B=b1−b2/Tr−b3/T 2
r−b4/T 3
r(2.84)
C=c1−c2/Tr+c3/T 3
r(2.85)
D=d1+d2/Tr(2.86)
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)Z(r)−Z(0),(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. Vrand Z(0) are calculated for the simple fluid at the fluid Trand Pr. using the equation
2.78, and with the constants for the simple fluid, as shown in the table 3, (H−H0)/RTc
is calculated. This term is (H−H0)/RTc(0). in this calculation, Zin the equation 2.78
is Z(0).
2. The step 1 is repeated, using the same Trand 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
(H−H0)/RTc(r). In this step, Zin 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
(H−H0)/RTc=(H−H0)/RTc(0) +
ω
ω(r)(H−H0)/RTc(r)−(H−H0)/RTc(0),(2.88)
where ω(r)= 0,3978.
Table 3: Constants for the Lee-Kesler model
Constant Simple Fluid Reference Fluid
b10.1181193 0.2026579
b20.265728 0.331511
b30.154790 0.027655
b40.030323 0.203488
c10.0236744 0.0313385
c20.0186984 0.0503618
c30.0 0.016901
c40.042724 0.041577
d1×1040155488 0.48736
d2×1040.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:
c=sK
ρ,(2.89)
where:
cSpeed of sound (m/s)
KBulk 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 TRANSPORT PROPERTIES
µ=∂T
∂P H
,(2.90)
The JT coefficient is calculated rigorously by the PR and SRK equations of state, while the
Goldzberg correlation is used for all other models,
µ=0.0048823Tpc 18/T 2
pr −1
PpcCpγ,(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
PC
Z[1+(1−Tr)2/7]
RA ,(3.1)
where:
VsSaturated molar volume (m /mol)
TcCritical temperature (K)
PcCritical pressure (Pa)
TrReduced temperature
ZRA Rackett constant of the component (or the mixture)
RIdeal Gas constant (8,314 J/[mol.K])
If T>Tcm,the Rackett
method does not provide a
value for Vsand, in this case,
DWSIM uses the
EOS-generated compressibility
factor to calculate the density
of the liquid phase.
For mixtures, the equation 3.1 becomes
Vs=RXxiTci
PciZ[1+(1−Tr)2/7]
RA ,(3.2)
with Tr=T/Tcm, and
Tcm=XXφiφjTcij ,(3.3)
φi=xiVci
PxiVci
,(3.4)
Tcij =
8VciVcj1/2
V1/3
ci+V1/3
cj3
TciTcj1/2,(3.5)
DWSIM - Technical Manual 18
3.1 Density 3 TRANSPORT PROPERTIES
where:
xiMolar fraction
VciCritical 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,
V=Vs1−(0.0861488 + 0.0344483ω) ln β+P
β+Pvp ,(3.7)
with
β
P=−1−9.070217 (1 −Tr)1/3+ 62.45326 (1 −Tr)2/3−135.1102 (1 −Tr) +
+ exp 4.79594 + 0.250047ω+ 1.14188ω2(1 −Tr)4/3,(3.8)
where:
VCompressed liquid volume (m /mol)
PPressure (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 )
VSpecific volume of the fluid (m /mol)
MM Liquid phase molecular volume (kg/kmol)
For the Ideal Gas Property
Package, the compressibility
factor is considered to be equal
to 1.
Vapor Phase Vapor phase density is calculated from the compressiblity factor generated by
the EOS model, according with the following equation:
ρ=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 TRANSPORT PROPERTIES
PPressure (Pa)
ZVapor phase compressibility factor
RIdeal Gas constant (8,314 J/[mol.K])
TTemperature (K)
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
i
xiln ηi!,(3.12)
where ηiis 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 +Cln T+DT E,(3.13)
where A, B, C, D and Eare 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, .807T0,618
r−0.357 exp(−0.449Tr)+0.34 exp(−4.058Tr)+0.018(3.14)
ξ= 0,176 Tc
MM3P4
c1/6
,(3.15)
where
ηViscosity (µP )
Tc, PcComponent (or mixture) critical properties
DWSIM - Technical Manual 20
3.3 Surface Tension 3 TRANSPORT PROPERTIES
TrReduced temperature, T/Tc
MM Molecular weight (kg/kmol)
In the second step, the experimental or calculated viscosity with the Lucas method is cor-
rected to take into account the effect of pressure, by the Jossi-Stiel-Thodos method [12],
"(η−η0)Tc
MM3P4
c1/6
+ 1#1/4
= 1.023 + 0.23364ρr+
+ 0.58533ρ2
r−0.40758ρ3
r+ 0.093324ρ4
r,(3.16)
where
η, η0Corrected viscosity / Lucas method calculated viscosity (µP )
Tc, PcComponent critical properties
ρrReduced 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],
σ
P2/3
cT1/3
c
= (0.132αc−0.279) (1 −Tr)11/9(3.17)
αc= 0.9076 1 + Tbr ln(Pc/1.01325)
1−Tbr ,(3.18)
where
σSurface tension (N/m)
TcCritical temperature (K)
PcCritical pressure (Pa)
Tbr Reduced normal boiling point, Tb/Tc
3.4 Isothermal Compressibility
Isothermal compressiblity of a given phase is calculated following the thermodynamic defi-
nition:
β=−1
V
∂V
∂P (3.19)
DWSIM - Technical Manual 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=1
β(3.20)
DWSIM - Technical Manual 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
T1/6
r
(4.1)
A=A∗T0.38
b
MMβTγ
c
,(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=PPφiφjλij (4.3)
λij = 2(λ−1
i+λ−1
j)−1(4.4)
φi=xiVci
PxiVci
,(4.5)
where
λLliquid 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η∗
MM 1.32 Cv−3R
2,(4.6)
where
λVvapor phase thermal conductivity (W/[m.K])
Cvconstant volume heat capacity (J/[mol.K])
λ∗and η∗are defined by:
λ∗=λ0H(4.7)
H=16.04E−3
MM/1000 1/2
f1/2/h2/3(4.8)
DWSIM - Technical Manual 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 Tr62, T +=Tr. If Tr>2, T += 2.
h=Vc
99.2φ(4.14)
η∗=η0HMM/1000
16.04E−3(4.15)
η0= 10−7
9
X
n=1
CnT(n−4)/3
0(4.16)
T0=T/f (4.17)
DWSIM - Technical Manual 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,
ln γ∗,m
±=ν+
νln γ∗,m
++ν−
νln γ∗,m
−(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:
Φ = −ln ai
MsPion mion
(5.2)
5.3 Freezing point depression
The Schrder and van Laar equation is used:
ln ai
(1 −(Tm,i/T )) =4mhi
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 (IFC-
67). 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 Formula-
tion 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 suf-
ficient 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 de-
crease. 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 in-
volving 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 Soave-
Redlich-Kwong.
6.4 FPROPS
FPROPS is a free open-source C-based library for high-accuracy evaluation of thermody-
namic 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 con-
version 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,band care the stoichiometric coefficients of reactants and product, respectively. A
is the limiting reactant and Bis in excess. The amount of each component at the end of the
reaction can then be calculated from the following stoichiometric relationships:
NA=NA0−NA0XA(7.2)
NB=NB0−b
aNA0XA(7.3)
NC=NC0+c
a(NA0XA)(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 XAis 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 REACTIONS
K=
n
Y
j=1
(qj)νj,(7.5)
where Kis the equilibrium constant, qis the basis of components (partial pressure in the vapor
phase or activity in the liquid phase) νis the stoichiometric coefficient of component jand nis
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
i
νij ξi,(7.6)
whereξiis the coordinate of the reaction iand νij is the stoichiometric coefficient of the j
component at reaction i. Defining the molar fraction of the component ias xj=nj/nt, where
ntis the total number of mols, including inerts, whe have the following expression for each
reaction i:
fi(ξ) = X
i
ln(xi)−ln(Ki)=0,(7.7)
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 Acomponent can be defined as
rA=k[A][B]−k0[C][D](7.9)
where
DWSIM - Technical Manual 31
8 PROPERTY ESTIMATION METHODS
k=Aexp (−E/RT )(7.10)
k0=A0exp (−E0/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
FA=FA0+
V
ˆrAdV, (7.12)
where FAis the molar flow of the Acomponent and Vis the reactor volume.
8 Property Estimation Methods
8.1 Petroleum Fractions
8.1.1 Molecular weight
Riazi and Al Sahhaf method [16]
MM =1
0.01964(6.97996 −ln(1080 −Tb)3/2
,(8.1)
where
MM Molecular weight (kg/kmol)
TbBoiling 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−4Tb−7.78712SG +
+2.08476 ×10−3TbSG)]T1.26007
bSG4.98308 (8.2)
Winn [17]
MM = 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.26007d154.98308 (8.4)
Lee-Kesler[17]
t1=−12272.6 + 9486.4d15 + (8.3741 −5.9917d15)P EM e (8.5)
t2= (1 −0.77084d15 −0.02058d152)×
×(0.7465 −222.466/P EM e)×107/P EM e (8.6)
t3= (1 −0.80882d15 −0.02226d152)×
×(0.3228 −17.335/P EM e)×1012/P EMe3(8.7)
MM =t1+t2+t3(8.8)
Farah
MM = exp(6.8117 + 1.3372A−3.6283B)(8.9)
MM = exp(4.0397 + 0.1362A−0.3406B−0.9988d15 + 0.0039P EM e),(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(8.12)
ln Pc= 5.689 −0.0566/SG −(0.43639 + 4.1216/SG + 0.21343/SG2)×
×10−3Tb+ (0.47579 + 1.182/SG + 0.15302/SG2)×10−6×T2
b−
−(2.4505 + 9.9099/SG2)×10−10 ×T3
b,(8.13)
DWSIM - Technical Manual 33
8.1 Petroleum Fractions 8 PROPERTY ESTIMATION METHODS
where
TbNBP (K)
TcCritical temperature (K)
PcCritical 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 EMe))(8.19)
Riazi-Daubert[17]
Tc= 9.5233 exp(−0.0009314P EM e −0.544442d15 + 0.00064791 ×P EM e ×d15) ×
×P EM e0.81067d150.53691 (8.20)
Pc= 31958000000 exp(−0.008505P EMe −4.8014d15 + 0.005749 ×P EMe ×d15) ×
×P EM e−0.4844d154.0846 (8.21)
Riazi[17]
Tc= 35.9413 exp(−0.00069P EM e −1.4442d15 + 0.000491 ×P EM e ×d15) ×
×P EM e0.7293d151.2771 (8.22)
8.1.4 Acentric Factor
Lee-Kesler method [16]
ω=−ln Pc
1.10325 −5.92714 + 6.09648/Tbr + 1.28862 ln Tbr −0.169347T6
br
15.2518 −15.6875/Tbr −13.472 ln Tbr + 0.43577T6
br
(8.23)
Korsten[17]
ω= 0.5899 ×((P EM V/Tc)1.3)/(1 −(P EM V/Tc)1.3)×log(Pc/101325) −1(8.24)
DWSIM - Technical Manual 34
8.1 Petroleum Fractions 8 PROPERTY ESTIMATION METHODS
8.1.5 Vapor Pressure
Lee-Kesler method[16]
ln Ppv
r= 5.92714 −6.09648/Tbr −1.28862 ln Tbr + 0.169347T6
br +(8.25)
+ω(15.2518 −15.6875/Tbr −13.4721 ln Tbr + 0.43577T6
br),
where
Ppv
rReduced vapor pressure, Ppv/Pc
Tbr Reduced NBP, Tb/Tc
ωAcentric factor
8.1.6 Viscosity
Letsou-Stiel [12]
η=ξ0+ξ1
ξ(8.26)
ξ0= 2.648 −3.725Tr+ 1.309T2
r(8.27)
ξ1= 7.425 −13.39Tr+ 5.933T2
r(8.28)
ξ= 176 Tc
MM3P4
c1/6
(8.29)
where
ηViscosity (Pa.s)
PcCritical pressure (bar)
TrReduced temperature, T/Tc
MM Molecular weight (kg/kmol)
Abbott[17]
t1= 4.39371 −1.94733Kw + 0.12769Kw2+ 0.00032629AP I2−0.0118246KwAP I +
+(0.171617Kw2+ 10.9943AP I + 0.0950663AP I2−0.869218KwAP I (8.30)
log v100 =t1
AP I + 50.3642 −4.78231Kw ,(8.31)
t2=−0.463634 −0.166532AP I + 0.000513447AP I2−0.00848995AP IKw +
+(0.080325Kw + 1.24899AP I + 0.19768AP I2(8.32)
log v210 =t2
AP I + 26.786 −2.6296Kw ,(8.33)
where
DWSIM - Technical Manual 35
8.2 Hypothetical Components 8 PROPERTY ESTIMATION METHODS
v100 Viscosity at 100 F (cSt)
v210 Viscosity at 210 F (cSt)
KwWatson characterization factor
AP I Oil API degree
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 TcJoback [12]
Critical pressure PcJoback [12]
Critical volume VcJoback [12]
Normal boiling point TbJoback [12]
Vapor pressure Ppv Lee-Kesler (Eq. 8.25)
Acentric factor ωLee-Kesler (Eq. 8.23)
Vaporization enthalpy ∆Hvap Vetere [12]
Ideal gas heat capacity Cgi
pHarrison-Seaton [18]
Ideal gas enthalpy of formation ∆H298
fMarrero-Gani [19]
DWSIM - Technical Manual 36
9 OTHER PROPERTIES
9 Other Properties
9.1 True Critical Point
The Gibbs criteria for the true critical point of a mixture of ncomponents may be expressed
of various forms, but the most convenient when using a pressure explicit cubic equation of state
is
L=
A11 A12 . . . A1n
A21 A22
.
.
.
An1. . . . . . Ann
= 0 (9.1)
M=
A11 A12 . . . A1n
A21 A22
.
.
.
An−1,1. . . . . . An−1,n
∂L
∂n1. . . . . . ∂L
∂nn
= 0,(9.2)
where
A12 =∂2A
∂n1∂n2T,V
(9.3)
All the Aterms in the equations 9.1 and 9.2 are the second derivatives of the total Helmholtz
energy Awith respect to mols and constant Tand 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,
fi
w=fH
w,(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:
µi
w=µH
w(9.5)
remembering that
DWSIM - Technical Manual 37
9.2 Natural Gas Hydrates 9 OTHER PROPERTIES
fi=xiPexp((µ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
νmln(1 −X
j
θmj ),(9.7)
where µβ
wis the chemical potential of water in the empty hydrate lattice (something like an
”ideal” chemical potential) and νmis the number of mcavities by water molecule in the lattice.
The fraction of cavities m-type cavities occupied by the gaseous component lis
θml = (Cmlfi)/((1 + X
j
Cmj fj)),(9.8)
where Cmj is the Langmuir constant and fi 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
C(T)=4π/kT ˆ∞
0
exp[(−w(r))/kT ]r2dr, (9.9)
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/(R11r)(δ10 +a/Rδ11)−σ6/(R5r)(δ4+a/Rδ5)],(9.10)
δN= [(1 −r/R −a/R)−N)−(1 + r/R −a/R)−N)]/N, (9.11)
where Nis equal to 4, 5, 10 or 11; zand Rare, respectively, the coordination number and
the cavity cell radius.
Supported hydrate formers CH4, C2H6, C3H8, iC4H10, H2S, N2 and CO2.
Model applicability range Temperature: 211 to 303 K; Pressure: 1 to 600 atm.
DWSIM - Technical Manual 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
fH
w= exp(Aβ
gln T+ (Bβ
g)/T + 2,7789 + Dβ
gT)×
exp Vβ
w[P−exp(Aβ
gln T+ (Bβ
g)/T + 2,7789 + Dβ
gT)]/RT ×
exp[X
m
νmln(1 −X
j
(Cmlfl)/(1 + XCmj fj))] (9.12)
The A,Band Dconstants are specific for each hydrate former and represent the vapor
pressure of the component in the empty hydrate lattice. Vβ
wrepresents the basic hydrate molar
volume (without the presence of guests) and the Langmuir constant (C) is calculated by the
following equation:
C(T) = 4π/kT ˆR−a
0
exp[(−w(r))/kT ]r2dr (9.13)
In the Klauda e Sandler method the spherically symmetric Kihara potential is also used,
w(r)=2ze[σ12/(R11r)(δ10 +a/Rδ11)−σ6/(R5r)(δ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 approx-
imation 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
fH
i=f0
i(1 −θi)α,(9.17)
where αdepends on the structure of the hydrate formed, and
f0
i=f0(P)f0(T)f0(xwγw),(9.18)
f0(P) = exp(βP/T ),(9.19)
f0(T) = A0exp(B0/(T−C0)),(9.20)
f0(xwγw)=(xwγw)(−1/λ2),(9.21)
where βand λ2depend on the structure of the hydrate formed and A’,B’ and C’ depends on
the hydrate former. xwand γmare, 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) = Xexp(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.8MeABP )+2.468SG −0.0005704(1.8M eABP )×SG)×
×(1.8MeABP )0.0572SG−0.72 (9.23)
r=1+2I
1−I1/2
(9.24)
where
rRefraction Index
SG Specific Gravity
MeABP 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
P F Flash Point (K)
t10AST M ASTM D93 10% vaporized temperature (K)
9.3.3 Pour Point
API Procedure 2B8.1
PFL = [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×MeABP (9.27)
where
P C Freezing Point (K)
KAPI characterization factor (API)
9.3.5 Cloud Point
API Procedure 2B12.1
P N =h10(−7.41+5.49 log(1.8MeABP )−0.712×(1.8MeABP )0.315 −0.133SG)i/1.8(9.28)
where
P N Cloud Point (K)
9.3.6 Cetane Index
API Procedure 2B13.1
IC = 415.26 −7.673AP I + 0.186 ×(1.8MeABP −458.67) + 3.503AP I ×
×log(1.8MeABP −458.67) −193.816 ×log(1.8MeABP −458.67) (9.29)
DWSIM - Technical Manual 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
VL1/2
(9.30)
where
∆HvMolar Heat of Vaporization
VLLiquid Molar Volume at 20 C
DWSIM - Technical Manual 42
References References
References
[1] J. Smith, Intro to Chemical Engineering Thermodynamics. New York: McGraw-Hill Com-
panies, 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,” Chem-
ical 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 three-
parameter 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 proper-
ties,” 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 Engi-
neering Journal, vol. 71, no. 2, pp. 145–151, December 1998, url.
DWSIM - Technical Manual 44

Navigation menu