Chapter 14 : Topics In Phase Equilibria¶
Example 14.1 page no : 243¶
In [2]:
import math
from numpy import array,round
# Variables
Tc = array([126.2, 190.6]);
T = 200.; #[K]
Tr = T/Tc;
Pc = array([34, 45.99]);
P = 30.; #[Bar]
Pr = P/Pc;
R = 83.14;
bi = round(0.08664*R*Tc/Pc,3);
ai = round(((0.42748*(R**2)*(Tc**2)*(Tr**(-0.5))/Pc)/(10**5)),3); #10**-5ai
y = array([0.4, 0.6]);
a = round(((y[0]**2)*ai[0])+(2*y[0]*y[1]*math.sqrt((ai[0]*ai[1])))+((y[1]**2)*ai[1]),2); #10**-5a
b = round((y[0]*bi[0])+(y[1]*bi[1]),3);
q = (a*10**5)/(b*R*T);
Beta = 0.051612;
#Z = 1+Beta-(q*(Z-Beta)/(Z*(Z+Beta)))
Z = 1; #initial
A = Z;
for i in range(11):
B = 1+Beta-((q*Beta)*(A-Beta)/(A*(A+Beta)));
if((B-A) == 0.0001):
break;
A = B;
Z = round(B,5);
I = math.log((Z+Beta)/Z)
a1 = round((2*y[0]*ai[0])+(2*y[1]*math.sqrt((ai[0]*ai[1]))-a),2); #10**-5a1
a2 = round((2*y[1]*ai[1])+(2*y[0]*math.sqrt((ai[1]*ai[0]))-a),2); #10**-5a2
b1 = bi[0];
b2 = bi[1];
q1 = round(q*(((a1+a)/a)-(b1/b)),5);
q2 = round(q*(((a2+a)/a)-(b2/b)),5);
ln_si1 = round((((b1/b)*(Z-1))-(math.log(Z-Beta))-(q1*I)),5); #ln si1
ln_si2 = round((((b2/b)*(Z-1))-(math.log(Z-Beta))-(q2*I)),5); #ln si2
si1 = round(math.exp(ln_si1),5);
si2 = round(math.exp(ln_si2),5);
q = [q1,q2];
ln_si = [ln_si1,ln_si2];
si = [si1,si2];
a_ = [a1,a2];
b_ = [b1,b2];
q_ = [q1,q2];
print 'a = ',a,'bar cm**6 mol**-2'
print 'b = ',b,'cm**3 mol6-1'
print 'q = ',q
Ans = [Tc,Tr,Pc,bi,ai,a_,b_,q_,ln_si,si];
print ' Tc Tr Pc bi ai a b q ln si si\n',Ans
Example 14.5 page no : 244¶
In [24]:
%matplotlib inline
from matplotlib.pyplot import plot,suptitle,xlabel,ylabel,subplot
from numpy import array,linspace,log,zeros
import math
def RF(A,B,K):
#By Regula Falsi Method
i = 1;
while(i == 1):
a = A;
F_a = math.log((1-a)/a)+(2*a*K)-K;
b = B;
F_b = math.log((1-b)/b)+(2*b*K)-K;
x1 = ((a*F_b)-(b*F_a))/(F_b-F_a);
F_x1 = math.log((1-x1)/x1)+(2*x1*K)-K;
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,2);
b_a = round(b,2);
a_a = round(a,2);
if(x1_a == b_a):
return round(x1,4);
elif(x1_a == a_a):
return round(x1,4);
#G_E/RT = Ax1x2 (A)
#ln V1 = A*x2**2 = A(1-x1)**2
#ln V2 = A*x1**2
#A[(1-x1_a)**2 - (1-x1_b)**2] = ln(x1_b/x1_a) (B)
#A[(x1_a)**2-(x1_b)**2] = ln((1-x1_b)/(1-x1_a)) (C)
#x1_b = 1-x1_a (D)
#A(1 - 2*x1) = ln((1-x1)/x1) (E)
#A = a/T + b - clnT (F)
#H_E = R(a + cT)x1x2 (G)
#Cp_E = (dH_E/dT) = Rcx1x2 (H)
T = linspace(250.1,450,20);
A = (-975*(T**-1)) + 22.4 - (3*log(T));
subplot(3,2,1)
plot(T,A)
x = [250.001, 450];
y = [2 ,2];
plot(x,y,'b--')
x = [272.9, 272.9];
y = [1.9001, 2.01];
plot(x,y,'r--')
x = [391.2, 391.2];
y = [1.9001, 2.01];
plot(x,y,'g--')
suptitle('(a)A vs T')
xlabel('T(K)')
ylabel('A')
subplot(3,2,2)
T = linspace(272.9,391.2,100);
K = round((-975*(T**-1))+22.4+(-3*log(T)),4);
root = zeros((100,1));
for z in range(13):
root[z] = RF(0.4,0.49,K[z]);
for z in range(13,80):
root[z] = RF(0.01,0.49,K[z]);
for z in range(80,100):
root[z] = RF(0.4,0.49,K[z]);
x1 = root;
plot(x1,T)
x = [0 ,0.55];
y = [272.9, 272.9];
plot(x,y)
x = [0 ,0.55];
y = [391.2, 391.2];
plot(x,y)
suptitle('(b)T vs x1')
xlabel('x1')
ylabel('T(K)')
root = 1-x1;
x1 = root;
plot(x1,T)#,rect = [0,250,1,450])
x = [0.5 ,0.51];
y = [272.9, 272.9];
plot(x,y)
x = [0.5, 0.51];
y = [391.2, 391.2];
plot(x,y)
#xset('window',1)
T = linspace(250.1,450,20);
A = (-540*(T**-1)) + 21.1 - (3*log(T));
subplot(3,2,3)
plot(T,A)
x = [250.001, 450];
y = [2 ,2];
plot(x,y,'b--')
x = [346, 346];
y = [1.51, 2.2];
plot(x,y,'r--')
suptitle('(a)A vs T')
xlabel('T(K)')
ylabel('A')
subplot(3,2,4)
T = linspace(250,346,100);
K = round((-540*(T**-1))+21.1+(-3*log(T)),4);
root = zeros(100);
for z in range(100):
root[z] = RF(0.1,0.49,K[z]);
x1 = root;
plot(x1,T)#,rect = [0,250,1,450])
x = [0 ,0.55];
y = [346, 346];
plot(x,y)#,style = 3)
suptitle('(b)T vs x1')
xlabel('x1')
ylabel('T(K)')
root = 1-x1;
x1 = root;
plot(x1,T)#,rect = [0,250,1,450])
x = [0.49, 0.51];
y = [345.4, 345.4];
plot(x,y)
#xset('window',2)
T = linspace(250.1,450,20);
A = (-1500*(T**-1)) + 23.9 - (3*log(T));
subplot(3,2,5)
plot(T,A)
x = [250.001, 450];
y = [2 ,2];
plot(x,y,'b--')
x = [339.7, 339.7];
y = [1.35, 2.2];
plot(x,y,'r--')
suptitle('(a)A vs T')
xlabel('T(K)')
ylabel('A')
subplot(3,2,6)
T = linspace(339.7,450,100);
K = round((-1500*(T**-1))+23.9+(-3*log(T)),4);
root = zeros(100);
for z in range(100):
root[z] = RF(0.1,0.49,K[z]);
x1 = root;
plot(x1,T)#,rect = [0,300,1,480])
x = [0, 0.55];
y = [339.7, 339.7];
plot(x,y)#,style = 3)
suptitle('(b)T vs x1')
xlabel('x1')
ylabel('T(K)')
root = 1-x1;
x1 = root;
plot(x1,T)#,rect = [0,300,1,480])
x = [0.49, 0.51];
y = [339.7, 339.7];
plot(x,y)
Out[24]:
Example 14.8 page no : 247¶
In [25]:
%matplotlib inline
from numpy import linspace,array
from matplotlib.pyplot import plot,suptitle,xlabel,ylabel
import math
# Variables
A12 = 4.62424;
A21 = 3.35629;
alpha12 = 3.78608;
alpha21 = 1.81775;
B11 = -996;
B22 = -1245;
B12 = -567;
P1_sat = 103.264; #[kPa]
P2_sat = 5.633; #[kPa]
T = 308.15; #[K]
R = 8314;
#G_E/RT = T1 = A21*x1 + A12*x2 - Q
#V1 = math.exp[(x2**2)*[A12 + (2*(A21-A12)*x1) - Q - (x1*(dQ/dx1))]]
#V2 = math.exp[(x1**2)*[A21 + (2*(A12-A21)*x2) - Q + (x2*(dQ/dx1))]]
#Q = (alpha12*x1*alpha21*x2)/((alpha12*x1) + (alpha21*x2))
#dQ/dx1 = dQ_x1 = (alpha12*alpha21*(alpha21*x2**2 - alpha12*x1**2))/((alpha12*x1 + alpha21*x2)**2)
#P = (x1*V1*P1_sat)/si1 + (x2*V2*P2_sat)/si2
#d12 = 2B12-B11-B22
#si1 = math.exp[((B11*(P-P1_sat)) + (P*y2**2*d12)))/RT]
#si2 = math.exp[((B22*(P-P2_sat)) + (P*y1**2*d12)))/RT]
#y1 = (x1*V1*P1_sat)/(si1*P)
#y2 = (x2*V2*P2_sat)/(si2*P)
#BUBL P
x1 = linspace(0.01,1,100)
x2 = 1-x1;
P = zeros(100)
y= zeros(100)
for i in range(99):
si1 = 1; #Assumed
si2 = 1; #Assumed
dP = 100;
while(dP>0.0001):
Q = round(((alpha12*x1[i]*alpha21*x2[i])/((alpha12*x1[i]) + (alpha21*x2[i]))),4);
dQ_x1 = round((alpha12*alpha21*((alpha21*((x2[i])**2)) - (alpha12*((x1[i])**2))))/(((alpha12*x1[i]) + (alpha21*x2[i]))**2),4);
V1 = round(math.exp((x2[i]**2)*(A12 + (2*(A21-A12)*x1[i]) - Q - (x1[i]*dQ_x1))),4);
V2 = round(math.exp((x1[i]**2)*(A21 + (2*(A12-A21)*x2[i]) - Q + (x2[i]*dQ_x1))),4);
Pi = round((((x1[i]*V1*P1_sat)/si1) + ((x2[i]*V2*P2_sat)/si2)),4);
y1 = round((x1[i]*V1*P1_sat)/(si1*Pi),4);
y2 = round((x2[i]*V2*P2_sat)/(si2*Pi),4);
d12 = (2*B12)-B11-B22;
si1 = round(math.exp(((B11*(Pi-P1_sat))+(Pi*(y2**2)*d12))/(R*T)),4);
si2 = round(math.exp(((B22*(Pi-P2_sat))+(Pi*(y1**2)*d12))/(R*T)),4);
Pf = round(((x1[i]*V1*P1_sat)/si1) + ((x2[i]*V2*P2_sat)/si2),4);
dP = abs(Pf-Pi);
P[i] = Pf;
y[i] = y1;
for i in range(99):
if(P[i]>104.61):
P[i] = 0;
x1[99]= 1
y[99] = 1;
P[99] = P1_sat;
subplot(1,2,1)
P_ = [5.633, 104.6];
x = [0, 0.0117];
plot(x,P_)
P_ = [104.6, 104.6];
x = [0, 0.02];
plot(x,P_,'r')
P_ = [5.633, 5.633];
y1 = [0 ,0.02];
plot(y1,P_,'b')
P_ = [104.6, 120];
xa = [0.0117, 0.0117];
plot(xa,P_,'g--')
legend('P vs x1')
suptitle('P-x-y')
xlabel('x1,y1')
ylabel('P/kPa')
P_ = [100, 120];
y1 = [0.02, 0.02];
plot(y1,P_,'w')
subplot(1,2,2)
P_ = [104.6, 104.6];
x = [0.943, 0.96];
plot(x,P_,'r')
P_ = [104.3, 104.6];
y1 = [0.946, 0.946];
plot(y1,P_,'b--')
P_ = [104.6 ,104.8];
xb = [0.95, 0.95];
plot(xb,P_,'g--')
plot(x1,P)
plot(y,P)
legend('P* = 104.6kPa')
suptitle('P-x-y(Ether Rich Region)')
xlabel('x1,y1')
ylabel('P/kPa')
Out[25]:
Example 14.9 page no : 248¶
In [26]:
%matplotlib inline
from matplotlib.pyplot import plot,suptitle,xlabel,ylabel,subplot
from numpy import linspace,array
import math
subplot(1,2,1)
# Variables
P = 101.33; #[kPa]
T = array([333.15, 343.15, 348.15, 353.15, 353.25, 363.15, 373.15]);
P1_sat = array([52.22, 73.47, 86.40, 101.05, 101.33, 136.14, 180.04]);
P2_sat = array([19.92, 31.16, 38.55, 47.36, 47.56, 70.11, 101.33]);
P_ = P1_sat+P2_sat;
plot(P_,T,'b');
Ans = [T,P1_sat,P2_sat,P_];
print ' T P1_sat P2_sat P1_sat+P2_sat'
print Ans
#Hence P lies between T = 333.15 and 343.15
P = [101.33, 101.33];
Tp = [330.15, 343.15];
plot(P,Tp,'r--')
#Thus by interpolation
P = [50.14, 104.63];
T_ = [342.15, 342.15];
plot(P,T_,'g--');
legend('P* vs T')
suptitle('T vs P1_sat + P2_sat')
xlabel('P1_sat + P2_sat(kPa)')
ylabel('T(K)')
T_ = 342.15; #[K]
subplot(1,2,2)
plot(P1_sat,T,'b')
#Hence P1_sat @ T = T_ = 342.15 by interpolation
P = [41 ,70];
Tp = [342.15, 342.15];
plot(P,Tp,'r--')
#Thus by interpolation P1_sat = 71.3kPa
P = [71.3, 71.3];
T_ = [330.15, 342.15];
plot(P,T_,'g--');
legend('P1_sat vs T')
suptitle('T vs P1_sat')
xlabel('P1_sat(kPa)')
ylabel('T(K)')
P_sat = 71.3;
T[0] = 342.15;
P = 101.33; #[kPa]
P1_sat[0] = P_sat;
P2_sat[0] = P-P_sat;
y_I = approx(P1_sat/P,3);
y_II = approx(1-(P2_sat/P),3);
for i in range(7):
if(y_I[i]>1):
y_I[i] = 0
elif(y_II[i]>1):
y_II[i] = 0
Ans = [T,y_I,y_II];
# Results
print ' T y1_I y1_II',Ans
Example 14.10 page no : 250¶
In [27]:
%matplotlib inline
from matplotlib.pyplot import plot,suptitle,xlabel,ylabel,subplot
from numpy import linspace,array
import math
subplot(2,1,1)
# Variables
m = 4.7087;
b = 2.1941;
t = 0.3984;
# Calculations and Results
P = linspace(0,40,10);
N = (m*P)/((b+(P**t))**(1/t));
plot(P,N)
m = 0.6206;
b = 1.5454;
t = 1;
n = (m*P)/((b+(P**t))**(1/t));
plot(P,n,'b--')
legend('Toth Equation')
suptitle('Adsorption Isotherm(n vs P)')
xlabel('P(kPa)')
ylabel('n(mol/kg)')
subplot(2,1,2)
C0 = 0.4016;
C1 = -0.6471;
C2 = 0.4567;
C3 = -0.12;
n = linspace(0,1.6,20);
K = C0+(C1*n)+(C2*(n**2))+(C3*(n**3));
plot(n,K)
n = linspace(0,0.5,20);
K = C0+(C1*n);
plot(n,K,'b--')
legend('Cubic Polynomial fit')
suptitle('n/P vs n for Ethylene')
xlabel('n(mol/kg)')
ylabel('n/P(mol/kg/kPa)')
show()
