# Chapter 4: Multi-Phase Problems

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
att=at(t);
if nargin==1
p=(exp(Aabc(:,1)-Aabc(:,2)./(att+Aabc(:,3))))';
else
p=(exp(Aabc(j,1)-Aabc(j,2)./(att+Aabc(j,3))))';
end```

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)
p=[p;(exp(Aabc(:,1)-Aabc(:,2)./(at(t(k))+Aabc(:,3))))'];
end
else
for k=1:length(t)
p=[p;(exp(Aabc(j,1)-Aabc(j,2)./(at(t(k))+Aabc(j,3))))'];
end
end

```

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

```Here are your compounds' names:
water
benzene
>> 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)
Ks=vappr(t)/p;

```

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:
benzene
toluene
>> 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:
benzene
toluene
o-xylene
>> 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.
Ex. Tdeg='C';
atk(298.15)
will give 25.
revised help by TYLC
```

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

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

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
teql(13.684)

water
benzene
>> teql(13.679)
ans =
325.3016
300.0000

```

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:
benzene
toluene
o-xylene
>> [t,ys]=bubpt(101,[.4 .4 .2])
t =
371.3326
ys =
0.6733    0.2775    0.0492
>> Tdeg
Tdeg =
K
>> 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
ts=350:2:420;
zs=[.4 .4 .2];
[f1,f2]=fbpdp(zs,101,ts);

```

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
problem.
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
vofs=0:.05:1;
ks=kval(375,101)
zs=[.4 .4 .2];
[f1,f2,f3]=tryf(zs,ks,vofs);

```

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

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 =
126.1392
>> dewpt(101.325,[1 1 1]/3)
ans =
130.3916

```

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
quitting.
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]