CENG 301 Material and Energy Balances

Chapter 4: Multi-Phase Problems

On this page:

Other sections in Chapter 4:

4.1: Phase Equilibria Between Vapor and Liquids

Problems involving a vapor and its liquid in equilibrium under ideal conditions may be analyzed with vapor pressures determined from the Antoine equation. If a system consists of two phases in equilibrium with the composition of one phase known, then according to the phase rule (p. 400 in Reklaitis) it has a one degree of freedom. Thus knowing the temperature of a vapor-liquid mixture allows us to find the vapor pressure of the system. Many pure compounds have vapor pressures that can be fit quite well with the Antoine equation:

	ln(pi*) = Ai - (Bi/(T+Ci))
	pi* 	is the vapor pressure of compound i,
	T	is the temperature
	Ai, Bi, Ci	are the Antoine coefficients found in the appendices
      in the text or at the end of this guide or in our computer 
      data files.


The partial pressure of compound i over an ideal liquid with mol fraction of that compound xi is:

pi = xipi*

If the vapor is also ideal and the vapor's total pressure is P and the mol fraction of i in the vapor is yi then:

pi = yiP

4.1.1 Vapor Pressure of Pure Compounds

The Vapor Pressure Function: vappr

The function vappr calculates the vapor pressure from the Antoine coefficients for compounds with data in the CENG database. Since the data in the CENG database gives pressures in kPa these are the units given by vappr. It uses the global variable Aabc that is set when the function setname is executed. The vapor pressure function vappr lists as:

function p=vappr(t,j)
% vappr: calculates vapor pressure using Antoine data
% function p=vappr(t,j)
% Argument List:
% Argument   Gives
%    t      temperature in Tdeg.
%    j      indices of compounds in cnms that will
%            have their vapor pressures determined.
%           If this argument is omitted, all vapor pressures
%             will be found.
% Aabc must be set for the compounds with setnme.
% Ex: start301 to specify compounds
%   Be sure to specify you do want to use property data
%     Tdeg='C';
%     vappr(300)
%   will give the vapor pressures of all the compounds at 300C.
% The vapor pressure is given in kPa for the Antoine Data in sdata.
global Aabc
if nargin==1

There is also a function called vapprvt that find the vapor pressures of the compounds at a vector of temperatures. It lists as:

function p=vapprvt(t,j)
% vapprvt: calculates Antoine vapor pressure for a vector of temperatures
% function p=vapprvt(t,j)
% Argument List:
% Argument   Gives
%    t      vector of temperatures in Tdeg.
%    j      indices of compounds in cnms that will
%            have their vapor pressures determined.
%           If this argument is omitted, all vapor pressures
%             will be found.
% Aabc must be set for the compounds with setnme.
% Ex: start301 with mass and energy data loaded
%   Specify the compounds you want.
%     Tdeg='C';
%     vapprvt([100 300])
%   will give the vapor pressures of all the compounds at 100 & 300C.
% The vapor pressure is given in kPa for the Antoine Data in the CENG database.
global Aabc
if nargin==1
   for k=1:length(t)
   for k=1:length(t)


We can find the vapor pressures of water and/or benzene after start301.

Here are your compounds' names:
>> vappr(300)
ans =
    3.5482   13.6790   <-- vapor pressure in kPa
>> vappr(300,2) 
ans =
   13.6790            <-- Getting the vapor pressure of C6H6


4.1.2 Bubble and Dew Points

The bubble point temperature is defined as the lowest temperature that produces a vapor bubble when a liquid mixture is slowly heated at constant pressure. In effect, we then know at the bubble point the total pressure and composition of a liquid in equilibrium with a vapor of unknown composition. The vapor composition and temperature must then be found. In an ideal system with all vapor pressures given by the Antoine equation the mol fractions in the vapor and liquid are related by:

yiP = xipi*(T)

The total pressure is P and the temperature is T. Since the sum of the mol fractions yi must be 1, we have a nonlinear equation to solve for the temperature T:

Once the bubble point temperature: TBP, is found, the original (unsummed) equation gives (with TBP used to find the vapor pressures) the vapor composition.

yi = xipi*(TBP)/P

We similarly define the dew point temperature as the temperature at which the first liquid drop would form when the temperature of a mixture of vapors is slowly decreased at a specified constant pressure. This time we know the vapor composition yi and the total pressure P. A nonlinear equation to find the dew point temperature can be found by solving the relation between yi and xi for xi and summing the result over all the compounds.

The Distribution Coefficient Function: kval

In both the bubble and dew point calculations it is conventional to define the ratio yi/xi as the distribution coefficient Ki. In our ideal system this coefficient is given by:

Ki(T,P) = yi/xi = pi*(T)/P

Note that the coefficient depends on the temperature T, pressure P and the compound index i. To find the dew point we then need to solve for T such that:

Once the dew point is found the liquid mol fractions may be calculated with:

xi = yi / Ki(TDP,P)

The program kval has as its first argument the temperature in the units of Tdeg and as its second argument the pressure in kPa. Kval returns the Ki for any set of compounds as shown below:

function Ks=kval(t,p)
% kval - Finds the distribution coefficient for a set of compounds
% function Ks=kval(t,p)
% Argument List:
% Argument    Gives
%    t      temperature in the units given by Tdeg.
%             (must be a scalar.)
%    p      pressure in kPa (must be a scalar).
% The vector of ideal Ks=(y1/x1, y2/x2,...) is returned.
% Ex: start301
%     Give the file name with property data for benzene and toluene
%    Tdeg='C'; 
%    kval(25,10)


Here is a session using kval to find the coefficients for benzene and toluene at 25°C and 10 kPa.

Here are your compounds' names:
>> Tdeg='C'; 
>> kval(25,10)
ans =
    1.2574    0.3791

Note that the temperature must be given in the units designated by Tdeg. Kval may be used to find the bubble point of a mixture with the composition (on a molar basis): 40% benzene, 40% toluene and 20%o-xylene at 101 kPa by:

Here are your compounds' names:
>> Tdeg='C';
>> ks=kval(25,101) 
ks =
    0.1245    0.0375    0.0088
>> xs=[.4 .4 .2]; 
>> sum(xs.*ks)		<-- 25°C is too low
ans =
    0.0666		<-- The yi's sum to .0666
>> sum(xs.*kval(100,101)) 
ans =
    1.0536		<-- 100°C is too high
>> ssec2([25 100],'1-sum([.4 .4 .2].*kval(x,101))') 
ans =
   98.1826		<-- ssec2 finds TBP as 98.18°
>> x=ans; 
>> [.4 .4 .2].*kval(x,101)                                    
ans =
    0.6733    0.2775    0.0492 		<-- Here are the yi
>> sum(ans) 
ans =
    1.0000		<-- They sum to 1


Bubble Points with: bubpt

We will now combine the steps associated with the bubble point determination in a single function. In order to get an initial range to look in for the bubble point the program uses the function teql. Teql inverts the Antoine Equation to find the temperatures for a given pressure of each compound. It gives its results in K, but these may be converted to the units of Tdeg by the function atk:

>> help atk
  atk   - Converts the given temperature to Kelvin
  function t=atk(tk)
  Argument List: tk is the temp. in K.
  The function returns the temp in Tdeg.
  See also AT
  Ex. Tdeg='C';
     will give 25.
 revised help by TYLC

Atk is the inverse of at as may be seen in:

>> Tdeg='C'; 
>> atk(300)
ans =
>> Tdeg='K'; 
>> atk(300)
ans =

Here are the comments in teql together with a demonstration of it for the same pressures found previously with vappr.

>> help teql
  teql: find equilibrium temperature for given pressure
  function t=teql(p,j)
  Argument List:
  Argument   Gives
     p      pressure in kPa
     j      indices of compounds in cnms that you want
            the temperature for.  If you omit this, all
            the compounds will have equil. temps. found.
   t is given in Tdeg units and is the equilibrium temperature at
 	the pressure specified
  Ex: start301
      Give the file name with property data for water and benzene

Here are your compounds' names:
>> teql(13.679) 
ans =


The maximum and minimum of the temperatures found by teql should always straddle the bubble point for an ideal mixture. Finally here is our bubpt function.

The procedure followed in the code is almost identical to the terminal calculation demonstrated earlier.

  1. The range of possible temperatures is found with teql.
  2. The data to be used is put into the vector cs.
  3. The ssec2 function is called to find the bubble point.
  4. The vapor compositions are found with kval.

Here is a demonstration that may be compared with the terminal calculations shown previously.

Here are your compounds' names:
>> [t,ys]=bubpt(101,[.4 .4 .2])
t =
ys =
    0.6733    0.2775    0.0492 
>> Tdeg 
Tdeg =
>> kval(t,101) 
ans =
    1.6833    0.6936    0.2460


The temperature returned by bubpt is in the units of Tdeg. In this case, the bubble point for our mixture is 371.332K and the vapor in equilibrium with the liquid has the mol fraction as shown in the following table:

      Compound      xi         yi       Ki
      C6H6         0.4       0.6733   1.6833
      Toluene      0.4       0.2775   0.6936
      o-xylene     0.2       0.0492   0.2460


The Dew Point function: dewpt

A function that finds the dew point in an analogous manner to bubpt is called dewpt. Here are the comments about it:

>> help dewpt
   dewpt		- Determines the dew point of a mixture
  function [t,xs]=dewpt(p,ys,erc)
  Argument List:
  Argument    Gives
     p       pressure in kPa
     ys      vector of vapor mol fractions.
    erc      error tolerance
  Returns     Giving
     t       Dew point temp. in Tdeg.
     xs      vector of liquid mol fractions in equil.
             with the vapor at the dew point.
  Ex: start301   to specify our compounds 
      [t,xs]=dewpt(101,[.4 .4 .2])
     (will give the dew point at 101kPa for the vapor mol fractions:
      Compound   vapor mol frac.
      first          .4
      second         .4
      third          .2


The trial and error solutions performed in bubpt and dewpt normally converge quite rapidly as long as you have a reasonable range of temperatures to search for. The function fbpdp illustrates the behavior of the functions that have zeros when we find the bubble point and dew point. Here are the comments in the function:

>> help fbpdp
  fbpdp - Shows behavior of functions in dew,bubble point calculations
  function [f1,f2]=fbpdp(zs,p,ts)
  Argument List:
  Argument    Gives
     zs      Liquid or vapor mol fractions.
     p       pressure in kPa
     ts      vector of temperature in Tdeg.
  Two functions used to solve for the bubble and dew points 

     of a mixture with composition zs.
  Vectors of fk(t) are returned where:
     f1 is 1-sum(zs*kval(t,p))
     f2 is 1-sum(zs/kval(t,p))
  Example: start301   to specify your compounds
           zs=[.4 .4 .2];


The following session shows the result of the suggested example:

 Here are your compound names:
 Compound #    Name
     1     benzene
     2     toluene
     3     o-xylene
Enter the number of streams: 0
>> teql(101)
ans =
  353.5035   <-- Thus the interval 350 to 420 spans 
  383.6375   <-- the possible range for the bubble 
  417.4291   <-- and dew points of the mixture. 
>> ts=350:2:420; 
>> zs=[.4 .4 .2]; 
>> [f1,f2]=fbpdp(zs,101,ts); 
>> plot(ts,[f1;f2]) 
>> grid
>> xlabel('Temperature (K)')
>> ylabel('Dew & Bubble Functions')
>> title('Functions used to Find TDP and TBP')
>> gtext('1-sum(xs*ks)')
>> gtext('1-sum(ys/ks)')


Here is the plot generated:

4.1.3 Flashing a Liquid Mixture

We may also analyze the behavior of a flash unit in which the exit liquid and vapor streams are assumed to be in equilibrium. The unit has one stream going into the unit and two streams leaving it. The unit makes use of differences in vapor pressures of the compounds in a mixture at a given temperature to achieve a separation. The isothermal flash equation is obtained from material balances around the unit:

     F = V + L
     Fzi = Vyi + Lxi
     where F is the feed rate in molar units
           V is the vapor rate
           L is the liquid rate
           zi, yi and xi are mol fractions of i in the feed,
           vapor and liquid streams respectively.


Assuming the liquid and vapor to be in equilibrium:

yi = Ki(T,P)xi

we may solve for the liquid compositions:

xi = Fzi / (VKi+L)

The vapor mol fractions and the liquid mol fractions must both sum to 1, so:

Thus, for a feed stream with mol fractions zi the equation that determines the vapor to liquid ratio is:

Ki for each compound can be obtained from the Antoine equation (as in kval ) at the specified temperature and pressure. After the ratio V/F is found, the composition of the product streams is found by:

        L/F = 1 - (V/F)
        xi = zi / ((V/F)Ki + (L/F))
        yi = Kixi


The function tryf was written to compare three functions that one might use to solve for V/F in the flash calculations. Its comments read:

>> help tryf 
function [f1,f2,f3]=tryf(zs,ks,vofs)
Three functions that could be used to solve the flash
Argument list
Argument   Gives
   zs    vector of feed mol fractions.
   ks    vector of distribution coefficients.
   vofs  values of V/F to determine each sum at.
Vectors of fk(V/F) are returned where:
   f1 is 1-sum(xs)
   f2 is 1-sum(ys)
   f3 is sum(ys-xs)
Example: start301: specify the names: benzene, toluene, o-xylene
         zs=[.4 .4 .2];


Following the example shown in these comments, we find:

 Here are your compound names:
 Compound #    Name
     1     benzene
     2     toluene
     3     o-xylene
Enter the number of streams: 0
>> vofs=0:.05:1; 
>> ks=kval(375,101) 
ks =
    1.8613    0.7754    0.2791
>> zs=[.4 .4 .2]; 
>> [f1,f2,f3]=tryf(zs,ks,vofs); 
>> plot(vofs,[f1;f2;f3])
>> xlabel('V/F') 
>> ylabel('fk') 
>> title('Three Functions for Flash Calculations') 
>> grid  
>> text(.7,-.11,'1-sum(xs)') 
>> text(.6,.08,'1-sum(ys)') 
>> text(.4,-.15,'sum(ys-xs)')
>> print


Here is the plot we get:

The following session shows how to find the vapor fraction V/F from a flash unit used to separate a mixture of benzene, toluene and o-xylene. The feed to the unit is given as: 30% benzene, 30% toluene and 40% o-xylene. The unit is operated at 112°C and 101kPa.

 Here are your compound names:
 Compound #    Name
     1     benzene
     2     toluene
     3     o-xylene
Enter the number of streams: 0
>> Tdeg='C';
>> ks=kval(112,101) 
ks =
    2.4284    1.0436    0.3903
>> fx='sum([.3 .3 .4].*(1-c)./(1+x.*(c-[1 1 1])))'; 
>> x=0; 
>> c=ks; 
>> eval(fx) 
ans =
   -0.1977                        <-- A negative result for x=0.
>> x=1; 
>> eval(fx)
ans =
    0.4358                        <-- A positive value for x=1.
>> ssec2([0 1],fx,ks)             <-- Using ssec2 to find
ans =
    0.3278                        <-- V/F = 0.3278
>> xs=[.3 .3 .4]./(1+ans*(ks-1)) 
xs =
    0.2043    0.2958    0.4999    <-- The liquid composition
>> ks.*xs                     
ans =
    0.4962    0.3087    0.1951    <-- The vapor composition


The following table summarizes the results of this session:

    T=112°C, P=101kPa, V/F=0.3278
                      Mol Fractions
    Compound   Feed       Liquid       Vapor       Ki
     Benzene    0.3       0.2043      0.4962     2.4284
     Toluene    0.3       0.2958      0.3087     1.0436
     o-Xylene   0.4       0.4999      0.1951     0.3903


The Flash function: flash

The flash function does all the operations shown in the last demonstration automatically. The solution of problem 7.8 below demonstrates its use:

Problem 7.8) One atmosphere is:


  1.000000    atm   =  101.3250    kPa


 Here are your compound names:
 Compound #    Name
     1     toluene
     2     p-xylene
     3     ethyl benzene
enter the number of streams: 0         


The bubble point was found with the bubpt function, and dew point with dewpt:

>> Tdeg='C'; 
>> bubpt(101.325,[1 1 1]/3)
ans =
>> dewpt(101.325,[1 1 1]/3)
ans =


Thus the bubble point is 126.14°C and the dew point is 130.39°C. Since 400K (126.85°C) is between the bubble point and the dew point, we would expect a mixture of vapor and liquid at that temperature. The program flash was used to find the vapor fraction and the composition of the liquid and vapor by the same procedure shown above. Here are the comments that tell how to use it:

>> help flash
  flash		- Performs a flash calculation with a table of results
  function flash(t,p,zs,erc,nmax)
  Argument List
  Argument   Gives
     t      temperature in Tdeg.
     p      pressure in kPa.
     zs     vector of feed mol fractions.
     erc    gives ssec2 the difference between successive values of x
            allowed when convergence is deemed to be achieved.
     nmax   is the maximum number of modified retrys by fzero before 
   If erc and nmax are not given, then the program uses:
      erc=.0001 and nmax=10
  Example:               to set up the property data file for
                         benzene, toluene, o-xylene
          >> start301  to specify the compounds 
          >> flash(375,101,[.4 .4 .2])
   will produce a table giving the flash V/F ratio
   as well as the liquid and vapor compositions.


Here is flash in use:

>> flash(126.85,101.325,[1 1 1]/3)
Temperature  Pressure    V/F
  126.8500  101.3250    0.1403
Compound        Feed      Liquid    Vapor     K=y/x
toluene         0.33333   0.30914   0.48157   1.55779
p-xylene        0.33333   0.34658   0.25216   0.72757
ethyl benzene   0.33333   0.34428   0.26627   0.77340
>> flash(125,101.325,[1 1 1]/3)
There is no V/F for this flash.
>> flash(130.5,101.325,[1 1 1]/3)
There is no V/F for this flash.

Note that temperatures below the bubble point or above the dew point result in a message that no value of V/F can be found.

[Go to next section: 4.2]
[ Go to CENG301 Homepage | CENG301 Notes Table of Contents]

Last modified July 28, 1998.