# Chapter 11 : Solution Thermodynamics Theroy¶

## Example 11.4 page no : 168¶

In [4]:
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  =  ')

(poly1d([      1,    -400,  -10800, 1008000,       0]), 'H1  =  ')
(poly1d([    1,  -640, 24000,     0,     0]), 'H2  =  ')
('J/mol', 420, 'H1_inf  =  ')
('J/mol', 640, 'H2_inf  =  ')


## Example 11.5 page no : 170¶

In [10]:
%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')

Out[10]:
<matplotlib.text.Text at 0x10f99d990>

## Example 11.7 page no : 171¶

In [11]:
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

Fugacity Coefficients are  0.8324 0.9511
Second Viral coefficient is  -72.14
Compressibility Factor is  0.87


## Example 11.8 page no : 176¶

In [12]:
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'

Fugacity coefficient is  0.638
fugacity is  44.7 bar


## Example 11.9 page no : 179¶

In [15]:
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

   ij     Tcij     Pcij   Vcij     Zcij      Wij     Trij      B0       B1       Bij
[array([11, 22, 12]), array([ 535.5,  591.8,  563. ]), array([ 41.5,  41.1,  41.3]), array([267, 316, 291]), array([ 0.249,  0.264,  0.257]), array([ 0.323,  0.262,  0.293]), array([ 0.603,  0.546,  0.574]), array([-0.865, -1.028, -0.943]), array([-1.3  , -2.045, -1.632]), array([-1378., -1872., -1611.])]
Fugacity Coefficients are  0.985 0.987


## Example 11.10 page no : 181¶

In [16]:
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'

G_E  =   344.4 J/mol
S_E  =   1.492 J/mol/K
H_E  =   826.4 J/mol
`