Chapter 10 : Vapor Liquid Equillibrium Introduction

Example 10.1 page no : 148

In [9]:
%matplotlib inline

from matplotlib.pyplot import plot,subplot,suptitle,xlabel,ylabel
from numpy import array,linspace,exp
import math 

#Antoinie Equations
#ln P1_sat = 14.2724-(2945.47/(T-49.15))  [KPa]
#ln P2_sat = 14.2043-(2972.64/(T-64.15))  [KPa]
#(a) Graph Showing P vs x1 and P vs y1 for T = 348.15K
T = 348.15;			#[K]
#using BUBL P calculations

#Calculation of P1_sat and P2_sat at T = 348.15K
P1_sat = round(math.exp(14.2724-(2945.47/(T-49.15))),2)  			#KPa
P2_sat = round(math.exp(14.2043-(2972.64/(T-64.15))),2)  			#KPa

#Using Eqn P = P2_sat+(P1_sat-P2_sat)x1

x = linspace(0,1,6)
P = round(P2_sat+((P1_sat-P2_sat)*x),2);
y = round(x*P1_sat/P,4);

print ('Explanations Of graph')

Ans = [x,y,P];
print '   x1      y1       P/PKa',Ans

y1 = 0.6;
y2 = 1-y1;
P_dew = round(1/((y1/P1_sat)+(y2/P2_sat)),2)
x1 = round(y1*P_dew/P1_sat,4)

# Plotting the graph
T = 348.15;			#[K]
P1_sat = round(exp(14.2724-(2945.47/(T-49.15))),2)  			#KPa
P2_sat = round(exp(14.2043-(2972.64/(T-64.15))),2)  			#KPa

x = linspace(0,1,6)
P = round(P2_sat+((P1_sat-P2_sat)*x),2);
y = round(x*P1_sat/P,4);

plot(x,P,'g-')  			#P vs x1
plot(y,P,'b-')  			#P vs y1
x = [0,0.1];
P = [P2_sat,P2_sat];
plot(x,P,'--')  			#P2_sat
x = [0.9,1];
P = [P1_sat,P1_sat];
plot(x,P,'r--') 			#P1_sat

x1 = 0.6;
P_b = round(P2_sat+((P1_sat-P2_sat)*x1),2);
y1 = round(x1*P1_sat/P_b,4);
x = [x1,y1];
P = [P_b,P_b];
plot(x,P,'bo-')  			#b--b'

y1 = 0.6;
y2 = 1-y1;
P_c = round(1/((y1/P1_sat)+(y2/P2_sat)),2)
x1 = round(y1*P_c/P1_sat,4)

x = [x1,y1];
P = [P_c,P_c];
plot(x,P,'ro-')  			#c'--c

P = [(P_b+10),P_b,P_c,(P_c-10)];
x = [0.6,0.6,0.6,0.6];
plot(x,P,'go-')  			#a--b--c--d--0.6

P = [(P_c-10),30];
x = [0.6,0.6];
plot(x,P,'yo--')

P = [110,80];
x = [0.6,0.6];
plot(x,P,'w')
suptitle('(a)T/t = 348.15K')
xlabel('x1,y1')
ylabel('P/Kpa')

print ("This is the liquid-phase composition at point c''")

#(b) Graph showing (t vs x1) and (t vs y1) for a pressure of 70KPa
#Example 10.2(b)
P = 70;			#[KPa]

T1_sat = round(2945.47/(14.2724-math.log(P))+49.15,2);
T2_sat = round(2972.64/(14.2043-math.log(P))+64.15,2);

T = array([T1_sat,347.15,351.15,355.15,359.15,T2_sat]);

P1_sat = round(exp(14.2724-(2945.47/(T-49.15))),2);  			#KPa
P2_sat = round(exp(14.2043-(2972.64/(T-64.15))),2);  			#KPa

x = round((P-P2_sat)/(P1_sat-P2_sat),3);
y = round((x*P1_sat)/P,3);

Ans = [x,y,T];
print '   x1       y1       T/t(K/C`)',Ans

#at x1 = 0.6;
x1_b = 0.6;
x2_b = 1-x1_b;

T_a = 347.15;			# Intermediate Temperature  (Point a in graph)
P1_sat_a = round(exp(14.2724-(2945.47/(T_a-49.15))),2);  			#KPa
P2_sat_a = round(exp(14.2043-(2972.64/(T_a-64.15))),2);  			#KPa
alpha = P1_sat_a/P2_sat_a;  			#Initial
a = T_a;
i = -1;
while(i == -1):
    P2_sat_b = P/((x1_b*alpha)+x2_b);
    b = round(2972.64/(14.2043-math.log(P2_sat_b))+64.15,2);
    dT = abs(a-b);
    if(dT == 0):
        i = 0;
        T_b = b;
    alpha = math.exp(0.0681-(2945.47/(b-49.15))+(2972.64/(b-64.15)));  			#Eqn C
    a = b;

P1_sat_b = round(math.exp(14.2724-(2945.47/(T_b-49.15))),2);  			#KPa
y1_b = round((x1_b*P1_sat_b)/P,4);  			#b`

print 'Hence by iteration Temp(Temp at b) at x1 = 0.6 is ',T_b,'K'
print 'Hence by iteration P1_sat at x1 = 0.6 is ',P1_sat_b,'KPa'
print 'Composition of Vapor(b`) at x1 = 0.6',y1_b


#At y1 = 0.6
y1_c = 0.6;
y2_c = 1-y1_c;
T_d = 355.15;			# Intermediate Temperature  (Point a in graph)
P1_sat_d = round(math.exp(14.2724-(2945.47/(T_d-49.15))),2);  			#KPa
P2_sat_d = round(math.exp(14.2043-(2972.64/(T_d-64.15))),2);  			#KPa
alpha = P1_sat_d/P2_sat_d;  			#Initial
d = T_d;
i = -1;
while(i == -1):
    P1_sat_c = P*(y1_c+(y2_c*alpha));
    c = round(2945.47/(14.2724-math.log(P1_sat_c))+49.15,2);
    dT = abs(d-c);
    if(dT == 0):
        i = 0;
        T_c = c;
    alpha = math.exp(0.0681-(2945.47/(c-49.15))+(2972.64/(c-64.15)));  			#Eqn C
    d = c;

P1_sat_c = round(math.exp(14.2724-(2945.47/(T_c-49.15))),2);  			#KPa
x1_c = round((y1_c*P)/P1_sat_c,4);  			#c`

print 'Hence by iteration Temp(Temp at b) at y1 = 0.6 is ',T_c,'K'
print 'Hence by iteration P1_sat at y1 = 0.6 is ',P1_sat_c,'KPa'
print 'Composition of liqiud(c`) at y1 = 0.6',x1_c

#Graph
T = linspace(T1_sat,T2_sat,10); 

P1_sat = round(exp(14.2724-(2945.47/(T-49.15))),2);  			#KPa
P2_sat = round(exp(14.2043-(2972.64/(T-64.15))),2);  			#KPa

x = round((P-P2_sat)/(P1_sat-P2_sat),3);
y = round((x*P1_sat)/P,3);

plot(x,T,'g-');    
plot(y,T,'b-');  

xsat = [0,0.1];
T2sat = [T2_sat,T2_sat];
plot(xsat,T2sat,'--')  			#T2_sat

xsat = [0.9,1];
T1sat = [T1_sat,T1_sat];
plot(xsat,T1sat,'r--') 			#T1_sat

Tcc = [T_c,T_c];
xc = [x1_c,y1_c];
plot(xc,Tcc,'ro-')  			#c--c'

Tbb = [T_b,T_b];
xb = [x1_b,y1_b];
plot(xb,Tbb,'bo-')  			#b--b'

Tabcd = [T_d,T_c,T_b,T_a];
xabcd = [0.6,0.6,0.6,0.6];
plot(xabcd,Tabcd,'go-')  			#a--b--c--d--0.6

Tao = [T_a,340];
xao = [0.6,0.6];
plot(xao,Tao,'yo--')
suptitle('(b)P = 70KPa')
xlabel('x1,y1')
ylabel('T(K)')
Explanations Of graph
   x1      y1       P/PKa [array([ 0. ,  0.2,  0.4,  0.6,  0.8,  1. ]), array([ 0.    ,  0.3313,  0.5692,  0.7483,  0.888 ,  1.    ]), array([ 41.98,  50.23,  58.47,  66.72,  74.96,  83.21])]
This is the liquid-phase composition at point c''
   x1       y1       T/t(K/C`) [array([ 1.   ,  0.738,  0.516,  0.318,  0.142,  0.   ]), array([ 1.   ,  0.849,  0.676,  0.474,  0.239,  0.   ]), array([ 342.99,  347.15,  351.15,  355.15,  359.15,  362.73])]
Hence by iteration Temp(Temp at b) at x1 = 0.6 is  349.57 K
Hence by iteration P1_sat at x1 = 0.6 is  87.17 KPa
Composition of Vapor(b`) at x1 = 0.6 0.7472
Hence by iteration Temp(Temp at b) at y1 = 0.6 is  352.73 K
Hence by iteration P1_sat at y1 = 0.6 is  96.54 KPa
Composition of liqiud(c`) at y1 = 0.6 0.4351
Out[9]:
<matplotlib.text.Text at 0x108748c50>

Example 10.2 page no : 149

In [8]:
# Variables
H = 990.;			#[Bar] Henry's Law const
T = 283.15;			#[K]
P2_sat = 0.01227;			#[Bar]  from Steam Tables
x1 = 0.01;			#Assumed
x2 = 1-x1;
y1 = 1;

# Calculations
P = round((x1*H)+(x2*P2_sat),3);
x1 = round((y1*P)/H,4);
x2 = 1-x1;
y2 = round((x2*P2_sat)/P,4);
y1 = 1-y2;

# Results
print 'Composition in liquid Phase',x1
print 'Composition in vapor Phase',y1
print 'Pressure Exerted on Can',P,'Bar'
print ('Hence Vapor phase chosen is nearly pure')
Composition in liquid Phase 0.01
Composition in vapor Phase 0.9988
Pressure Exerted on Can 9.912 Bar
Hence Vapor phase chosen is nearly pure

Example 10.3 page no : 150

In [7]:
import math 

#Equations to be Used
# ln v1 = A*(x2**2)  ln v2 = A*(x1**2)  Where A = 2.771-0.00523T

#Antoine Equations
#ln P1_sat  =  16.59158-(3643.31/(T-33.424)) 
#ln P2_sat  =  14.25326-(2665.54/(T-53.424))
#P  =  E(xi * Vi * Pi_sat)  E--Summation        Eqn 10.6
#P  =  1/E(yi / (vi * Pi_sat))  E--Summation    Eqn 10.7

#(a) Calculate P and (yi) , for T = 318.15K and x1 = 0.25

# Variables
T = 318.15;			#[K]  Given
x1 = 0.25; 			#Given
x2 = 1-x1;

# Calculations and Results
P1_sat  =  round(math.exp(16.59158-(3643.31/(T-33.424))),2);			#[KPa]
P2_sat  =  round(math.exp(14.25326-(2665.54/(T-53.424))),2);			#[KPa]
A = round(2.771-(0.00523*T),3);
v1 = round(math.exp(A*(x2**2)),3);
v2 = round(math.exp(A*(x1**2)),3);

#Form Eqn(10.6)
P_a = round((x1*v1*P1_sat)+(x2*v2*P2_sat),2);			#[KPa]
y1_a = round((x1*v1*P1_sat)/P_a,3);
y2_a = round((x2*v2*P2_sat)/P_a,3);

print ('(a)P and [yi] for T = 318.15K and x1 = 0.25')
print ('BUBL P calculations')
print 'P  =  ',P_a,'KPa'
print 'y1  =  ',y1_a
print 'y2  =  ',y2_a

#(b) Calculate P and (xi) , for T = 318.15K and y1 = 0.60

#DEW P calculation
y1 = 0.6;
y2 = 1-y1;
T = 318.15;			#[K]
P1_sat  =  round(math.exp(16.59158-(3643.31/(T-33.424))),2);			#[KPa]
P2_sat  =  round(math.exp(14.25326-(2665.54/(T-53.424))),2);			#[KPa]
A = round(2.771-(0.00523*T),3);
v1 = 0.1;			#Assumed
v2 = 0.1;			#Assumed
a1 = v1;
a2 = v2;
i = -1;
while(i == -1):
    P = round(1/((y1/(a1*P1_sat))+(y2/(a2*P2_sat))),2); 
    x1 = round(y1*P/(a1*P1_sat),4);
    x2 = 1-x1;
    b1 = round(math.exp(A*(x2**2)),4);
    b2 = round(math.exp(A*(x1**2)),4);
    dt = abs(b1-a1);
    if(dt == 0):
        i = 0;
        v1 = b1;
        v2 = b2;
        break;
    a1 = b1;
    a2 = b2;

x1_b = x1;
x2_b = 1-x1_b;
P_b = P;
v1_b = v1;
v2_b = v2;
print ('(b)P and [xi] for T = 318.15K and y1 = 0.60')
print ('DEW P calculations')
print 'P  =  ',P_b,'kPa'
print 'x1  =  ',x1_b
print 'x2  =  ',x2_b

#(c) Calculate T and (yi) for P  =  101.33 KPa and x1 = 0.85

#BUBL T calculation
P = 101.33;
x1 = 0.85;
x2 = 1-x1;
T1_sat = round((3643.31/(16.59158-math.log(P)))+33.424,2);
T2_sat = round((2665.54/(14.25326-math.log(P)))+53.424,2);
T = (x1*T1_sat)+(x2*T2_sat);
a = T;			#Initial
i = -1;
while(i == -1):
    A = round(2.771-(0.00523*a),4);
    v1 = round(math.exp(A*(x2**2)),4);
    v2 = round(math.exp(A*(x1**2)),4);
    P1_sat  =  round(math.exp(16.59158-(3643.31/(a-33.424))),2);			#[KPa]
    P2_sat  =  round(math.exp(14.25326-(2665.54/(a-53.424))),2);			#[KPa]
    alpha = P1_sat/P2_sat;  
    P1_sat = round(P/((x1*v1)+(x2*v2/alpha)),2);
    b = round((3643.31/(16.59158-math.log(P1_sat)))+33.424,2);
    dt = abs(b-a);
    if(dt == 0):
        i = 0;
        T = b;
        break;
    a = b;

T_c = T;
y1_c = round(x1*v1*P1_sat/P,3);
y2_c = 1-y1_c;
print ('(c)T and [yi] for P = 101.33kPa and x1 = 0.')
print ('BUBL T calculations')
print 'Temperature  =  ',T_c,'K'
print 'y1  =  ',y1_c
print 'y2  =  ',y2_c

#(d) Calculate T and (xi) for P  =  101.3 KPa and y1 = 0.4
P = 101.3;
y1 = 0.4;
y2 = 1-y1;
T1_sat = round((3643.31/(16.59158-math.log(P)))+33.424,2);
T2_sat = round((2665.54/(14.25326-math.log(P)))+53.424,2);
T = (y1*T1_sat)+(y2*T2_sat);
v1 = 1;  			#Initially
v2 = 1;  			#Initially
a = T;  			#Initial
i = -1;
while(i == -1):
    A = round(2.771-(0.00523*a),4);
    P1_sat  =  round(math.exp(16.59158-(3643.31/(a-33.424))),2);			#[KPa]
    P2_sat  =  round(math.exp(14.25326-(2665.54/(a-53.424))),2);			#[KPa]
    alpha = P1_sat/P2_sat;
    x1 = round((y1*P)/(v1*P1_sat),4);
    x2 = 1-x1;
    v1 = round(math.exp(A*(x2**2)),4);
    v2 = round(math.exp(A*(x1**2)),4);
    P1_sat = P*((y1/v1)+(y2*alpha/v2));
    b = round((3643.31/(16.59158-math.log(P1_sat)))+33.424,2);
    dt = abs(a-b);
    if(dt == 0):
        T = a;
        i = 0;
        break;
    a = b;
T_d = T;
x1_d = x1;
x2_d = x2;
print ('(d)T and [xi] for P = 101.33kPa and y1 = 0.40')
print ('DEW T calculations')
print 'T  =  ',T,'K'
print 'x1  =  ',x1_d
print 'x2  =  ',x2_d

#(e) Taz , (xi_az) and (yi_az)  for T  =  318.15K 
T  =  318.15;
# Relative Volatility   alpha_12 = (y1/x1)/(y2/x2)
#At Azeotrope y1 = x1 and y2 = x2 and alpha_12 = 1
P1_sat  =  round(math.exp(16.59158-(3643.31/(T-33.424))),2);			#[KPa]
P2_sat  =  round(math.exp(14.25326-(2665.54/(T-53.424))),2);			#[KPa]
#From eqn (10.5) alpha_12 = (v1*P1_sat)/(v2*P2_sat)
A = round(2.771-(0.00523*T),4);

#When x1 = 0 v2 = 1 and v1 = math.exp(A)
alpha_12_x10 = P1_sat*math.exp(A)/P2_sat; 

#When x1 = 1 v1 = 1 and v2 = math.exp(A)
alpha_12_x11 = P1_sat/(P2_sat*math.exp(A));

#But this is not Azeotrope  (at Azeotrope alpha_12 = 1)

#v1_az/v2_az = (P2_sat/P1_sat) = K
K = P2_sat/P1_sat;

#ln(v1/v2) = ln(K) = A(1-(2*x1))
x1_az = round((A-math.log(K))/(2*A),3);
x2_az = 1-x1_az;
y1_az = x1_az;
y2_az = x2_az;
v1_az = round(math.exp(A*(x2_az**2)),3);
v2_az = round(math.exp(A*(x1_az**2)),3);
P_az = round(v1_az*P1_sat,2);

print ('Azeotropic Pressure and Azeotropic Composition for T  =  318.15K')
print 'Azeotropic Pressure  =  ',P_az,'KPa'
print 'x1_az',x1_az
print 'y1_az',y1_az
(a)P and [yi] for T = 318.15K and x1 = 0.25
BUBL P calculations
P  =   73.52 KPa
y1  =   0.282
y2  =   0.718
(b)P and [xi] for T = 318.15K and y1 = 0.60
DEW P calculations
P  =   62.89 kPa
x1  =   0.8168
x2  =   0.1832
(c)T and [yi] for P = 101.33kPa and x1 = 0.
BUBL T calculations
Temperature  =   331.2 K
y1  =   0.67
y2  =   0.33
(d)T and [xi] for P = 101.33kPa and y1 = 0.40
DEW T calculations
T  =   326.69 K
x1  =   0.4598
x2  =   0.5402
Azeotropic Pressure and Azeotropic Composition for T  =  318.15K
Azeotropic Pressure  =   73.71 KPa
x1_az 0.325
y1_az 0.325

Example 10.4 page no : 151

In [6]:
from numpy import array,round
import math 

# Variables
T = 283.15;			#[K]
#(a) Dew Point Pressure
#Species = [" Methane ";" Ethane ";" Propane "];
y = array([0.1,0.2,0.7]);

P1 = 6.9;			#[bar]
K1 = array([20,3.25,0.92]);
x1 = round(y/K1,3);

P2 = 10.34;			#[bar]
K2 = array([13.2,2.25,0.65]);
x2 = round(y/K2,3);

P3 = 8.7;			#[bar]
K3 = array([16,2.65,0.762]);
x3 = round(y/K3,3);

P = [P1,P2,P3];
x = [x1,x2,x3];
E1 = zeros(3);

# Calculations and Results
for i in range(3):
    for j in range(3):
        E1[i] = E1[i]+x[i][j];			#Summation

P_dew  =  8.7;
Ans = [[y,K1,x1,K2,x2,K3,x3],[1,0,E1[0],0,E1[1],0,E1[2]]];
print (    '               P = 6.9 bar      P = 10.34 bar     P = 8.7 bar')
print '   yi      Ki     yi/Ki     Ki     yi/Ki     Ki      yi/Ki',Ans
print ('Last Row Represents the summation')
print 'The dew Point Pressure',P_dew,'KPa'

T = 283.15;			#[K]
#(b) Bubble Point Pressure
#Species = [" Methane ";" Ethane ";" Propane "];
x = array([0.1,0.2,0.7]);

P1 = 26.2;			#[bar]
K1 = array([5.6,1.11,0.335]);
y1 = round(x*K1,3);

P2 = 27.6;			#[bar]
K2 = array([5.25,1.07,0.32]);
y2 = round(x*K2,3);

P3 = 26.54;			#[bar]
K3 = [5.49,1.1,0.33];
y3 = round(x*K3,3);

i = 1;
j = 1;
P = [P1,P2,P3];
y = [y1,y2,y3];
E2 = zeros(3);
for i in range(3):
    for j in range(3):
        E2[i] = E2[i] + y[i][j] 			#Summation

P_Bubble  =  26.54;
Ans = [[x,K1,y1,K2,y2,K3,y3],[1,0,E2[0],0,E2[1],0,E2[2]]]
print (    '               P = 26.2 bar      P = 27.6 bar     P = 26.54 bar')
print '   xi      Ki     xiKi     Ki     xiKi     Ki      xiKi',Ans
print ('Last Row Represents the summation')
print 'The Bubble Point Pressure',P_Bubble,'KPa'
               P = 6.9 bar      P = 10.34 bar     P = 8.7 bar
   yi      Ki     yi/Ki     Ki     yi/Ki     Ki      yi/Ki [[array([ 0.1,  0.2,  0.7]), array([ 20.  ,   3.25,   0.92]), array([ 0.005,  0.062,  0.761]), array([ 13.2 ,   2.25,   0.65]), array([ 0.008,  0.089,  1.077]), array([ 16.   ,   2.65 ,   0.762]), array([ 0.006,  0.075,  0.919])], [1, 0, 0.82800000000000007, 0, 1.1739999999999999, 0, 1.0]]
Last Row Represents the summation
The dew Point Pressure 8.7 KPa
               P = 26.2 bar      P = 27.6 bar     P = 26.54 bar
   xi      Ki     xiKi     Ki     xiKi     Ki      xiKi [[array([ 0.1,  0.2,  0.7]), array([ 5.6  ,  1.11 ,  0.335]), array([ 0.56 ,  0.222,  0.234]), array([ 5.25,  1.07,  0.32]), array([ 0.525,  0.214,  0.224]), [5.49, 1.1, 0.33], array([ 0.549,  0.22 ,  0.231])], [1, 0, 1.016, 0, 0.96299999999999997, 0, 1.0]]
Last Row Represents the summation
The Bubble Point Pressure 26.54 KPa

Example 10.5 page no : 152

In [5]:
import math 

# Variables
z1 = 0.45;
z2 = 0.35;
z3 = 0.2;
P = 110.;			#[KPa]
T = 353.15;			#[K]
P1_sat = 195.75;			#[KPa]
P2_sat = 97.84;			#[KPa]
P3_sat = 50.32;			#[KPa]

# Calculations
#BUBL Calculation
x1 = z1;
x2 = z2;
x3 = z3;
P_BUBL = (x1*P1_sat)+(x2*P2_sat)+(x3*P3_sat);

#DEW Calculation
y1 = z1;
y2 = z2;
y3 = z3;
P_Dew = 1/((y1/P1_sat)+(y2/P2_sat)+(y3/P3_sat));

#Since P_Bubl<P<P_dew
#Flash Calculation
K1 = P1_sat/P;
K2 = P2_sat/P;
K3 = P3_sat/P;

#Finding V from Eqn(10.17)
#E((zi*Ki)/(1+(V*(Ki-1)))) = 1

x = 0;
F_x = (((z1*K1)/(1+((K1-1)*x)))+((z2*K2)/(1+((K2-1)*x)))+((z3*K3)/(1+((K3-1)*x)))-1);
F_a = F_x;

x = 0.9;
F_x = (((z1*K1)/(1+((K1-1)*x)))+((z2*K2)/(1+((K2-1)*x)))+((z3*K3)/(1+((K3-1)*x)))-1);
F_b = F_x;
A = 0;
B = 0.9;
i = 1;
while(i == 1):
    a = A;
    F_a = (((z1*K1)/(1+((K1-1)*a)))+((z2*K2)/(1+((K2-1)*a)))+((z3*K3)/(1+((K3-1)*a)))-1);
    b = B;
    F_b = (((z1*K1)/(1+((K1-1)*b)))+((z2*K2)/(1+((K2-1)*b)))+((z3*K3)/(1+((K3-1)*b)))-1);
    x1 = ((a*F_b)-(b*F_a))/(F_b-F_a);
    F_x1 = (((z1*K1)/(1+((K1-1)*x1)))+((z2*K2)/(1+((K2-1)*x1)))+((z3*K3)/(1+((K3-1)*x1)))-1);

    if((F_a*F_x1)<0):
        flag = 1;
        A = a;
        B = x1;
    elif((F_x1*F_b)<0):
        flag = 2;
        A = x1;
        B = b;  

    x1_a = round(x1,1);
    b_a = round(b,1);
    a_a = round(a,1);
    if(x1_a == b_a):
        V = round(x1,4);
        i = 0;
        break; 
    elif(x1_a == a_a):
        root = round(x1,4);
        i = 0;
        break;   

# Results
print 'Hence By solving the polynomial V  =  ',V
L = 1-V;
#from eqn 10.16
#yi = (zi*Ki)/(1+(V*(Ki-1)))
y1 = round((z1*K1)/(1+((K1-1)*V)),4);
y2 = round((z2*K2)/(1+((K2-1)*V)),4);
y3 = round((z3*K3)/(1+((K3-1)*V)),4);
x1 = round(y1/K1,4);
x2 = round(y2/K2,4);
x3 = round(y3/K3,4);
y = [y1, y2, y3];
x = [x1, x2, x3];
print 'Moles Of liquid',L
print 'Moles Of vapor',V
print 'Mole fraction Of liquid',x
print 'Mole fraction Of vapor',y

# note : answers might be different because of rounding off error.
Hence By solving the polynomial V  =   0.8789
Moles Of liquid 0.1211
Moles Of vapor 0.8789
Mole fraction Of liquid [0.26700000000000002, 0.38769999999999999, 0.38229999999999997]
Mole fraction Of vapor [0.47520000000000001, 0.3448, 0.1749]

Example 10.6 page no : 153

In [2]:
# Program To Find the Composition of Vapor and Liquid Phases

from numpy import array,zeros,linspace,exp,round
import math 

# Variables
z = array([0.1,0.2,0.7]);
K = array([10,1.76,0.52]);

y_01 = zeros((10,3))
Sum = zeros(10)
V_01 = linspace(0.1,1,10)

# Calculations
for i in range(10):
    for j in range(3):
        y_01[i,:] = round((z*K)/(1+(V_01[i]*(K-1))),3);
    Sum[i] = sum(y_01[i,:]);

for i in range(10):
    if(Sum[i+1]<1):
        a1 = V_01[i];
        i = -1;
        break;

V_001 = linspace(a1,a1+0.1,10)

y_001 = zeros((10,3))
for i in range(10):
    for j in range(3):
        y_001[i,:] = round((z*K)/(1+(V_001[i]*(K-1))),3);
    Sum[i] = sum(y_001[i,:]);

for i in range(10):
    if(Sum[i]<1):
        a01 = V_001[i];
        i = -1;
        break;

V_0001 = linspace(a01,a01+0.01,10)

y_0001 = zeros((10,3))
for i in range(10):
    for j in range(3):
        y_0001[i,:] = round((z*K)/(1+(V_0001[i]*(K-1))),3);
    Sum[i] = sum(y_0001[i,:]);

for i in range(10):
    if(Sum[i]<1):
        a001 = V_0001[i];
        i = -1;
        break;

V = a001;
y_02 = round((z*K)/(1+(a1*(K-1))),3);
y_03 = round((z*K)/(1+((a1+0.1)*(K-1))),3);
y_0273 = round((z*K)/(1+(V*(K-1))),3);
x_0273 = round(y_0273/K,3);

Ans = [[z,K,y_02,y_03,y_0273,x_0273],[1,0,sum(y_02),sum(y_03),sum(y_0273),approx(sum(x_0273),2)]];

# Results
print '   z      K      y,V = 0.2  y,V = 0.3  y,V = 0.273 x,V = 0.273',Ans
print ('NOTE : Last Row represents the summation')
print ('Hence for V  =  0.273 E(yi)  =  0 and E(xi)  =  0')
print 'Fraction of Vapor (V)  =  ',V
   z      K      y,V = 0.2  y,V = 0.3  y,V = 0.273 x,V = 0.273 [[array([ 0.1,  0.2,  0.7]), array([ 10.  ,   1.76,   0.52]), array([ 0.357,  0.306,  0.403]), array([ 0.27 ,  0.287,  0.425]), array([ 0.286,  0.291,  0.42 ]), array([ 0.029,  0.165,  0.808])], [1, 0, 1.0660000000000001, 0.98199999999999998, 0.99699999999999989, 1.002]]
NOTE : Last Row represents the summation
Hence for V  =  0.273 E(yi)  =  0 and E(xi)  =  0
Fraction of Vapor (V)  =   0.277777777778