import math
# Variables,
T = 300; #[K] - Temperature
b = 100; #[cal/mol]
R = 1.987; #[cal/mol*K] - Universal gas constant
# R*T*math.log(Y_1) = b*x_2**(2)
# R*T*math.log(Y_2) = b*x_1**(2)
#For equimolar mixture
x_1 = 0.5; #Mole fraction of component 1
x_2 = 0.5; #Mole fraction of component 2
# Calculations and Results
#The excess Gibbs free energy is given by
# G_excess = R*T*(x_1*math.log(Y_1) + x_2*math.log(Y_2)) = b*x_1*x_2**(2) + b*x_2*x_1**(2) = b*x_1*(x_1 + x_2) = b*x_1*x_2
G_excess = b*x_1*x_2; #[cal/mol]
#The ideal Gibbs free energy change of mixing is given by,
delta_G_id_mix = R*T*(x_1*math.log(x_1)+x_2*math.log(x_2)); #[cal/mol]
#The Gibbs free energy of mixing is given by
delta_G_mix = delta_G_id_mix + G_excess; #[cal/mol]
# delta_S_mix = delta_S_id_mix = - R*sum(x_i*math.log(x_i))
#delta_G_mix = delta_H_mix - T*delta_S_mix = delta_H_mix + R*T*(x_1*math.log(x_1)+x_2*math.log(x_2))
delta_H_mix = b*x_1*x_2; #[cal/mol]
print "The value of Gibbs free energy change for equimolar mixture formation is %f cal/mol"%(delta_G_mix);
print "The value of enthalpy change for equimolar mixture formation is %f cal/mol"%(delta_H_mix);
#Work required for separation of mixture into pure components is
W = delta_G_mix;
print "The least amount of work required for separation at 300 K is %f cal/mol"%(W);
import math
from numpy import zeros
from scipy.stats import linregress
# Variables,
T = 25 +173.15; #[K] - Temperature
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=[8.0640,7.6260,7.2780,7.2370,5.9770,4.5434,3.4019,2.8023,1.7694,1.3780,1.2302,1.1372,1.0478,1.0145,1.0070];
y_2=[1.0037,1.0099,1.0102,1.0047,1.0203,1.0399,1.1051,1.1695,1.4462,1.8520,2.2334,2.6886,3.7489,4.8960,5.6040];
x_2 = zeros(15); # x_2 = (1 - x_1)
G_RT = zeros(15); # G_RT = G_excess/(R*T)
x1x2_G_RT = zeros(15); # x1x2_G_RT = (x_1*x_2/(G_excess/(R*T)))
G_RTx1x2 = zeros(15); # G_RTx1x1 = G_excess/(R*T*x_1*x_2)
# Calculations and Results
for i in range(15):
x_2[i]=(1-x_1[i]);
G_RT[i]=(x_1[i]*math.log(y_1[i]))+(x_2[i]*math.log(y_2[i]));
x1x2_G_RT[i]=(x_1[i]*x_2[i])/G_RT[i];
G_RTx1x2[i]=1/x1x2_G_RT[i];
slop,intr,sig,d,e=linregress(x_1,x1x2_G_RT);
A = 1/intr;
B = 1/(slop+(1/A));
print " The value of van Laar parameters are A = %f and B = %f"%(A,B);
# Now from Margules equation
# G_RTx1x2 = G_excess/(R*T*x_1*x_2) = B1*x_1 + A1*x_1 = A1 + (B1 - A1)*x_1
#slope = (B1 - A1) and intercept = A1
# Again employing the concept of linear regression of the data ( x_1 , G_RTx1x2 ) to find the value of intercept and slope of the above equation
#Let slope = slop1 and intercept = intr1
slop1,intr1,sig1,d,e=linregress(x_1,G_RTx1x2);
A1 = intr1;
B1 = slop1 + A1;
print " The value of Margules parameters are A = %f and B = %f"%(A1,B1);
print " Because of the higher value of correlation factor for Van Laar model it fits the data better."
# Variables,
T = 60 + 273.15; #[K] - Temperature
R = 1.987; #[cal/mol*K] - Universal gas constant
#component 1 = acetone
#component 2 = water
x_1 = 0.3; # Mole fraction of component 1
x_2 = 1 - x_1; #Mole fraction of component 2
V_mol_1 = 74.05; #[cm**(3)/mol] - Molar volume of pure component 1
V_mol_2 = 18.07; #[cm**(3)/mol] - Molar volume of pure component 2
#for Wilson equation
a_12 = 291.27; #[cal/mol]
a_21 = 1448.01; #[cal/mol]
#For NRTL
b_12 = 631.05; #[cal/mol]
b_21 = 1197.41; #[cal/mol]
alpha = 0.5343;
# Calculations and Results
#Froom wilson equation
A_12=(V_mol_2/V_mol_1)*(math.exp(-a_12/(R*T)));
A_21 = (V_mol_1/V_mol_2)*(math.exp(-a_21/(R*T)));
Beta = A_12/(x_1+x_2*A_12) - A_21/(x_2+x_1*A_21);
#math.log(Y1) = -math.log(x_1 + x_2*A_12) + x_2*Beta;
Y1 = math.exp(-math.log(x_1+x_2*A_12)+x_2*Beta);
#similarly for Y2
Y2 = math.exp(-math.log(x_2+x_1*A_21)-x_1*Beta);
print "The value of activity coefficients for Wilson equation are Y1 = %f \t and \t Y2 = %f"%(Y1,Y2);
#From NRTL equation,
t_12 = b_12/(R*T);
t_21 = b_21/(R*T);
G_12 = math.exp(-alpha*t_12);
G_21 = math.exp(-alpha*t_21);
#math.log(Y1) = x_1**(2)*[t_12*(G_12/(x_1+x_2*G_12))**(2) + (G_12*t_12)/((G_12/(x_1+x_2*G_12))**(2))]
Y1_prime = math.exp(x_2**(2)*(t_21*(G_21/(x_1+x_2*G_21))**(2)+(G_12*t_12)/(((x_2+x_1*G_12))**(2))));
#Similarly for Y2
Y2_prime = math.exp(x_1**(2)*(t_12*(G_12/(x_2+x_1*G_12))**(2)+(G_21*t_21)/(((x_1+x_2*G_21))**(2))));
print "The value of activity coefficients for NRTL equation are Y1 = %f \t and \t Y2 = %f"%(Y1_prime,Y2_prime);
# Variables
T = 307; #[K]
x_1 = 0.047;
x_2 = 1 - x_1;
# The subgroups in the two components are
# Acetone (1) : 1 CH3, 1 CH3CO
# n-pentane (2) : 2 CH3, 3 CH2
#Group volume (Rk) and surface area (Qk) parameters of the subgroup are
k_CH3 = 1;
k_CH2 = 2;
k_CH3CO = 19;
Rk_CH3 = 0.9011;
Rk_CH2 = 0.6744;
Rk_CH3CO = 1.6724;
Qk_CH3 = 0.848;
Qk_CH2 = 0.540;
Qk_CH3CO = 1.4880;
# Interaction parameters of the subgroups in kelvin (K) are
a_1_1 = 0;
a_1_2 = 0;
a_1_19 = 476.40;
a_2_1 = 0;
a_2_2 = 0;
a_2_19 = 476.40;
a_19_1 = 26.76;
a_19_2 = 26.76;
a_19_19 = 0;
# Calculations
r_1 = 1*Rk_CH3 + 1*Rk_CH3CO;
r_2 = 2*Rk_CH3 + 3*Rk_CH2;
q_1 = 1*Qk_CH3 + 1*Qk_CH3CO;
q_2 = 2*Qk_CH3 + 3*Qk_CH2;
J_1 = r_1/(r_1*x_1+r_2*x_2);
J_2 = r_2/(r_1*x_1+r_2*x_2);
L_1 = q_1/(q_1*x_1+q_2*x_2);
L_2 = q_2/(q_1*x_1+q_2*x_2);
t_1_1 = math.exp(-a_1_1/T);
t_1_2 = math.exp(-a_1_2/T);
t_1_19 = math.exp(-a_1_19/T);
t_2_1 = math.exp(-a_2_1/T);
t_2_2 = math.exp(-a_2_2/T);
t_2_19 = math.exp(-a_2_19/T);
t_19_1 = math.exp(-a_19_1/T);
t_19_2 = math.exp(-a_19_2/T);
t_19_19 = math.exp(-a_19_19/T);
e_1_1 = 1*Qk_CH3/q_1;
e_2_1 = 0;
e_19_1 = (1*Qk_CH3CO/q_1);
e_1_2 = 2*Qk_CH3/q_2;
e_2_2 = 3*Qk_CH2/q_2;
e_19_2 = 0;
B_1_1 = e_1_1*t_1_1 + e_2_1*t_2_1 + e_19_1*t_19_1;
B_1_2 = e_1_1*t_1_2 + e_2_1*t_2_2 + e_19_1*t_19_2;
B_1_19 = e_1_1*t_1_19 + e_2_1*t_2_19 + e_19_1*t_19_19;
B_2_1 = e_1_2*t_1_1 + e_2_2*t_2_1 + e_19_2*t_19_1;
B_2_2 = e_1_2*t_1_2 + e_2_2*t_2_2 + e_19_2*t_19_2;
B_2_19 = e_1_2*t_1_19 + e_2_2*t_2_19 + e_19_2*t_19_19;
theta_1 = (x_1*q_1*e_1_1 + x_2*q_2*e_1_2)/(x_1*q_1 + x_2*q_2);
theta_2 = (x_1*q_1*e_2_1 + x_2*q_2*e_2_2)/(x_1*q_1 + x_2*q_2);
theta_19 = (x_1*q_1*e_19_1 + x_2*q_2*e_19_2)/(x_1*q_1 + x_2*q_2);
s_1 = theta_1*t_1_1 + theta_2*t_2_1 + theta_19*t_19_1;
s_2 = theta_1*t_1_2 + theta_2*t_2_2 + theta_19*t_19_2;
s_19 = theta_1*t_1_19 + theta_2*t_2_19 + theta_19*t_19_19;
# math.log(Y1_C) = 1 - J_1 + math.log(J_1) - 5*q_1*(1- (J_1/L_1) + math.log(J_1/L_1))
# math.log(Y2_C) = 1 - J_2 + math.log(J_2) - 5*q_2*(1- (J_2/L_2) + math.log(J_2/L_2))
Y1_C = math.exp(1 - J_1 + math.log(J_1) - 5*q_1*(1- (J_1/L_1) + math.log(J_1/L_1)));
Y2_C = math.exp(1 - J_2 + math.log(J_2) - 5*q_2*(1- (J_2/L_2) + math.log(J_2/L_2)));
# For species 1
summation_theta_k_1 = theta_1*(B_1_1/s_1) + theta_2*(B_1_2/s_2) + theta_19*(B_1_19/s_19);
summation_e_ki_1 = e_1_1*math.log(B_1_1/s_1) + e_2_1*math.log(B_1_2/s_2) + e_19_1*math.log(B_1_19/s_19);
# For species 2
summation_theta_k_2 = theta_1*(B_2_1/s_1) + theta_2*(B_2_2/s_2) + theta_19*(B_2_19/s_19);
summation_e_ki_2 = e_1_2*math.log(B_2_1/s_1) + e_2_2*math.log(B_2_2/s_2) + e_19_2*math.log(B_2_19/s_19);
# math.log(Y1_R) = q_1(1 - summation_theta_k_1 + summation_e_ki_1)
# math.log(Y2_R) = q_2(1 - summation_theta_k_2 + summation_e_ki_2)
Y1_R = math.exp(q_1*(1 - summation_theta_k_1 + summation_e_ki_1));
Y2_R = math.exp(q_2*(1 - summation_theta_k_2 + summation_e_ki_2));
# math.log(Y1) = math.log(Y1_C) + math.log(Y1_R)
# math.log(Y2) = math.log(Y2_C) + math.log(Y2_R)
Y1 = math.exp(math.log(Y1_C) + math.log(Y1_R));
Y2 = math.exp(math.log(Y2_C) + math.log(Y2_R));
# Results
print " The activity coefficients are Y1 = %f and Y2 = %f"%(Y1,Y2);
# Variables,
T = 25 + 273.15; #[K] - Temperature
R = 1.987; #[cal/mol*K] - Universal gas constant
#component 1 = chloroform
#component 2 = carbon tetrachloride
x_1 = 0.5; #Mole fraction of component 1 #Equimolar mixture
x_2 = 0.5; #Mole fraction of component 2
V_mol_1 = 81; #[cm**(3)/mol] - Molar volume of pure component 1
V_mol_2 = 97; #[cm**(3)/mol] - Molar volume of pure component 2
del_1 = 9.2; #[(cal/cm**(3))**(1/2)] - Mole fraction of component 1
del_2 = 8.6; #[(cal/cm**(3))**(1/2)] - Mole fraction of component 2
# Calculations
#Scatchard - Hilderbrand model
phi_1 = (x_1*V_mol_1)/(x_1*V_mol_1+x_2*V_mol_2);
phi_2 = (x_2*V_mol_2)/(x_1*V_mol_1+x_2*V_mol_2);
#math.log(Y1) = (V_mol_1/(R*T))*phi_1**(2)*(del_1-del_2)**(2)
Y1 = math.exp((V_mol_1/(R*T))*(phi_1**(2))*((del_1-del_2)**(2)));
#Similarly, for Y2
Y2 = math.exp((V_mol_2/(R*T))*(phi_2**(2))*((del_1-del_2)**(2)));
# Results
print "The value of activity coefficients for Scatchard-Hilderbrand model are Y1 = %f \t and \t Y2 = %f"%(Y1,Y2);
# Variables,
T = 25 + 273.15; #[K] - Temperature
mol_HCl = 0.001; #[mol/kg] - Molality of HCl
A = 0.510; #[(kg/mol)**(1/2)]
Z_positive = 1; #Stoichiometric coefficient of 'H' ion
Z_negative = -1; #Stoichiometric coefficient of 'Cl' ion
m_H_positive = mol_HCl; #
m_Cl_negative = mol_HCl;
# Calculations and Results
# I = 1/2*[((Z_positive)**(2))*m_H_positive + ((Z_negative)**(2))*m_Cl_negative]
I = 1./2*(((Z_positive)**(2))*m_H_positive + ((Z_negative)**(2))*m_Cl_negative);
#Using Debye-Huckel limiting law wee get,
# math.log(Y1) = -A*(abs(Z_positive*Z_negative))*(I**(1/2)))
Y = 10**(-A*(abs(Z_positive*Z_negative))*(I**(1./2)));
print "The mean activity coefficient at 25 C using Debye-Huckel limiting law is Y = %f"%(Y);
#Using Debye-Huckel extended model we get
#math.log(Y_prime) = (-A*(abs(Z_positive*Z_negative))*(I**(1/2)))/(1 + (I**(1/2)));
Y_prime = 10**((-A*(abs(Z_positive*Z_negative))*(I**(1./2)))/(1 + (I**(1./2))));
print "The mean activity coefficient at 25 C using Debye-Huckel extended model is Y = %f"%(Y_prime);
# Variables,
T = 25 + 273.15; #[K] - Temperature
mol_CaCl2 = 0.001; #[mol/kg] - Molality of HCl
A = 0.510; #[(kg/mol)**(1/2)]
Z_positive = 2; #Stoichiometric coefficient of 'Ca' ion
Z_negative = -1; #Stoichiometric coefficient of 'Cl' ion
m_Ca_positive = mol_CaCl2;
m_Cl_negative = 2*mol_CaCl2;
# Calculations
# I = 1/2*[((Z_positive)**(2))*m_Ca_positive + ((Z_negative)**(2))*m_Cl_negative]
I = 1./2*(((Z_positive)**(2))*m_Ca_positive + ((Z_negative)**(2))*m_Cl_negative);
#Using Debye-Huckel limiting law wee get,
# math.log(Y1) = -A*(abs(Z_positive*Z_negative))*(I**(1/2)))
Y = 10**(-A*(abs(Z_positive*Z_negative))*(I**(1./2)));
# Results
print "The mean activity coefficient at 25 C using Debye-Huckel limiting law is Y = %f"%(Y);
# Variables,
T = 50 + 273.15; #[K] - Temperature
R=8.314; #[J/mol*K] - Universal gas constant
x_1 = 0.3; # Mole fraction of component 1
x_2 = (1-x_1); # Mole fraction of component 2
#Increment of 1% means Y2 = 1.01*Y1
# Calculations
#Excess volume of the mixture is given by,
V_excess = 4*x_1*x_2; #[cm**(3)/mol]
#Amd threfore
V_1_excess = 4*x_2*x_2*10**(-6); #[m**(3)/mol] - Exces volume of component 1
V_2_excess = 4*x_1*x_1*10**(-6); #[m**(3)/mol] - Exces volume of component 2
#We have from equation 14.89 of the book,
#V_i_excess/(R*T) = (del_math.log(Y_i)/del_P)_T,x
#Rearranging above equation
#d(math.log(Y1)) = (V1_excess/(R*T))dP
#Integrating the above equation at constant 'T' and 'x' in the limit from 'Y1' to '1.01*Y1' in the LHS and from 'P' to 'P+delta_P' in the RHS
#On simplification we get
#math.log(1.01*Y1/Y1) = (V_1_exces/(R*T))*delta_P
delta_P = math.log(1.01)/(V_1_excess/(R*T)); #[N/m**(2)]
delta_P = delta_P*10**(-6); #[MPa]
# Results
print "The required pressure to increase the activity coefficient by 1%% is %f MPa"%(delta_P);
# Variables
T = 293.15; #[K] - Temperature
R=8.314; #[J/mol*K] - Universal gas constant
A = 1280; #[J/mol]
#(dA/dT)_P,x = del_A (say)
dA_dT = -7.0; #[J/mol-K]
#For equilomar mixture,
x_1 = 0.5; # Mole fraction of component 1
x_2 = 0.5; # Mole fraction of component 2
# Calculations
#math.log(Y1) = (A/(R*T))*x_2**(2)
#math.log(Y2) = (A/(R*T))*x_1**(2)
Y1 = math.exp((A/(R*T))*x_2**(2));
Y2 = math.exp((A/(R*T))*x_1**(2));
#G_excess/(R*T*) = x_1*math.log(Y1) + x_2*math.log(Y2) = (A/(R*T))*x_1*x_2
G_excess = A*x_1*x_2; #[J/mol] - Excess Gibbs free energy
#H_excess/(R*T**(2)) = -[d/dT(G_excess/(R*T))]_P,x
#H_excess/(R*T**(2)) = -((x_1*x_2)/R)*[d/dT(A/T)]_P,x
#On simplification and putting the values we get
H_excess = A*x_1*x_2 - T*dA_dT*x_1*x_2; #[J/mol] - Excess enthalpy
#Now excess entropy is given by
S_excess = (H_excess - G_excess)/T; #[J/mol-K] - Excess entropy
# Results
print "For equimolar mixture";
print "Excess Gibbs free energy G_excess) is %f J/mol"%(G_excess);
print "Excess enthalpy H_excess) is %f J/mol"%(H_excess);
print "Excess entropy S_excess) is %f J/mol"%(S_excess);
print "The value of activity coefficient Y1 is %f"%(Y1);
print "The value of activity coefficient Y2 is %f"%(Y2);
# Given
T = 60 + 273.15; #[K] - Temperature
R = 8.314; #[J/mol*K] - Universal gas constant
# math.log(Y1_inf) = math.log(Y2_inf) = 0.15 + 10/T
# From van L# Variablesr equation
# A = math.log(Y1_inf) and B = math.log(Y2_inf) and since it is given that math.log(Y1_inf) = math.log(Y2_inf), therefore A = B
#(x_1*x_2)/(G_excess/R*T) = x_1/B + x_2/A = X_1/A + x_2/A = 1/A
# G_excess/(R*T) = A*x_1*x_2
# For equilomar mixture,
x_1 = 0.5; # Mole fraction of component 1
x_2 = 0.5; # Mole fraction of component 2
# Calculations and Results
# Expression for A can be written as
# A = 0.15 + 10/T, where T is in C. Therefore
A = 0.15 + 10/(T - 273.15);
# Differentiating it with respect to temprature we get
dA_dT = - 10/((T-273.15)**(2));
# The excess Gibbs free energy can be calculated as
G_excess = A*x_1*x_2*(R*T); #[J/mol]
# The ideal Gibbs free energy change can be calculated as
delta_G_id_mix = R*T*(x_1*math.log(x_1) + x_2*math.log(x_2)); #[J/mol]
# Finally we have,
delta_G_mix = G_excess + delta_G_id_mix; #[J/mol]
print "The Gibbs free energy change of mixing for equimolar mixture is %f J/mol"%(delta_G_mix);
# Now let us determine the excess enthalpy. We know that
# H_excess/(R*T**(2)) = -[d/dT(G_excess/R*T)]_P,x = - x_1*x_2*[d/dT(A)]_P,x
# Therefore at 'T' = 60 C the excess enthalpy is given by
H_excess = -R*(T**(2))*x_1*x_2*dA_dT; #[J/mol]
delta_H_id_mix = 0; #[J/mol] - Enthalpy change of mixing for ideal solution is zero.
#Thus enthalpy change of mixing for an equimolar mixture at 333.15 K is given by
delta_H_mix = delta_H_id_mix + H_excess; #[J/mol]
print "The enthalpy change of mixing for equimolar mixture is %f J/mol"%(delta_H_mix);