# Chapter 16 : Other Phase Equilibria¶

### Example 16.1 Page Number : 564¶

In :

import math

# Variables
T = 0 + 273.15;			#[K] - Temperature
P = 20*10**(5);			#[Pa] - Pressure
R = 8.314;			#[J/mol*K] - Universal gas constant

#componenet 1 : methane (1)
#componenet 2 : methanol (2)

H_constant = 1022;			#[bar] - Henry's law constant
H_constant = H_constant*10**(5);			#[Pa]

# The second virial coefficients are
B_11 = -53.9;			#[cm**(3)/mol]
B_11 = B_11*10**(-6);			#[m**(3)/mol]
B_12 = -166;			#[cm**(3)/mol]
B_12 = B_12*10**(-6);			#[m**(3)/mol]
B_22 = -4068;			#[cm**(3)/mol]
B_22 = B_22*10**(-6);			#[m**(3)/mol]

den_meth = 0.8102;			#[g/cm**(3)] - Density of methanol at 0 C
Mol_wt_meth = 32.04;			# Molecular weight of methanol
P_2_sat = 0.0401;			#[bar] - Vapour pressure of methanol at 0 C

# Calculations
#The molar volume of methanol can be calculated as
V_2_liq = (1/(den_meth/Mol_wt_meth))*10**(-6);			#[m**(3)/mol]

#The phase equilibrium equation of the components at high pressure
#y1*phi_1*P = x_1*H_1
#y2*phi_2*P = x_2*H_2

#Since methane follows Henry's law therefore methanol follows the lewis-Rnadall rule
#f_2 is the fugacity of the compressed liquid which is calculated using
#f_2 = f_2_sat*math.exp[V_2_liq*(P - P_sat_2)/(R*T)]
#where f_2_sat can be calculated using virial equation
# math.log(phi_2_sat) = math.log(f_2_sat/P_2_sat) = (B_22*P_2_sat)/(R*T)

f_2_sat = P_2_sat*math.exp((B_22*P_2_sat*10**(5))/(R*T));			#[bar]

#Putting the value of 'f_2_sat' in the math.expression of f_2 , we get
f_2 = f_2_sat*math.exp(V_2_liq*(P - P_2_sat*10**(5))/(R*T));			#[bar]

#Now let us calculate the fugacity coefficients of the species in the vapour mixture
del_12 = 2*B_12 - B_11 - B_22;			#[m**(3)/mol]

#math.log(phi_1) = (P/(R*T))*(B_11 + y2**(2)*del_12)
#math.log(phi_2) = (P/(R*T))*(B_22 + y1**(2)*del_12)

#The calculation procedure is to assume a value of y1, calculate 'phi_1' and 'phi_2' and calculate 'x_1' and 'x_2' from the phase equilibrium equations and see whether x_1 + x_2 = 1,if not then another value of y1 is assumed.

y2 = 0.1;
error=10;

while(error>0.001):
y1=1-y2;
phi_1 = math.exp((P/(R*T))*((B_11 + y2**(2)*del_12)));
phi_2 = math.exp((P/(R*T))*((B_22 + y1**(2)*del_12)));
x_1 = (y1*phi_1*P)/H_constant;
x_2 = (y2*phi_2*P)/(f_2*10**(5));
x = x_1 + x_2;
error=abs(1-x);
y2=y2 - 0.000001;

# Results
print " The solubility of methane in methanol is given by x1 = %f"%(x_1);

 The solubility of methane in methanol is given by x1 = 0.018614


### Example 16.2 Page Number : 566¶

In :

import math
from scipy.integrate import quad

# Variables
x_C2H6_1 = 0.33*10**(-4);			# Solubility of ethane in water at 25 C and 1 bar

#componenet 1 : ethane (1)
#componenet 2 : water (2)

# Z = 1 - 7.63*10**(3)*P - 7.22*10**(-5)*P**(2)

#The phase equilibrium equation of ethane is
#f_1_V = x_1*H_1
#since vapour is pure gas,  f_1_V = x_1*H_1   or, phi_1*P = x_1*H_1,  where 'phi_1' is fugacity coefficient of pure ethane
# math.log(phi) = integral('Z-1)/P) from limit '0' to 'P'

P1 = 0;
P2 = 1;
P3 = 35;

# Calculations
def f51(P):
return (1-7.63*10**(-3)*P-7.22*10**(-5)*P**(2)-1)/P

intgral =  quad(f51,P1,P2)

phi_1_1 = math.exp(intgral);			# - Fugacity coefficient of ethane at 1 bar
f_1_1 = phi_1_1*P2;			#[bar] - Fugacity of ethane at 1 bar

#Similarly

def f52(P):
return (1-7.63*10**(-3)*P-7.22*10**(-5)*P**(2)-1)/P

intgral_1 =  quad(f52,P1,P3)

phi_1_35 = math.exp(intgral_1);			# Fugacity coefficient of ethane at 35 bar
f_1_35 = phi_1_35*P3;			#[bar] - Fugacity of ethane at 35 bar

# At ethane pressure of 1 bar ,  x_C2H6_1*H_1 = phi_1_1
H_1 = phi_1_1/x_C2H6_1;			#[bar] - Henry's constant

# At ethane pressure of 35 bar ,  x_C2H6_35*H_1 = phi_1_35
x_C2H6_35 = f_1_35/H_1;			# Solubility of ethane at 35 bar pressure

# Results
print "The solubility of ethane at 35 bar is given by   x_C2H6 = %e "%(x_C2H6_35);

The solubility of ethane at 35 bar is given by   x_C2H6 = 8.525648e-04


### Example 16.4 Page Number : 571¶

In :

from scipy.optimize import fsolve
import math

# Variables
T = 200;			#[K]
R = 8.314;			#[J/mol*K] - universal gas constant
# G_E = A*x_1*x_2
A = 4000;			#[J/mol]
x_1 = 0.6;			# Mle fraction of feed composition

# Calculations and Results
# Since A is given to be independent of temperature
UCST = A/(2*R);			#[K] - Upper critical solution temperature
print " The UCST of the system is %f K"%(UCST);

# Since the given temperature is less than UCST therefore two phase can get formed at the given temperature.

# x1_alpha = 1 - x1_beta
# We know that,  x1_alpha*Y_1_alpha = x2_alpha*Y_2_alpha
# x1_alpha*math.exp[(A/(R*T))*(x2_alpha)**(2)] = (1 - x1_alpha)*math.exp[(A/(R*T))*(x1_alpha)**(2)]
# where use has beeen made of the fact that x1_alpha = 1 - x1_beta and x2_beta = 1 - x1_beta = x1_alpha .Taking math.logarithm of both side we get
# math.log(x1_alpha) + (A/(R*T))*(1 - x1_alpha)**(2) = math.log(1 - x1_alpha) + (A/(R*T))*x1_alpha**(2)
# math.log(x1_alpha/(1-x1_alpha)) = (A/(R*T))*(2*x1_alpha - 1)

def f(x1_alpha):
return  math.log(x1_alpha/(1-x1_alpha)) - (A/(R*T))*(2*x1_alpha - 1)
x1_alpha = fsolve(f,0.1)
x1_beta = fsolve(f,0.9)
# Because of symmetry 1 - x1_beta = x1_alpha

# It can be seen that the equation,  math.log(x1/(1-x1)) = (A/(R*T))*(2*x1 - 1) has two roots.

print " The composition of two liquid phases in equilibrium is given by x1_alpha = %f and x1_beta = %f"%(x1_alpha,x1_beta);

z1 = 0.6;
# z1 = x1_alpha*alpha + x1_beta*beta
# beta = 1 - alpha
alpha = (z1-x1_beta)/(x1_alpha-x1_beta);			#[mol]
Beta = 1 - alpha;			#[mol]
print " The relative amount of phases is given by alpha = %f mol and beta = %f mol"%(alpha,Beta);

# the relative amounts of the phases changes with the  feed composition

#math.log(x1/(1-x1)) = (A/(R*T))*(2*x1 - 1)

print " math.logx1/1-x1 = A/R*T*2*x1 - 1";
print " If the above equation has two real roots of x1 one for phase alpha and the other for phase beta then two liquid phases get formed";
print " and if it has no  real roots then a homogeneous liquid mixture is obtained";

 The UCST of the system is 240.558095 K
The composition of two liquid phases in equilibrium is given by x1_alpha = 0.169102 and x1_beta = 0.830898
The relative amount of phases is given by alpha = 0.348896 mol and beta = 0.651104 mol
math.logx1/1-x1 = A/R*T*2*x1 - 1
If the above equation has two real roots of x1 one for phase alpha and the other for phase beta then two liquid phases get formed
and if it has no  real roots then a homogeneous liquid mixture is obtained


### Example 16.5 Page Number : 573¶

In :

# Variables
T = 300;			#[K]
R = 8.314;			#[J/mol*K] - universal gas constant
A = 7000;			#[J/mol]

# math.log(x_1/(1-x_1)) = (A/(R*T))*(2*x_1-1)

# Calculations
def f(x_1):
return math.log(x_1/(1-x_1))-((A/(R*T))*(2*x_1-1))

x1_alpha=fsolve(f,0.1)

x1_beta=1-x1_alpha;

# Results
print "The equilibrium compositin of the two liquid phase system is given by x1_alpha \t = %f  x1_beta \t = %f"%(x1_alpha,x1_beta);

The equilibrium compositin of the two liquid phase system is given by x1_alpha 	 = 0.091897  x1_beta 	 = 0.908103


### Example 16.7 Page Number : 579¶

In :

R = 8.314;			#[J/mol*K] - Universal gas constant
M_wt_meth = 32;			# Molecular weight of methanol
M_wt_water = 18;			# Molecular weight of water
m_meth = 0.01;			#[g] - Mass of methanol added per cm**(3) of solution

den_sol = 1;			#[g/cm**(3)]

# Calculations
#The mole fraction of solute is given by
#x_2 = (moles of solute in cm**(3) of solution)/(moles of solute + moles of water) in 1 cm**(3) of solution
x_2 = (m_meth/M_wt_meth)/((m_meth/M_wt_meth)+((1-m_meth)/M_wt_water));

#We know that heat of fusion of water is
H_fus = -80;			#[cal/g] - Enthalpy change of fusion at 0 C
H_fus = H_fus*4.186*M_wt_water;			#[J/mol]

#Therefore freezing point depression is given by
# T - T_m = (R*(T**(2))*x_2)/H_fus
T_f = 273.15;			#[K] - Freezing point of water
delta_T_f = (R*(T_f**(2))*x_2)/H_fus;			#[K]

# Results
print "The depression in freezing point is given by  delta_T = %f K"%(delta_T_f);

The depression in freezing point is given by  delta_T = -0.581403 K


### Example 16.8 Page Number : 580¶

In :

# Variables
R = 8.314;			#[J/mol*K] - universal gas constant
T_f = 273.15;			#[K] - Freezing point of water
m_water = 100;			#[g] - Mass of water
m_NaCl = 3.5;			#[g] - Mass of NaCl
M_wt_water = 18.015;			# Molecular weight of water
M_wt_NaCl = 58.5;			# Molecular weight of NaCl
mol_water = m_water/M_wt_water;			#[mol] - Moles of water
mol_NaCl = m_NaCl/M_wt_NaCl;			#[mol] - Moles of NaCl

H_fus = -80;			#[cal/g] - Enthalpy change of fusion at 0 C
H_fus = H_fus*4.186*M_wt_water;			#[J/mol]

# Calculations
#Mole fraction of the solute (NaCl) is given by
x_2 = mol_NaCl/(mol_NaCl+mol_water);

x_2_act = 2*x_2;			# Actual mole fraction

#Now depression in freezing point is given by
# T - T_m = (R*(T**(2))*x_2_act)/H_fus
delta_T_f = (R*(T_f**(2))*x_2_act)/H_fus;			#[C]

#Thus freezing point of seawater = depression in freezing point

# Results
print "The freezing point of seawater is  %f C"%(delta_T_f);

The freezing point of seawater is  -2.192853 C


### Example 16.10 Page Number : 583¶

In :

# Variables
R = 8.314;			#[J/mol*K] - universal gas constant
T_b = 373.15;			#[K] - Boiling point of water
m_water = 100.;			#[g] - Mass of water
m_C12H22 = 5.;			#[g] - Mass of glucise (C12H22)
M_wt_water = 18.015;			# Molecular weight of water
M_wt_C12H22 = 342.30;			# Molecular weight of C12H22
mol_water = m_water/M_wt_water;			#[mol] - Moles of water
mol_C12H22 = m_C12H22/M_wt_C12H22;			#[mol] - Moles of C12H22

# Calculations
H_vap = 540.;			#[cal/g] - Enthalpy change of vaporisation
H_vap = H_vap*4.186*M_wt_water;			#[J/mol]

#Mole fraction of the solute (C12H22) is given by
x_2 = mol_C12H22/(mol_C12H22+mol_water);

#The boiling point elevation is given by
# T - T_b = (R*T_b**(2)*x_2**(2))/H_vap**(2)

delta_T_b = (R*T_b**(2)*x_2)/(H_vap);

# Results
print "The elevation in boiling point is given by  delta_T = %f C"%(delta_T_b)

The elevation in boiling point is given by  delta_T = 0.074611 C


### Example 16.11 Page Number : 584¶

In :

# Variables
R = 8.314;			#[J/mol*K] - Universal gas constant
T = 25 + 273.15;			#[K] - Surrounding temperature
den_water = 1000.;			#[kg/m**(3)] - Density of water
m_water = 100.;			#[g] - Mass of water
m_C12H22 = 5.;			#[g] - Mass of glucise (C12H22)
M_wt_water = 18.015;			# Molecular weight of water
M_wt_C12H22 = 342.30;			# Molecular weight of C12H22
mol_water = m_water/M_wt_water;			#[mol] - Moles of water
mol_C12H22 = m_C12H22/M_wt_C12H22;			#[mol] - Moles of C12H22

# Calculations
#Mole fraction of the water is given by
x_1 = mol_water/(mol_C12H22+mol_water);

#Molar volume of water can be calculated as
V_l_water = (1./den_water)*M_wt_water*10**(-3);			#[m**(3)/mol]

#The osmotic pressure is given by
pi = -(R*T*math.log(x_1))/V_l_water;			#[N/m**(2)]
pi = pi*10**(-5);			#[bar]

# Results
print "The osmotic pressure of the mixture is %f bar"%(pi);

The osmotic pressure of the mixture is 3.616073 bar


### Example 16.12 Page Number : 585¶

In :
# Variables
R = 8.314;			#[J/mol*K] - universal gas constant
T = 25 + 273.15;			#[K] - Surrounding temperature
den_water = 1000.;			#[kg/m**(3)] - Density of water
m_water = 100.;			#[g] - Mass of water
m_NaCl = 3.5;			#[g] - Mass of NaCl
M_wt_water = 18.015;			# Molecular weight of water
M_wt_NaCl = 58.5;			# Molecular weight of NaCl
mol_water = m_water/M_wt_water;			#[mol] - Moles of water
mol_NaCl = m_NaCl/M_wt_NaCl;			#[mol] - Moles of NaCl

# Calculations
H_fus = -80.;			#[cal/g] - Enthalpy change of fusion at 0 C
H_fus = H_fus*4.186*M_wt_water;			#[J/mol]

#Mole fraction of the solute (NaCl) is given by
x_2 = mol_NaCl/(mol_NaCl+mol_water);

x_2_act = 2*x_2;			# Actual mole fraction

x_1 = 1 - x_2_act;

#Molar volume of water can be calculated as
V_l_water = (1/den_water)*M_wt_water*10**(-3);			#[m**(3)/mol]

#The osmotic pressure is given by
pi = -(R*T*math.log(x_1))/V_l_water;			#[N/m**(2)]
pi = pi*10**(-5);			#[bar]

# Results
print "The minimum pressure to desalinate sea water is %f bar"%(pi);

The minimum pressure to desalinate sea water is 29.662232 bar


### Example 16.13 Page Number : 586¶

In :

# Variables
R = 8.314;			#[J/mol*K] - universal gas constant
T = 173.15;			#[K] - Surrounding temperature
P = 60;			#[bar]
P = P*10**(5);			#[Pa]

#componenet 1 : CO2 (1)
#componenet 2 : H2 (2)
P_1_sat = 0.1392;			#[bar] - Vapour pressre of pure solid CO2
P_1_sat = P_1_sat*10**(5);			#[bar]
V_s_1 = 27.6;			#[cm**(3)/mol] - Molar volume of solid CO2
V_s_1 = V_s_1*10**(-6);			#[m**(3)/mol]
n_1 = 0.01;			#[mol] - Initial number of moles of CO2
n_2 = 0.99;			#[mol] - Initial number of moles of H2

#Let us determine the fugacity of solid CO2 (1) at 60 C and 173.15 K
# f_1 = f_1_sat*math.exp(V_s_1*(P-P_1_sat)/(R*T))

# Calculations
#Since vapour pressure of pure solid CO2 is very small, therefore
f_1_sat = P_1_sat;
f_1 = f_1_sat*math.exp(V_s_1*(P-P_1_sat)/(R*T));

#Since gas phase is ideal therefore
# y1*P = f_1
y1 = f_1/P;

#Number of moles of CO2 in vapour phase at equilibrium can be calculated as
#y1 = (n_1_eq/(n_1_eq + n_2)). Therefore
n_1_eq = n_2*y1/(1-y1);

#Therefore moles of CO2 precipitated is
n_ppt = n_1 - n_1_eq;			#[mol]

# Results
print "The moles of CO2 precipitated is %f mol"%(n_ppt);

The moles of CO2 precipitated is 0.007417 mol


### Example 16.14 Page Number : 586¶

In :

# Variables
R = 8.314;			#[J/mol*K] - universal gas constant
T = 350.;			#[K] - Surrounding temperature

#componenet 1 : organic solid (1)
#componenet 2 : CO2 (2)

P_1_sat = 133.3;			#[N/m**(2)] - Vapour pressre of organic solid
V_s_1 = 200.;			#[cm**(3)/mol] - Molar volume of organic solid
V_s_1 = V_s_1*10.**(-6);			#[m**(3)/mol]

#/At 350 K, the values of the coefficients
B_11 = -500.;			#[cm**(3)/mol]
B_22 = -85.;			#[cm****(3)/mol]
B_12 = -430.;			#[cm**(3)/mol]

# Calculations
#From phase equilibrium equation of component 1, we get
# y1*P*phi_1 = f_1
# f_1 = f_1_sat*math.exp(V_s_1*(P-P_1_sat)/(R*T))

#Since vapour pressure of organic solid is very small, therefore
f_1_sat = P_1_sat;

# Now let us determine the fugacity coefficient of organic solid in the vapour mixture.
# math.log(phi_1) = (P/(R*T))*(B_11 + y2**(2)*del_12)
del_12 = (2*B_12 - B_11 - B_22)*10**(-6);			#[m**(3)/mol]

#It is given that the partial pressure of component 1 in the vapour mixture is 1333 N?m**(2)
# y1*P = 1333 N/m**(2) or, y1 = 1333/P
# y2 = 1- 1333/P
# math.log(phi_1) = (P/(R*T))*(B_11 + (1- 1333/P)**(2)*del_12)

#The phase equilibrium equation becomes
# y1*P*phi_1 = f_1_sat*math.exp(V_s_1*(P-P_1_sat)/(R*T))
#Taking math.log on both side we have
# math.log(y1*P) + math.log(phi_1) = math.log(f_1_sat) + (V_s_1*(P-P_1_sat)/(R*T))
# (V_s_1*(P-P_1_sat)/(R*T)) - math.log(phi_1) = math.log(1333/133.3) = math.log(10)

#substituting for math.log(phi_1) from previous into the above equation we get
# (V_s_1*(P-P_1_sat)/(R*T)) - (P/(R*T))*(B_11 + (1- 1333/P)**(2)*del_12) - math.log(10) = 0
# On simplification we get,
# 975*P**(2) - 6.7*10**(9)*P + 4.89*10**(8) = 0
# Solving the above qyadratic equation using shreedharcharya rule

P3 = (6.7*10**(9) + ((-6.7*10**(9))**(2)-4*975*4.98*10**(8))**(1./2))/(2*975);			#[Pa]
P4 = (6.7*10**(9) - ((-6.7*10**(9))**(2)-4*975*4.98*10**(8))**(1./2))/(2*975);			#[Pa]
# The second value is not possible, therefore pressure of the system is P3
P3 = P3*10**(-5);			#[bar]

# Results
print " The total pressure of the system is %f bar"%(P3);

 The total pressure of the system is 68.717948 bar