A new numerical procedure to determine the VLE curve

doi:10.1016/S0097-8485(00)00112-1

Computers & Chemistry

Volume 25, Issue 5, September 2001, Pages 483-488

https://doi.org/10.1016/S0097-8485(00)00112-1 Get rights and content

Abstract

A fast and accurate (to any required precision) numerical method is presented to calculate the vapour-liquid equilibrium for any subcritical temperature from a given equation of state. The calculations are made using the popular program mathematica, and we think that this method could be very useful in chemical physics and chemical engineering applications, as well as for teaching purposes.

Introduction

The knowledge of phase transitions of pure or mixed fluids is essential in many practical applications. Engineers and applied physicists and chemists commonly use an equation of state (EOS) to evaluate the vapour-liquid equilibrium (VLE) properties. However, this calculation is no easy task numerically, and the shortcut techniques that have been developed are not applicable to every EOS (Wisniak et al., 1998).

As is well known, to determine the VLE curve of a pure fluid, it is necessary to solve the equations $P ρ_{V},T =P ρ_{L},T μ ρ_{V},T =μ ρ_{L},T$ where the pressure, P, and the chemical potential, μ, are analytical functions of temperature, T, and density, ρ.

One way to obtain the VLE properties from a given EOS is through the equality of the areas for each temperature. This is equivalent to the equality of the chemical potentials corresponding to the gas and liquid phases at a given temperature. In practice however, the accurate determination of the position of the straight section of the isotherm for which the pressures and areas must be equal is no trivial problem, because the isotherm has a steep slope for densities (molar volumes) corresponding to the liquid phase. A small deviation in the equality of the areas involves small deviations in the densities ρ_L of the liquid phase and ρ_V of the gas phase. Because of the steepness of the slope of the isotherm in the liquid region, the deviation in the determination of P_L is extraordinarily sensitive to the accuracy of the determination of the liquid density, ρ_L, and hence, so is the accuracy in the determination of the mechanical equilibrium between the two phases.

This problem may be avoided by solving simultaneously the system of equations (1), to which end one would need analytical expressions not only for the EOS but also for the chemical potential. Whereas there have been numerous EOSs in the literature, analytical expressions for the chemical potential are more difficult to find (Mulero et al., 1999).

Various EOSs have been constructed on the basis of intensive computer simulation data for Lennard–Jones (LJ) fluids (Johnson et al., 1993, Kolafa and Nezbeda, 1994, Cuadros et al., 1996). But the complexity of the corresponding $P=P ρ,T$ and $μ=μ ρ,T$ expressions makes it very difficult to solve system (1) analytically. Instead, it is necessary to use numerical techniques. Thus, to determine each pair of points $ρ_{V},ρ_{L}$ , one must solve system (1) for a given temperature, with the inherent associated problem of choosing the initial value, especially in the neighbourhood of the critical point.

We present here an easier and faster method of solving system (1). This method yields practically the entire VLE curve, without having to solve the system point by point. It is based on numerically solving a system of nonlinear differential equations using the program mathematica (Wolfram, 1991). This program offers to the user an interactive tool for numerical calculations and a versatile programming language for fast and accurate solutions to the most diverse problems. For these reasons, it is becoming ever more commonly used by scientists, engineers, and higher education students.

The main goal of the present paper is to describe a new procedure to calculate the coexistence curve between the liquid and gas phases of a LJ system. The method could be extrapolated to any system for which one knows the EOS and the expression $μ=μ ρ,T .$

Section snippets

Theoretical basis of the method

For any T below the critical temperature T_C, the functions $ρ_{V} =ρ_{V} T$ and $ρ_{L} =ρ_{L} T$ satisfy Eq. (1). Since P and μ are known analytically, one then differentiates Eq. (1) with respect to T: $d P ρ_{V} T,T d T = d P ρ_{L} T,T d T d μ ρ_{V} T,T d T = d μ ρ_{L} T,T d T$ and hence $∂P ∂T ρ_{V} T,T + ∂P ∂ρ ρ_{V} T,T ρ ̇_{V} = ∂P ∂T ρ_{L} T,T + ∂P ∂ρ ρ_{L} T,T ρ ̇_{L} ∂μ ∂T ρ_{V} T,T + ∂μ ∂ρ ρ_{V} T,T ρ ̇_{V} = ∂μ ∂T ρ_{L} T,T + ∂μ ∂ρ ρ_{L} T,T ρ ̇_{L}$ where $ρ ̇_{V}$ and $ρ ̇_{L}$ denote the derivatives of ρ_V and ρ_L with respect to T.

From Eq. (2), one obtains the system of two first-order nonlinear differential equations $ρ ̇_{V} =f ρ_{V},ρ_{L},T ρ ̇_{L} =f ρ_{L},ρ$

Algorithm

1.
Initiate T₀.
2.
Define the functions $P ρ,T$ and $μ ρ,T .$
3.
Search for a numerical solution, using Newton's method, to the simultaneous equations $P ρ_{V},T_{0} =P ρ_{L},T_{0}$ and $μ ρ_{V},T_{0} =μ ρ_{L},T_{0}$
4.
Assign the solution's values to ρ_V₀ and ρ_L₀.
5.
Find a numerical solution to the ordinary differential equations $d P ρ_{V} T,T d T = d P ρ_{L} T,T d T,$ $d P μ_{V} T,T d T = d P μ_{L} T,T d T,$ $ρ_{V} T_{0} =ρ_{V_{0}}$ and $ρ_{L} T_{0} =ρ_{L_{0}}$ , for the functions ρ_V and ρ_L with the independent variable T in the range 0.7 to 1.35.

As we noted above, to solve the system of two nonlinear differential equations

Applications to calculating the VLE curve

For our calculations, we chose three EOS models: Johnson et al., 1993, Kolafa and Nezbeda, 1994, and Cuadros et al., 1996. All the model EOSs will be expressed in reduced LJ units.

Conclusions

1.
We have described a new approach to calculating the VLE curve, system (1), which is mathematically equivalent to the classical requirement of solving the system that results from differentiating with respect to T.
2.
If the EOS is sufficiently regular, with the imposition of determined initial conditions $T_{0},ρ_{V} T_{0},ρ_{L} T_{0}$ , the solution of the system (2) exists and is unique.
3.
The complexity of the EOS expressions makes it practically impossible to obtain analytical formulae for ρ_V and ρ_L We therefore

Acknowledgements

The authors express their gratitude to the Consejerı́a de Educación y Juventud of the Junta de Extremadura and the Fondo Social Europeo for the financial aid granted through the Project IPR98B004. Also, M.I. Parra would like to express her gratitude to the Socrates–Erasmus program for its support through a grant, as well as for facilitating her academic visit to the Institute of Mathematics of the Politecnika Zielonogórska.

References (10)

J. Kolafa et al.
Fluid Phase Equilibria
(1994)
J. Wisniak et al.
Chem. Eng. Sci.
(1998)
N.F. Carnahan et al.
J. Chem. Phys.
(1969)
F. Cuadros et al.
Phase Transitions
(1996)
P. Hartmann
Ordinary Differential Equations
(1964)

There are more references available in the full text version of this article.

Cited by (7)

Simple approach to approximate predictions of the vapor-liquid equilibrium curve near the critical point and its application to Lennard-Jones fluids
2012, Physics Letters, Section A: General, Atomic and Solid State Physics
Citation Excerpt :
The result of this procedure was to remove a small amount of data ‘scatter’ of graphs illustrating the percentage deviation of the values being compared (see Figs. 2, 4). Argon and methane are the subject of numerous scientific studies [16–20], providing coverage of both theoretical and experimental studies. It has been payed much attention to the thermodynamic properties, due to widespread use of this element both in industry and science.
We propose a simple analytical form of the vapor–liquid equilibrium curve near the critical point for Lennard-Jones fluids. Coexistence densities curves and vapor pressure have been determined using the Van der Waals and Dieterici equation of state. In described method the Bernoulli differential equations, critical exponent theory and some type of Maxwellʼs criterion have been used. Presented approach has not yet been used to determine analytical form of phase curves as done in this Letter. Lennard-Jones fluids have been considered for analysis. Comparison with experimental data is done. The accuracy of the method is described.
Many-body and quantum effects in the surface tension and surface energy of liquid neon and argon using the Fowler's approximation
2012, Chemical Physics
Citation Excerpt :
Although the RDF expression of Cuadros et al. [47] is a complex analytical model, but it is (such as the RDF equation of Morsali et al. [40]) in good agreement with the real RDF obtained in computer simulations [9,47]. Therefore, the difference between the calculated values in this work and those of Mulero et al. [9] can be due to the difference between the experimental densities [46] and the density equations of Okrasinski et al. [48,49]. In the other words, the Okrasinski et al. expressions produce the density values somewhat bigger than those of Ref. [46].
The Fowler’s expression for calculation of the reduced surface tension and surface energy has been used with Lennard–Jones (LJ) and two-body Hartree–Fock dispersion (HFD)-like potentials for neon and argon, respectively. The required radial distribution functions (RDFs) have been used from two recently determined expressions in the literature and a new equation proposed in this work. Quantum corrections for neon system have been considered using the Feynman–Hibbs (FH) and Wigner–Kirkwood (WK) approaches. To take many-body forces into account for argon system, the simple three-body potentials of Wang and Sadus (2006) [33] and Hauschild and Prausnitz (1993) [30] used with the HFD-like potential without requiring an expensive three-body calculation. The results show that the quantum and three-body effects improve the prediction of the surface tension of liquid neon and argon using the Fowler’s expression.
Modeling and simulation of supercritical extraction columns using computer algebra
2005, Journal of Supercritical Fluids
A procedure for the simulation of extraction columns is developed using a program, Thermath, based on the Mathematica^® computer algebra (CA) software. The procedure requires the simultaneous solution, using the Newton–Raphson method, of mass and energy balances, phase equilibrium equations, and equations needed to match the number of degrees of freedom at each stage. Starting from the equipment model equations, analytical derivatives with respect to the problem unknowns were obtained using CA, and automatically implemented in FORTRAN using Thermath. Similarly, fugacity coefficients and residual enthalpies were derived from the Peng–Robinson and the Soave–Redlich–Kwong equations of state, using either the classical or the Huron–Vidal mixing rules, and also implemented using Thermath. The code was used to simulate the separation of products from an enzymatic reaction using supercritical carbon dioxide. Extraction columns for the systems carbon dioxide–ethanol–water and carbon dioxide–limonene–linalool were also simulated. The results show trends compatible with those that should be expected from supercritical extraction columns.
The slowly varying density approximation for the surface tension and surface energy of Lennard-Jones fluids revisited
2003, Journal of Physics Condensed Matter
Fowler's approximation for the surface tension and surface energy of Lennard-Jones fluids revisited
2003, Journal of Physics Condensed Matter
Calculation of the surface tension and the surface energy of Lennard-Jones fluids from the radial distribution function in the interface zone
2006, Molecular Physics

View all citing articles on Scopus

View full text

A new numerical procedure to determine the VLE curve

Abstract

Introduction

Section snippets

Theoretical basis of the method

Algorithm

Applications to calculating the VLE curve

Conclusions

Acknowledgements

Fluid Phase Equilibria

Chem. Eng. Sci.

J. Chem. Phys.

Phase Transitions

Ordinary Differential Equations