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

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,
• Raoult's law - partial pressure of a component is equal to the liquid phase mole fraction times the saturation pressure,

The last two assumptions lead to a relationship

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

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

In this we solve for

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

% 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

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 then substituting into the second equation

% 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 and a vapor exit stream with flowrate and a liquid exit stream with molar flow . At steady state

with component material balances

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

or

where 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 , and to define a vapor fraction phi. The second equation becomes

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

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

Alternatively, one can compute the vapor phase mole fractions which must also sum to one

Which of these equations to use to solve for ? 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

This is the equation we solve to find .

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