Documentation

Vapor pressure of a pure substance is calculated usually with the Antoine-formula, which is an empirical equation with three material-specific constants.

where:

• p_sat.,i is the vapor pressure (bar)
• T is the temperature
• A,B,C are the Antoine constants

For the majority of substances in the database of VLE - Calc .com these constants are stored. Compounds for which the Antoine constans are missing, vapor pressure is calculated by the Clausius-Clapeyron equation:

where:

• p_sat.,i is the vapor pressure (bar)
• ΔH_V,i is the heat of vaporization (ethalpy of vaporization) (J/mol)
• R is the universal gas constant (8.314 J/mol⋅K)
• T is the temperature (Kelvin)
• T_BP is the normal boiling point (at 1 atmosphere) (Kelvin)

In case of a multicomponent system, and when the temperature and pressure are not too close to the critical state (the vapor phase can be considered as ideal gas), one can use Raoult's law, modified by the activity coefficient:

where:

• p_mixture is the vapor pressure of the mixture (bar)
• γ_i(gamma) is the activity coefficient of component i (see in chapter 4.4)
• p_sat.,i is the vapor pressure of pure component i (bar)
• x_i is the liquid phase mole fraction of component i
• y_i is the vapor phase mole fraction of component i

There are substances whose molecules form associates in the vapor phase, and therefore erroneous vapor composition results using the modified Raoult's law. A typical example is the pair of acetic acid and water: the modified Raoult's law alone results in a minimal boiling azeotrope, although in reality this system is zeotropic. In this case, the association must be taken into account even at low pressures.

The method of J. Gmehling (see chapter 7) was adopted, whose essence is considering the association as an equilibrium chemical reaction.

To this equilibrium reaction we assign an equilibrium constant that depends only on the temperature:

K_D = p_i,D / p²_i,M

where p_i,M is the partial pressure of the monomers, p_i,D is the partial pressure of the dimers.

We apply the following simplifications:

• we only deal with the association of carboxylic acids
• only dimerization is assumed
• cross-dimerization (association of different compounds) is not taken into account
• the same equilibrium constant is used for all substances containing carboxylic acid groups (no substance-specific equilibrium constants are used)
• (The liquid phase association is included in the activity coefficient.)

The calculation of vapor phase dimerization goes according to the following:

• At first the equilibrium constant of the dimerization is calculated at the given temperature T: ln K_D = A_D + (B_D / T)
  (A_D and B_D are constants.)
  
  
• Then the partial pressures of the monomers and dimers have to be calculated:
  
  p_sat.,i,M = (−1 + √1 + 4⋅K_D⋅p_sat.,i) / 2⋅K_D
  p_i,M = x_i ⋅ γ_i ⋅ p_sat.,i,M
  p_i,D = K_D ⋅ p²_i,M
  
  
• Finally it should be considered that the dimer species counts as 1 molecule in the vapor pressure, but it counts as 2 molecules in the vapor composition:
  
  p_i = p_i,M + p_i,D
  y_i = (p_i,M + 2⋅p_i,D) / p_mixture
  

Model calculations of multicomponent systems usually need some method to take into account the interaction of different type materials. This program uses the so-called activity coefficient. (This is a kind of measure of interaction strongness between different molecules in a liquid solution, the "extent of non-ideality".)

VLE - Calc .com currently uses the UNIQUAC and the UNIFAC methods. The former is used for specific substances and works with interaction parameters derived from specific experimental measurements. The latter is a group-contribution method, so it can be used - theoretically - for any molecule, without having any interaction parameters based on experiments. Very simplified, UNIFAC works in a way, that molecules are defined as a set of different functional groups, and the activity coefficient is calculated from the sum of the molecular interactions between each group-pairs (using rather complex formulas).

In case of some common compound pairs the predictive power (accuracy) of UNIFAC seemed to be quite poor, so for some compound pairs UNIQUAC interaction parameters have been defined and recorded in the VLE-Calc.com database, derived from experimental data sets. When a user deals with a multicomponent system, in which all component-pairs' UNIQUAC parameters are available, the program automatically uses them instead of UNIFAC.

Two versions of UNIFAC is used by this app: one version for equation of state (EOS) usage, another version if no EOS is used.

In principle UNIFAC enables combining substances arbitrarily, unfortunately in practice it is not so simple. Certain combinations are allowed, while others are not. It depends on whether the interaction parameters are available for all functional group pairs occurring in the applied compounds, or not. Because if not, the program rather disallow the selection of such substances at the same time, in order to prevent giving erroneous results due to missing group interaction parameters.

This type of blocking is implemented in a way, that before displaying the list of substances (for compound selection) the program checks all the compounds which can be paired with the already selected one(s) and which not. Those which cannot be paired do not appear in the compound selector list, so they cannot be selected.

At normal conditions (meaning: common solvents, from room temperature up to about 100°C, at normal or reduced pressure) Raoult's law (in chapter 4.2) is fully suitable for vapor pressure calculations (naturally modified with the activity coefficient). (The vapor phase can be considered as ideal gas.)

But when the system's temperature or/and pressure is somewhat close to the critical point (or the mixture contains supercritical compontents), the vapor phase cannot be considered as ideal gas anymore. In such case an equation of state (EOS) should be used, which can describe the properties of the mixture realistic, even near to the critical state.

VLE - Calc .com uses currently the PSRK (Predictive Soave-Redlich-Kwong) equation of state. It's advantage is that it uses the UNIFAC model which could be fitted excellent to the activity coefficient part (chapter 4.4).

Several specific data are needed for EOS calculations (e.g. the critical constants of the material and others). When opening the compound selector list, the program checks all substances that the EOS-related data are available, or not. Because if not, the program either disallow the selection of such substances (by hiding them in the list), or disables the use of EOS and only enables only "normal" conditions (i.e. temperatures/pressures far from the critical state) for the VLE calculations.

Calculating the vapor pressure for a given temperature is relatively easy according to section 4.2 (Raoult's-law). The solution of its reverse in case of a single material is also relatively simple (see: both the Antoine and the Clausius-Clapeyron equations can be inverted to express ecplicitly the T). But for a multicomponent system it is not trivial how to get the boiling temperature T of p_mix.

Literature recommends in general iterative methods, just like the quasi-Newton secant method, which have been used by early versions of VLE-Calc.com. Now, the App uses a method wich is based on the fact that the vapor pressure function of ln(p) − 1/T is quasi-linear (even in case of a multi-component system), altough not in a very wide temperature range. The algorithm:

• first we calcualte the boiling temparature of all of the pure components at the give pressure
• then we calculate the vapor pressure of the mixture at the lowest and at the highest temperature
• based on the 2-2 pressure-temperature data we now fit a linear function on ln(p) − 1/T (slope and intersect)
• based on the fitted linear function one could calculate the searched temperature for the given temperature - this is the bubble point of the mixture at the given pressure.

Note: this result is relatively accurate. Error is about 0.1°C. For less demanding applications one can stop at this point. But this result can be made much more accurate at relatively low cost, with several orders of magitude.


• For the third temperature we got, we calculate again the vapor pressure of the mixture. (This should be very close to the given pressure. If not, something were wrong along the calculations.)
• Now we may fit a cubic function (instead of a linear one) based on the calculated 3-3 ln(p) and 1/T values, and we get the bubble temperature of given pressure from this.
• (When vapor composition is necessary, the one should calculate one more time the mixture vapor pressure at the final bubble temperature, and the composition obtains as a by-product of the results.)

This method is computationally cheap (thus fast), mixture vapor pressure should be calculated only four-times in order to get the most accurate result, while the secant method usually needs more.

Introduction

Various liquid materials are either miscible, or not. (But that is certain.) Well, it is well-known that hydrocarbons and lots of other organic stuffs do not mix with water, but there are plenty of cases where it is not trivial at all whether the mixture splits to multiple liquid phases, or not.

In principle, liquid phase split can be predicted from the curvature of the Gibbs free energy of mixing function (g^E = ΔG_mix.(x,T)/RT): when the function is convex from below at a composition range, i.e. ∂²g^E(x,T)/∂x² < 0, then the mixture will split to multiple liquid phases. This range is the unstable zone (see as green in the picture).

In case of two components the composition of the phases can be predicted by aligning a tangent line according to the figure (see as red line). In equilibrium the activity (a_i=γ_i⋅x_i) of all components are equal in all phases, so a computationally equivalent solution is to find the compositions with isoactivity (a_i,phase-1 = a_i,phase-2).

In practice, there is some trouble with this. One can notice in the figure that the g^E function above the red tangent line is not everywhere convex from below (see as blue in the picture). This is the metastable zone. If the gross composition of a mixture falls into the metastable zone, testing the sign of the second derivate will not show that the mixture is biphasic.

In case of three or even more components the situation becomes more complicated, because the concepts (zones) described above should be interpreted and computationally handled on two or more dimensional surfaces. Practice of this is very not trivial, several article were published in this topic in the recent decades. Former versions of VLE - Calc .com used to use methods based on second derivate tests and isoactivity searches described above, thus the algorithms could not detect heterogenity in all cases. Another critical problem was that these algorithms produced sometimes erroneous equilibrium concentrations, the tie lines on the LLE-diagrams did not lined up consistently, even crossing tie-lines occurred. (The latter was not due to a bug in the program code, but due to properties of the thermodynamic model - see later in this chapter).

Liquid stability test

Due to the problems described in the introduction, the method of Wasylkievich et al were adopted to current version of VLE - Calc .com. This one is a very sophisticated approach that seemes to be reliable - though it was extremely difficult to reproduce and implement here.

The essence of the method is that the local minima (or inflexion points) of the so-called Tangent-Plane-Distance-Function (TPDF) have to be located, and the equilibral compositions of the liquid phases should be searched starting from the found local minima. The method has several advantage:

• indicates heterogenity even in the metastable zone
• can be used for any number of components (in theory)
• it provides very good starting points for finding the equilibrium compositions
• indicates when the found equilibral composition are false
• indicates when the system splits into more than 2 phases (in case of more than 2 components)

At a given T temperature and z gross mixture composition:

TPDF(x) = g^E(x) − L(x, z)

where L(x, z) is a function tangent to g^E at z gross mixture composition.

Is TPDF(x, z) negative in any x, then the liquid mixture is unstable, and will split to multiple phases. (See with blue on the TPDF figure)

Two-component systems

In case of a binary system the local minima of the TPDF (as well as inflexions and local maximum) can be found at the zeros of the 1^st derivates. In current implementation these roots are searched by a quasi-Newton method. Not all roots will be searched, but the two probably closest to the minima of the TPDF.

In this method we exploit some useful properties of the binary g^E function:

• there is no always local maximum, but if it exist, there is always only one.
• the slope of the function from the borders of the composition space into inward direction is always negative

Thus even if a local maximum exist, starting a search from the borders one would find always local minuma or inflexion point. The 1^st derivate of the TPDF can be expressed by an approximate analytical derivative, and this expression should be made equal to zero:

∂TPDF(x_1,2) / ∂x₁ ≅  (ln x₁ + ln γ₁(x) −  ln z₁ − ln γ₁(z)) −  (ln x₂ + ln γ₂(x) −  ln z₂ − ln γ₂(z)) = 0

The algorithm:

• Compute the activity coefficients at gross composition z, and store them (they are constants during the processγ_1,2(z)).
• Start "guessed x"with point x₁=0.
• Compute the activity coefficients at current x (γ_1,2(x)).
• Using the analytical derivate described above one can express the formulae below, from which the new x can be obtained:
  
  c₁ = ln γ₁(x) −  ln z₁ − ln γ₁(z)
  c₂ = ln γ₂(x) −  ln z₂ − ln γ₂(z)
  x_1,NEW = exp (−c₁ +  ln x₂ + c₂)
  x_2,NEW = 1 − x_1,NEW
  
  
• Calculate the activity coefficients again with the new x value, and iterate the process until convergence. The resultant x is the closest minimum or inflexion to x₁=0.
• Repeat the process described above starting from x₁=1. From this process can be obtained the second local minimum or inflexion - if it exist. Namely, if there is only one minimum and no inflexion then both process gets to the same result, hence the system is homogenous.

Check the sign of the TPDF at the resultant x values. If it is negative, the system is unstable (heterogenous, biphasic). Now proceed searching for equilibral compositions of the phases, using these x values as starting points (this process is the "LLE-flash"). In general the Rachford-Rice equation is used for that (see below).

Systems with three or more components

As mentioned in the introduction, the situation in case of a multicomponent system is much more complicated as in binary case. The method of S. Wasylkievich et al handles the situation with high reliability, in case of any number of components. However the process is computationally quite intensive. I would not go into the formulas here anymore, everything is described in detail in the cited article.

The essence of the method is that in the case of n component the a TPDF is a (n−1)-dimensional surface, and the local minima of it has to be find. This surface is "crumpled cloth"-like, it has local minima (full red dots on the figure), local maximum can occur on it, and there can be saddle points (empty red dots on the figure). These points are called the stationary points of TPDF. In these points all the partial derivatives of the TPDF, by mole fractions of all i components, are zero ( ∂ TPDF(x) / ∂x_i = 0). These stationary points are connected by "ridges" and "valleys", and according to the method, theese points are searched by tracking these ridges and valleys, along the whole (n−1)-dimensional composition space.

Once the minima are found, similarly to the binary case, the sign of TPDF at the minimums are checked, whether it is negative or positive. If any of the minima is negative, the system is unstable (heterogenous). Then proceed searching for the equilibral compositions of the phases, using the x values of the local minima as starting points (this process is the "LLE-flash"), and with solving the Rachford-Rice equation is used for that (see below).

Unfortunately, in multicomponent case, the first obtained phase compositions are not surely correct - it has to be checked. It can happen that not only one, but more isoactivity point-pairs can be found for a given z gross composition. From these solutions, the one with the lowest total Gibbs free energy of the system is correct.

So the obtained phase compositions have to be tested for stability, according to the ridge tracking technique described above (let z = x_{one-of-the-phases}). If x_{one-of-the-phases} is stable (i.e. no negative local minimum found), then the found equilibral phase compositions are correct. If the stability test results unstable (negative local minimum are found), then the found equilibral phase compositions are incorrect, and the LLE-flash process has to be repeated using the mole fractions of the local minima obtained in the last stability test, as initial values.

From the number of local minima, it can be assumed how many liquid phases are, but there is no guarantee that it is actually as much. E.g. if there are three minima, it is worth assuming this, and then look for three equilibrium phases based on the minima. If there were no three phases, it would be revealed either in the Rachford-Rice calculations (the system automatically reduces to 2-phase) or in the stability test of the found phase compositions.

The Rachford-Rice equation

Once it is known that the system is heterogenous, it is time to calculate the accurate composition of the coexisting luqid phases. In equilibrium the activity of the components (a_i=γ_i⋅x_i) is equal in all phases. At the same time it has to be also take into account that the material balance of each component must be in order, i.e. the total amount of each component in the separate phases must be equal to that of in the gross mixture (z). Rachford and Rice derived a very useful little equation for this task, solving it leads to results that meet both expectations at the same time:

• F(α) is the objective function, which has to be equal to zero
• α is a number between 0 and 1, which measures the molar amount of the 1^st separate phase, from 1 mole gross mixture. (The molar amount of the 2^nd phase is logically 1−α.)
• z_i is the mole fraction of component i in the gross mixture
• K_i is the distribution ratio of component i: K_i = x_i,phase 1/x_i,phase 2 = γ_i,phase 2/γ_i,phase 1

The Rachford-Rice equation can be solved by iterative way:

• In the first step calculate the distribution ratios K_i using the assumed phase compositions and activity coefficients originating from the liquid stability test.
• Determine α by the Newton-Raphson method. Substitue everything into F(α) and into it's first derivate according to α, and the next approximation of α can be obtained by the formula: α_new = α_previous − F(α_previous) / F'(α_previous). Repeat with the new α and iterate until convergence.
• Calculate the phase compositions with the formulae:
  x_i,phase 1 = z_i / (1 + α(K_i − 1))
  x_i,phase 2 = x_i,phase 1 ⋅ K_i
• Normalize the obtained mole fractions (for all phases): x_i = x_i / ∑ⁿ_j=1 x_j
• Calculate the activity coefficients using the new phase compositions, and calculate again all K_i, and iterate the whole process until convergence.

Boiling point / vapor pressure map by contour lines

Although it is not the most useful representation, but this type is the oldest triangle diagram in this app, since it is relatively simple and can be obtained quickly. Basically it represents the boiling point (bubble point, in isobaric case) or the vapor pressure (in isothermal case) of a ternary mixture, using contour lines. It's disadvantage is that the vapor composition (as additional information) does not appear in any form in this mode of representation.

At first, the triangular composition space is divided into many small triangles with a resolution of 20x20, and a VLE calculation is performed at each corner point of each small triangle. The minimum and maximum temperature or vapor pressure values are determined from all of the VLE results, and 8 values are selected between them, for the 8 contour lines. These are the "constant values", since these values are constant along the contour lines. With each constant value, all the small triangles are checked wether this value appears in it, or not. If a constant falls between the lowest and highest temperature/pressure values at the corners of the small triangle, then the contour line of the constant value passes through the given small triangle. If so, the position of the crossing section in the triangle is calculated by linear interpolation, and stored.

Plotting all of the crossing sections the contour lines outlined.

Residue curve map

For distillation problems, the residue curve is a much more useful representation than the boiling point map. The residue curve lies along the path as the composition of the bottom product changes during a batch distillation process. Logically, it moves from the more volatile regions to the less volatile regions, and the direction of the curve at each point is determined by the equilibrium x − y vector, characteristic of that point.

The shape of residue curve diagrams can be very diverse. In fact, the azeotropic points of a given ternary system, basically determine the shape of the figure. The azeotropic points (together with the corner points of the entire triangular diagram, which represent the pure components) divide the total composition space into smaller areas. Simple batch distillation processes can operate generally in these areas, in which the starting composition falls, so with the help of the residue curve set, it can be easily determined visually whether a given separation task is possible or not.

At first the azeotropic points of the system are located, then the corner points of the whole triangular diagram are added to the set. These points are called nodes. The set of nodes is arranged in triangles by Delaunay triangulation. Within each small triangle, residue curves are generated starting from one or more points.

During generation of a residue curve, the area of the small triangle is scanned from the starting point first in the x → y direction and then in the x ← y direction until a node is reached. Then another starting point is chosen and scan again.

Liquid miscibility phase-diagram

Liquid-liquid equilibrium diagrams have numerous types, these are not discussed here just the method of figure making.

At first the algorithm tests for liquid stability the binary subsystems representing the edges of the triangle diagram. The phase compositions of the miscibility gaps of the heterogenous binaries are stored for drawing on the one hand, on the other hand they are used for define starting points in the middle of the gaps.

Then it starts stepping from the starting points in small steps towards the inside of the triangle and determines the ternary tie-lines or conodes. The direction is determined by the slope of the tie-line (the direction is always perpendicular to the last tie-line), and the step size depends on the length of the tie-line - the shorter the tie-line, the smaller the step size. (These rules are basically for aesthetic purposes.)

If the algorithm finds a three-phase region, it continues to stepping from the middle of the other two edges of the triangle formed by the compositions of the three phases. Finally, plots all the found tie-lines.

In case of ternary vapor-liquid equilibrium calculations, the liquid-liquid equilibrium tie-lines are precalculated, stored, and arranged in small and thin triangles. Before starting the VLE calculation, the program checks if the given liquid composition falls into any of the LLE triangles. If so, it determines the composition of the separated liquid phases in a simplified procedure and performs VLE calculations with these compositions. (This precalculation and store-the-LLE-map mess is necessary because the liquid stability test is computationally very intensive. Less stability testing is required for an LLE diagram than for a batch distillation.)

During a single stage continuous distillation (flash distillation) process the starting liquid mixture (feed stream) is evaporated partially. The composition of the head- and bottom product depends on the evaporated proportion. In case of total evaporation (no bottom product) the composition of the distillate is the same as that of the starting mixture. In case of no evaporation (only bottom product), logically the bottom product is the same as the feed. The intermediate (real) situation is calculated using the Rachford-Rice equation (just like in case of liquid-liquid equilibrium).

where

• F(α) is the objective function, which has to be equal to zero
• α is the evaporated proportion, a number between 0 and 1, the molar amount of the vapor phase from 1 mol feed mixture. (The molar amount of the bottom product
• z_i is the mole fraction of component i in the feed stream
• K_i is the distribution ratio of component i between the vapor and liquid phase: K_i = y_i/x_i

The Rachford-Rice equation can be solved by iterative way:

• The evaporated proportion α is given, constant.
• At first perform a bubble point calculation using the z feed composition.
• Calculate distribution ratios (as astarting approximation) using the resultant y values: K_i = y_i/z_i
• Calculate the composition of vapor and liquid phase:
  x_i = z_i / (1 + α(K_i − 1))
  y_i = x_i ⋅ K_i
• Normalize the mol fractions:
  x_i = x_i / ∑ⁿ_j=1 x_j
  y_i = y_i / ∑ⁿ_j=1 y_j
• Perform another bubble point calculation using the resultant x values, and then recalculate the distribution ratios with the new y values (K_i = y_i/x_i), and iterate until convergence.

Batch distillation is a timely process, the composition of both the boiling liquid and the distillate changes over time. To model this, the process is divided into n equal parts. VLE-Calc.com always uses a resolution of n = 200, even if a lower resolution is given in the form, and only shows as much of the results in the table as specified. The program always performs the calculations in molar units.

The initial x composition of the mixture is given. The amount of the starting mixture is always n = 200 mol as described above (this made the code simpler), which is of course converted at the end according to the amount (and unit) specified by the user. A bubble point calculation is performed from x, and with the obtained vapor composition y one can calculate the amount of all components from the first evaporated dn = 1 mol:

n_new = n_old − dn
x_i,new = (x_i,old ⋅ n_old − y_i ⋅ dn) / n_new

Store the new boiling liquid composition and calculate from this again the vapor composition, and so on, iterate until the liquid runs out completely.