MATLAB Examples

Antoine Class for Vapor-Liquid Equilibrium Calculations

The Antoine class provides a computational framework to help doing simple vapor-liquid equilibrium calculations in Matlab. These notes introduce Antoine's equation, then shows how the Antoine class can be used to do typical vapor-liquid equilibrium calculations.

Contents

These are a few functions and variables useful in formatting the output of subsequenct calculations.

clear all;
clc;

% Utility functions for formatting computational output.
sep = '------------------------------------------------------------------';
problem = @(s) disp(sprintf('\n%s\nProblem: %s\n%s',sep,s,sep));
answer  = @(str,val,units) disp(sprintf('%s = %6.3g [%s]',str,val,units));

Review: Saturation Pressure of a Pure Component

Antoine's equation is a representation of the vapor-liquid equilibrium for a pure component. A common form of the equation is

$$\log_{10}P_{Sat} = A - \frac{B}{T+C}$$

where coefficients A,B, and C are determined from experimental data. Values of these coefficients for a large number of chemical species can be found from many sources, including http://webbook.nist.gov/chemistry. Commonly used units are [mmHg] for pressure and [deg C] for temperature.

Here we demonstrate use of Antoine's equation to plot the saturation pressure of pure water as a function of temperature.

% For water between 60 and 150 degrees C

A = 7.96681;
B = 1668.21;
C = 228;

% Temperature range

T = 60:150;

% Antoine's equation

log10P = A - B./(T+C);

% Computing P

Psat = 10.^log10P;

% Construct the desired plot

plot(T,Psat);
xlabel('Temperature [C]');
ylabel('Pressure [mmHg]');
title('Saturation Pressure from Antoine Equation');

Antoine Class: Object Oriented Vapor/Liquid Calculations

Let's now introduce the Antoine class object. Recent versions of Matlab include methods for objected oriented programming. These methods encapsulate data and functions for in convenient tools for engineering calculations. Objects provide an elegant way to implement engineering calculations once you get used to the concepts and syntax,

The Antoine Class Dataset

The Antoine class includes parameters for a small number of compounds. These can be accessed as indicated below.

% Load available data into a structure p

p = Antoine.data;

% Show list of available compounds

problem('Antoine Class Dataset');
disp(p);

% Display Antoine data for one member of the dataset

problem('Antoine Data for Ammonia');
disp(p.nh3)
disp(sep);
------------------------------------------------------------------
Problem: Antoine Class Dataset
------------------------------------------------------------------
     h2o: [1x1 Antoine]
      n2: [1x1 Antoine]
      o2: [1x1 Antoine]
     nh3: [1x1 Antoine]
    meoh: [1x1 Antoine]
    etoh: [1x1 Antoine]
     bnz: [1x1 Antoine]
     tol: [1x1 Antoine]
     CH4: [1x1 Antoine]
    C2H6: [1x1 Antoine]
    C3H8: [1x1 Antoine]
     nC4: [1x1 Antoine]
     nC5: [1x1 Antoine]
     nC6: [1x1 Antoine]
     nC7: [1x1 Antoine]
     nC8: [1x1 Antoine]


------------------------------------------------------------------
Problem: Antoine Data for Ammonia
------------------------------------------------------------------
  Antoine

  Properties:
    Tmin: -83
    Tmax: 60
       A: 7.3605
       B: 926.1320
       C: 240.1700
    Pmin: 29.3732
    Pmax: 1.8843e+04


------------------------------------------------------------------

Creating Objects for Compounds not in the Antoine Class Dataset

Because only a small number of compounds are included in the dataset included with the Antoine class, users will generally need to create Antoine objects using data gleaned from other sources. The syntax for creating an Antoine object is

     x = Antoine(Tmin, Tmax, A, B, C)

where [Tmin,Tmax] is the temperature range, and A, B, C are the Antoine coefficients.

The following codes create two Antoine objects corresponding to benzene and toluene, respectively. Data comes from Table B.4 of Murphy's textbook.

bnz = Antoine(  8, 113, 6.90656, 1211.033, 220.79)
tol = Antoine(  6, 137, 6.95464, 1344.8,   219.48)
bnz = 

  Antoine

  Properties:
    Tmin: 8
    Tmax: 113
       A: 6.9066
       B: 1.2110e+03
       C: 220.7900
    Pmin: 41.0537
    Pmax: 1.8986e+03



tol = 

  Antoine

  Properties:
    Tmin: 6
    Tmax: 137
       A: 6.9546
       B: 1.3448e+03
       C: 219.4800
    Pmin: 9.7831
    Pmax: 1.5212e+03


Piecewise specification of Antoine Parameters

Sometimes the Antoine parameters are given in 'piecewise' fashion over multiple temperature intervals. The Antoine class will handle these cases provided the temperature intervals are (1) contiguous and (2) non-overlapping. An example

h2o  = Antoine([  0,   60,  8.10785, 1750.286, 235.0; ...
                 60,  150,  7.96681, 1668.210, 228.0])

h2o.plot;
h2o = 

  Antoine

  Properties:
    Tmin: [2x1 double]
    Tmax: [2x1 double]
       A: [2x1 double]
       B: [2x1 double]
       C: [2x1 double]
    Pmin: [2x1 double]
    Pmax: [2x1 double]


Searching the NIST Chemistry Webbook

The NIST Chemistry Webbook (http://webbook.nist.gov/chemistry/) is an excellent source of thermodynamic data. The Antoine.search method streamlines searching the webbook, and the Antoine.SI creates an Antoine object where parameters are in SI units (K, bar) rather than the more traditional (C, mmHg).

% Search for Antoine parameters for propane

Antoine.search('butane');

% Create Antoine object for butane given parameters in SI units

nC4  = Antoine.SI([ 135.42, 212.89, 4.70812, 1200.475, -13.013; ...
                    212.89, 272.66, 3.85002,  909.650, -36.146; ...
                    272.66, 425,    4.35576, 1175.581,  -2.071]);

Psat: Method for Computing Saturation Pressure

% The Antoine class includes methods for several standard computations.
% This demonstrate calculation of saturation pressure for benzene and
% toluene using the |Psat| method.

T = 70:110;
plot(T, bnz.Psat(T), T, tol.Psat(T));

legend('Benzene','Tolune','Location','NW');
xlabel('Temperature [deg C]');
ylabel('Pressure [mmHg]');
title('Saturation Pressure');

Tsat: Method for Computing Saturation Temperature

The Tsat method computes that saturation temperature as a function of pressure. This is a convenient way to compute normal boiling points.

% Retrieve Antoine parameters for water
p = Antoine.data;
h2o = p.h2o;

% Use Tsat method to report boiling point of water
problem('Normal Boiling Points');
answer('Normal Boiling Point of Water',p.h2o.Tsat(760),'deg C');
disp(sep);

% Use Tsat method to report normal boiling points for all compounds in the
% Antoine class data set.
disp('Normal Boiling Points of Species in Antoine''s Database');
f = fields(p);
for i = 1:length(f)
    disp(sprintf('%6s  %7.2f  [deg C]',f{i},p.(f{i}).Tsat(760)))
end
disp(sep)
------------------------------------------------------------------
Problem: Normal Boiling Points
------------------------------------------------------------------
Normal Boiling Point of Water =    100 [deg C]
------------------------------------------------------------------
Normal Boiling Points of Species in Antoine's Database
   h2o   100.00  [deg C]
    n2  -195.42  [deg C]
    o2  -182.85  [deg C]
   nh3   -33.43  [deg C]
  meoh    64.71  [deg C]
  etoh    78.33  [deg C]
   bnz    80.03  [deg C]
   tol   110.63  [deg C]
   CH4  -161.45  [deg C]
  C2H6   -88.67  [deg C]
  C3H8   -42.73  [deg C]
   nC4    -0.38  [deg C]
   nC5    35.06  [deg C]
   nC6    68.73  [deg C]
   nC7    98.43  [deg C]
   nC8   125.68  [deg C]
------------------------------------------------------------------

plot: Plotting Saturation Pressure versus Temperature

The Antoine class includes a simple plotting function for saturation pressure.

p = Antoine.data;

clf;
hold on;
p.n2.plot([],'b');
p.o2.plot([],'r');
hold off;
legend('Nitrogen','Oxygen','Location','NW')
grid

Application: Normal Boiling Point of a Pure Component

The boiling point of a pure component is the temperature at which the saturation pressure is equal to the ambient pressure. A simple way to find the boiling point is to use Matlab's fzero function to solve Antoine's equation for temperature as a function of pressure. To use fzero, we need a function which evaluates to 0 for the desired temperature and pressure. Here we use an anonymous function for this purpose (if you haven't encountered these before, anonymous functions are useful Matlab technique for representing and manipulating functions that can be represented by a single line of Matlab code).

% Antoine's equation in a form f(T,P) = 0

f = @(P,T) log10(P) - A + B/(T+C);

% The normal boiling point of water is the solution to f(T,760) = 0

Tb = fzero(@(T)f(760,T),[60,150]);

problem('Normal Boiling Point of Water');
answer('Normal boiling point of water (fzero)      ',Tb,'deg C');

% The |Tsat| method provides a direct method for computing the saturation
% temperature as a function of pressure.

p = Antoine.data;
answer('Normal boiling point of water (Tsat method)',p.h2o.Tsat(760),'deg C');
disp(sep);
------------------------------------------------------------------
Problem: Normal Boiling Point of Water
------------------------------------------------------------------
Normal boiling point of water (fzero)       =    100 [deg C]
Normal boiling point of water (Tsat method) =    100 [deg C]
------------------------------------------------------------------

Application: What is the boiling point of water in Littleton, Colorado?

The average barometric pressure in Littleton, Colorado, is 24.63 in of Hg (compared to sea level 29.92 in Hg). What is the boiling point of water at this reduced pressure?

% Convert inHg to mmHg (25.4 mm/in)

P = 24.63*25.4;

% Boiling point of water

problem('Boiling Point of Water in Littleton, Colorado');
answer('Atmospheric Pressure',P,'mm Hg');
answer('       Boiling Point', p.h2o.Tsat(P),'deg C');
disp(sep);
------------------------------------------------------------------
Problem: Boiling Point of Water in Littleton, Colorado
------------------------------------------------------------------
Atmospheric Pressure =    626 [mm Hg]
       Boiling Point =   94.6 [deg C]
------------------------------------------------------------------

Application: Bubble Point Calculations

Given a liquid phase mixture, the bubble point is the temperature at which the first bubble of vapor appears, and the composition of that bubble. To estimate the bubble point temperature and composition, we invoke three assumptions:

  • Antoine's equation - provides an estimate of the pure component saturation pressure.
  • Ideal gas - partial pressure of a component is equal to the vapor phase mole fraction times total pressure, $P_i = y_i P$
  • Raoult's law - partial pressure of a component is equal to the liquid phase mole fraction times the saturation pressure, $P_i = x_i P^{Sat}_i$

The last two assumptions lead to a relationship

$$ y_i P = x_i P^{Sat}_i(T)$$

which is the key to these vapor/liquid equilibrium problems. For bubble points, we know the pressure $P$ and liquid phase mole fractions $x_i$. What's left to do is adjust the temperature until the vapor phase mole fractions add to one.

$$ \sum_i x_i\frac{P^{Sat}_i(T)}{P} = 1 $$

Here we show how to formulate and solve for the bubble point of 50/50 mol% of benzene and toluene at 1 atm.

% Use z to denote feed composition.

P = 760;
z_bnz = 0.5;
z_tol = 0.5;

% Using Raoult's law, solve for vapor phase mole fractions has a function
% of temperature

y_bnz = @(T)z_bnz*bnz.Psat(T)/P;
y_tol = @(T)z_tol*tol.Psat(T)/P;

% The bubble point is the temperature at which the vapor phase mole
% fractions sum to one.

Tbub = fzero(@(T)y_bnz(T) + y_tol(T) - 1, 100);

problem('Bubble Point Calculation')
answer('z_bnz',z_bnz,'mole fraction');
answer('z_tol',z_tol,'mole fraction');
disp(sep);
answer('Bubble Point Temperature',Tbub,'deg C');
answer('Mole Fraction Benzene',y_bnz(Tbub),'mole fraction');
answer('Mole Fraction Toluene',y_tol(Tbub),'mole fraction');

% A graphical interpretation of the solution process

T = 70:110;
plot(T,y_bnz(T)+y_tol(T),T,y_bnz(T),T,y_tol(T));
xlabel('Temperature [deg C]');
ylabel('Mole Fraction');
title('Benzene and Toluene vapor phase mole fractions');
legend('y_{bnz}+y_{tol}','y_{bnz}','y_{tol}','Location','NW');
ax = axis;
hold on;
plot([Tbub Tbub ax(1)],[ax(3) 1 1],'--');
plot(Tbub,1,'.','Markersize',15);
text(Tbub,1.05,sprintf('   T_{bubble} = %5.2f deg C',Tbub), ...
    'HorizontalAlignment','center');
plot(Tbub,y_bnz(Tbub),'g.','Markersize',15);
plot([Tbub,ax(1)],y_bnz(Tbub)*[1 1],'g--');
text(Tbub,y_bnz(Tbub),sprintf('   y_{bnz} = %5.3f',y_bnz(Tbub)));
plot(Tbub,y_tol(Tbub),'r.','Markersize',15);
plot([Tbub,ax(1)],y_tol(Tbub)*[1 1],'r--');
text(Tbub,y_tol(Tbub),sprintf('   y_{tol} = %5.3f',y_tol(Tbub)));
hold off;
------------------------------------------------------------------
Problem: Bubble Point Calculation
------------------------------------------------------------------
z_bnz =    0.5 [mole fraction]
z_tol =    0.5 [mole fraction]
------------------------------------------------------------------
Bubble Point Temperature =   92.1 [deg C]
Mole Fraction Benzene =  0.714 [mole fraction]
Mole Fraction Toluene =  0.286 [mole fraction]

Application: Dew Point Calculation

Given a vapor mixture at a fixed pressure, the dew point is the temperature at which the first drop of condensate appears. As for the bubble point computation, the key relationship is

$$y_i P = x_i P^{Sat}_i(T)$$

In this we solve for $x_i$

$$x_i = y_i\frac{P}{P^{Sat}_i(T)}$$

The we find the temperature for which the liquid phase mole fractions sum to one,

$$ \sum_i y_i\frac{P}{P^{Sat}_i(T)} =  1$$

% Use z to denote feed composition.

P = 760;
z_bnz = 0.5;
z_tol = 0.5;

% Solve for liquid phase mole fractions as a function of temperature

x_bnz = @(T)z_bnz*P./bnz.Psat(T);
x_tol = @(T)z_tol*P./tol.Psat(T);

% The bubble point is the temperature at which the vapor phase mole
% fractions sum to one.

Tdew = fzero(@(T)x_bnz(T) + x_tol(T) - 1, 100);

problem('Dew Point Calculation')
answer('z_bnz',z_bnz,'mole fraction');
answer('z_tol',z_tol,'mole fraction');
disp(sep);
answer('Dew Point Temperature',Tdew,'deg C');
answer('Mole Fraction Benzene',x_bnz(Tdew),'mole fraction');
answer('Mole Fraction Benzene',x_tol(Tdew),'mole fraction');
disp(sep);

% A graphical interpretation of the solution process

T = 70:110;
plot(T,x_bnz(T)+x_tol(T),T,x_bnz(T),T,x_tol(T));
xlabel('Temperature [deg C]');
ylabel('Mole Fraction');
title('Benzene and Toluene liquid phase mole fractions');
legend('x_{bnz}+x_{tol}','x_{bnz}','x_{tol}','Location','NE');
ax = axis;
hold on;
plot([Tdew Tdew ax(1)],[ax(3) 1 1],'--');
plot(Tdew,1,'.','Markersize',15);
text(Tdew,1.05,sprintf('   T_{dew} = %5.2f deg C',Tdew), ...
    'HorizontalAlignment','center');
plot(Tdew,x_bnz(Tdew),'g.','Markersize',15);
plot([Tdew,ax(1)],x_bnz(Tdew)*[1 1],'g--');
text(Tdew,x_bnz(Tdew),sprintf('   x_{bnz} = %5.3f',x_bnz(Tdew)));
plot(Tdew,x_tol(Tdew),'r.','Markersize',15);
plot([Tdew,ax(1)],x_tol(Tdew)*[1 1],'r--');
text(Tdew,x_tol(Tdew),sprintf('   x_{tol} = %5.3f',x_tol(Tdew)));
hold off;
------------------------------------------------------------------
Problem: Dew Point Calculation
------------------------------------------------------------------
z_bnz =    0.5 [mole fraction]
z_tol =    0.5 [mole fraction]
------------------------------------------------------------------
Dew Point Temperature =   98.8 [deg C]
Mole Fraction Benzene =   0.29 [mole fraction]
Mole Fraction Benzene =   0.71 [mole fraction]
------------------------------------------------------------------

Application: Txy Diagram

Vapor/liquid equilibrium behavior of a binary mixture can be conveniently summarized in a Txy diagram. The Txy diagram depicts the solution to the equations

$$ x_A + x_B = 1$$

$$ x_AP^{Sat}_A(T) +  x_BP^{Sat}_B(T) = P$$

The second of these equations comes from setting the partial pressures of the indivdual components equal to the total pressure. Solving the first equation for $x_B = 1 - x_A$ then substituting into the second equation

$$ x_A = \frac{P-P^{Sat}_B(T)}{P^{Sat}_A(T)-P^{Sat}_B(T)} $$

$$ y_A = x_A \frac{P^{Sat}_A(T)}{P} $$

% Fix total pressure

P = 760;

% Find the boiling points of the pure components. The Antoine class has a
% method Tsat for doing this calculation.

Tbnz = bnz.Tsat(P);
Ttol = tol.Tsat(P);

% A vector of temperatures at which to compute vapor and liquid phase
% compositions

T = Tbnz:((Ttol-Tbnz)/100):Ttol;

% Compute the liquid and vapor phase mole fractions for benzene

x = (P-tol.Psat(T))./(bnz.Psat(T)-tol.Psat(T));
y = x.*bnz.Psat(T)/P;

% Plot the results

plot(x,T,y,T);
title('Txy diagram: benzene/toluene');
xlabel('Mole fraction benzene');
ylabel('Temperature [deg C]');
legend('Dew Point','Bubble Point');

% Show where the previous bubble and dew point calculations fit on the Txy
% diagram.

hold on;
plot(x_bnz(Tdew),Tdew,'.','Markersize',15);
plot(y_bnz(Tbub),Tbub,'g.','Markersize',15);
plot([0.5 0.5],[Tbnz Ttol],'k--');
plot([0.5 x_bnz(Tdew)],[Tdew Tdew],'b--');
plot([0.5 y_bnz(Tbub)],[Tbub Tbub],'g--');
hold off;

Application: xy Diagram

The xy diagram provides much of the same information as the Txy diagram in format where the vapor phase composition is plotted as a function of the liquid phase composition. Temperature becomes a implicit parameter.

% Reuse data for the Txy diagram to construct the xy diagram

plot(x,y)
xlabel('x_{bnz}');
ylabel('y_{bnz}');
axis('equal')
axis('square');

% Annotate xy diagram with temperature

hold on;
for k = 1:10:length(T)
    plot(x(k),y(k),'.','Markersize',5);
    text(x(k),y(k),sprintf('   %4.1f',T(k)));
end
hold off

% Show where previous bubble point and dew point calculations fit on this
% plot

hold on
plot([0.5 0.5 0],[0 y_bnz(Tbub) y_bnz(Tbub)],'g--');
plot([0 x_bnz(Tdew) x_bnz(Tdew)],[0.5 0.5 0],'g--');
hold off

Application: Isothermal Flash

The isothermal flash is a calculation coupling material balances with vapor/liquid equilibrium. Consider a flash unit with a molar feedrate $n_{Feed}$ and a vapor exit stream with flowrate $n_{Vap}$ and a liquid exit stream with molar flow $n_{Liq}$. At steady state

$$n_{Feed} = n_{Liq} + n_{Vap} $$

with component material balances

$$ z_i n_{Feed} = x_i n_{Liq} + y_i n_{Vap} $$

where $z_i$ denotes the mole fraction of component i in the feed stream. From Raoult's law

$$ y_i P = x_i P^{Sat}_i(T) $$

or

$$ y_i = K_i x_i $$

where $K_i = P^{Sat}_i(T)/P$ are the 'K-values'. Given the temperature, pressure, molar feedrate and feed compositions, the task is to compute the flowrates and composition of the exit streams. Traditionally, the first step in this calculation is to divide the first two equations by $n_{Feed}$, and to define a vapor fraction phi. The second equation becomes

$$ z_i  = (1-\phi) x_i   + \phi y_i $$

Using Raoult's law to eliminate $y_i$, solving each component balance for $x_i$ gives

$$ x_i = \frac{z_i}{1 + \phi(K_i-1)} $$

The liquid phase mole fractions must sum to one. Therefore, for a fixed temperature $T$, we need to find $\phi$ such that

$$ \sum+i \frac{z_i}{1 + \phi(K_i-1)} = 1 $$

Alternatively, one can compute the vapor phase mole fractions $y_i = K_ix_i$ which must also sum to one

$$ \sum_i \frac{K_iz_i}{1 + \phi(K_i-1)} = 1$$

Which of these equations to use to solve for $\phi$? While either can be made to work, it turns out that the difference between these equations has better numerical properties. The resulting equation is called the Rachford-Rice equation

$$ \sum_i \frac{(K_i-1)z_i}{1 + \phi(K_i-1)} = 0 $$

This is the equation we solve to find $\phi$.

% Feed and operating conditions
z_bnz = 0.5;
z_tol = 0.5;

% Operating Conditions
P = 760;
T = 95;

% K-values
K_bnz = bnz.Psat(T)/P;
K_tol = tol.Psat(T)/P;

% Rachford-Rice equation
r = @(phi) z_bnz*(K_bnz-1)./(1+phi*(K_bnz-1)) + ...
           z_tol*(K_tol-1)./(1+phi*(K_tol-1));

% Solve for vapor phase fraction
phi = fzero(r,[0,1]);

% Show solution
problem('Isothermal Flash')
answer('Flash Temperature',T,'deg C');
answer('Flash Pressure',P,'mmHg');
answer('z_bnz',z_bnz,'mole fraction');
answer('z_tol',z_tol,'mole fraction');
disp(sep);
answer('Vapor Phase Fraction',phi,'fraction');
answer('y_bnz',K_bnz*z_bnz/(1+phi*(K_bnz-1)),'mole fraction');
answer('y_tol',K_tol*z_tol/(1+phi*(K_tol-1)),'mole fraction');
disp(sep);
answer('Liquid Phase Fraction',1-phi,'fraction');
answer('x_bnz',z_bnz/(1+phi*(K_bnz-1)),'mole fraction');
answer('x_tol',z_tol/(1+phi*(K_tol-1)),'mole fraction');
disp(sep);
------------------------------------------------------------------
Problem: Isothermal Flash
------------------------------------------------------------------
Flash Temperature =     95 [deg C]
Flash Pressure =    760 [mmHg]
z_bnz =    0.5 [mole fraction]
z_tol =    0.5 [mole fraction]
------------------------------------------------------------------
Vapor Phase Fraction =  0.436 [fraction]
y_bnz =  0.625 [mole fraction]
y_tol =  0.375 [mole fraction]
------------------------------------------------------------------
Liquid Phase Fraction =  0.564 [fraction]
x_bnz =  0.403 [mole fraction]
x_tol =  0.597 [mole fraction]
------------------------------------------------------------------