Thermodynamic Models

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.

The following sections describe the calculation methods used in DWSIM for the physical and chemical description of the elements of a simulation.

Thermodynamic Properties

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:

`f_{i}^{L}=f_{i}^{V}`

The fugacity of a component in a mixture depends on temperature, pressure and composition. in order to relate `f_{i}^{V}` with temperature, pressure and molar fraction, we define the fugacity coefficient,

`\phi_{i}=\frac{f_{i}^{V}}{y_{i}P},`

which can be calculated from PVT data, commonly obtained from an equation of state. For a mixture of ideal gases, `\phi_{i}=1`.

The fugacity of the i component in the liquid phase is related to the composition of that phase by the activity coefficient `\gamma_{i}` , which by itself is related to `x_{i}` and standard-state fugacity `f_{i}^{0}` by

`\gamma_{i}=\frac{f_{i}^{L}}{x_{i}f_{i}^{0}}.`

The standard state fugacity `f_{i}^{0}` 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

`\phi_{i}=\frac{f_{i}^{L}}{x_{i}P},`

with the fugacity coefficient `\phi_{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\phi_{i}^{L} = \int_{V^{L}}^{\infty}[(\frac{\partial P}{\partial n_{i}}\)_{T,V,n_{j}}-\frac{RT}{V}\]dV-RT\ln Z^{L},`

` RT\ln\phi_{i}^{V}=\int_{V^{V}}^{\infty}[(\frac{\partial P}{\partial(n_{i})}\)_{T,V,n_{j}}-\frac{RT}{V}\]dV-RT\ln Z^{V},`

where the compressibility factor `Z` is given by

`Z^{L}=\frac{PV^{L}}{RT}`

`Z^{V}=\frac{PV^{V}}{RT}`

Fugacity Coefficient calculation models

Peng-Robinson Equation of State

The Peng-Robinson equation 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=\frac{RT}{(V-b)}-\frac{a(T)}{V(V+b)+b(V-b)}`

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\omega-0.26992\omega^{2})(1-T_{r}^{(1/2)})]^{2}0.45724(R^{2}T_{c}^{2})/P_{c}`

`b=0.07780(RT_{c})/P_{c}`

where

`\omega` acentric factor

`T_{c}` critical temperature

`P_{c}` critical pressure

`T_{r}` reduced temperature, T/Tc

For mixtures, equation [eq:PR] can be used, replacing a and b by mixture-representative values. a and b mixture values are normally given by the basic mixing rule,

`a_{m}=\sum_{i}\sum_{j}x_{i}x_{j}\sqrt{(a_{i}a_{j})}(1-k_{ij})`

`b_{m}=\sum_{i}x_{i}b_{i}`

where

`x_{i,j}` molar fraction of the i or j component in the phase (liquid or vapor)

`a_{i,j}` i or j component a constant

`b_{i,j}` i or j component b constant

`k_{ij}` binary interaction parameter which characterizes the i-j pair

The fugacity coefficient obtained with the Peng-Robinson EOS in given by

`\lnfrac{f_{i}}{x_{i}P}=\frac{b_{i}}{b_{m}}(Z-1)-\ln(Z-B)-\frac{A}{2\sqrt{2}B}(\frac{\sum_{k}x_{k}a_{ki}}{a_{m}}-\frac{b_{i}}{b_{m}})\ln(\frac{Z+2,414B}{Z-0,414B}),`

where Z in the phase compressibility factor (liquid or vapor) and can be obtained from the equation,

`Z^{3}-(1-B)Z^{2}+(A-3B^{2}-2B)Z-(AB-B^{2}-2B)=0,`

`A=\frac{a_{m}P}{R^{2}T^{2}}`

`B=\frac{b_{m}P}{RT}`

`Z=\frac{PV}{RT}`

Soave-Redlich-Kwong Equation of State

The Soave-Redlich-Kwong Equation [20] is also a cubic equation of state in volume,

`P=\frac{RT}{(V-b)}-\frac{a(T)}{V(V+b)},`

The a and b parameters are given by:

`a(T)=[1+(0.48+1.574\omega-0.176\omega^{2})(1-T_{r}^{(1/2)})]^{2}0.42747(R^{2}T_{c}^{2})/P_{c}`

`b=0.08664(RT_{c})/P_{c}`

The equations [eq:mixrule1] and [eq:mixrule2] are used to calculate mixture parameters. Fugacity is calculated by

`\ln\dfrac{f_{i}}{x_{i}P}=\frac{b_{i}}{b_{m}}(Z-1)-\ln(Z-B)-\frac{A}{B}(\frac{\sum_{k}x_{k}a_{ki}}{a_{m}}-\frac{b_{i}}{b_{m}})\ln(\frac{Z+B}{Z})`

The phase compressibility factor Z is obtained from

`Z^{3}-Z^{2}+(A-B-B^{2})Z-AB=0,`

`A=\frac{a_{m}P}{R^{2}T^{2}}`

`B=\frac{b_{m}P}{RT}`

`Z=\frac{PV}{RT}`

The cubic-Z equations, 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 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`, 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,

`v_{L}=v_{L}^{EOS}-\sum x_{i}c_{i}`

Volume translation can be applied to any two-constant cubic equation, thereby eliminating the volumetric defficiency suffered by all two-constant equations.

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:

`\kappa=\kappa_{0}+\kappa_{1}(1+T_{r}^{0.5})(0.7-T_{r})`

`\kappa_{0}=0.378893+1.4897153\,\omega-0.17131848\,\omega^{2}+0.0196554\,\omega^{3}`

where `\kappa_{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

A subsequent modification published in 1986 (PRSV2) [21] further improved the model's accuracy by introducing two additional pure component parameters to the previous attraction term modification.

The modification is:

`\kappa=\kappa_{0}+[\kappa_{1}+\kappa_{2}(\kappa_{3}-T_{r})(1-T_{r}^{0}.5)](1+T_{r}^{0.5})(0.7-T_{r})`

`\kappa_{0}=0.378893+1.4897153\,\omega-0.17131848\,\omega^{2}+0.0196554\,\omega^{3}`

where `\kappa_{1}`, `\kappa_{2}`, and `\kappa_{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.

Chao-Seader and Grayson-Streed models

Chao-Seader ([2]) and Grayson-Streed ([4]) are older, semi-empirical models. The Grayson-Streed 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 10342 kPa (1500 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 hydrogen content.

Calculation models for the liquid phase activity coefficient

The activity coefficient \gamma 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

`\gamma_{i}=(\frac{\partialnG^{E}/{RT}}(\partialn_{i}))_{P,T,n_{j\neq i}}`

where `G^{E}` 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, `\gamma_{i} = 1`. Expressions for `G^{E}/{RT}` provide values for the activity coefficients.

UNIQUAC and UNIFAC models

The UNIQUAC equation considers `g\equiv G^{E}/{RT}` formed by 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\equiv g^{C}+g^{R}`

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,

`g^{C}=\sum_{i}x_{i}\ln\phi_{i}/x_{i}+5\sum_{i}q_{i}x_{i}\ln\theta_{i}/\phi_{i}`

and

`g^{R}=-\sum_{i}q_{i}x_{i}\ln(\sum_{j}\theta_{j}\tau_{j}i)`

where

`\phi_{i}\equiv(x_{i}r_{i})/(\sum_{j}x_{j}r_{j})`v

and

`\theta_{i}\equiv(x_{i}q_{i})/(\sum_{j}x_{j}q_{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 `\tau_{ij}\neq\tau_{ji}`. When `i=j`, `\tau_{ii}=\tau_{jj}=1`.

In these equations, `r_{i}` (a relative molecular volume) and `q_{i}` (a relative molecular surface area) are pure species parameters. The influence of temperature in `g` enters by means of the `\tau_{ij}` parameters, which are temperature-dependent:

`\tau_{ij}=\exp(u_{ij}-u_{jj})/{RT}`

This way, the UNIQUAC parameters are values of `(u_{ij}-u_{jj})`.

An expression for `\gamma_{i}` is found through the application of the following relation:

`\ln\gamma_{i}=[\partialnG^{E}/{RT}/(\partialn_{i})]_{(P,T,n_{j\neq i})}`

The result is represented by the following equations:

`\ln\gamma_{i}=\ln\gamma_{i}^{C}+\ln\gamma_{i}^{R}`

`\ln\gamma_{i}^{C}=1-J_{i}+\ln J_{i}-5q_{i}(1-J_{i}/L_{i}+\ln J_{i}/L_{i})`

`\ln\gamma_{i}^{R}=q_{i}(1-\ln s_{i}-\sum_{j}\theta_{j}\tau_{ij}/s_{j})`

where

`J_{i}=r_{i}/(\sum_{j}r_{j}x_{j})`

`L=q_{i}/(\sum_{j}q_{j}x_{j})`

`s_{i}=\sum_{l}\theta_{l}\tau_{li}`

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 `\tau_{ij}=1` for `i=j`. The parameters values `(u_{ij}-u_{jj})` 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. When applied to a solution of groups, the equations are written in the form:

`\ln\gamma_{i}^{C}=1-J_{i}+\ln J_{i}-5q_{i}(1-J_{i}/L_{i}+\ln J_{i}/L_{i})`

`\ln\gamma_{i}^{R}=q_{i}(1-\sum_{k}(\theta_{k}\beta_{ik}/s_{k})-e_{ki}ln\beta_{ik}/s_{k})`

Furthermore, the following definitions apply: `r_{i}=\sum_{k}\nu_{k}^{(i)}R_{k}`

`q_{i}=\sum_{k}\nu_{k}^{(i)}Q_{k}`

`e_{ki}=(\nu_{k}^{(i)}Q_{k})/q_{i}`

`\beta_{ik}=\sum_{m}e_{mk}\tau_{mk}`

`\theta_{k}=(\sum_{i}x_{i}q_{i}e_{ki})/(\sum_{i}x_{j}q_{j})`

`s_{k}=\sum_{m}\theta_{m}\tau_{mk}`

`s_{i}=\sum_{l}\theta_{l}\tau_{li}`

`\tau_{mk}=\exp(-a_{mk})/T`

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 parameter `\nu_{k}^{(i)}` is the number of the k subgroup in a molecule of the i species. The subgroup parameter values `R_{k}` and `Q_{k}` and the interaction parameters `-a_{mk}` 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 `\gamma-\phi` 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:

`\ln\gamma_{i}^{C}=1-J_{i}+\ln J_{i}-5q_{i}(1-J_{i}/L_{i}+\ln J_{i}/L_{i})`

`J_{i}=r_{i}^{3/4}/(\sum_{j}r_{j}^{3/4}x_{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 `R_{k}` and `Q_{k}` 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:

`\tau_{mk}=\exp(-a_{mk}^{(0)}+a_{mk}^{(1)}T+a_{mk}^{(2)}T^{2})/T`

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 [7] 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 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 [16] 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\gamma_{i}=\frac{\underset{j=1}{\overset{n}{\sum}}\tau_{ji}x_{j}G_{ji}}{\underset{k=1}{\overset{n}{\sum}}x_{k}G_{ki}}+\underset{j=1}{\overset{n}{\sum}}\frac{x_{j}G_{ij}}{\underset{k=1}{\overset{n}{\sum}}x_{k}G_{kj}}(\tau_{ij}-\frac{\underset{m=1}{\overset{n}{\sum}}\tau_{mj}x_{m}G_{mj}}{\underset{k=1}{\overset{n}{\sum}}x_{k}G_{kj}}),`

`G_{ij}=exp(-\tau_{ij}\alpha_{ij}),`

`\tau_{ij}=a_{ij}/{RT},`

where

`\gamma_{i}` Activity coefficient of component i

`x_{i}` Molar fraction of component i

`a_{ij}` Interaction parameter between i-j `(a_{ij}\neq a_{ji})` (cal/mol)

`T` Temperature (K)

`\alpha_{ij}` non-randomness parameter for the i-j pair `(\alpha_{ij}=\alpha_{ji})`

The significance of `G_{ij}` is similar to `\Lambda_{ij}` from Wilson's equation, that is, they are characteristic 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.

Enthalpy, Entropy and Heat Capacities

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 [15],

`\frac{H-H^{id}}{RT}=X`;`\frac{S-S^{id}}{R}=Y`

values for X and Y are calculated by the PR and SRK EOS, according to the table:

	`\frac{H-H^{id}}{RT}`	`\frac{S-S^{id}}{R}`
PR	`Z-1-\frac{1}{2^{1.5}bRT}[a-T\frac{da}{dT}]\ln[\frac{V+2.414b}{V+0.414b}]`	`\ln(Z-B)-\ln\frac{P}{P^{0}}-\frac{A}{2^{1.5}bRT}[\frac{T}{a}\frac{da}{dT}]\ln[\frac{V+2.414b}{V+0.414b}]`
SRK	`Z-1-\frac{1}{bRT}[a-T\frac{da}{dT}]\ln[1+\frac{b}{V}]`	`\ln(Z-B)-\ln\frac{P}{P^{0}}-\frac{A}{B}[\frac{T}{a}\frac{da}{dT}]\ln[1+\frac{B}{Z}]`

`P_{o}= 101325 Pa`.

`H^{id}` 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.

Heat capacities are obtained directly from the EOS, by using the following thermodynamic relations:

`C_{p}-C_{p}^{id}=T\intop_{\infty}^{V}(\frac{\partial^{2}P}{\partial T^{2}})dV-\frac{T(\partial P/\partial T)_{V}^{2}}{(\partial P/\partial V)_{T}}-R`

`C_{p}-C_{v}=-T\frac{(\frac{\partial P}{\partial T})_{V}^{2}}{(\frac{\partial P}{\partial V})_{T}}`

Lee-Kesler

Enthalpies, entropies and heat capacities are calculated by the Lee-Kesler model [9] through the following equations:

`\frac{H-H^{id}}{RT_{c}}=T_{r}(Z-1-\frac{b_{2}+2b_{3}/T_{r}+3b_{4}/T_{r}^{2}}{T_{r}V_{r}}-\frac{c_{2}-3c_{3}/T_{r}^{2}}{2T_{r}V_{r}^{2}}+\frac{d_{2}}{5T_{r}V_{r}^{2}}+3E)`

`\frac{S-S^{id}}{R}+\ln(\frac{P}{P_{0}})=\ln Z-\frac{b_{2}+b_{3}/T_{r}^{2}+2b_{4}/T_{r}^{3}}{V_{r}}-\frac{c_{1}-2c_{3}/T_{r}^{3}}{2V_{r}^{2}}+\frac{d_{1}}{5V_{r}^{5}}+2E`

`\frac{C_{v}-C_{v}^{id}}{R}=\frac{2(b_{3}+3b_{4}/T_{r})}{T_{r}^{2}V_{r}}-\frac{3c_{3}}{T_{r}^{3}V_{r}^{2}}-6E`

`\frac{C_{p}-C_{p}^{id}}{R}=\frac{C_{v}-C_{v}^{id}}{R}-1-T_{r}\frac{(\frac{\partial P_{r}}{\partial T_{r}^ {}})_{V_{r}}^{2}}{(\frac{\partial P_{r}}{\partial V_{r}})_{T_{r}}}`

`E=\frac{c_{4}}{2T_{r}^{3}\gamma}[\beta+1-(\beta+1+\frac{\gamma}{V_{r}^{2}})\exp(-\frac{\gamma}{V_{r}^{2}})]`

An iterative method is required to calculate `V_{r}`. 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=\frac{P_{r}V_{r}}{T_{r}}=1+\frac{B}{V_{r}}+\frac{C}{V_{r}^{2}}+\frac{D}{V_{r}^{5}}+\frac{c_{4}}{T_{r}^{3}V_{r}^{2}}(\beta+\frac{\gamma}{V_{r}^{2}})\exp(-\frac{\gamma}{V_{r}^{2}})`

`B=b_{1}-b_{2}/T_{r}-b_{3}/T_{r}^{2}-b_{4}/T_{r}^{3}`

`C=c_{1}-c_{2}/T_{r}+c_{3}/T_{r}^{3}`

`D=d_{1}+d_{2}/T_{r}`

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)}+\frac{\omega}{\omega^{(r)}}(Z^{(r)}-Z^{(0)}),`

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. `V_{r}` and `Z^{(0)}` are calculated for the simple fluid at the fluid `T_{r}` and `P_{r}`. using the equation [eq:LKH], and with the constants for the simple fluid, as shown in the table , `(H-H^{0})/{RT}_{c}` is calculated. This term is `[(H-H^{0})/{RT}_{c}]^{(0)}` . in this calculation, Z in the equation is `Z^{(0)}`.

2. Step 1 is repeated, using the same `T_{r}` and `P_{r}`, but using the constants for the reference fluid as shown in table. With these values, the equation allows the calculation of `[(H-H^{0})/{RT}_{c}]^{(r)}`. In this step, `Z` in the equation is `Z^{(r)}`.

3. Finally, one determines the residual enthalpy for the fluid of interest by

`[(H-H^{0})/{RT}_{c}] = [(H-H^{0})/{RT}_{c}]^{(0)}+\frac{\omega}{\omega^{(r)}}([(H-H^{0})/{RT}_{c}]^{(r)}-[(H-H^{0})/{RT}_{c}]^{(0)}),`

where `\omega^{(r)}=0.3978`.

Constants for the Lee-Kesler model

Constant	Simple Fluid	Reference Fluid
`b_{1}`	0.1181193	0.2026579
`b_{2}`	0.265728	0.331511
`b_{3}`	0.154790	0.027655
`b_{4}`	0.030323	0.203488
`c_{1}`	0.0236744	0.0313385
`c_{2}`	0.0186984	0.0503618
`c_{3}`	0.0	0.016901
`c_{4}`	0.042724	0.041577
`d_{1}\times10^{4}`	0155488	0.48736
`d_{2}\times10^{4}`	0.623689	0.0740336
`\beta`	0.65392	1.226
`\gamma`	0.060167	0.03754

Speed of Sound

The speed of sound in a given phase is calculated by the following equations:

`c=\sqrt{\frac{K}{\rho}},`

where:

`c` Speed of sound (m/s)

`K` Bulk Modulus (Pa)

`\rho` Phase Density (kg/m?)

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:

`\mu=(\frac{\partial T}{\partial P})_{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,

`\mu=\frac{0.0048823T_{pc}(18/T_{pr}^{2}-1)}{P_{pc}C_{p}\gamma},`

for gases, and

`\mu=-\frac{1}{\rho C_{p}},`

for liquids.

Transport Properties

Density

Liquid Phase

Liquid phase density is calculated with the Rackett equation for non-EOS models when experimental data is not available [15],

`V_{s}=\frac{RT_{C}}{P_{C}}Z_{RA}^{[1+(1-T_{r})^{2/7}]},`

where:

`V_{s}` Saturated molar volume (m?/mol)

`T_{c}` Critical temperature (K)

`P_{c}` Critical pressure (Pa)

`T_{r}` Reduced temperature

`Z_{RA}` Rackett constant of the component/mixture

`R` Ideal Gas constant (8.314 J/[mol.K])

If `T\geq T_{cm}`, the Rackett method does not provide a value for `V_{s}` and, in this case, DWSIM uses the EOS-generated compressibility factor to calculate the density of the liquid phase.

For mixtures, the equation becomes

`V_{s}=R(\sum\frac{x_{i}T_{c_{i}}}{P_{c_{i}}})Z_{RA}^{[1+(1-T_{r})^{2/7}]},`

with `T_{r}=T/T_{cm}`, and

`T_{c_{m}}=\sum\sum\phi_{i}\phi_{j}T_{c_{ij}},`

`\phi_{i}=\frac{x_{i}V_{c_{i}}}{\sum x_{i}V_{c_{i}}},`

`T_{c_{ij}}=[\frac{8(V_{c_{i}}V_{c_{j}})^{1/2}}{(V_{c_{i}}^{1/3}+V_{c_{j}}^{1/3})^{3}}](T_{c_{i}}T_{c_{j}})^{1/2},`

where:

`x_{i}` Molar fraction

`V_{c_{i}}` Critical volume (m?/mol)

If `Z_{RA}` isn't available, it is calculated from the component acentric factor,

`Z_{RA}=0.2956-0.08775\omega,`

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=V_{s}[1-(0.0861488+0.0344483\omega)\ln\frac{\beta+P}{\beta+P_{vp}}],`

with

`\frac{\beta}{P} = -1-9.070217(1-T_{r})^{1/3}+62.45326(1-T_{r})^{2/3}-135.1102(1-T_{r})+ +\exp(4.79594+0.250047\omega+1.14188\omega^{2})(1-T_{r})^{4/3},`

where:

`V` Compressed liquid volume (m?/mol)

`P` Pressure (Pa)

`P_{vp}` Vapor pressure / Bubble point pressure (Pa)

Finally, density is calculated from the molar volume by the following relation:

`\rho=\frac{MM}{1000V},`

where:

`\rho` Density (kg/m?)

`V` Specific volume of the fluid (m?/mol)

`MM` Liquid phase molecular volume (kg/kmol)

Vapor Phase

For the Ideal Gas Property Package, the compressibility factor is considered to be equal to 1.

Vapor phase density is calculated from the compressiblity factor generated by the EOS model, according with the following equation:

`\rho=\frac{MM\,P}{1000ZRT},`

where:

`\rho` Density (kg/m?)

`MM` Molecular weight of the vapor phase (kg/kmol)

`P` Pressure (Pa)

`Z` Vapor phase compressibility factor

`R` Ideal Gas constant (8.314 J/[mol.K])

`T` Temperature (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:

`\rho_{m}=f_{l}\rho_{l}+f_{v}\rho_{v},`

where:

`\rho_{m,l,v}` Density of the mixture / liquid phase / vapor phase (kg/m?)

`f_{l,v}` Volumetric fraction of the liquid phase / vapor phase (kg/kmol)

Viscosity

Liquid Phase

When experimental data is not available, liquid phase viscosity is calculated from

`\eta_{L}=\exp(\sum_{i}x_{i}\ln\eta_{i}),`

where `\eta_{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

`\eta=\exp(A+B/T+C\ln T+DT^{E}),`

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 [15],

`\eta\xi=[0,.807T_{r}^{0,618}-0.357\exp(-0.449T_{r})+0.34\exp(-4.058T_{r})+0.018]`

`\xi=0,176(\frac{T_{c}}{MM^{3}P_{c}^{4}})^{1/6},`

where

`\eta` Viscosity `(\mu P)`

`T_{c},P_{c}` Component (or mixture) critical properties

`T_{r}` Reduced temperature, `T/T_{c}`

`MM` Molecular weight (kg/kmol)

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 [15],

`[(\eta-\eta_{0})(\frac{T_{c}}{MM^{3}P_{c}^{4}})^{1/6}+1]^{1/4} = 1.023+0.23364\rho_{r}+ + 0.58533\rho_{r}^{2}-0.40758\rho_{r}^{3}+0.093324\rho_{r}^{4},`

where

`\eta,\eta_{0}` Corrected viscosity / Lucas method calculated viscosity `(\mu P)`

`T_{c},P_{c}` Component critical properties

`\rho_{r}` Reduced density, `\rho/\rho_{c}=V/V_{c}`

`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.

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 [15],

`\frac{\sigma}{P_{c}^{2/3}T_{c}^{1/3}}=(0.132\alpha_{c}-0.279)(1-T_{r})^{11/9}`

`\alpha_{c}=0.9076[1+\frac{T_{br}\ln(P_{c}/1.01325)}{1-T_{br}}],`

where

`\sigma` Surface tension (N/m)

`T_{c}` Critical temperature (K)

`P_{c}` Critical pressure (Pa)

`T_{br}` Reduced normal boiling point, `T_{b}/T_{c}`

Isothermal Compressibility

Isothermal compressiblity of a given phase is calculated following the thermodynamic definition:

`\beta=-\frac{1}{V}\frac{\partial V}{\partial P}`

The above expression is calculated rigorously by the PR and SRK equations of state. For the other models, a numerical derivative approximation is used.

Bulk Modulus

The Bulk Modulus of a phase is defined as the inverse of the isothermal compressibility:

`K=\frac{1}{\beta}`

Thermal Properties

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 [15],

`\lambda_{i} = \frac{A(1-T_{r})^{0.38}}{T_{r}^{1/6}}`

`A = \frac{A^{*}T_{b}^{0.38}}{MM^{\beta}T_{c}^{\gamma}},`

where `A^{*}`,`\alpha`,`\beta` and `\gamma` 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 [15],

`\lambda_{L} =\sum\sum\phi_{i}\phi_{j}\lambda_{ij} `

`\lambda_{ij} = 2(\lambda_{i}^{-1}+\lambda_{j}^{-1})^{-1}`

`\phi_{i}=\frac{x_{i}V_{c_{i}}}{\sum x_{i}V_{c_{i}}},`

where

`\lambda_{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 [15],

`\lambda_{V}=\lambda^{*}+\frac{1000\eta^{*}}{MM}1.32(C_{v}-\frac{3R}{2}),`

where

`\lambda_{V}` vapor phase thermal conductivity (W/[m.K])

`C_{v}` constant volume heat capacity (J/[mol.K])

`\lambda^{*}` and `\eta^{*}` are defined by:

`\lambda^{*}=\lambda_{0}H`

`H=(\frac{16.04E-3}{MM/1000})^{1/2}f^{1/2}/h^{2/3}`

`\lambda_{0}=1944\eta_{0}`

`f=\frac{T_{0}\theta}{190.4}`

`h=\frac{V_{c}}{99.2}\phi`

`\theta=1+(\omega-0.011)(0.56553-0.86276\ln T^{+}-0.69852/T^{+}`

`\phi=[1+(\omega-0.011)(0.38650-1.1617\ln T^{+})]0.288/Z_{c}`

If `T_{r}\leq 2, T^{+}=T_{r}.` If `T_{r}>2,T^{+}=2.`

`h=\frac{V_{c}}{99.2}\phi`

`\eta^{*}=\eta_{0}H\frac{MM/1000}{16.04E-3}`

`\eta_{0}=10^{-7}\sum_{n=1}^{9}C_{n}T_{0}^{(n-4)/3}`

`T_{0}=T/f`

Specialized Models / Property Packages

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.

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 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

Property Estimation Methods

Petroleum Fractions

Molecular weight

Riazi and Al Sahhaf method [18]

`MM=[\frac{1}{0.01964}(6.97996-\ln(1080-T_{b})]^{3/2},`

where

`MM` Molecular weight (kg/kmol)

`T_{b}` 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\times10^{-4}T_{b}-7.78712SG+ +2.08476\times10^{-3}T_{b}SG)]T_{b}^{1.26007}SG^{4.98308}`

Winn [17]

`MM=0.00005805PEMe^{2.3776}/d15^{0.9371},`

where

`PEMe` Mean Boiling Point (K)

`d15` Specific Gravity @ 60 ?F

Riazi[17]

`MM = 42.965\exp(0.0002097PEMe-7.78d15+0.00208476\times PEMe\times d15)\times \times PEMe^{1.26007}d15^{4.98308}`

Lee-Kesler[17]

`t_{1} = -12272.6+9486.4d15+(8.3741-5.9917d15)PEMe`

`t_{2} = (1-0.77084d15-0.02058d15^{2})(0.7465-222.466/PEMe)10^{7}/PEMe`

`t_{3} = (1-0.80882d15-0.02226d15^{2})(0.3228-17.335/PEMe)\times10^{12}/PEMe^{3}`

`MM = t_{1}+t_{2}+t_{3}`

Farah

`MM = \exp(6.8117+1.3372A-3.6283B)`

`MM = \exp(4.0397+0.1362A-0.3406B-0.9988d15+0.0039PEMe),`

where

`A,B` Walther-ASTM equation parameters for viscosity calculation

Specific Gravity

Riazi e Al Sahhaf [18]

`SG=1.07-\exp(3.56073-2.93886MM^{0.1}),`

where

`SG` Specific Gravity

`MM` Molecular weight (kg/kmol)

Critical Properties

Lee-Kesler [18]

`T_{c}=189.8+450.6SG+(0.4244+0.1174SG)T_{b}+(0.1441-1.0069SG)10^{5}/T_{b}`

`\ln P_{c} = 5.689-0.0566/SG-(0.43639+4.1216/SG+0.21343/SG^{2})10^{-3}T_{b}+(0.47579+1.182/SG+0.15302/SG^{2})10^{-6}T_{b}^{2}- -(2.4505+9.9099/SG^{2})\times10^{-10}T_{b}^{3},`

where

`T_{b}` NBP (K)

`T_{c}` Critical temperature (K)

`P_{c}` Critical pressure (bar)

Farah

`T_{c} = 731.968+291.952A-704.998B`

`T_{c} = 104.0061+38.75A-41.6097B+0.7831PEMe`

`T_{c} = 196.793+90.205A-221.051B+309.534d15+0.524PEMe`

`P_{c} = \exp(20.0056-9.8758\ln(A)+12.2326\ln(B))`

`P_{c} = \exp(11.2037-0.5484A+1.9242B+510.1272/PEMe)`

`P_{c} = \exp(28.7605+0.7158\ln(A)-0.2796\ln(B)+2.3129\ln(d15)-2.4027\ln(PEMe))`

Riazi-Daubert[17]

`T_{c} = 9.5233\exp(-0.0009314PEMe-0.544442d15+0.00064791\times PEMe\times d15)\times \times PEMe^{0.81067}d15^{0.53691}`

`P_{c} = 31958000000\exp(-0.008505PEMe-4.8014d15+0.005749\times PEMe\times d15)\times \times PEMe^{-0.4844}d15^{4.0846}`

Riazi[17]

`T_{c} = 35.9413\exp(-0.00069PEMe-1.4442d15+0.000491\times PEMe\times d15)\times \times PEMe^{0.7293}d15^{1.2771}`

Acentric Factor

Lee-Kesler method [18]

`\omega=\frac{-\ln\frac{P_{c}}{1.10325}-5.92714+6.09648/T_{br}+1.28862\ln T_{br}-0.169347T_{br}^{6}}{15.2518-15.6875/T_{br}-13.472\ln T_{br}+0.43577T_{br}^{6}}`

Korsten[17]

`\omega=0.5899\times((PEMV/T_{c})^{1.3})/(1-(PEMV/T_{c})^{1.3})\times\log(P_{c}/101325)-1`

Vapor Pressure

Lee-Kesler method[18]

`\ln P_{r}^{pv} = 5.92714-6.09648/T_{br}-1.28862\ln T_{br}+0.169347T_{br}^{6}+ +\omega(15.2518-15.6875/T_{br}-13.4721\ln T_{br}+0.43577T_{br}^{6}),`

where

`P_{r}^{pv}` Reduced vapor pressure, `P^{pv}/P_{c}`

`T_{br}` Reduced NBP, `T_{b}/T_{c}`

`\omega` Acentric factor

Viscosity

Letsou-Stiel [15]

`\eta = \frac{\xi_{0}+\xi_{1}}{\xi}`

`\xi_{0} = 2.648-3.725T_{r}+1.309T_{r}^{2}`

`\xi_{1} = 7.425-13.39T_{r}+5.933T_{r}^{2}`

`\xi = 176(\frac{T_{c}}{MM^{3}{P}_{c}^{4}})^{1/6}`

where

`\eta `Viscosity (Pa.s)

`P_{c} `Critical pressure (bar)

`T_{r}` Reduced temperature, `T/T_{c}`

`MM` Molecular weight (kg/kmol)

Abbott [17]

`t_{1} = 4.39371-1.94733Kw+0.12769Kw^{2}+0.00032629API^{2}-0.0118246KwAPI+ +(0.171617Kw^{2}+10.9943API+0.0950663API^{2}-0.869218KwAPI`

`\log v_{100} = \frac{t_{1}}{API+50.3642-4.78231Kw},`

`t_{2} = -0.463634-0.166532API+0.000513447API^{2}-0.00848995APIKw+ +(0.080325Kw+1.24899API+0.19768API^{2}`

`\log v_{210} = \frac{t_{2}}{API+26.786-2.6296Kw},`

where

`v_{100}` Viscosity at 100 ?F (cSt)

`v_{210}` Viscosity at 210 ?F (cSt)

`K_{w}` Watson characterization factor

API Oil API degree

User-Created Compounds

The majority of properties of the user-created compounds is calculated, when necessary, using the group contribution methods, with the UNIFAC structure of the hypo as the basis of calculation. The table below lists the properties and their calculation methods.

Hypo calculation methods

Property	Symbol	Method
Critical temperature	`T_{c}`	Joback [15]
Critical pressure	`P_{c}`	Joback [15]
Critical volume	`V_{c}`	Joback [15]
Normal boiling point	`T_{b}`	Joback [15]
Vapor pressure	`P^{pv}`	Lee-Kesler
Acentric factor	`\omega`	Lee-Kesler
Vaporization enthalpy	`\Delta H_{vap}`	Vetere [15]
Ideal gas heat capacity	`C_{p}^{gi} `	Harrison-Seaton [5]
Ideal gas enthalpy of formation	`\Delta H_{f}^{298}`	Marrero-Gani [10]

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

`\delta=(\frac{\Delta H_{v}-RT}{V_{L}})^{1/2}`

where

`\Delta H_{v}` Molar Heat of Vaporization

`V_{L}` Liquid Molar Volume at 20 ?C

References

[1] Bell, Ian H. and Wronski, Jorrit and Quoilin, Sylvain and Lemort, Vincent, "Pure and Pseudo-pure Fluid Thermophysical Property Evaluation and the Open-Source Thermophysical Property Library CoolProp", Industrial and Engineering Chemistry Research (2014), 2498--2508.

[2] K. C. Chao and J. D. Seader, "A General Correlation of Vapor-Liquid Equilibria in Hydrocarbon Mixtures", AICHE Journal (1961), 599-605.

[3] Chen, Guang-Jin and Guo, Tian-Min, "A new approach to gas hydrate modelling", Chemical Engineering Journal (1998), 145--151.

[4] H. G. Grayson and C. W. Streed, "Vapor-Liquid Equilibria for High Temperature, High Pressure Hydrogen-Hydrocarbon Systems", N/A (1963).

[5] Harrison, B. Keith and Seaton, William H., "Solution to missing group problem for estimation of ideal gas heat capacities", Industrial and Engineering Chemistry Research (1988), 1536-1540.

[6] Heidemann, Robert A. and Khalil, Ahmed M., "The calculation of critical points", AIChE Journal (1980), 769--779.

[7] Jeong Won Kang and Vladimir Diky and Michael Frenkel, "New modified {UNIFAC} parameters using critically evaluated phase equilibrium data", Fluid Phase Equilibria (2015), 128 - 141.

[8] Klauda, J.B. and Sandler, S.I., "A Fugacity Model for Gas Hydrate Phase Equilibria", Industrial and Engineering Chemistry Research (2000), 3377-3386.

[9] Lee, Byung Ik and Kesler, Michael G., "A generalized thermodynamic correlation based on three-parameter corresponding states", AIChE Journal (1975), 510--527.

[10] Marrero, Jorge and Gani, Rafiqul, "Group-contribution based estimation of pure component properties", Fluid Phase Equilibria (2001), 183--208.

[11] Michelsen, Michael and Mollerup, Jorgen, "Thermodynamic Models: Fundamentals and Computational Aspects", Tie-Line Publications (2007).

[12] Andr? Mohs and J?rgen Gmehling, "A revised {LIQUAC} and {LIFAC} model (LIQUAC*/LIFAC*) for the prediction of properties of electrolyte containing solutions", Fluid Phase Equilibria (2013), 311 - 322.

[13] Parrish, W. R. and Prausnitz, J. M., "Dissociation Pressures of Gas Hydrates Formed by Gas Mixtures", Industrial and Engineering Chemistry Process Design and Development (1972), 26-35.

[14] Peng, Ding-Yu and Robinson, Donald B., "A New Two-Constant Equation of State", Industrial and Engineering Chemistry Fundamentals (1976), 59-64.

[15] Reid, Robert, "The Properties of Gases and Liquids", McGraw-Hill (1987).

[16] Renon, H. and J. M. Prausnitz, "Local Compositions in Thermodynamic Excess Functions for Liquid Mixtures", AICHE Journal (1968), 135-144.

[17] Riazi, M.-R., "Characterization and properties of petroleum fractions", ASTM (2005).

[18] Riazi, M. R. and Al-Adwani, H. A. and Bishara, A., "The impact of characterization methods on properties of reservoir fluids and crude oils: options and restrictions", Journal of Petroleum Science and Engineering (2004), 195--207.

[19] Smith, Joseph, "Intro to Chemical Engineering Thermodynamics", McGraw-Hill Companies (1996).

[20] Soave, Giorgio, "Equilibrium constants from a modified Redlich-Kwong equation of state", Chemical Engineering Science (1972), 1197--1203.

[21] Stryjek, R. and Vera, J. H., "PRSV2: A cubic equation of state for accurate vapor?liquid equilibria calculations", The Canadian Journal of Chemical Engineering (1986), 820--826.

[22] Thomsen, K., "Aqueous electrolytes: model parameters and process simulation" (1997).

[23] Curtis H. Whitson and Michael R. Brule, "Phase Behavior (SPE Monograph Series Vol. 20)", Society of Petroleum Engineers (2000).