# Chapter 14 : Activity Coefficients Models for Liquid Mixtures¶

### Example 14.4 Page Number : 461¶

In :

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);

The value of Gibbs free energy change for equimolar mixture formation is -388.185034 cal/mol
The value of enthalpy change for equimolar mixture formation is 25.000000 cal/mol
The least amount of work required for separation at 300 K is -388.185034 cal/mol


### Example 14.10 Page Number : 466¶

In :

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."

 The value of van Laar parameters are A = 2.271326 and B = 1.781420
The value of Margules parameters are A = 2.293274 and B = 1.746000
Because of the higher value of correlation factor for Van Laar model it fits the data better.


### Example 14.11 Page Number : 470¶

In :

# 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);

The value of activity coefficients for Wilson equation are Y1 = 2.172258 	 and 	 Y2 = 1.254121
The value of activity coefficients for NRTL equation are Y1 = 2.143030 	 and 	 Y2 = 1.262506


### Example 14.12 Page Number : 474¶

In :

# 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);

 The activity coefficients are Y1 = 4.992034  and  Y2 = 1.005260


### Example 14.13 Page Number : 481¶

In :

# 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);

The value of activity coefficients for Scatchard-Hilderbrand model are Y1 = 1.010245 	 and 	 Y2 = 1.017658


### Example 14.14 Page Number : 485¶

In :

# 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);

The mean activity coefficient at 25 C using Debye-Huckel limiting law is Y = 0.963546
The mean activity coefficient at 25 C using Debye-Huckel extended model is Y = 0.964643


### Example 14.15 Page Number : 485¶

In :

# 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);

The mean activity coefficient at 25 C using Debye-Huckel limiting law is Y = 0.879290


### Example 14.17 Page Number : 488¶

In :

# 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);

The required pressure to increase the activity coefficient by 1% is 13.639411 MPa


### Example 14.21 Page Number : 490¶

In :

# 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);

For equimolar mixture
Excess Gibbs free energy G_excess) is 320.000000 J/mol
Excess enthalpy H_excess) is 833.012500 J/mol
Excess entropy S_excess) is 1.750000 J/mol
The value of activity coefficient Y1 is 1.140305
The value of activity coefficient Y2 is 1.140305


### Example 14.22 Page Number : 491¶

In :

# 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);

The Gibbs free energy change of mixing for equimolar mixture is -1700.608815 J/mol
The enthalpy change of mixing for equimolar mixture is 640.806876 J/mol