# 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
a  =   17.56 bar cm**6 mol**-2
b  =   28.607 cm**3 mol6-1
q  =   [2.3919299999999999, 4.5587999999999997]
Tc       Tr        Pc        bi         ai       a         b         q        ln si       si
[array([ 126.2,  190.6]), array([ 1.58478605,  1.04931794]), array([ 34.  ,  45.99]), array([ 26.737,  29.853]), array([ 10.995,  22.786]), [10.23, 22.449999999999999], [26.736999999999998, 29.853000000000002], [2.3919299999999999, 4.5587999999999997], [-0.056640000000000003, -0.19971], [0.94493000000000005, 0.81896999999999998]]

## 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]:
[<matplotlib.lines.Line2D at 0x1134d2050>]

## 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]:
<matplotlib.text.Text at 0x113348dd0>

## 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
T       P1_sat    P2_sat   P1_sat+P2_sat
[array([ 333.15,  343.15,  348.15,  353.15,  353.25,  363.15,  373.15]), array([  52.22,   73.47,   86.4 ,  101.05,  101.33,  136.14,  180.04]), array([  19.92,   31.16,   38.55,   47.36,   47.56,   70.11,  101.33]), array([  72.14,  104.63,  124.95,  148.41,  148.89,  206.25,  281.37])]
T      y1_I      y1_II [array([ 342.15,  343.15,  348.15,  353.15,  353.25,  363.15,  373.15]), array([ 0.70364157,  0.72505675,  0.85265963,  0.99723675,  1.        ,
0.        ,  0.        ]), array([ 0.70364157,  0.69248988,  0.61955985,  0.5326162 ,  0.53064246,
0.30810224,  0.        ])]

## 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')
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()