%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)')
# 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')
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
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'
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.
# 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