import math
# Variables
T = 310.93; #[K] - Temperature
P = 2.76*10**(6); #[Pa] - Pressure
y1 = 0.8942; #[mol] - mole fraction of component 1
y2 = 1 - y1; #[mol] - mole fraction of component 2
R=8.314; #[J/mol*K] - Universal gas constant
#For component 1 (methane)
Tc_1 = 190.6; #[K] - Critical temperature
Pc_1 = 45.99*10**(5); #[N/m**(2)] - Critical pressure
Vc_1 = 98.6; #[cm**(3)/mol] - Critical molar volume
Zc_1 = 0.286; # - Critical compressibility factor
w_1 = 0.012; # - Critical acentric factor
#Similarly for component 2 (n-Butane)
Tc_2 = 425.1; #[K]
Pc_2 = 37.96*10**(5); #[N/m**(2)]
Vc_2 = 255; #[cm**(3)/mol]
Zc_2=0.274;
w_2=0.2;
# Calculations
#For component 1
Tr_1 = T/Tc_1; #Reduced temperature
#At reduced temperature
B1_0 = 0.083-(0.422/(Tr_1)**(1.6));
B1_1 = 0.139-(0.172/(Tr_1)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0+(w*B_1)
B_11 = ((B1_0+(w_1*B1_1))*(R*Tc_1))/Pc_1; #[m**(3)/mol]
#Similarly for component 2
Tr_2 = T/Tc_2; #Reduced temperature
#At reduced temperature Tr_2,
B2_0 = 0.083-(0.422/(Tr_2)**(1.6));
B2_1 = 0.139-(0.172/(Tr_2)**(4.2));
B_22 = ((B2_0+(w_2*B2_1))*(R*Tc_2))/Pc_2; #[m**(3)/mol]
#For cross coeffcient
Tc_12 = (Tc_1*Tc_2)**(1./2); #[K]
w_12 = (w_1+w_2)/2;
Zc_12 = (Zc_1+Zc_2)/2;
Vc_12 = (((Vc_1)**(1./3)+(Vc_2)**(1./3))/2)**(3); #[cm**(3)/mol]
Vc_12 = Vc_12*10**(-6); #[m**(3)/mol]
Pc_12 = (Zc_12*R*Tc_12)/Vc_12; #[N/m**(2)]
# Variables, Z = 1 + (B*P)/(R*T)
#Now we have,(B_12*Pc_12)/(R*Tc_12) = B_0 + (w_12*B_1)
#where B_0 and B_1 are to be evaluated at Tr_12
Tr_12 = T/Tc_12;
#At reduced temperature Tr_12
B_0 = 0.083-(0.422/(Tr_12)**(1.6));
B_1 = 0.139-(0.172/(Tr_12)**(4.2));
B_12 = ((B_0+(w_12*B_1))*R*Tc_12)/Pc_12; #[m**(3)/mol]
#For the mixture
B = y1**(2)*B_11+2*y1*y2*B_12+y2**(2)*B_22; #[m**(3)/mol]
#Now del_12 can be calculated as,
del_12 = 2*B_12 - B_11 - B_22; #[m**(3)/mol]
#We have the relation, math.log(phi_1) = (P/(R*T))*(B_11 + y2**(2)*del_12), therefore
phi_1 = math.exp((P/(R*T))*(B_11 + y2**(2)*del_12));
#Similarly for component 2
phi_2 = math.exp((P/(R*T))*(B_22 + y1**(2)*del_12));
# Results
print " The value of fugacity coefficient of component 1 phi_1) is %f"%(phi_1);
print " The value of fugacity coefficient of component 2 phi_2) is %f"%(phi_2);
#Finally fugacity coefficient of the mixture is given by
#math.log(phi) = y1*math.log(phi_1) + y2*math.log(phi_2);
phi = math.exp(y1*math.log(phi_1) + y2*math.log(phi_2));
print " The value of fugacity coefficient of the mixture phi) is %f "%(phi);
#The fugacity coefficient of the mixture can also be obtained using
#math.log(phi) = (B*P)/(R*T)
from scipy.optimize import fsolve
import math
# Variables
T = 460.; #[K] - Temperature
P = 40.*10**(5); #[Pa] - Pressure
R=8.314; #[J/mol*K] - Universal gas constant
# component 1 = nitrogen
# component 2 = n-Butane
y1 = 0.4974; # Mole percent of nitrogen
y2 = 0.5026; # Mole percent of n-Butane
Tc_nit = 126.2; #[K]
Pc_nit = 34.00*10**(5); #[Pa]
Tc_but = 425.1; #[K]
Pc_but = 37.96*10**(5); #[Pa]
# (1). van der Walls equation of state
# The fugacity coefficient of component 1 in a binary mixture following van der Walls equation of state is given by,
# math.log(phi_1) = b_1/(V-b) - math.log(Z-B) -2*(y1*a_11 + y2*a_12)/(R*T*V)
# and for component 2 is given by,
# math.log(phi_2) = b_2/(V-b) - math.log(Z-B) -2*(y1*a_12 + y2*a_22)/(R*T*V)
# Where B = (P*b)/(R*T)
# Calculations
# For componenet 1 (nitrogen)
a_1 = (27*R**(2)*Tc_nit**(2))/(64*Pc_nit); #[Pa-m**(6)/mol**(2)]
b_1 = (R*Tc_nit)/(8*Pc_nit); #[m**(3)/mol]
# Similarly for componenet 2 (n-Butane)
a_2 = (27*R**(2)*Tc_but**(2))/(64*Pc_but); #[Pa-m**(6)/mol**(2)]
b_2 = (R*Tc_but)/(8*Pc_but); #[m**(3)/mol]
# Here
a_11 = a_1;
a_22 = a_2;
# For cross coefficient
a_12 = (a_1*a_2)**(1./2); #[Pa-m**(6)/mol**(2)]
# For the mixture
a = y1**(2)*a_11 + y2**(2)*a_22 + 2*y1*y2*a_12; #[Pa-m**(6)/mol**(2)]
b = y1*b_1 + y2*b_2; #[m**(3)/mol]
# The cubic form of the van der Walls equation of state is given by,
# V**(3) - (b+(R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0
# Substituting the value and solving for V, we get
# Solving the cubic equation
def f(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V_1=fsolve(f,-1)
V_2=fsolve(f,0)
V_3=fsolve(f,1)
# The molar volume V = V_3, the other two roots are imaginary
V = V_3; #[m**(3)/mol]
# The comprssibility factor of the mixture is
Z = (P*V)/(R*T);
# And B can also be calculated as
B = (P*b)/(R*T);
# The fugacity coefficient of component 1 in the mixture is
phi_1 = math.exp(b_1/(V-b) - math.log(Z-B) -2*(y1*a_11 + y2*a_12)/(R*T*V));
# Similarly fugacity coefficient of component 2 in the mixture is
phi_2 = math.exp(b_2/(V-b) - math.log(Z-B) -2*(y1*a_12 + y2*a_22)/(R*T*V));
# The fugacity coefficient of the mixture is given by,
# math.log(phi) = y1*math.log(phi_1) + y2*math.log(phi_2)
phi = math.exp(y1*math.log(phi_1) + y2*math.log(phi_2));
# Also the fugacity coefficient of the mixture following van der Walls equation of state is given by,
# math.log(phi) = b/(V-b) - math.log(Z-B) -2*a/(R*T*V)
phi_dash = math.exp(b/(V-b) - math.log(Z-B) -2*a/(R*T*V));
# The result is same as obtained above
# Results
print " 1van der Walls equation of state";
print " The value of fugacity coefficient of component 1 nitrogen) is %f"%(phi_1);
print " The value of fugacity coefficient of component 2 n-butane) is %f"%(phi_2);
print " The value of fugacity coefficient of the mixture is %f"%(phi);
print " Also the fugacity coefficient of the mixture from van der Walls equation of state is %f which is same as above)"%(phi_dash);
# (2). Redlich-Kwong equation of state
# For component 1,
a_1_prime = (0.42748*R**(2)*Tc_nit**(2.5))/Pc_nit; #[Pa-m**(6)/mol**(2)]
b_1_prime = (0.08664*R*Tc_nit)/Pc_nit; #[m**(3)/mol]
#similarly for component 2,
a_2_prime = (0.42748*R**(2)*Tc_but**(2.5))/Pc_but; #[Pa-m**(6)/mol**(2)]
b_2_prime = (0.08664*R*Tc_but)/Pc_but; #[m**(3)/mol]
# For cross coefficient
a_12_prime = (a_1_prime*a_2_prime)**(1./2); #[Pa-m**(6)/mol**(2)]
# For the mixture
a_prime = y1**(2)*a_1_prime + y2**(2)*a_2_prime +2*y1*y2*a_12_prime; #[Pa-m**(6)/mol**(2)]
b_prime = y1*b_1_prime +y2*b_2_prime; #[m**(3)/mol]
#The cubic form of Redlich Kwong equation of state is given by,
# V**(3)-((R*T)/P)*V**(2)-((b**(2))+((b*R*T)/P)-(a/(T**(1/2)*P))*V-(a*b)/(T**(1/2)*P)=0
# Solving the cubic equation
def f1(V):
return V**(3)-((R*T)/P)*V**(2)-((b_prime**(2))+((b_prime*R*T)/P)-(a_prime/(T**(1./2)*P)))*V-(a_prime*b_prime)/(T**(1./2)*P)
V_4=fsolve(f1,1)
V_5=fsolve(f1,0)
V_6=fsolve(f1,-1)
# The molar volume V = V_4, the other two roots are imaginary
V_prime = V_4; #[m**(3)/mol]
# The compressibility factor of the mixture is
Z_prime = (P*V_prime)/(R*T);
# And B can also be calculated as
B_prime = (P*b_prime)/(R*T);
# The fugacity coefficient of component 1 in the binary mixture is given by
# math.log(phi_1) = b_1/(V-b) - math.log(Z-B) + ((a*b_1)/((b**(2)*R*T**(3/2))))*(math.log((V+b)/V)-(b/(V+b)))-(2*(y1*a_1+y2*a_12)/(R*T**(3/2)b))*(math.log(V+b)/b)
phi_1_prime = math.exp((b_1_prime/(V_prime-b_prime))-math.log(Z_prime-B_prime)+((a_prime*b_1_prime)/((b_prime**(2))*R*(T**(3./2))))*(math.log((V_prime+b_prime)/V_prime)-(b_prime/(V_prime+b_prime)))-(2*(y1*a_1_prime+y2*a_12_prime)/(R*(T**(3./2))*b_prime))*(math.log((V_prime+b_prime)/V_prime)));
# Similarly fugacity coefficient of component 2 in the mixture is
phi_2_prime = math.exp((b_2_prime/(V_prime-b_prime))-math.log(Z_prime-B_prime)+((a_prime*b_2_prime)/((b_prime**(2))*R*(T**(3./2))))*(math.log((V_prime+b_prime)/V_prime)-(b_prime/(V_prime+b_prime)))-(2*(y1*a_12_prime+y2*a_2_prime)/(R*(T**(3./2))*b_prime))*(math.log((V_prime+b_prime)/V_prime)));
# The fugacity coefficient of the mixture is given by,
# math.log(phi) = y1*math.log(phi_1) + y2*math.log(phi_2)
phi_prime = math.exp(y1*math.log(phi_1_prime) + y2*math.log(phi_2_prime));
# Also the fugacity coefficient for the mixture following Redlich-kwong equation of state is also given by
# math.log(phi) = b/(V-b) - math.log(Z-B) - (a/(R*T**(3/2)))*(1/(V+b)+(1/b)*math.log((V+b)/V))
phi_prime_dash = math.exp(b_prime/(V_prime-b_prime) - math.log(Z_prime-B_prime) - (a_prime/(R*T**(3./2)))*(1/(V_prime+b_prime)+(1./b_prime)*math.log((V_prime+b_prime)/ V_prime)));
print " \nRedlich-Kwong equation of state";
print " The value of fugacity coefficient of component 1 nitrogen) is %f"%(phi_1_prime);
print " The value of fugacity coefficient of component 2 n-butane) is %f"%(phi_2_prime);
print " The value of fugacity coefficient of the mixture is %f"%(phi_prime);
print " Also the fugacity coefficient for the mixture from Redlich-kwong equation of state is %f which is same as above)"%(phi_prime_dash);