Chapter 15 : Vapour Liquid Equilibrium

Example 15.1 Page Number : 503

In [1]:
 
# Variables
T = 90+ 273.15;			#[K] - Temperature
P = 1;			#[atm] - Pressure
x_1 = 0.5748;			# Equilibrium composition of liquid phase
y_1 = 0.7725;			# Equilibrium composition of vapour phase

			# We start with 1 mol of mixture of composition z_1 = 0.6, from marterial balance we get
			# (L + V)*z_1 = L*x_1 + V*y_1
			# Since total number of moles is 1, we get
			# z_1 = L*x_1 + (1-L)*y_1

# Calculations and Results
z_1_1 = 0.6;			# - Mole fraction of benzene
L_1 = (z_1_1 - y_1)/(x_1 - y_1);
V_1 = 1 - L_1;

print " For z_1 = 0.6";
print " The moles in the liquid phase is %f mol"%(L_1);
print " The moles in the vapour phase is %f mol"%(V_1);

z_1_2 = 0.7;			# - Mole fraction of benzene
L_2 = (z_1_2 - y_1)/(x_1 - y_1);
V_2 = 1 - L_2;

print " For z_1 = 0.7";
print " The moles in the liquid phase is %f mol"%(L_2);
print " The moles in the vapour phase is %f mol"%(V_2);


			# For z = 0.8
			# The feed remains vapour and the liquid is not formed at all as it is outside the two-phase region (x_1 = 0.5748 and y_1 = 0.7725).
print " For z_1 = 0.8";

 For z_1 = 0.6
 The moles in the liquid phase is 0.872534 mol
 The moles in the vapour phase is 0.127466 mol
 For z_1 = 0.7
 The moles in the liquid phase is 0.366717 mol
 The moles in the vapour phase is 0.633283 mol
 For z_1 = 0.8

Example 15.2 Page Number : 515

In [2]:
 
import math

# Variables
			# math.log(P_1_sat) = 14.39155 - 2795.817/(t + 230.002)
			# math.log(P_2_sat) = 16.26205 - 3799.887/(t + 226.346)

			#(a)
x_1_a =0.43;			# Equilibrium composition of liquid phase
t_a = 76;			#[C] - Temperature
x_2_a = 1 - x_1_a;

# Calculations and Results
			# Since liquid phase composition is given we use the relation 
			# P = x_1*P_1_sat + x_2*P_2_sat
			# At t = 76 C
P_1_sat_a = math.exp(14.39155 - 2795.817/(t_a + 230.002));
P_2_sat_a = math.exp(16.26205 - 3799.887/(t_a + 226.346));
			# Therefore total pressure is
P_a = x_1_a*P_1_sat_a + x_2_a*P_2_sat_a;			#[kPa]
y_1_a = (x_1_a*P_1_sat_a)/(P_a);
y_2_a = (x_2_a*P_2_sat_a)/(P_a);

print "a).The system pressure is = %f kPa"%(P_a);
print "    The vapour phase composition is y_1 = %f"%(y_1_a);

			#(b)
y_1_b = 0.43;			# Equilibrium composition of vapour phase
y_2_b = 1 - y_1_b;
t_b = 76;			#[C] - Temperature

P_1_sat_b = math.exp(14.39155 - 2795.817/(t_b + 230.002));
P_2_sat_b = math.exp(16.26205 - 3799.887/(t_b + 226.346));

			# Since vapour phase composition is given ,the system pressure is given by
			# 1/P = y_1/P_1_sat + y_2/P_2_sat
P_b = 1/(y_1_b/P_1_sat_b + y_2_b/P_2_sat_b);

x_1_b = (y_1_b*P_b)/P_1_sat_b;
x_2_b = (y_2_b*P_b)/P_2_sat_b;

print "b).The system pressure is P = %f kPa"%(P_b);
print "    The liquid phase composition is x_1 = %f"%(x_1_b);

			#(c)
x_1_c = 0.32;			# Equilibrium composition of liquid phase
x_2_c = 1 - x_1_c;
P_c = 101.33;			#[kPa] - Pressure of the system

			# We have,  P = x_1*P_1_sat + x_2*P_2_sat
t_1_sat = 2795.817/(14.39155 - math.log(P_c)) - 230.002;
t_2_sat = 3799.887/(16.26205 - math.log(P_c)) - 226.346;
t = x_1_c*t_1_sat + x_2_c*t_2_sat;

error = 10;
while(error>0.1):
    P_1_sat = math.exp(14.39155 - 2795.817/(t + 230.002));
    P_2_sat = math.exp(16.26205 - 3799.887/(t + 226.346));
    P = x_1_c*P_1_sat + x_2_c*P_2_sat;
    error=abs(P - P_c);
    t = t - 0.1;

P_1_sat_c = math.exp(14.39155 - 2795.817/(t + 230.002));
P_2_sat_c = math.exp(16.26205 - 3799.887/(t + 226.346));

y_1_c = (x_1_c*P_1_sat_c)/(P_c);
y_2_c = 1 - y_1_c;

print "c).The system temperature is t = %f C"%(t);
print "    The vapour phase composition is y_1 = %f"%(y_1_c);

			#(d)
y_1_d = 0.57;			# Vapour phase composition
y_2_d = 1 - y_1_d;
P_d = 101.33;			#[kPa] - Pressure of the system

			# Since vapour phase composition is given, we can use the relation
			# 1/P = y_1/P_1_sat + y_2/P_2_sat
t_1_sat_d = 2795.817/(14.39155 - math.log(P_d)) - 230.002;
t_2_sat_d = 3799.887/(16.26205 - math.log(P_d)) - 226.346;
t_d = y_1_d*t_1_sat_d + y_2_d*t_2_sat_d;

fault = 10;
while(fault>0.1):
    P_1_sat_prime = math.exp(14.39155 - 2795.817/(t_d + 230.002));
    P_2_sat_prime = math.exp(16.26205 - 3799.887/(t_d + 226.346));
    P_prime = 1/(y_1_d/P_1_sat_prime + y_2_d/P_2_sat_prime);
    fault=abs(P_prime - P_d);
    t_d = t_d + 0.01;

P_1_sat_d = math.exp(14.39155 - 2795.817/(t_d + 230.002));
P_2_sat_d = math.exp(16.26205 - 3799.887/(t_d + 226.346));

x_1_d = (y_1_d*P_d)/(P_1_sat_d);
x_2_d = 1 - x_1_d;

print "d).The system temperature is t = %f C"%(t_d);
print "    The liquid phase composition is x_1 = %f"%(x_1_d);

a).The system pressure is = 105.268363 kPa
    The vapour phase composition is y_1 = 0.782290
b).The system pressure is P = 60.894254 kPa
    The liquid phase composition is x_1 = 0.136725
c).The system temperature is t = 79.943314 C
    The vapour phase composition is y_1 = 0.679347
d).The system temperature is t = 84.600005 C
    The liquid phase composition is x_1 = 0.234934

Example 15.3 Page Number : 516

In [3]:
 
x_1_a = 0.25;			# Equilibrium composition of liquid phase
x_2_a = 0.35;
x_3_a = 1 - x_1_a - x_2_a;
t_a = 80;			#[C] - Temperature

# Calculations and Results
			# At t = 80 C
P_1_sat_a = math.exp(14.3916 - 2795.82/(t_a + 230.00));
P_2_sat_a = math.exp(14.2724 - 2945.47/(t_a + 224.00));
P_3_sat_a = math.exp(14.2043 - 2972.64/(t_a + 209.00));

			# Since liquid phase composition is given we use the relation 
P_a = x_1_a*P_1_sat_a + x_2_a*P_2_sat_a + x_3_a*P_3_sat_a;			#[kPa]

y_1_a = (x_1_a*P_1_sat_a)/(P_a);
y_2_a = (x_2_a*P_2_sat_a)/(P_a);
y_3_a = (x_3_a*P_3_sat_a)/(P_a);

print "a).The system pressure is P = %f kPa"%(P_a);
print "    The vapour phase composition is given by y_1 = %f  y_2 = %f and y_3 = %f "%(y_1_a,y_2_a,y_3_a);

			#(2)
y_1_b = 0.50;			# Equilibrium composition of liquid phase
y_2_b = 0.30;
y_3_b = 1 - y_1_a - y_2_a;
t_b = 85;			#[C] - Temperature

			# At t = 80 C
P_1_sat_b = math.exp(14.3916 - 2795.82/(t_b + 230.00));
P_2_sat_b = math.exp(14.2724 - 2945.47/(t_b + 224.00));
P_3_sat_b = math.exp(14.2043 - 2972.64/(t_b + 209.00));

			# Since vapour phase composition is given we use the relation 
P_b = 1/(y_1_b/P_1_sat_b + y_2_b/P_2_sat_b + y_3_b/P_3_sat_b);			#[kPa]

			# Therefore we have
x_1_b = (y_1_b*P_b)/P_1_sat_b;
x_2_b = (y_2_b*P_b)/P_2_sat_b;
x_3_b = (y_3_b*P_b)/P_3_sat_b;

print "b).The system pressure is P = %f kPa"%(P_b);
print "    The liquid phase composition is given by x_1 = %f x_2 = %f and x_3 = %f "%(x_1_b,x_2_b,x_3_b);

			#(c)
x_1_c = 0.30;			# Equilibrium composition of liquid phase
x_2_c = 0.45;
x_3_c = 1 - x_1_c - x_2_c;
P_c = 80;			#[kPa] - Pressure of the system

			# We have,  P = x_1*P_1_sat + x_2*P_2_sat + x_3*P_3_sat
t_1_sat = 2795.82/(14.3916 - math.log(P_c)) - 230.00;			#[C]
t_2_sat = 2945.47/(14.2724 - math.log(P_c)) - 224.00;			#[C]
t_3_sat = 2972.64/(14.2043 - math.log(P_c)) - 209.00;			#[C]
t = x_1_c*t_1_sat + x_2_c*t_2_sat + x_3_c*t_3_sat;

error = 10;
while(error>0.5):
    P_1_sat = math.exp(14.3916 - 2795.82/(t + 230.00));
    P_2_sat = math.exp(14.2724 - 2945.47/(t + 224.00));
    P_3_sat = math.exp(14.2043 - 2972.64/(t + 209.00));
    P = x_1_c*P_1_sat + x_2_c*P_2_sat + x_3_c*P_3_sat;
    error=abs(P - P_c);
    t = t - 0.2;

P_1_sat_c = math.exp(14.3916 - 2795.82/(t + 230.00));
P_2_sat_c = math.exp(14.2724 - 2945.47/(t + 224.00));
P_3_sat_c = math.exp(14.2043 - 2972.64/(t + 209.00));

y_1_c = (x_1_c*P_1_sat_c)/(P_c);
y_2_c = (x_2_c*P_2_sat_c)/(P_c);
y_3_c = 1 - y_1_c - y_2_c;

print "c).The system temperature is t = %f C"%(t);
print "    The vapour phase composition is y_1 = %f y_2 %f and y_3 = %f"%(y_1_c,y_2_c,y_3_c);

			#(d)
y_1_d = 0.6;			# Vapour phase composition
y_2_d = 0.2;
y_3_d = 1 - y_1_d - y_2_d;
P_d = 90;			#[kPa] - Pressure of the system

			# Since vapour phase composition is given, we can use the relation
			# 1/P = y_1/P_1_sat + y_2/P_2_sat + y_3/P_3_sat
t_1_sat_d = 2795.82/(14.3916 - math.log(P_d)) - 230.00;
t_2_sat_d = 2945.47/(14.2724 - math.log(P_d)) - 224.00;
t_3_sat_d = 2972.64/(14.2043 - math.log(P_d)) - 209.00;
t_d = y_1_d*t_1_sat_d + y_2_d*t_2_sat_d + y_3_d*t_3_sat_d;

fault = 10;
while(fault>0.5):
    P_1_sat_prime = math.exp(14.3916 - 2795.82/(t_d + 230.00));
    P_2_sat_prime = math.exp(14.2724 - 2945.47/(t_d + 224.00));
    P_3_sat_prime = math.exp(14.2043 - 2972.64/(t_d + 209.00));
    P_prime = 1/(y_1_d/P_1_sat_prime + y_2_d/P_2_sat_prime + y_3_d/P_3_sat_prime);
    fault=abs(P_prime - P_d);
    t_d = t_d + 0.2;

P_1_sat_d = math.exp(14.3916 - 2795.82/(t_d + 230.00));
P_2_sat_d = math.exp(14.2724 - 2945.47/(t_d + 224.00));
P_3_sat_d = math.exp(14.2043 - 2972.64/(t_d + 209.00));

x_1_d = (y_1_d*P_d)/(P_1_sat_d);
x_2_d = (y_2_d*P_d)/(P_2_sat_d);
x_3_d = 1 - x_1_d - x_2_d;

print "d).The system temperature is t = %f C"%(t_d);
print "    The liquid phase composition is x_1 = %f x_2 = %f and x_3 = %f"%(x_1_d,x_2_d,x_3_d);

a).The system pressure is P = 108.239332 kPa
    The vapour phase composition is given by y_1 = 0.497672  y_2 = 0.316379 and y_3 = 0.185948 
b).The system pressure is P = 129.287998 kPa
    The liquid phase composition is given by x_1 = 0.259997 x_2 = 0.338895 and x_3 = 0.401108 
c).The system temperature is t = 67.020387 C
    The vapour phase composition is y_1 = 0.544826 y_2 0.357251 and y_3 = 0.097922
d).The system temperature is t = 73.127683 C
    The liquid phase composition is x_1 = 0.307471 x_2 = 0.230183 and x_3 = 0.462346

Example 15.4 Page Number : 519

In [4]:
 
import math

# Variables
T = 120;			#[C] - Temperature
P_1 = 1;			#[atm] - Initial pressure
P_1 = P_1*101.325;			#[kPa]
R = 8.314;			#[J/mol*K] - Universal gas constant

y_1 = 0.3;			# Mole fraction of propane
y_2 = 0.5;			# Mole fraction of butane
y_3 = 0.2;			# Mole fraction of hexane

			# math.log(P_1_sat) = 13.7713 - 1892.47/(t + 248.82)
			# math.log(P_2_sat) = 13.7224 - 2151.63/(t + 236.91)
			# math.log(P_3_sat) = 13.8216 - 2697.55/(t + 224.37)

# Calculations and Results
			#(a)
P_1_sat = math.exp(13.7713 - 1892.47/(T + 248.82));
P_2_sat = math.exp(13.7224 - 2151.63/(T + 236.91));
P_3_sat = math.exp(13.8216 - 2697.55/(T + 224.37));

			# Since vapour phase composition is given we can use the relation,
P_2 = 1/(y_1/P_1_sat + y_2/P_2_sat + y_3/P_3_sat);			#[kPa]

print " a).The pressure of the mixture when first drop condenses is given by P = %f kPa"%(P_2);

			#(b)
x_1 = (y_1*P_2)/P_1_sat;
x_2 = (y_2*P_2)/P_2_sat;
x_3 = (y_3*P_2)/P_3_sat;

print " b).The liquid phase composition is given by x_1 propane = %f x_2 butane = %f and x_3 hexane = %f "%(x_1,x_2,x_3);

			# (c)
			# Work transfer per mol is given by
			# W = integral(P*dV) = integral((R*T/V)*dV) = R*T*math.log(V_2/V_1) = R*T*math.log(P_1/P_2)
w = R*(T+273.15)*math.log(P_1/P_2);			#[J/mol]
			# For 100 mol
W = w*100;			#[J]
W = W*10**(-3);			#[kJ]
print " c).The work required for 100 mol of mixture handeled is %f kJ"%(W);

 a).The pressure of the mixture when first drop condenses is given by P = 1278.060855 kPa
 b).The liquid phase composition is given by x_1 propane = 0.067811 x_2 butane = 0.291139 and x_3 hexane = 0.641049 
 c).The work required for 100 mol of mixture handeled is -828.526086 kJ

Example 15.5 Page Number : 520

In [3]:
 
import math
# Variables
T = 27;			#[C] - Temperature

			# math.log(P_1_sat) = 13.8216 - 2697.55/(t + 224.37)
			# math.log(P_2_sat) = 13.8587 - 2911.32/(t + 216.64)

			#(a)
x_1_a = 0.2;
x_2_a = 1 - x_1_a;

# Calculations and Results
			# At t = 27 C
P_1_sat = math.exp(13.8216 - 2697.55/(T + 224.37));			#[kPa]
P_2_sat = math.exp(13.8587 - 2911.32/(T + 216.64));			#[kPa]
P_a = x_1_a*P_1_sat + x_2_a*P_2_sat;			#[kPa]

y_1_a = x_1_a*P_1_sat/P_a;
y_2_a = x_2_a*P_2_sat/P_a;

print "a).The total pressure is P = %f kPa"%(P_a);
print "    The vapour phase composition is given by y_1 = %f and y_2 = %f"%(y_1_a,y_2_a);

			#(b)
y_1_b = 0.2;
y_2_b = 1 - y_1_b;
			# Since vapour phase composition is given we can use the relation 
P_b = 1/(y_1_b/P_1_sat + y_2_b/P_2_sat);			#[kPa]

			# Therefore we have
x_1_b = (y_1_b*P_b)/P_1_sat;
x_2_b = (y_2_b*P_b)/P_2_sat;

print "b).The total pressure is, P = %f kPa"%(P_b);
print "    The liquid phase composition is given by x_1 = %f and x_2 = %f"%(x_1_b,x_2_b);

			#(c)
P_c = 30;			#[kPa] - Total pressure
x_1_c = 0.2;
x_2_c = 1 - x_1_c;

			# We have,  P = x_1*P_1_sat + x_2*P_2_sat
t_1_sat = 2697.55/(13.8216 - math.log(P_c)) - 224.37;
t_2_sat = 2911.32/(13.8587 - math.log(P_c)) - 216.64;
t = x_1_c*t_1_sat + x_2_c*t_2_sat;

fault = 10;
while(fault>0.3):
    P_1_sat = math.exp(13.8216 - 2697.55/(t + 224.37));
    P_2_sat = math.exp(13.8587 - 2911.32/(t + 216.64));
    P = x_1_c*P_1_sat + x_2_c*P_2_sat;
    fault = abs(P - P_c);
    t = t - 0.1;

P_1_sat_c = math.exp(13.8216 - 2697.55/(t + 224.37));
P_2_sat_c = math.exp(13.8587 - 2911.32/(t + 216.64));

y_1_c = (x_1_c*P_1_sat_c)/(P_c);
y_2_c = 1 - y_1_c;

print "c).The system temperature is t = %f C"%(t);
print "    The vapour phase composition is y_1 = %f and y_2 = %f "%(y_1_c,y_2_c);

a).The total pressure is P = 9.795726 kPa
    The vapour phase composition is given by y_1 = 0.448801 and y_2 = 0.551199
b).The total pressure is, P = 7.835131 kPa
    The liquid phase composition is given by x_1 = 0.071288 and x_2 = 0.928712
c).The system temperature is t = 53.904663 C
    The vapour phase composition is y_1 = 0.413592 and y_2 = 0.586408

Example 15.6 Page Number : 521

In [5]:
 
# Variables
P = 90;			#[kPa] - Pressure
R = 8.314;			#[J/mol*K] - Universal gas constant

			# math.log(P_sat) = A - B/(t + C)

			# For benzene
A_1 = 13.8594;
B_1 = 2773.78;
C_1 = 220.07;
			# For ethyl benzene
A_2 = 14.0045;
B_2 = 3279.47;
C_2 = 213.20;

x_1 = 0.5;			# Equimolar mixture
x_2 = 0.5;

			# The bubble point temperature equation is
			# P = x_1*P_1_sat + x_2*P_2_sat

# Calculations and Results
t_1_sat = B_1/(A_1 - math.log(P)) - C_1;
t_2_sat = B_2/(A_2 - math.log(P)) - C_2;
t = x_1*t_1_sat + x_2*t_2_sat;

fault = 10;
while(fault>0.3):
    P_1_sat = math.exp(A_1 - B_1/(t + C_1));
    P_2_sat = math.exp(A_2 - B_2/(t + C_2));
    P_net = x_1*P_1_sat + x_2*P_2_sat;
    fault=abs(P_net - P);
    t = t - 0.1;

print " The bubble point temperature is %f C"%(t);

			# Now let us determine dew point temperature for y_1 = 0.5, and P = 90 kPa
y_1 = 0.5;			# Equimolar mixture
y_2 = 0.5;

			# 1/P = y_1/P_1_sat + y_2/P_2_sat
			# Let us statrt with t = 104.07
t_old = 104.07;
error = 10;
while(error>0.3):
    P_1_sat_prime = math.exp(A_1 - B_1/(t_old + C_1));
    P_2_sat_prime = math.exp(A_2 - B_2/(t_old + C_2));
    P_net_prime = 1/(y_1/P_1_sat_prime + y_2/P_2_sat_prime);
    error=abs(P_net_prime - P);
    t_old = t_old + 0.1;

print " The dew point temperature is %f C"%(t_old);

 The bubble point temperature is 93.862003 C
 The dew point temperature is 114.470000 C

Example 15.7 Page Number : 522

In [6]:
 
# Variables
P = 1;			#[bar] - Pressure
P = P*10**(2);			#[kPa]

			# math.log(P_1_sat) = 13.8594 - 2773.78/(t + 220.07)
			# math.log(P_2_sat) = 14.0098 - 3103.01/(t + 219.79)

			# The bubble point equation is
			# P = x_1*P_1_sat + x_2*P_2_sat;

t_1_sat = 2773.78/(13.8594 - math.log(P)) - 220.07;
t_2_sat = 3103.01/(14.0098 - math.log(P)) - 219.79;

			# For x_1 = 0.1
			# t = x_1_1*t_1_sat + x_2_1*t_2_sat;
x_1 = [0.1,0.5,0.9];
x_2 = [0,0,0]

# Calculations and Results
for i in range(3):
    x_2[i] = 1 - x_1[i];
    t = x_1[i]*t_1_sat + x_2[i]*t_2_sat;
    fault = 10;
    while(fault>0.3):
        P_1_sat = math.exp(13.8594 - 2773.78/(t + 220.07));
        P_2_sat = math.exp(14.0098 - 3103.01/(t + 219.79));
        P_net = x_1[i]*P_1_sat + x_2[i]*P_2_sat;
        fault=abs(P_net - P);
        t = t - 0.1;
    print " The bubble point temperature for x_1 = %f) is %f C"%(x_1[i],t);



			# Now let us determine dew point temperature
			# 1/P = y_1/P_1_sat + y_2/P_2_sat

y_1 = [0.1,0.5,0.9];
y_2 = [0,0,0]
for j in range(3):
    y_2[j] = 1 - y_1[j];
    t_prime = y_1[j]*t_1_sat + y_2[j]*t_2_sat;
    error = 10;
    while(error>0.3):
        P_1_sat = math.exp(13.8594 - 2773.78/(t_prime + 220.07));
        P_2_sat = math.exp(14.0098 - 3103.01/(t_prime + 219.79));
        P_net = 1/(y_1[j]/P_1_sat + y_2[j]/P_2_sat);
        error=abs(P_net - P);
        t_prime = t_prime + 0.1;
    
    print " The dew point temperature for y_1 = %f is %f C"%(y_1[j],t_prime);


			#Therefore
print " The temperature range for z_1 = 0.1 which two phase exist is 105.61 C to 108.11 C";
print " The temperature range for z_1 = 0.5 which two phase exist is 91.61 C to 98.40 C";
print " The temperature range for z_1 = 0.9 which two phase exist is 81.71 C to 84.51 C";

 The bubble point temperature for x_1 = 0.100000) is 105.605549 C
 The bubble point temperature for x_1 = 0.500000) is 91.607993 C
 The bubble point temperature for x_1 = 0.900000) is 81.710437 C
 The dew point temperature for y_1 = 0.100000 is 108.105549 C
 The dew point temperature for y_1 = 0.500000 is 98.407993 C
 The dew point temperature for y_1 = 0.900000 is 84.510437 C
 The temperature range for z_1 = 0.1 which two phase exist is 105.61 C to 108.11 C
 The temperature range for z_1 = 0.5 which two phase exist is 91.61 C to 98.40 C
 The temperature range for z_1 = 0.9 which two phase exist is 81.71 C to 84.51 C

Example 15.8 Page Number : 524

In [7]:
 
from scipy.optimize import fsolve 
import math 

# Variables
x_1 = 0.20;
x_2 = 0.45;
x_3 = 0.35;
P = 10;			#[atm]
P = P*101325*10**(-3);			#[kPa]

			# math.log(P_1_sat) = 13.7713 - 1892.47/(t + 248.82)
			# math.log(P_2_sat) = 13.7224 - 2151.63/(t + 236.91)
			# math.log(P_3_sat) = 13.8183 - 2477.07/(t + 233.21)

			#(a)
			# The bubble point equation is
			# P = x_1*P_1_sat + x_2*P_2_sat + x_3*P_3_sat

# Calculations and Results
t_1_sat = 1892.47/(13.7713 - math.log(P)) - 248.82;
t_2_sat = 2151.63/(13.7224 - math.log(P)) - 236.91;
t_3_sat = 2477.07/(13.8183 - math.log(P)) - 233.21;
t = x_1*t_1_sat + x_2*t_2_sat + x_3*t_3_sat;

fault = 10;
while(fault>0.1):
    P_1_sat = math.exp(13.7713 - 1892.47/(t + 248.82));
    P_2_sat = math.exp(13.7224 - 2151.63/(t + 236.91));
    P_3_sat = math.exp(13.8183 - 2477.07/(t + 233.21));
    P_net = x_1*P_1_sat + x_2*P_2_sat + x_3*P_3_sat;
    fault=abs(P_net - P);
    t = t - 0.003;

BPT = t;
print " a).The bubble point temperature is %f C"%(BPT);

			# (b)
			# Now let us determine dew point temperature for y_1 = 0.5, and P = 90 kPa
y_1 = 0.20;
y_2 = 0.45;
y_3 = 0.35;

			# 1/P = y_1/P_1_sat + y_2/P_2_sat + y_3/P_3_sat

t_old = 90;			#[C]
error = 10;
while(error>0.1):
    P_1_sat_prime = math.exp(13.7713 - 1892.47/(t_old + 248.82));
    P_2_sat_prime = math.exp(13.7224 - 2151.63/(t_old + 236.91));
    P_3_sat_prime = math.exp(13.8183 - 2477.07/(t_old + 233.21));
    P_net_prime = 1/(y_1/P_1_sat_prime + y_2/P_2_sat_prime + y_3/P_3_sat_prime);
    error=abs(P_net_prime - P);
    t_old = t_old + 0.003;

DPT = t_old;
print " b).The dew point temperature is %f C"%(DPT);

			# (c)
			# Therefore at 82 C two phase exists
			# At 82 C and P = 1013.25 kPa pressure
T_c = 82;			#[C]
P_c = 1013.25;			#[kPa]
z_1 = 0.20;
z_2 = 0.45;
z_3 = 0.35;

P_1_sat_c = math.exp(13.7713 - 1892.47/(T_c + 248.82));
P_2_sat_c = math.exp(13.7224 - 2151.63/(T_c + 236.91));
P_3_sat_c = math.exp(13.8183 - 2477.07/(T_c + 233.21));

K_1 = P_1_sat_c/P_c;
K_2 = P_2_sat_c/P_c;
K_3 = P_3_sat_c/P_c;

			# We have to find such a V that the following equation is satisfied.
			# summation(y_i) = K_i*z_i/(1-V+V*K_i) = 1
			# K_1*z_1/(1-V+V*K_1) + K_2*z_2/(1-V+V*K_2) + K_3*z_3/(1-V+V*K_3) = 1

def f1(V): 
	 return  K_1*z_1/(1-V+V*K_1) + K_2*z_2/(1-V+V*K_2) + K_3*z_3/(1-V+V*K_3)-1
V = fsolve(f1,0.4)

			# Therefore now we have
y_1_c = K_1*z_1/(1-V+V*K_1);
y_2_c = K_2*z_2/(1-V+V*K_2);
y_3_c = K_3*z_3/(1-V+V*K_3);
x_1_c = y_1_c/K_1;
x_2_c = y_2_c/K_2;
x_3_c = y_3_c/K_3;

print " c).The proportion of vapour is given by V = %f"%(V);
print "     The composition of vapour foemed is given by y_1 = %f y_2 = %f and y_3 = %f "%(y_1_c,y_2_c,y_3_c);
print "     The composition of liquid formed is given by x_1 = %f x_2 = %f and x_3 = %f "%(x_1_c,x_2_c,x_3_c);

 a).The bubble point temperature is 71.833954 C
 b).The dew point temperature is 97.254000 C
 c).The proportion of vapour is given by V = 0.332143
     The composition of vapour foemed is given by y_1 = 0.365018 y_2 = 0.466574 and y_3 = 0.168408 
     The composition of liquid formed is given by x_1 = 0.117932 x_2 = 0.441757 and x_3 = 0.440311

Example 15.9 Page Number : 526

In [8]:
 
from scipy.optimize import fsolve 
import math 

# Variables
T = 27;			#[C] - Temperature
z_1 = 0.4;
z_2 = 0.3;
z_3 = 0.3;

			# math.log(P_sat) = A - B/(t + C)

			# For propane
A_1 = 13.7713;
B_1 = 1892.47;
C_1 = 248.82;
			# For i-butane
A_2 = 13.4331;
B_2 = 1989.35;
C_2 = 236.84;
			# For n-butane
A_3 = 13.7224;
B_3 = 2151.63;
C_3 = 236.91;

			#(a)
			
y_1 = z_1;
y_2 = z_2;
y_3 = z_3;

# Calculations and Results
			# At 27 C,
P_1_sat = math.exp(A_1 - B_1/(T + C_1));
P_2_sat = math.exp(A_2 - B_2/(T + C_2));
P_3_sat = math.exp(A_3 - B_3/(T + C_3));

			# The dew point pressure is given by
P_1 = 1/(y_1/P_1_sat + y_2/P_2_sat + y_3/P_3_sat); 


x_1 = z_1;
x_2 = z_2;
x_3 = z_3;

			# The bubble point pressure is given by
P_2 = x_1*P_1_sat + x_2*P_2_sat + x_3*P_3_sat;

print " a).The two phase region exists between %f and %f kPa"%(P_1,P_2);

			#(b)
			# The mean of the two-phase pressure range is given by
P_mean = (P_1 + P_2)/2;

			# Now let us calculate the K values of the components
K_1 = P_1_sat/P_mean;
K_2 = P_2_sat/P_mean;
K_3 = P_3_sat/P_mean;

			# summation of y_i = 1, gives
			# (K_1*z_1)/(1-V-K_1*V) + (K_2*z_2)/(1-V-K_2*V) + (K_3*z_3)/(1-V-K_3*V) = 1
			# Solving we get
def f(V): 
	 return (K_1*z_1)/(1-V+K_1*V) + (K_2*z_2)/(1-V+K_2*V) + (K_3*z_3)/(1-V+K_3*V)-1
V = fsolve(f,0.1)

y_1_prime = (z_1*K_1)/(1-V+K_1*V);

print " b).The mole fraction of propane in vapour phase is %f whereas in the feed is %f and fraction of vapour in the system is %f"%(y_1_prime,z_1,V);

 a).The two phase region exists between 421.886904 and 588.370027 kPa
 b).The mole fraction of propane in vapour phase is 0.559489 whereas in the feed is 0.400000 and fraction of vapour in the system is 0.425313

Example 15.10 Page Number : 527

In [9]:
 
# Variables
T = 50;			#[C] - Temperature
P = 64;			#[kPa] - Pressure
z_1 = 0.7;
z_2 = 0.3;

# math.log(P_sat) = A - B/(t + C)

# For acetone
A_1 = 14.37824;
B_1 = 2787.498;
C_1 = 229.664;
# For acetonitrile
A_2 = 14.88567;
B_2 = 3413.099;
C_2 = 250.523;

# Calculations
# At 50 C,
P_1_sat = math.exp(A_1 - B_1/(T + C_1));			#[kPa]
P_2_sat = math.exp(A_2 - B_2/(T + C_2));			#[kPa]

# Now let us calculate the K values of the components
K_1 = P_1_sat/P;
K_2 = P_2_sat/P;

# summation of y_i = 1, gives
# (K_1*z_1)/(1-V-K_1*V) + (K_2*z_2)/(1-V-K_2*V) = 1
# Solving we get
def f(V): 
    return (K_1*z_1)/(1-V+K_1*V) + (K_2*z_2)/(1-V+K_2*V) -1
V = fsolve(f,0.1)
L = 1 - V;
# Therefore
y_1 = (K_1*z_1)/(1-V+K_1*V);
y_2 = (K_2*z_2)/(1-V+K_2*V);

x_1 = y_1/K_1;
x_2 = y_2/K_2;

# Results
print " The value of V = %f"%(V);
print " The value of L = %f"%(L);
print " The liquid phase composition is x_1 = %f and x_2 = %f"%(x_1,x_2);
print " The vapour phase composition is y_1 = %f and y_2 = %f"%(y_1,y_2);

 The value of V = 0.450368
 The value of L = 0.549632
 The liquid phase composition is x_1 = 0.619963 and x_2 = 0.380037
 The vapour phase composition is y_1 = 0.797678 and y_2 = 0.202322

Example 15.11 Page Number : 528

In [10]:
 
# Variables
P = 12.25*101325*10**(-3);			#[kPa]
z_1 = 0.8;
z_2 = 1 - z_1;
V = 0.4;
# math.log(P_1_sat) = 13.7713 - 2892.47/(T + 248.82)
# math.log(P_2_sat) = 13.7224 - 2151.63/(T + 236.91)

# P_1_sat = math.exp(13.7713 - 21892.47/(T + 248.82));
# P_2_sat = math.exp(13.7224 - 2151.63/(T + 236.91));

# Let the total mixture be 1 mol
# We have to assume a temperature such that,
# y_1 + y_2 = (K_1*z_1)/(1-V-K_1*V) + (K_2*z_2)/(1-V-K_2*V) = 1




x_1 = z_1;
x_2 = z_2;

# The bubble point pressure is given by
# P = x_1*(math.exp(13.7713 - 21892.47/(T + 248.82))) + x_2*(math.exp(13.7224 - 2151.63/(T + 236.91)));

# Calculations
def f(T): 
	 return x_1*(math.exp(13.7713 - 1892.47/(T + 248.82))) + x_2*(math.exp(13.7224 - 2151.63/(T + 236.91))) - P
T_1 = fsolve(f,0.1)
BPT = T_1;


y_1 = z_1;
y_2 = z_2;

# The dew point equation is given by
# 1/P = y_1/P_1_sat + y_2/P_2_sat
def f1(T): 
	 return 1/(y_1/(math.exp(13.7713 - 1892.47/(T + 248.82))) + y_2/(math.exp(13.7224 - 2151.63/(T + 236.91)))) - P
T_2 = fsolve(f1,0.1)
DPT = T_2;

T = 47;			#[C]
error = 10;
while(error>0.001):
    P_1_sat = math.exp(13.7713 - 1892.47/(T + 248.82));
    P_2_sat = math.exp(13.7224 - 2151.63/(T + 236.91));
    K_1 = P_1_sat/P;
    K_2 = P_2_sat/P;
    y1 = (K_1*z_1)/(1-V+K_1*V);
    y2 = (K_2*z_2)/(1-V+K_2*V);
    y = y1 + y2;
    error=abs(y - 1);
    T = T - 0.0001;

# Results
print " The temperature when 40 mol  of mixture is in the vapour is %f C"%(T);

 The temperature when 40 mol  of mixture is in the vapour is 45.333000 C

Example 15.12 Page Number : 529

In [11]:
 
# Variables
T = 105;			#[C]
P = 1.5;			#[atm]
P = P*101325*10**(-3);			#[kPa]
z = [0.4,0.3667,0.2333];			# Feed composition
x = [0.3,0.3,0.4];			# Equilibrium liquid phase composition
y = [0.45,0.40,0.15];			# Equilibrium vapour phase composition

# From the material balance equation of component 1, we get
# (L + V)*z_1 = L*x_1 + V*y_1

# Since total moles are one, therefore  L + V = 1 and thus
# z_1 = L*x_1 + (1-L)*y_1
# Calculations and Results
for i in range(3):
    L = (z[i] - y[i])/(x[i] - y[i]);
    V = 1 - L;
    print " The number of moles in liquid phase z = %f is given by L = %f"%(z[i],L);
    print " The number of moles in vapour phase z = %f is given by V = %f"%(z[i],V);

 The number of moles in liquid phase z = 0.400000 is given by L = 0.333333
 The number of moles in vapour phase z = 0.400000 is given by V = 0.666667
 The number of moles in liquid phase z = 0.366700 is given by L = 0.333000
 The number of moles in vapour phase z = 0.366700 is given by V = 0.667000
 The number of moles in liquid phase z = 0.233300 is given by L = 0.333200
 The number of moles in vapour phase z = 0.233300 is given by V = 0.666800

Example 15.13 Page Number : 530

In [13]:
 
# Variables
T = 90;			#[C]
P = 1;			#[atm]
P = P*101325*10**(-3);			#[kPa]
z_1 = [0.1,0.5,0.8];

# Calculations
# math.log(P_1_sat) = 13.8594 - 2773.78/(t + 220.07)
# math.log(P_2_sat) = 14.0098 - 3103.01/(t + 219.79)

#At T = 90 C
P_1_sat = math.exp(13.8594 - 2773.78/(T + 220.07));
P_2_sat = math.exp(14.0098 - 3103.01/(T + 219.79));
K_1 = P_1_sat/P;
K_2 = P_2_sat/P;

# For z_1 = 0.1
# y1 = (K_1*z_1(i))/(1-V+K_1*V);
# y2 = (K_2*z_2)/(1-V+K_2*V);

x_1 = (P - P_2_sat)/(P_1_sat - P_2_sat);
y_1 = (x_1*P_1_sat)/(P);

# For two phase region to exist at 90 C and 101.325 kPa,z_1 sholud lie between x_1 and y_1

# Results
print " For two phase region to exist at 90 C and 101.325 kPa, z_1 sholud lie between %f and %f"%(x_1,y_1);
print " For z_1 = 0.1 and z_1 = 0.5 only liquid phase exists V = 0."
print " For z_1 = 0.8 only vapour exists V = 1."

 For two phase region to exist at 90 C and 101.325 kPa, z_1 sholud lie between 0.574884 and 0.772458
 For z_1 = 0.1 and z_1 = 0.5 only liquid phase exists V = 0.
 For z_1 = 0.8 only vapour exists V = 1.

Example 15.14 Page Number : 531

In [14]:
 
from numpy import zeros

# Variables
T = 90;			#[C]
P = 1;			#[atm]
P = P*101325*10**(-3);			#[kPa]
z_1 = [0.1,0.5,0.8];

# math.log(P_1_sat) = 13.8594 - 2773.78/(t + 220.07)
# math.log(P_2_sat) = 14.0098 - 3103.01/(t + 219.79)

# Calculations and Results
#(a)
#At T = 90 C
P_1_sat = math.exp(13.8594 - 2773.78/(T + 220.07));
P_2_sat = math.exp(14.0098 - 3103.01/(T + 219.79));
K_1 = P_1_sat/P;

x_1 = [0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0];
P_prime = zeros(11);
x_2 = zeros(11);
y_1 = zeros(11);

print " a.";
print " x_1    \t\t  P   \t\t   y_1      ";

for i in range(11):
    x_2[i] = 1 - x_1[i];
    P_prime[i] = x_1[i]*P_1_sat + x_2[i]*P_2_sat;
    y_1[i] = (x_1[i]*P_1_sat)/P_prime[i];
    print " %f \t  %f  \t   %f  "%(x_1[i],P_prime[i],y_1[i]);

			#(b)
T_1_sat = 2773.78/(13.8594-math.log(P)) - 220.07;			#[C]
T_2_sat = 3103.01/(14.0098-math.log(P)) - 219.79;			#[C]

T_prime = [110.62,107,104,101,98,95,92,89,86,83,80.09];

P1_sat = zeros(11);
P2_sat = zeros(11);
x_1 = zeros(11);
y_1 = zeros(11);

print "  b.";
print " TC  \t\t  P_1_sat kPa  \t\t P_2_sat kPa \t\t  x_1  \t\t  y_1 ";

for j in range(11):
    P1_sat[j] = math.exp(13.8594 - 2773.78/(T_prime[j] + 220.07));
    P2_sat[j] = math.exp(14.0098 - 3103.01/(T_prime[j] + 219.79));
    x_1[j] = (P-P2_sat[j])/(P1_sat[j]-P2_sat[j]);
    y_1[j] = (x_1[j]*P1_sat[j])/P;
    print " %f \t    %f    \t   %f    \t      %f  \t       %f "%(T_prime[j],P1_sat[j],P2_sat[j],x_1[j],y_1[j]);

 a.
 x_1    		  P   		   y_1      
 0.000000 	  54.233834  	   0.000000  
 0.100000 	  62.425251  	   0.218098  
 0.200000 	  70.616668  	   0.385597  
 0.300000 	  78.808085  	   0.518277  
 0.400000 	  86.999502  	   0.625971  
 0.500000 	  95.190920  	   0.715131  
 0.600000 	  103.382337  	   0.790162  
 0.700000 	  111.573754  	   0.854176  
 0.800000 	  119.765171  	   0.909433  
 0.900000 	  127.956588  	   0.957615  
 1.000000 	  136.148005  	   1.000000  
  b.
 TC  		  P_1_sat kPa  		 P_2_sat kPa 		  x_1  		  y_1 
 110.620000 	    237.827187    	   101.332530    	      -0.000055  	       -0.000129 
 107.000000 	    216.742019    	   91.320455    	      0.079767  	       0.170629 
 104.000000 	    200.376847    	   83.629571    	      0.151570  	       0.299740 
 101.000000 	    184.975751    	   76.460482    	      0.229134  	       0.418300 
 98.000000 	    170.500977    	   69.787769    	      0.313139  	       0.526923 
 95.000000 	    156.915208    	   63.586616    	      0.404360  	       0.626206 
 92.000000 	    144.181607    	   57.832824    	      0.503680  	       0.716718 
 89.000000 	    132.263848    	   52.502830    	      0.612106  	       0.799008 
 86.000000 	    121.126156    	   47.573718    	      0.730789  	       0.873601 
 83.000000 	    110.733338    	   43.023234    	      0.861050  	       0.941001 
 80.090000 	    101.331285    	   38.950606    	      0.999899  	       0.999961

Example 15.15 Page Number : 533

In [15]:
 
# Variables
P_1_sat = 79.80;			#[kPa]
P_2_sat = 40.45;			#[kPa]

#(1)
T = 373.15;			#[K]
x_1 = 0.05;
x_2 = 1 - x_1;
Y1 = math.exp(0.95*x_2**(2));
Y2 = math.exp(0.95*x_1**(2));

# Calculations and Results
# The total pressure of the system is given by
P = x_1*Y1*P_1_sat + x_2*Y2*P_2_sat;			#[kPa]
y_1 = x_1*Y1*P_1_sat/P;
y_2 = x_2*Y2*P_2_sat/P;

print " 1).The first bubble is formed at %f kPa and the composition y_1 = %f"%(P,y_1);

#(2)
T = 373.15;			#[K]
y_1_prime = 0.05;
y_2_prime = 1 - y_1_prime;

# Let us assume a value of x_1,
x_1_prime = 0.0001;

error = 10;
while(error>0.001):
    x_2_prime = 1 - x_1_prime;
    Y1_prime = math.exp(0.95*x_2_prime**(2));
    Y2_prime = math.exp(0.95*x_1_prime**(2));
    P_prime = x_1_prime*Y1_prime*P_1_sat + x_2_prime*Y2_prime*P_2_sat;
    x_1 = (y_1_prime*P_prime)/(Y1_prime*P_1_sat);
    error=abs(x_1_prime - x_1);
    x_1_prime = x_1_prime + 0.00001;

P_2 = x_1_prime*Y1_prime*P_1_sat + x_2_prime*Y2_prime*P_2_sat;

print " 2).The first drop is formed at %f kPa and has the composition x_1 = %f"%(P_2,x_1_prime);

 1).The first bubble is formed at 47.923166 kPa and the composition y_1 = 0.196237
 2).The first drop is formed at 41.974204 kPa and has the composition x_1 = 0.009370

Example 15.16 Page Number : 534

In [16]:
# Variables
T = 78.15;			#[C]
P_1_sat = 755;			#[mm Hg]
P_2_sat = 329;			#[mm Hg]

z_1 = 0.3;
V = 0.5;

# math.log(Y1) = 0.845/(1 + 0.845*(x_1/x_2))**(2)
# math.log(Y2) = 1/(1 + 1.183*(x_2/x_1))**(2)

# A value of x_1 is to determined for which  V = 0.5
# Let us assume a value of x_1, say x_1 = 0.150
x_1 = 0.150;

# Calculations
error = 10;
while(error>0.001):
    x_2 = 1 - x_1;
    Y1 = math.exp(0.845/(1 + 0.845*(x_1/x_2))**(2));
    Y2 = math.exp(1/(1 + 1.183*(x_2/x_1))**(2));
    P = x_1*Y1*P_1_sat + x_2*Y2*P_2_sat;
    y_1 = (x_1*Y1*P_1_sat)/P;
    V_prime = (z_1 - x_1)/(y_1 - x_1);
    error=abs(V_prime - V);
    x_1 = x_1 + 0.00001;

P_prime = x_1*Y1*P_1_sat + x_2*Y2*P_2_sat;			#[mm hg]

# At x_1 ,  V = 0.5, 
# Therefore when the mixture is 50 % vaporized at 78.15 C the mole fraction of component 1 in the liquid phase is x_1 and the system pressure is P_prime

# Results
print " The required pressure is %f mm Hg"%(P_prime);
print " and the mole fraction of component 1 in the liquid phase for this pressure is x_1 = %f"%(x_1);

 The required pressure is 505.838606 mm Hg
 and the mole fraction of component 1 in the liquid phase for this pressure is x_1 = 0.157650

Example 15.17 Page Number : 536

In [17]:
 
from scipy.stats import linregress

# Variables
T = 25;			#[C] - Temperature
P = [118.05,124.95,137.90,145.00,172.90,207.70,227.70,237.85,253.90,259.40,261.10,262.00,258.70,252.00,243.80];			#[mm Hg]
x_1 = [0.0115,0.0160,0.0250,0.0300,0.0575,0.1125,0.1775,0.2330,0.4235,0.5760,0.6605,0.7390,0.8605,0.9250,0.9625];
y_1 = [0.1810,0.2250,0.3040,0.3450,0.4580,0.5670,0.6110,0.6325,0.6800,0.7050,0.7170,0.7390,0.8030,0.8580,0.9160];

# Pressure value for which x_1 = y_1 = 0, corresponds to P_2_sat,therefore
P_2_sat = 97.45;			#[mm Hg]
# Pressure value for which x_1 = y_1 = 1, corresponds to P_1_sat,therefore
P_1_sat = 230.40;			#[mm Hg]

x_2 = zeros(15);
y_2 = zeros(15);
Y1 = zeros(15);
Y2 = zeros(15);
GE_RT = zeros(15);
x1x2_GE_RT = zeros(15);

# Calculations
for i in range(15):
    x_2[i] = 1 - x_1[i];
    y_2[i] = 1 - y_1[i];
    Y1[i] = (y_1[i]*P[i])/(x_1[i]*P_1_sat);
    Y2[i] = (y_2[i]*P[i])/(x_2[i]*P_2_sat);
    GE_RT[i] = x_1[i]*math.log(Y1[i]) + x_2[i]*math.log(Y2[i]);			# G_E/(R*T)
    x1x2_GE_RT[i] = (x_1[i]*x_2[i])/GE_RT[i];


#[M,N,sig]=reglin(x_1,x1x2_GE_RT);
M,N,sig,d,e = linregress(x_1,x1x2_GE_RT);

# Linear regression between x_1 and x_1*x_2/(G_E/R*T) gives intercept = N and slope = M

# van Laar equation is x_1*x_2/(G_E/R*T) = 1/A + (1/B - 1/A)
# 1/A = N
A = 1/N;
B = 1/(M + 1/A);

# Results
print " The value of Van Laar coefficient A =  %f"%(A);
print " The value of Van Laar coefficient B =  %f"%(B);

 The value of Van Laar coefficient A =  2.270624
 The value of Van Laar coefficient B =  1.781821

Example 15.18 Page Number : 541

In [19]:
 
from scipy.optimize import fsolve 
import math 

# Variables
T = 343.15;			#[K] - Temperature
# At 343.15 K
# math.log(Y1) = 0.95*x_2**(2)
# math.log(Y2) = 0.95*x_1**(2)
P_1_sat = 79.80;			#[kPa]
P_2_sat = 40.50;			#[kPa]

# Calculations
# At x_1 = 0,
Y1_infinity = math.exp(0.95);
alpha_12_x0 = (Y1_infinity*P_1_sat)/(P_2_sat);
# At x_1 = 1,
Y2_infinity = math.exp(0.95);
alpha_12_x1 = (P_1_sat)/(Y2_infinity*P_2_sat);


# At azeotrope, Y1*P1_sat = Y2*P2_sat
# Y2/Y1 = P_1_sat/P_2_sat
# Taking math.logarithm of both sides we get
# math.log(Y2) - math.log(Y1) = math.log(P_1_sat/P_2_sat)
# 0.95*x_1**(2) - 0.95*x_2**(2) = math.log(P_1_sat/P_2_sat)
# Solving the above equation
def f(x_1): 
	 return 0.95*x_1**(2) - 0.95*(1-x_1)**(2) - math.log(P_1_sat/P_2_sat)
x_1 = fsolve(f,0.1)

# At x_1
x_2 = 1 - x_1;
Y1 = math.exp(0.95*x_2**(2));
Y2 = math.exp(0.95*x_1**(2));
P = x_1*Y1*P_1_sat + x_2*Y2*P_2_sat;			#[kPa] - Azeotrope pressure
y_1 = (x_1*Y1*P_1_sat)/P;

# Since x_1 = y_1, (almost equal) ,the above condition is of azeotrope formation



# Results
print " Since alpha_12_x=0 = %f and alpha_12_x=1 = %f "%(alpha_12_x0,alpha_12_x1);
print " and since alpha_12 is a continuous curve and in between a value of alpha_12 = 1, shall come and at this composition the azeotrope shall get formed."
print " The azeotrope composition is x_1 = y_1 = %f"%(x_1);
print " The azeotrope presssure is %f kPa"%(P);

 Since alpha_12_x=0 = 5.094806 and alpha_12_x=1 = 0.762023 
 and since alpha_12 is a continuous curve and in between a value of alpha_12 = 1, shall come and at this composition the azeotrope shall get formed.
 The azeotrope composition is x_1 = y_1 = 0.856959
 The azeotrope presssure is 81.366308 kPa

Example 15.19 Page Number : 541

In [20]:
 
# Variables
T = 45;			#[C] - Temperature

x_1 = [0.0455,0.0940,0.1829,0.2909,0.3980,0.5069,0.5458,0.5946,0.7206,0.8145,0.8972,0.9573];
y_1 = [0.1056,0.1818,0.2783,0.3607,0.4274,0.4885,0.5098,0.5375,0.6157,0.6913,0.7869,0.8916];
P = [31.957,33.553,35.285,36.457,36.996,37.068,36.978,36.778,35.792,34.372,32.331,30.038];

			# Pressure value for which x_1 = y_1 = 0, corresponds to P_2_sat,therefore
P_2_sat = 29.819;			#[kPa]
			# Pressure value for which x_1 = y_1 = 1, corresponds to P_1_sat,therefore
P_1_sat = 27.778;			#[kPa]

x_2 = zeros(12);
y_2 = zeros(12);
Y1 = zeros(12);
Y2 = zeros(12);
alpha_12 = zeros(12);
GE_RT = zeros(12);
x1x2_GE_RT = zeros(12);
# Calculations and Results
print " x_1  \t\t  y_1    \t   P   \t\t    Y1   \t\tY2  \t    alpha_12    \t  G_E/RT \t  x1*x2/G_E/RT";

for i in range(12):
    x_2[i] = 1 - x_1[i];
    y_2[i] = 1 - y_1[i];
    Y1[i] = (y_1[i]*P[i])/(x_1[i]*P_1_sat);
    Y2[i] = (y_2[i]*P[i])/(x_2[i]*P_2_sat);
    alpha_12[i] = (y_1[i]/x_1[i])/(y_2[i]/x_2[i]);
    GE_RT[i] = x_1[i]*math.log(Y1[i]) + x_2[i]*math.log(Y2[i]);			# G_E/(R*T)
    x1x2_GE_RT[i] = (x_1[i]*x_2[i])/GE_RT[i];
    print " %f\t  %f\t  %f \t     %f  \t   %f \t     %f\t   %f  \t%f"%(x_1[i],y_1[i],P[i],Y1[i],Y2[i],alpha_12[i],GE_RT[i],x1x2_GE_RT[i]);

			#[M,N,sig]=reglin(x_1,x1x2_GE_RT);
M,N,sig,d,e=linregress(x_1,x1x2_GE_RT);

			# Linear regression between x_1 and x_1*x_2/(G_E/R*T) gives intercept = N and slope = M

			# Now let us assume the system to follow van Laar  activity coefficient model. 
			# x_1*x_2/(G_E/(R*T)) = x_1/B + x_2/A =  x_1/B + (1 - x_1)/A  =  1/A + (1/B - 1/A)*x_1 = N + M*x_1

			# 1/A = N
A = 1/N;
			# (1/B - 1/A) = M
B = 1/(M + 1/A);

print ""
print " The value of van Laar parameters are, A = %f and B = %f "%(A,B);

Y1_infinity = math.exp(A);
Y2_infinity = math.exp(B);


			
			# Now let us calculate the azeotrope composition.
			# At azeotrope, Y1*P1_sat = Y2*P2_sat
			# math.log(Y1/Y2) = math.log(P_2_sat/P_1_sat)
			# From van Laar model we get
			# math.log(P_2_sat/P_1_sat) = (A*B**(2)*2*x_2**(2))/(A*x_1 + B*x_2)**(2) + (B*A**(2)*2*x_1**(2))/(A*x_1 + B*x_2)**(2)
			# Solving the above equation
def f(x_1): 
	 return  math.log(P_2_sat/P_1_sat) - (A*B**(2)*(1-x_1)**(2))/(A*x_1 + B*(1-x_1))**(2) + (B*A**(2)*x_1**(2))/(A*x_1 + B*(1-x_1))**(2)
x_1 = fsolve(f,0.1)

print " The azeotrope composition is given by x_1 = y_1 = %f"%(x_1);

 x_1  		  y_1    	   P   		    Y1   		Y2  	    alpha_12    	  G_E/RT 	  x1*x2/G_E/RT
 0.045500	  0.105600	  31.957000 	     2.670039  	   1.004220 	     2.476833	   0.048705  	0.891698
 0.094000	  0.181800	  33.553000 	     2.336127  	   1.016177 	     2.141582	   0.094298  	0.903137
 0.182900	  0.278300	  35.285000 	     1.932808  	   1.045150 	     1.722733	   0.156610  	0.954268
 0.290900	  0.360700	  36.457000 	     1.627355  	   1.102263 	     1.375325	   0.210697  	0.979023
 0.398000	  0.427400	  36.996000 	     1.430228  	   1.180094 	     1.129007	   0.242105  	0.989635
 0.506900	  0.488500	  37.068000 	     1.285998  	   1.289486 	     0.929034	   0.252871  	0.988458
 0.545800	  0.509800	  36.978000 	     1.243394  	   1.338371 	     0.865446	   0.251278  	0.986567
 0.594600	  0.537500	  36.778000 	     1.196853  	   1.407094 	     0.792366	   0.245302  	0.982671
 0.720600	  0.615700	  35.792000 	     1.100930  	   1.650961 	     0.621199	   0.209369  	0.961630
 0.814500	  0.691300	  34.372000 	     1.050218  	   1.918247 	     0.510015	   0.160745  	0.939933
 0.897200	  0.786900	  32.331000 	     1.020818  	   2.247586 	     0.423097	   0.101740  	0.906550
 0.957300	  0.891600	  30.038000 	     1.007145  	   2.557286 	     0.366877	   0.046909  	0.871410

 The value of van Laar parameters are, A = 1.049085 and B = 1.064514 
 The azeotrope composition is given by x_1 = y_1 = 0.468254

Example 15.20 Page Number : 541

In [21]:
 
# Variables
T = 25;			#[C] - Temperature
P_1_sat = 230.4;			#[mm Hg]
P_2_sat = 97.45;			#[mm Hg]
Y1_infinity = 8.6;
Y2_infinity = 6.6;

			# Assuming ideal vpour behaviour means that phi = 1 and since system pressure is low, therefore
			# f_i = P_i_sat 
			# Assuming the activity coefficients to follow van Laar model we get
A = math.log(Y1_infinity);
B = math.log(Y2_infinity);

			#math.log(Y1) = A/(1+ (A*x_1)/(B*x_2))**(2)
			# math.log(Y2) = B/(1+ (B*x_2)/(A*x_1))**(2)

x_1 = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9];

x_2 = zeros(9);
Y1 = zeros(9);
Y2 = zeros(9);
y1_P = zeros(9);
y2_P = zeros(9);
P = zeros(9);
y_1 = zeros(9);

print " a.";
print " x_1      \t\t    Y1    \t\t    Y2    \t\t    y1*P    \t\t    y2*P   \t\t   P    \t\t   y_1 ";
# Calculations and Results
for i in range(9):
    x_2[i] = 1 - x_1[i];
    Y1[i] = math.exp(A/(1+ (A*x_1[i])/(B*x_2[i]))**(2));
    Y2[i] = math.exp(B/(1+ (B*x_2[i])/(A*x_1[i]))**(2));
    y1_P[i] = x_1[i]*Y1[i]*P_1_sat;
    y2_P[i] = x_2[i]*Y2[i]*P_2_sat;
    P[i] = x_1[i]*Y1[i]*P_1_sat + x_2[i]*Y2[i]*P_2_sat;
    y_1[i] = (x_1[i]*Y1[i]*P_1_sat)/P[i];
    print " %f\t\t %f\t\t %f \t\t %f \t\t %f \t\t %f \t   %f"%(x_1[i],Y1[i],Y2[i],y1_P[i],y2_P[i],P[i],y_1[i]);

			#(b)
			# The total system pressure versus x_1 shows a maxima and thus azeotrope is formed by the VLE system
			# The maxima occurs in the range of x_1 = 0.6 to 0.8, so an  azeotrope is formed in this composition range

			# At the azeotrope  point, Y1*P1_sat = Y2*P2_sat
			# math.log(Y1) - math.log(Y2) = math.log(P_2_sat/P_1_sat)
			# On putting the values and then solving the above equation we get
def f(x_1): 
	 return  A/(1+1.14*x_1/(1-x_1))**(2)- B/(1+0.877*(1-x_1)/x_1)**(2) - math.log(P_2_sat/P_1_sat)
x_1_prime = fsolve(f,0.1)

			# At x_1
x_2_prime = 1 - x_1_prime;
Y1_prime = math.exp(A/(1+ (A*x_1_prime)/(B*x_2_prime))**(2));
Y2_prime = math.exp(B/(1+ (B*x_2_prime)/(A*x_1_prime))**(2));
P_prime = x_1_prime*Y1_prime*P_1_sat + x_2_prime*Y2_prime*P_2_sat;			#[kPa] - Azeotrope pressure
y_1_prime = (x_1_prime*Y1_prime*P_1_sat)/P_prime;

			# Since x_1 = y_1,azeotrope formation will take place
print " b";
print " The total system pressure versus x_1 shows a maxima and thus azeotrope is formed by the VLE system";
print " The azeotrope composition is x_1 = y_1 = %f"%(x_1_prime);
print " The azeotrope presssure is %f mm Hg"%(P_prime);

 a.
 x_1      		    Y1    		    Y2    		    y1*P    		    y2*P   		   P    		   y_1 
 0.100000		 5.446877		 1.024149 		 125.496041 		 89.822961 		 215.319002 	   0.582838
 0.200000		 3.680305		 1.097308 		 169.588473 		 85.546152 		 255.134625 	   0.664702
 0.300000		 2.640401		 1.225500 		 182.504521 		 83.597451 		 266.101972 	   0.685844
 0.400000		 2.002736		 1.421865 		 184.572186 		 83.136461 		 267.708646 	   0.689452
 0.500000		 1.599580		 1.708524 		 184.271623 		 83.247849 		 267.519472 	   0.688816
 0.600000		 1.340316		 2.120132 		 185.285311 		 82.642744 		 267.928055 	   0.691549
 0.700000		 1.174189		 2.709824 		 189.373156 		 79.221704 		 268.594861 	   0.705051
 0.800000		 1.072057		 3.558780 		 197.601501 		 69.360613 		 266.962114 	   0.740186
 0.900000		 1.017109		 4.791470 		 210.907696 		 46.692876 		 257.600573 	   0.818739
 b
 The total system pressure versus x_1 shows a maxima and thus azeotrope is formed by the VLE system
 The azeotrope composition is x_1 = y_1 = 0.706548
 The azeotrope presssure is 268.599631 mm Hg

Example 15.21 Page Number : 544

In [21]:
 
T = 50;			#[C]
			# At 50 C
P_1_sat = 0.67;			#[atm]
P_2_sat = 0.18;			#[atm]
Y1_infinity = 2.5;
Y2_infinity = 7.2;

# Calculations and Results
			#(1)
			# alpha_12 = (y_1/x_1)/(y_2/x_2) = (Y1*P_1_sat)/((Y2*P_2_sat))
			# At  x_1 tending to zero,
alpha_12_x0 = (Y1_infinity*P_1_sat)/(P_2_sat);
			# At  x_1 tending to 1,
alpha_12_x1 = (P_1_sat)/((Y2_infinity*P_2_sat));


print " 1).Since alpha_12_x=0) = %f and alpha_12_x=1) = %f "%(alpha_12_x0,alpha_12_x1);
print "     and since alpha_12 is a continuous curve and in between a value of alpha_12 = 1 shall come and at this composition azeotrope shall get formed."

			
print " 2).Since the activity coefficient values are greater than 1 therefore the deviations from Roults law is positive"
print "     and the azeotrope is maximum pressure or minimum boiling";

			#(3)
			# Let us assume the system to follow van Laar activity coefficient model
A = math.log(Y1_infinity);
B = math.log(Y2_infinity);

			# math.log(Y1) = A/(1+ (A*x_1)/(B*x_2))**(2)
			# math.log(Y2) = B/(1+ (B*x_2)/(A*x_1))**(2)

			# At the azeotrope  point, Y1*P1_sat = Y2*P2_sat
			# math.log(Y1) - math.log(Y2) = math.log(P_2_sat/P_2_sat)
			# On putting the values and then solving the above equation
def f(x_1): 
	 return  A/(1+ (A*x_1)/(B*(1-x_1)))**(2)- B/(1+ (B*(1-x_1))/(A*x_1))**(2) - math.log(P_2_sat/P_1_sat)
x_1 = fsolve(f,0.1)

			# At x_1
x_2 = 1 - x_1;
Y1 = math.exp(A/(1+ (A*x_1)/(B*x_2))**(2));
Y2 = math.exp(B/(1+ (B*x_2)/(A*x_1))**(2));
P = x_1*Y1*P_1_sat + x_2*Y2*P_2_sat;			#[kPa] - Azeotrope pressure
y_1 = (x_1*Y1*P_1_sat)/P;

			# Since x_1 = y_1,the azeotrope formation will take place

print " 3).The azeotrope composition is x_1 = y_1 = %f"%(x_1);
print "     The azeotrope presssure is %f atm"%(P);

 1).Since alpha_12_x=0) = 9.305556 and alpha_12_x=1) = 0.516975 
     and since alpha_12 is a continuous curve and in between a value of alpha_12 = 1 shall come and at this composition azeotrope shall get formed.
 2).Since the activity coefficient values are greater than 1 therefore the deviations from Roults law is positive
     and the azeotrope is maximum pressure or minimum boiling
 3).The azeotrope composition is x_1 = y_1 = 0.910173
     The azeotrope presssure is 0.689143 atm

Example 15.22 Page Number : 545

In [22]:
 
# Variables
T = 25;			#[C]
			# At 50 C
P_1_sat = 7.866;			#[kPa]
P_2_sat = 3.140;			#[kPa]

			# G_E/(R*T) = 1.4938*x_1*x_2/(1.54*x_1 + 0.97*x_2)

			# The excess Gibbs free energy math.expression can be written as
			# x_1*x_2/(G_E/(R*T)) = 1.54*x_1/1.4938 + 0.97*x_2/1.4938 = x_1/0.97 + x_2/1.54

			# Comparing with the van Laar math.expression
			# x_1*x_2/(G_E/(R*T)) = x_1/B + x_2/A,  we get
A = 1.54;
B = 0.97;

			# The activity coefficients are thus given by
			# math.log(Y1) = A/(1+ (A*x_1)/(B*x_2))**(2)
			# math.log(Y2) = B/(1+ (B*x_2)/(A*x_1))**(2)

x_1 = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95];

x_2 = zeros(10);
Y1 = zeros(10);
Y2 = zeros(10);
P = zeros(10);
y_1 = zeros(10);

print " x_1      \t\t    Y1    \t\t    Y2    \t\t P kPa    \t\t  y_1 ";

# Calculations and Results
for i in range(10):
    x_2[i] = 1 - x_1[i];
    Y1[i] = math.exp(A/(1+ (A*x_1[i])/(B*x_2[i]))**(2));
    Y2[i] = math.exp(B/(1+ (B*x_2[i])/(A*x_1[i]))**(2));
    P[i] = x_1[i]*Y1[i]*P_1_sat + x_2[i]*Y2[i]*P_2_sat;
    y_1[i] = (x_1[i]*Y1[i]*P_1_sat)/P[i];
    print " %f\t\t %f\t\t %f \t\t %f \t\t  %f "%(x_1[i],Y1[i],Y2[i],P[i],y_1[i]);


			# The azeotrope is formed near x_1 = 0.95 as in this region a maxima in pressure is obtained.

			# At the azeotrope  point, Y1*P1_sat = Y2*P2_sat
			# math.log(Y1) - math.log(Y2) = math.log(P_2_sat/P_2_sat)
			# On putting the values and then solving the above equation
def f(x_1): 
	 return  A/(1+ (A*x_1)/(B*(1-x_1)))**(2)- B/(1+ (B*(1-x_1))/(A*x_1))**(2) - math.log(P_2_sat/P_1_sat)
x_1_prime = fsolve(f,0.1)

			# At x_1
x_2_prime = 1 - x_1_prime;
Y1_prime = math.exp(A/(1+ (A*x_1_prime)/(B*x_2_prime))**(2));
Y2_prime = math.exp(B/(1+ (B*x_2_prime)/(A*x_1_prime))**(2));
P_prime = x_1_prime*Y1_prime*P_1_sat + x_2_prime*Y2_prime*P_2_sat;			#[kPa] - Azeotrope pressure
y_1_prime = (x_1_prime*Y1_prime*P_1_sat)/P_prime;

			# Since x_1_prime = y_1_prime,the azeotrope formation will take place

print " Part II ";
print " The azeotrope composition is x_1 = y_1 = %f"%(x_1_prime);
print " The azeotrope presssure is %f kPa "%(P_prime);

 x_1      		    Y1    		    Y2    		 P kPa    		  y_1 
 0.100000		 3.042798		 1.022050 		 5.281779 		  0.453155 
 0.200000		 2.201629		 1.081457 		 6.180223 		  0.560433 
 0.300000		 1.725242		 1.172375 		 6.648107 		  0.612389 
 0.400000		 1.438293		 1.292347 		 6.960227 		  0.650186 
 0.500000		 1.258593		 1.440723 		 7.211980 		  0.686364 
 0.600000		 1.144175		 1.617876 		 7.432102 		  0.726584 
 0.700000		 1.072060		 1.824770 		 7.621913 		  0.774475 
 0.800000		 1.028913		 2.062721 		 7.770131 		  0.833286 
 0.900000		 1.006610		 2.333241 		 7.858834 		  0.906775 
 0.950000		 1.001587		 2.481217 		 7.874109 		  0.950528 
 Part II 
 The azeotrope composition is x_1 = y_1 = 0.958680
 The azeotrope presssure is 7.874467 kPa

Example 15.23 Page Number : 547

In [23]:
 
# Variables
T = 58.7;			#[C]
P = 1;			#[atm]
P = P*101325*10**(-3);			#[kPa]

			# math.log(P_sat) = 16.6758 - 3674.49/(t + 226.45) - For ethyl alcohol
			# math.log(P_sat) = 13.8216 - 2697.55/(t + 224.37) - For hexane

			# Let us take hexane as (1) and ethanol as (2)
			# At 58.7 C
P_1_sat = math.exp(13.8216 - 2697.55/(T + 224.37));			#[kPa]
P_2_sat = math.exp(16.6758 - 3674.49/(T + 226.45));			#[kPa]

# Calculations and Results
Y1 = P/P_1_sat;
Y2 = P/P_2_sat;

x_2 = 0.332;			# Mol % of ethanol (given)
x_1 = 1 - x_2;			# Mol % of hehane

			# The van Laar parameters are given by
A = ((1 + (x_2*math.log(Y2))/(x_1*math.log(Y1)))**(2))*math.log(Y1);
B = ((1 + (x_1*math.log(Y1))/(x_2*math.log(Y2)))**(2))*math.log(Y2);

print " The value of van Laar parameters are A = %f and B = %f "%(A,B);

			# Now let us calvulate the distribution coefficient K
x_1_prime = 0.5;			#[given]
x_2_prime = 1 - x_1_prime;

			# The activity coefficients are thus given by
			# math.log(Y1) = A/(1+ (A*x_1)/(B*x_2))**(2)
			# math.log(Y2) = B/(1+ (B*x_2)/(A*x_1))**(2)

Y1_prime = math.exp(A/(1+ (A*x_1_prime)/(B*x_2_prime))**(2));
Y2_prime = math.exp(B/(1+ (B*x_2_prime)/(A*x_1_prime))**(2));
P_prime = x_1_prime*Y1_prime*P_1_sat + x_2_prime*Y2_prime*P_2_sat;

			# We have,  K_1 = y_1/x_1 = Y1*P_1_sat/P
K_1 = Y1_prime*P_1_sat/P_prime;

print " The distribution coefficient is given by  K_1  =  %f"%(K_1)

 The value of van Laar parameters are A = 1.669879 and B = 2.662289 
 The distribution coefficient is given by  K_1  =  1.352866