from numpy import poly1d,roots
#H = 400x1 + 600x2 + x1x2(40x1 + 20x2) Given
#Substituting x2 = 1-x1
H = poly1d([600, -180, 0, -20],'x1','c') #(A)
#K = dH/dx1
K = poly1d([-180, 0, -60],'x1','c')
# Calculations
#Using Eqn 11.15 H1 = H+x2*K
#substituting x2 = 1-x1
H1 = poly1d([420, 0, -60, 40],'x1','c') #(B)
#Similarly for H2
H2 = poly1d([600 ,0 ,0, 40],'x1','c') #(C)
#Now to calculate H1_inf and H2_inf
#x1 = 0 in (B)
H1_inf = 420; #[J/mol]
#x2 = 0 in (C) i.e. x1 = 1
H2_inf = 640; #[J/mol]
# Results
print (H1,'H1 = ');
print (H2,'H2 = ');
print ('J/mol',H1_inf,'H1_inf = ')
print ('J/mol',H2_inf,'H2_inf = ')
%matplotlib inline
from matplotlib.pyplot import plot,suptitle,xlabel,ylabel,subplot
from numpy import array,exp,round,linspace
import math
# Variables
#Using Eqn 11.30
#G = t(T)+RT ln f
#G` = t(T)+rt ln f`
#Hence ln(f/f`) = (1/RT)*(G-G`)
#G = H-TS
#G` = H-TS`
#ln (f/f`) = (n/R)*(((H-H`)/T)-(S-S`)) (A)
R = 8.314;
n = 18.015;
f_ = 1; #[kPa]
P_ = f_;
H_ = 3076.8; #[J/g]
S_ = 10.3450; #[J/g/K]
T = 573.15; #[K]
P = array([10, 50, 100, 200, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000, 6500, 7000, \
7600, 8000, 8400,0]);
H = array([3076.6, 3075.7, 3074.2, 3072.1, 3064.8, 3052.1, 3038.9, 3025 ,3010.4, 2995.1, 2979 ,2962, \
2944.2, 2925.5, 2905.8, 2885 ,2862.9, 2839.4, 2808.8, 2786.8, 2763.1,0]);
S = array([9.282, 8.5380, 8.2166, 7.8937, 7.4614, 7.1251, 6.9207, 6.7696, 6.6470, 6.5422, \
6.4491, 6.3642, 6.2852, 6.2105, 6.1388, 6.0692, 6.0008, 5.9327, 5.8503, 5.7942, 5.7366,0]);
# Calculations and Results
K = round(exp((n/R)*(((H-H_)/T)-(S-S_))),0);
f = K*f_;
P[21] = 8592.7; #[kPa]
P_sat = P[21]
f[21] = 6738.9;
f_sat = f[21]
si = round(f/P,4);
si_sat = si[21]
Vl = round(1.403*n,2) #[cm**3/mol]
#Umath.sing Eqn 11.41
P_new = linspace(8592.7,10000,10);
f_new = round(f_sat*exp(Vl*(P_new-P_sat)/(R*1000*T)),1);
si_new = f_new/P_new;
subplot(211)
plot(P/1000,f/1000,'b')
plot(P_new/1000,f_new/1000,'g')
dotsx = [0 ,P_sat/1000];
dotsy = [f_sat/1000, f_sat/1000];
plot(dotsx,dotsy,'b--')
dotsx = [0, f_sat/1000];
dotsy = [P_sat/1000 ,P_sat/1000];
plot(dotsy,dotsx,'g--')
dotsx = [11,8];
dotsx = [6,6];
plot(dotsx,dotsy,'w')
suptitle('f vs P')
xlabel('P X 10**-3 /kPa')
ylabel('f X 10**-3 /kPa')
subplot(222)
plot(P/1000,si,'r')
plot(P_new/1000,si_new,'g')
dotsx = [0, P_sat/1000];
dotsy = [si_sat ,si_sat];
plot(dotsx,dotsy,'g--')
dotsx = [0.55, si_sat];
dotsy = [P_sat/1000, P_sat/1000];
plot(dotsy,dotsx,'r--')
suptitle('si vs P')
xlabel('P X 10**-3 /kPa')
ylabel('si')
import math
# Variables
T = 200.; #[K]
P = 30.; #[bar]
R = 83.14;
x1 = 0.4; #[N2]
x2 = 1-x1; #[CH4]
B11 = -35.2; #[cm**3/mol]
B22 = -105.; #[cm**3/mol]
B12 = -59.8; #[cm**3/mol]
# Calculations
delta_12 = round((2*B12)-B11-B22,1);
si_1 = round(math.exp((P/(R*T))*(B11+(x2**2*delta_12))),4);
si_2 = round(math.exp((P/(R*T))*(B22+(x1**2*delta_12))),4);
B = round((x1**2*B11)+(2*x1*x2*B12)+(x2**2*B22),2);
Z = round(1+((B*P)/(R*T)),2);
# Results
print 'Fugacity Coefficients are ',si_2,si_1
print 'Second Viral coefficient is ',B
print 'Compressibility Factor is ',Z
import math
# Variables
T = 473.15; #[K]
P = 70.; #[bar]
Tc = 420.; #[K]
Pc = 40.43; #[bar]
omega = 0.191;
# Calculations
#By interpolation in Tables E.15 and E.16
si_0 = 0.627;
si_1 = 1.096;
#Using Eqn(11.64)
si = round(si_0*(si_1**omega),3);
f = round(si*P,1);
# Results
print 'Fugacity coefficient is ',si
print 'fugacity is ',f,'bar'
from numpy import array,exp,round
import math
# Variables
P = 25.; #[KPa]
T = 323.15; #[K]
R = 83.14;
x1 = 0.5;
x2 = 1-x1;
# Calculations
ij = array([11,22,12]);
Tc_ij = array([535.5,591.8,563.0]);
Pc_ij = array([41.5,41.1,41.3]);
Vc_ij = array([267,316,291]);
Zc_ij = round((Pc_ij*Vc_ij)/(R*Tc_ij),3);
omega_ij = array([0.323,0.262,0.293]);
Tr_ij = round(T/Tc_ij,3);
B0 = round(0.083-(0.422/(Tr_ij**1.6)),3)
B1 = round(0.139-(0.172/(Tr_ij**4.2)),3)
B_ij = round((R*Tc_ij/Pc_ij)*(B0+(omega_ij*B1)),0);
delta_12 = (2*B_ij[2])-B_ij[0]-B_ij[1];
R = 8314;
si_1 = round(exp((P/(R*T))*(B_ij[0]+(x2**2*delta_12))),3);
si_2 = round(exp((P/(R*T))*(B_ij[2]+(x1**2*delta_12))),3);
Ans = [ij,Tc_ij,Pc_ij,Vc_ij,Zc_ij,omega_ij,Tr_ij,B0,B1,B_ij];
# Results
print ' ij Tcij Pcij Vcij Zcij Wij Trij B0 B1 Bij'
print Ans
print 'Fugacity Coefficients are ',si_2,si_1
import math
# Variables
T0 = 298.15; #[K]
T = 323.15; #[K]
Cp_E = -2.86; #[J/mol/K]
Ho_E = 897.9; #[J/mol]
Go_E = 384.5; #[J/mol]
# Calculations
#(a) Derivations
#G_E = -a*(T ln T - T)+ bT + c
#S_E = a ln T - b
#H_E = aT + c
#Where
#a = Cp_E
#c = Ho_E-aT0
#b = ((Go_E+a(T ln T0 - T0)-c)/T0)
#(b)
a = Cp_E;
c = round(Ho_E-(a*T0),1);
b = round((Go_E+(a*((T0*math.log(T0))-T0)-c))/T0,4);
G_E = round((-a*(T*math.log(T)-T))+(b*T)+c,1);
S_E = round((a*math.log(T))-b,3);
H_E = round((a*T)+c,1);
# Results
print 'G_E = ',G_E,'J/mol'
print 'S_E = ',S_E,'J/mol/K'
print 'H_E = ',H_E,'J/mol'