import math
# Variables
T = 40 + 273.15; #[C] - Temperature
P_1 = 0; #[bar]
P_2 = 10; #[bar]
V_liq = 90.45; #[cm**(3)/mol]
V_liq = V_liq*10**(-6); #[m**(3)/mol]
P_sat = 4.287; #[bar]
# For butadiene
T_c = 425.0; #[K] - Critical temperature
P_c = 43.3; #[bar] - Critical pressure
P_c = P_c*10**(5); #[N/m**(2)]
w = 0.195; # Acentric factor
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations and Results
# Let us calculate second virial coefficient at 40 C
Tr = T/T_c; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0+(w*B_1)
B = ((B_0 + (w*B_1))*(R*T_c))/P_c; #[m**(3)/mol] - Second virial coefficient
# math.log(f/P) = (B*P)/(R*T)
# f = P*exp((B*P)/(R*T))
print " The table is as follows"
print " Pbar \t\t fbar \t\t phi";
from numpy import zeros
P = [1,2,3,4,4.287,5,6,8,10];
f = zeros(9);
phi = zeros(9);
for i in range(5):
f[i]=P[i]*(math.exp((B*P[i]*10**(5))/(R*T))); #[bar] # Pressure inside the exponential term has to be in N/m**(2)
phi[i]= (f[i]/P[i]);
print " %f \t %f \t\t\t %f"%(P[i],f[i],phi[i])
f_sat = f[4];
# From pressure of 4.287 bar onwards the valid equation to compute fugacity of compressed liquid is given by
# f = f_sat*exp[V_liq*(P-P_sat)/(R*T)]
for j in range(5,9):
f[j] = f_sat*math.exp((V_liq*(P[j]-P_sat)*10**(5))/(R*T)); #[bar] # Pressure inside the exponential term has to be in N/m**(2)
phi[j] = f[j]/P[j];
print " %f \t %f \t\t\t %f"%(P[j],f[j],phi[j]);
import math
from scipy.integrate import quad
# Variables
n = 100.; #[mol] - No of moles
T_1 = 600.; #[K] - Initial temperature
T_2 = 300.; #[K] - Final temperature
P_1 = 10.; #[atm] - Initial pressure
P_1 = P_1*101325; #[Pa]
P_2 = 5.; #[atm] - Final presssure
P_2 = P_2*101325; #[Pa]
Tc = 369.8; #[K] - Critical temperature
Pc = 42.48; #[bar] - Critical pressure
Pc = Pc*10**(5); #[Pa]
w = 0.152;
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations
# At 600 K
Tr = T_1/Tc; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0 + (w*B_1))*(R*Tc))/Pc; #[m**(3)/mol] - Second virial coefficient
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6); # (dB_0/dT)
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2); # (dB_1/dT)
dB_dT = ((R*Tc)/(Pc))*(dB0_dT + w*dB1_dT); # dB/dT
# Now let us calculate B and dB/dT at 300 K
Tr_prime = T_2/Tc; # Reduced temperature
B_0_prime = 0.083-(0.422/(Tr_prime)**(1.6));
B_1_prime = 0.139-(0.172/(Tr_prime)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B_prime = ((B_0_prime + (w*B_1_prime))*(R*Tc))/Pc; #[m**(3)/mol] - Second virial coefficient
dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_2**(-2.6); # (dB_0/dT)
dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_2**(-5.2); # (dB_1/dT)
dB_dT_prime = ((R*Tc)/(Pc))*(dB0_dT_prime + w*dB1_dT_prime); # dB/dT
# The change in enthalpy for ideal gas is given by
def f16(T):
return -0.966+7.279*10**(-2)*T-3.755*10**(-5)*T**(2)+7.58*10**(-9)*T**(3)
delta_H_ig = quad(f16,T_1,T_2)[0]
delta_H_ig = delta_H_ig*4.184; #[J/mol]
# We know that delta_H_ig = delta_U_ig + R*delta_T. Therefore change in internal energy is given by
delta_U_ig = delta_H_ig - R*(T_2 - T_1); #[J/mol]
# The change in entropy of ideal gas is given by
def f17(T):
return Cp_0/T
#delta_S_ig = quad(f17,T_1,T_2) - R*math.log(P_2/P_1)[0]
def f18(T):
return (-0.966+7.279*10**(-2)*T-3.755*10**(-5)*T**(2)+7.58*10**(-9)*T**(3))/T
delta_S_ig = quad(f18,T_1,T_2)[0]*4.184 - R*math.log(P_2/P_1)
# Now let us calculate the change in enthalpy of gas. We know that
# delta_H = delta_H_ig + delta_H_R
# delta_H_R = H_2_R - H_1_R
H_2_R = B_prime*P_2 - P_2*T_2*dB_dT_prime; # [J/mol]
H_1_R = B*P_1 - P_1*T_1*dB_dT; # [J/mol]
delta_H_R = H_2_R - H_1_R; # [J/mol]
delta_H = delta_H_ig + delta_H_R; # [J/mol]
# Let us calculate the residual entropy of gas
S_2_R = -P_2*dB_dT_prime; #[J/mol-K]
S_1_R = -P_1*dB_dT; #[J/mol-K]
delta_S = delta_S_ig + (S_2_R - S_1_R); #[J/mol-K]
# Let us calculate the residual internal energy of gas
U_2_R = -P_2*T_2*dB_dT_prime; #[J/mol-K]
U_1_R = -P_1*T_1*dB_dT; #[J/mol-K]
delta_U = delta_U_ig + (U_2_R - U_1_R); #[J/mol-K]
# For 100 mol sample,
delta_H_ig = delta_H_ig*n*10**(-3); #[kJ/mol]
delta_H = delta_H*n*10**(-3); #[kJ/mol]
delta_U_ig = delta_U_ig*n*10**(-3); #[kJ/mol]
delta_U = delta_U*n*10**(-3); #[kJ/mol]
delta_S_ig = delta_S_ig*n*10**(-3); #[kJ/mol]
delta_S = delta_S*n*10**(-3); #[kJ/mol]
# Results
print " The value of delta_H = %f kJ/mol"%(delta_H);
print " The value of delta_H_ig ideal gas)= %f J/mol"%(delta_H_ig);
print " The value of delta_U = %f kJ/mol"%(delta_U);
print " The value of delta_U_ig ideal gas) = %f kJ/mol"%(delta_U_ig);
print " The value of delta_S = %f kJ/mol"%(delta_S);
print " The value of delta_S_ig ideal gas) = %f kJ/mol"%(delta_S_ig);
# Variables
T = 35 + 273.15; #[K] - Temperature
P = 10; #[atm] - Pressure
P = P*101325; #[Pa]
# Methane obeys the equation of state
# Z = 1 + (P*B)/(R*T)
# Calculations
# At 35 C,
B = -50; #[cm**(3)/mol]
dB_dT = 1.0; #[cm**(3)/mol-K] - dB/dT
dB_dT = dB_dT*10**(-6); #[m**(3)/mol-K]
d2B_dT2 = -0.01; #[cm**(3)/mol-K**(2)] - d**2(B)/d(T**2)
d2B_dT2 = d2B_dT2*10**(-6); #[m**(3)/mol-K**(2)]
# Ideal gas molar heat capacity of methane is given by
# Cp_0 = 4.75 + 1.2*10**(-2)*T + 0.303*10**(-5)*T**(2) - 2.63*10**(-9)*T**(3)
# The molar heat capacity is given by
# Cp = Cp_0 + Cp_R
# For virial gas equation of state
Cp_R = -P*T*d2B_dT2; #[J/mol-K]
# thus heat capacity is given by
# Cp = a + b*T + c*T**(2) + d*T**(3) - P*T*d2B_dT2
# Putting the values, we get
Cp = (4.75 + 1.2*10**(-2)*T + 0.303*10**(-5)*T**(2) - 2.63*10**(-9)*T**(3))*4.184 - P*T*d2B_dT2; #[J/mol-K]
# Results
print " The molar heat capacity of methane is %f J/mol-K"%(Cp);
from scipy.optimize import fsolve
# Variables
T_1 = 360; #[K] - Initial temperature
P_1 = 10; #[bar] - Initial pressure
P_1 = P_1*10**(5); #[Pa]
Tc = 408.1; #[K] - Critical temperature
Pc = 36.48; #[bar] - Critical pressure
Pc = Pc*10**(5); #[Pa]
w = 0.181;
R = 8.314; #[J/mol*K] - Universal gas constant
Cv_0 = 106.0; #[J/mol-K]
# Calculations
# At 360 K
Tr = T_1/Tc; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0 + (w*B_1))*(R*Tc))/Pc; #[m**(3)/mol] - Second virial coefficient
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6); # (dB_0/dT)
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2); # (dB_1/dT)
dB_dT = ((R*Tc)/(Pc))*(dB0_dT + w*dB1_dT); # dB/dT
# At state 1
V_1 = B + (R*T_1)/P_1; #[m**(3)/mol] - Molar volume
# At state 1
V_2 = 10*V_1; #[m**(3)/mol] - Molar volume
T_2 = -(P_1*T_1*(dB_dT))/Cv_0 + T_1; #[K]
T = 350;
fault = 10;
def f(T):
return 106*(T-T_1)+972.72-((R*T**(2))/(V_2-B_prime))*dB_dT_prime
while(fault>0.007):
Tr_prime = T/Tc; # Reduced temperature
B_0_prime = 0.083-(0.422/(Tr_prime)**(1.6));
B_1_prime = 0.139-(0.172/(Tr_prime)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B_prime = ((B_0_prime + (w*B_1_prime))*(R*Tc))/Pc; #[m**(3)/mol] - Second virial coefficient
dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_2**(-2.6); # (dB_0/dT)
dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_2**(-5.2); # (dB_1/dT)
dB_dT_prime = ((R*Tc)/(Pc))*(dB0_dT_prime + w*dB1_dT_prime); # dB/dT
T_prime = fsolve(f,0.15)
fault=abs(T-T_prime);
T = T + 0.001;
# Results
print " The final temperature is %f K"%(T);
# Variables
T = 220 + 273.15; #[K] - Temperature
Tc = 562.2; #[K] - Critical temperature
Pc = 48.98; #[bar] - Critical pressure
Pc = Pc*10**(5); #[Pa]
w = 0.210;
R = 8.314; #[J/mol*K] - Universal gas constant
P_sat = 1912.86; #[kPa] - Saturation pressure at 220 C
P_sat = P_sat*10**(3); #[Pa]
Mol_wt = 78.114; #[g/mol] - Molecular weight of benzene
# Calculations and Results
#(1)
# At 220 C
Tr = T/Tc; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0 + (w*B_1))*(R*Tc))/Pc; #[m**(3)/mol] - Second virial coefficient
# We know that math.log(f/P) = (B*P)/(R*T)
# Thus at saturated conditions
# math.log(f_sat/P_sat) = B*P_sat/(R*T)
f_sat = P_sat*(math.exp((B*P_sat)/(R*T))); #[Pa]
f_sat = f_sat*10**(-3); #[kPa]
print " 1).The fugacity of liquid benzene is %f kPa"%(f_sat);
#(2)
P = 2014.7; # [psia] - Total gauge pressure
P = 138.94; # [bar]
P = P*10**(5); # [Pa]
den = 0.63; # [g/cm**(3)] - density of benzene
den = den*10**(3); # [kg/m**(3)]
# Therefore specific volume is
V = 1/den; #[m/**(3)/kg]
# Molar volume is given by
V = V*Mol_wt*10**(-3); #[m**(3)/mol]
# Thus fugacity at 220 C and pressure P is given by
f = f_sat*(math.exp((V*(P-P_sat))/(R*T)));
print " 2).The fugacity of liquid benzene is %f kPa"%(f);
# Variables
# C = -0.067 + 30.7/T
# D = 0.0012 - 0.416/T
T = 80 + 273.15; #[K]
P = 30; #[bar]
#P = P; #[N/m**(2)]
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations
# We have the relation derived in the book
# d(G/(R*T)) = (V/(R*T))dP - (H/(R*T**(2)))dT
# Writing the same equation for ideal gas and subtracting it from the above equation we get
# d(G_R/(R*T)) = (V_R/(R*T))dP - (H_R/(R*T**(2)))dT
# Therefore, H_R/(R*T**(2)) = -[del((G_R)/(R*T))/del(T)]_P
# Substituting the relation G_R/(R*T) = math.log(f/P), we get
# H_R/(R*T**(2)) = -[del(math.log(f/P))/del(T)]_P = -[del(-C*P - D*P**(2))/del(T)]_P
# H_R/R = - 30.7*P + 0.416*P**(2)
# Substituting the given conditions we get
H_R = R*(-30.7*P + 0.416*P**(2)); #[J/mol]
# Results
print " The molar enthalpy of the gas relative to that of the ideal gas at 80 C and 30 bar pressure is H_R = %f J/mol"%(H_R);
# Variables
# (1)
T = 311; #[K] - Temperature
R = 8.314; #[J/mol*K] - Universal gas constant
# Pressure in 'bar' is given below
P = [0.690,1.380,2.760,5.520,8.280,11.034,13.800,16.550];
# Molar volume in 'm**(3)/mol' is given below
V = [0.0373,0.0186,0.00923,0.00455,0.00298,0.00220,0.00175,0.00144];
# Z = 1 + (B/V) + (C/V**(2))
# (Z-1)*V = B + (C/V)
from numpy import zeros
Z=zeros(8);
k=zeros(8);
t=zeros(8);
# Calculations and Results
for i in range(8):
Z[i]=(P[i]*10**(5)*V[i])/(R*T);
k[i]=(Z[i]-1)*V[i];
t[i]=1./V[i];
from scipy.stats import linregress
from numpy import array
t = array(t)
k = array(k)
C,B,c,d,e = linregress(t.T,k.T)
#[C,B,sig]=reglin(t',k');
#From the regression, we get intercept = B and slope = C,and thus,
print " 1).The second virial coefficient of CO2 is given by B = %e m**3)/mol"%(B);
print " The thied virial coefficient of CO2 is given by C = %e m**6)/mol**2)"%(C);
# (2)
P_final = 13.8; #[bar]
def f40(P):
return V-(R*T)/P
# We know that R*T*math.log(f/P) = quad(f40,0,P)[0]
# Therefore we have to plot V - (R*T)/P versus P and calculate the area beneath the curve from 0 to 13.8 bar
# For this we need the value of the term V - (R*T)/P at P = 0. At low pressure the virial equation becomes
# Z = 1 + (B/V)
#and V - (R*T)/P = (Z*R*T)/P - (R*T)/P = (1 + (B/V))*((R*T)/P) - (R*T)/P = (B*R*T)/(P*V) = (B/Z)
# Thus lim P tending to zero (V - (R*T)/P) = B ( as P tend to zero, Z tend to 1 )
P_prime = [0.000,0.690,1.380,2.760,5.520,8.280,11.034,13.800];
V_prime = [0.000,0.0373,0.0186,0.00923,0.00455,0.00298,0.00220,0.00175];
summation = 0;
x=zeros(8);
y=zeros(8);
z=zeros(8);
for j in range(2,8):
x[j]=V_prime[j]-(R*T)/(P_prime[j]*10**(5)); #[m**(3)/mol]
y[j]=(x[j] + x[j-1])/2;
z[j]=y[j]*((P_prime[j]-P_prime[j-1]))*10**(5);
summation = summation + z[j] ; #[J/mol]
summation = summation + 2*z[1] - z[1]; # Because in the above calculation,in order to calculate the average a summation of z(2) is not included,only half of it gets added
def f41(P):
return V - (R*T)/P
# Now we have, area = quad(f41,0,13.8*10**(5))[0]
# R*T*math.log(f/P) = summation
f = P_final*(math.exp(summation/(R*T))); #[bar]
print " 2).The fugacity of steam at 311 K and 13.8 bar pressure is %f bar"%(f);
# Variables
T = 0 + 273.15; #[K] - Temperature
R = 8.314; #[J/mol*K] - Universal gas constant
# Pressure in 'atm' is given below
P = [100,200,300,400,500,600,700,800,900,1000];
# The compressibility factor values are
Z = [1.069,1.138,1.209,1.283,1.356,1.431,1.504,1.577,1.649,1.720];
# Z = 1 + (B/V) + (C/V**(2))
# (Z-1)*V = B + (C/V)
from numpy import zeros
V = zeros(10);
k = zeros(10);
t = zeros(10);
# Calculations and Results
for i in range(10):
V[i]=Z[i]*R*T/(P[i]*101325); #[m**(3)/mol]
k[i]=(Z[i]-1)*V[i];
t[i]=1/V[i];
from scipy.stats import linregress
C,B,c,d,e = linregress(t,k)
#[C,B,sig]=reglin(t,k);
#From the regression, we get intercept = B and slope = C,and thus,
print " 1).The second virial coefficient of H2 is given by B = %e m**3)/mol"%(B);
print " The thied virial coefficient of H2 is given by C = %e m**6)/mol**2)"%(C);
# (2)
# We know that, limit P tending to zero (V-(R*T)/P) = B, therfore P = 0, V-(R*T)/P = B
def f42(P):
return V-(R*T)/P #and determine the integral integrate((V-(R*T)/P)
# Now let us tabulate V-(R*T)/P and determine the integral quad(f42,0,1000)[0]
P_prime = [0,100,200,300,400,500,600,700,800,900,1000];
Z_prime = [0,1.069,1.138,1.209,1.283,1.356,1.431,1.504,1.577,1.649,1.720];
summation = 0;
V_prime = zeros(11);
x = zeros(11);
y = zeros(11);
z = zeros(11);
for j in range(1,11):
V_prime[j]=Z_prime[j]*R*T/(P_prime[j]*101325); #[m**(3)/mol]
x[j]=V_prime[j]-(R*T)/(P_prime[j]*101325);
y[j]=(x[j] + x[j-1])/2;
z[j]=y[j]*((P_prime[j]-P_prime[j-1]))*101325;
summation = summation + z[j] ; #[J/mol]
summation = summation + 2*z[1] - z[1];
# Now we have
# R*T*math.log(f/P) = summation
P_dash = 1000; #[atm] - Pressure at which fugacity is to be calculated
T_dash = 273.15; #[K] - Temperature at which fugacity is to be calculated
f = P_dash*math.exp(summation/(R*T_dash)); #[atm]
print " 2).The fugacity of H2 at 0 C and 1000 atm pressure is f = %f atm"%(f);
from scipy.optimize import fsolve
import math
from scipy.integrate import quad
# Variables
P_1 = 1*10.**(6.); #[Pa] - Initial pressure
T_1 = 200. + 273.15; #[K] - Initial temperature
P_2 = 8*10.**(6); #[Pa] - Final pressure
Tc = 647.1; #[K] - Critical temperature of water
Pc = 220.55; #[bar] - Critical pressure of water
Pc = Pc*10**(5); #[Pa]
w = 0.345;
R = 8.314; #[J/mol*K] - Universal gas constant
# For the virial gas the following are the relations for residual enthalpy and entropy
# H_R = B*P - P*T*(dB/dT)
# S_R = -P*(dB/dT)
# Where, (dB/dT) = ((R*Tc)/Pc)*((dB_0/dT) + w*(dB_1/dT))
# dB0_dT = 0.422*1.6*Tc**(1.6)*T**(-2.6); # (dB_0/dT)
# dB1_dT = 0.172*4.2*Tc**(4.2)*T**(-5.2); # (dB_1/dT)
# (1)
Cp_0 = 29.114; #[J/mol-K] - Specific heat capacity at constant pressure
# For the isentropic process entropy change is zero, thus
# delta_S = Cp_0*math.log(T_2/T_1) - P_2*(dB/dT)_2 + P_1*(dB/dT)_1 = 0
# Calculations
# At state 1,
Tr_1 = T_1/Tc;
B0_1 = 0.083 - 0.422/(Tr_1**(1.6));
B1_1 = 0.139 - 0.172/(Tr_1**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B_1 = ((B0_1 + w*B1_1)*(R*Tc))/Pc; # [m**(3)/mol] - Second virial coefficient at state 1
dB0_dT_1 = 0.422*1.6*Tc**(1.6)*T_1**(-2.6); # (dB_0/dT)
dB1_dT_1 = 0.172*4.2*Tc**(4.2)*T_1**(-5.2); # (dB_1/dT)
dB_dT_1 = ((R*Tc)/Pc)*((dB0_dT_1) + w*(dB1_dT_1)); # (dB/dT)_1
# Now let us assume the exit temperature to be 870 K, at this temperature
# T_2 = 870; #[K] -
# At this temperature
# delta_S = Cp_0*math.log(T_2/T_1) - P_2*(dB/dT)_2 + P_1*(dB/dT)_1 =
T_2 = 860.; #[K] - Exit temperature
# Therefore at state 2, we have
Tr_2 = T_2/Tc;
B0_2 = 0.083 - 0.422/(Tr_2**(1.6));
B1_2 = 0.139 - 0.172/(Tr_2**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B_2 = ((B0_2 + w*B1_2)*(R*Tc))/Pc; # [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT_2 = 0.422*1.6*Tc**(1.6)*T_2**(-2.6); # (dB_0/dT)
dB1_dT_2 = 0.172*4.2*Tc**(4.2)*T_2**(-5.2); # (dB_1/dT)
dB_dT_2 = ((R*Tc)/Pc)*((dB0_dT_2) + w*(dB1_dT_2)); # (dB/dT)_2
delta_H_s = Cp_0*(T_2 - T_1) + B_2*P_2 -P_2*T_2*(dB_dT_2) - B_1*P_1 + P_1*T_1*(dB_dT_1); #[J/mol] - Enthalpy change
# As no heat exchange is assumed to take place with the surroundings,work transfer is given by
W_1 = - delta_H_s; # [J/mol]
# Results
print " 1).The exit temperature is %f K"%(T_2);
print " The required amount of work is %f J/mol"%(W_1);
# (2)
eff = 0.8; # Adiabatic efficiency
delta_H_a = delta_H_s/0.8; # Actual enthalpy change
# Now for calculating the value of T_exit
# delta_H_a = Cp_0*(T_exit - T_1) + B*P_2 -P_2*T_exit*(dB_dT) - B_1*P_1 + P_1*T_1*(dB_dT_1)
# On simplification we get
# 29.114*(T_2 - T_1)*B_2*8*10**(6)-8*10**(6)*T_2*(dB/dT)_2 = 12643.77
# Let us assume a temperature of say
T = 900.; #[K]
fault=10.;
def f(T_exit):
return delta_H_a - Cp_0*(T_exit - T_1) + B*P_2 -P_2*T_exit*(dB_dT) - B_1*P_1 + P_1*T_1*(dB_dT_1)
while(fault>0.3):
Tr = T/Tc;
B0 = 0.083 - 0.422/(Tr**(1.6));
B1 = 0.139 - 0.172/(Tr**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B = ((B0 + w*B1)*(R*Tc))/Pc; # [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT = 0.422*1.6*Tc**(1.6)*T**(-2.6); # (dB_0/dT)
dB1_dT = 0.172*4.2*Tc**(4.2)*T**(-5.2); # (dB_1/dT)
dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT)); # (dB/dT)_1
T_exit = fsolve(f,900)
fault=abs(T-T_exit);
T = T + 0.2;
Texit = T;
W_2 = - delta_H_a; # [J/mol]
print " 2).The exit temperature is %f K"%(Texit);
print " The required amount of work is %f J/mol"%(W_2);
def f20(T):
return Cp_0/T
T_prime = 700.; #[K]
fault1=10.;
def f1(T_out):
return 32.2168*math.log(T_out) + 0.1922*10**(-2)*T_out + 0.5274*10**(-5) \
*T_2**(2) - 1.1976*10**(-9)*T_out**(3)-8*10**(6)*dB_dT_prime -216.64
while(fault1>0.5):
Tr_prime = T_prime/Tc;
B0_prime = 0.083 - 0.422/(Tr_prime**(1.6));
B1_prime = 0.139 - 0.172/(Tr_prime**(4.2));
B_prime = ((B0_prime + w*B1_prime)*(R*Tc))/Pc; # [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_prime**(-2.6); # (dB_0/dT)
dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_prime**(-5.2); # (dB_1/dT)
dB_dT_prime = ((R*Tc)/Pc)*((dB0_dT_prime) + w*(dB1_dT_prime)); # (dB/dT)_1
T_out = fsolve(f1,10)
fault1=abs(T_prime-T_out);
T_prime = T_prime + 0.5;
T_out = T_prime;
# Now we have to calculate enthalpy change as W = -delta_H
def f21(T):
return (7.7 + 0.04594*10**(-2.)*T + 0.2521*10**(-5.)*T**(2.) - 0.8587*10.**(-9.)*T**(3.))*4.184
delta_H_3 = quad(f21,T_1,T_out)[0] + B_prime*P_2 - P_2*T_out*dB_dT_prime - B_1*P_1 + P_1*T_1*dB_dT_1
print T_1,T_out
W_3 = - delta_H_3; # [J/mol]
print " 3).The exit temperature is %f K"%(T_out);
print " The required amount of work is %f J/mol"%(W_3);
#(4)
n = 0.8; # Adiabatic efficiency
delta_H_a_4 = delta_H_3/n; #[J/mol]
W_4 = -delta_H_a_4; #[J/mol]
T_prime1 = 700.; #[K]
fault2=10.;
def f2(T_2):
return 7.7*4.184*(T_2-T_1) + ((0.04594*4.184*10**(-2))/2)* \
(T_2**(2)-T_1**(2)) + ((0.2521*4.184*10**(-5))/3)*(T_2**(3)-T_1**(3)) - \
((0.8587*4.184*10**(-9))/4)*(T_2**(4)-T_1**(4)) + B_prime1*8*10**(6) \
- 8*10**(6)*T_2*dB_dT_prime1 + 191.7 + 496.81 - delta_H_a_4
while(fault2>0.5):
Tr_prime1 = T_prime1/Tc;
B0_prime1 = 0.083 - 0.422/(Tr_prime1**(1.6));
B1_prime1 = 0.139 - 0.172/(Tr_prime1**(4.2));
# (B*Pc)/(R*Tc) = B0 + w*B1
B_prime1 = ((B0_prime1 + w*B1_prime1)*(R*Tc))/Pc; # [m**(3)/mol] - Second virial coefficient at state 2
dB0_dT_prime1 = 0.422*1.6*Tc**(1.6)*T_prime1**(-2.6); # (dB_0/dT)
dB1_dT_prime1 = 0.172*4.2*Tc**(4.2)*T_prime1**(-5.2); # (dB_1/dT)
dB_dT_prime1 = ((R*Tc)/Pc)*((dB0_dT_prime1) + w*(dB1_dT_prime1)); # (dB/dT)_1
T_out1 = fsolve(f2,100)
fault2=abs(T_prime1-T_out1);
T_prime1 = T_prime1 + 0.5;
T_out1 = T_prime1;
print " 4).The exit temperature is %f K"%(T_out1);
print " The required amount of work is %f J/mol"%(W_4);
# answers are varies because of rounding error. please check it manually.
# Variables
Vol = 0.15; #[m**(3)] - Volume of the cylinder
P_1 = 100.; #[bar] - Initial pressure
P_1 = P_1*10**(5); #[Pa]
T_1 = 170.; #[K] - Initial temperature
n_withdrawn = 500.; #[mol] - Withdrawn number of moles
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations and Results
#(1)
Y = 1.4; # Coefficient of adiabatic expansion
n_total = (P_1*Vol)/(R*T_1); #[mol] - Total number of moles
n_2 = n_total - n_withdrawn; #[mol] - Left number of moles
V_1 = Vol/n_total; #[m**(3)/mol] - Molar volume at initial state.
# At final state
V_2 = Vol/n_2; #[m**(3)/mol] - Molar volume at final state
# During duscharging P_1*V_1**(Y) = P_2*V_2**(Y), therefore
P_2_1 = P_1*((V_1/V_2)**(Y)); #[Pa] - Final pressure
P_2_1 = P_2_1*10**(-5); #[bar]
T_2_1 = ((P_2_1*10**(5))*V_2)/R; #[K] - Final temperature
print " 1).The final temperature %f K"%(T_2_1);
print " The final pressure %f bar"%(P_2_1);
def f53(T):
return Cp_0/T
T_prime = 150; #[K]
error = 10;
while(error>1):
f_T = 18.886*math.log(T_prime) + 4.2*10**(-3)*T_prime - 92.4;
f_T_dash = 18.886/T_prime + 4.2*10**(-3);
T_new = T_prime - (f_T/f_T_dash);
error=abs(T_prime - T_new);
T_prime = T_new;
T_2_2 = T_prime; #[K] - Final temperature
P_2_2 = ((n_2*R*T_2_2)/Vol)*10**(-5); #[bar] - Final pressure
print " 2).The final temperature %f K"%(T_2_2);
print " The final pressure %f bar"%(P_2_2);
#(3)
Tc = 126.2; #[K] - Critical temperature of nitrogen
Pc = 34.0; #[bar] - Critical pressure of nitrogen
Pc = Pc*10**(5); #[Pa]
w = 0.038; # Acentric factor
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6); # (dB_0/dT) at state 1
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2); # (dB_1/dT) at state 1
dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT)); # (dB/dT) at state 1
# The residual entropy at the initial state is given by
S_R_1 = -P_1*(dB_dT); #[J/mol-K]
# Now let us calculate molar volume at initial state
Tr = T_1/Tc; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0+(w*B_1))*(R*Tc))/Pc; #[m**(3)/mol]
V_1_3 = B + (R*T_1)/P_1; #[m**(3)/mol]
# Therefore number of moles in the initial state is
n_1_3 = Vol/V_1_3; #[mol]
# Therefore final number of moles is
n_2_3 = n_1_3 - n_withdrawn;
# Therefore molar volume at final state is
V_2_3 = Vol/n_2_3; #[m**(3)/mol]
# Now let us determine the relation between pressure and temperature in the final state
# P_2_3 = (R*T_2_3)/(V_2_3 - B_2)
#delta_S = 0, thus delta_S_ig + delta_S_R = 0
delta_S_R = - S_R_1;
def f54(T):
return Cp_0/T
T_2_prime = 135; #[K]
delta = 0.1;
error = 10;
while(error>0.01):
T_r = T_2_prime/Tc; # Reduced temperature
B_0_3 = 0.083-(0.422/(T_r)**(1.6));
B_1_3 = 0.139-(0.172/(T_r)**(4.2));
B_3 = ((B_0_3+(w*B_1_3))*(R*Tc))/Pc; #[m**(3)/mol]
dB0_dT_3 = 0.422*1.6*Tc**(1.6)*T_2_prime**(-2.6); # (dB_0/dT)
dB1_dT_3 = 0.172*4.2*Tc**(4.2)*T_2_prime**(-5.2); # (dB_1/dT)
dB_dT_3 = ((R*Tc)/Pc)*((dB0_dT_3) + w*(dB1_dT_3)); # (dB/dT)
P_2_3 = (R*T_2_prime)/(V_2_3 - B_3);
delta_S = 27.2*(math.log(T_2_prime/T_1)) + 4.2*10**(-3)*(T_2_prime - T_1) - R*(math.log(P_2_3/P_1)) - P_2_3*(dB_dT_3) + delta_S_R;
T_new = T_2_prime + delta;
error=abs(delta_S);
T_2_prime = T_new;
T_2_3 = T_2_prime; #[K] - Final temperature
# Therefore at T_2_3
P_2_3 = P_2_3*10**(-5); #[bar] - Final pressure
print " 3).The final temperature %f K"%(T_2_3);
print " The final pressure %f bar"%(P_2_3);
# Variables
P_1 = 80.; #[bar] - Initial pressure
P_1 = P_1*10.**(5); #[Pa]
T_1 = 300. + 273.15; #[T] - Initial temperature
P_2 = 40.; #[bar] - Final pressure
P_2 = P_2*10**(5); #[Pa]
T_2 = 300. + 273.15; #[K] - Final temperature
T_0 = 25. + 273.15; #[K] - Surrounding temperature
P_0 = 1.; #[atm] - Surrounding pressure
P_0 = P_0*101325; #[Pa]
Tc = 647.1; #[K]
Pc = 220.55; #[bar]
Pc = Pc*10.**(5); #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc); #[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc); #[m**(3)/mol]
def f(V):
return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1
V_1_1=fsolve(f,0.1)
V_1_2=fsolve(f,10)
V_1_2=fsolve(f,100)
# The largest root is considered because of vapour
V_1 = V_1_1;
U_R_1 = -a/V_1; #[J/mol] - Internal energy
H_R_1 = P_1*V_1 - R*T_1 - a/V_1; #[J/mol] - Enthalpy
S_R_1 = R*math.log((P_1*(V_1-b))/(R*T_1));
def f1(V):
return V**(3)-(b+(R*T_2)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2
V_2_1 = fsolve(f1,0.1)
V_2_2 = fsolve(f1,10)
V_2_3 = fsolve(f1,100)
V_2 = V_2_1;
U_R_2 = -a/V_2; #[J/mol] - Internal energy
H_R_2 = P_2*V_2 - R*T_2 - a/V_2; #[J/mol] - Enthalpy
S_R_2 = R*math.log((P_2*(V_2-b))/(R*T_2));
delta_U_R = U_R_2 - U_R_1; #
delta_H_R = H_R_2 - H_R_1; #
delta_S_R = S_R_2 - S_R_1; #
delta_U_ig = 0; #[J/mol] - As temperature is constant
delta_H_ig = 0; #[J/mol] - As temperature is constant
delta_S_ig = - R*math.log(P_2/P_1); # [J/mol-K]
delta_U = delta_U_R + delta_U_ig; #[J/mol]
delta_H = delta_H_R + delta_H_ig; #[J/mol]
delta_S = delta_S_R + delta_S_ig; #[J/mol-K]
# Change in exergy is given by
# delta_phi = phi_1 - phi_2 = U_1 - U_2 + P_0*(V_1 - _V_2) - T_0*(S_1 - S_2)
delta_phi = - delta_U + P_0*(V_1 - V_2) - T_0*(-delta_S); #[J/mol]
# Results
print " The change in internal energy is %f J/mol"%(delta_U);
print " The change in enthalpy is %f J/mol"%(delta_H);
print " The change in entropy is %f J/mol-K"%(delta_S);
print " The change in exergy is %f J/mol"%(delta_phi);
# Variables
T_1 = 500.; #[K] - Initial temperature
P_1 = 30.; #[atm] - Initial pressure
P_1 = P_1*101325; #[Pa]
P_2 = 1; #[atm] - Final pressure
P_2 = P_2*101325; #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
# For chlorine
Tc = 417.2; #[K] - Critical temperature
Pc = 77.10; #[bar] - Critical pressure
Pc = Pc*10**(5); #[Pa]
#Redlich Kwong equation of state,
a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc; # [Pa*m**(6)*K**(1/2)/mol]
b = (0.08664*R*Tc)/Pc; # [m**(3)/mol]
# Calculations
def f1(V):
return V**(3)-((R*T_1)/P_1)*V**(2)-((b**(2))+((b*R*T_1)/P_1)-(a/(T_1**(1./2)*P_1)))*V-(a*b)/(T_1**(1./2)*P_1)
V_1=fsolve(f1,1)
V_2=fsolve(f1,10)
V_3=fsolve(f1,100)
V = V_1; #[m**(3)/mol]
Z = (P_1*V_1)/(R*T_1); #compressibility factor
# The residual enthalpy at state 1 is given by
H_R_1 = (Z-1)*R*T_1 + ((3*a)/(2*b*T_1**(1./2)))*(math.log(V/(V+b))); #[J/mol]
H_R_2 = 0; # Residual enthalpy at state 2
delta_H_R = H_R_2 - H_R_1; #[J/mol] - Residual enthalpy change
# and since isothermal conditions are maintained, therfore
delta_H_ig = 0; # Enthalpy change under ideal condition
delta_H = delta_H_R + delta_H_ig; #[J/mol]
# Results
print " The change in enthalpy is given by delta_H = %f J/mol"%(delta_H);
# Variables
Vol_1 = 0.1; #[m**(3)] - Initial volume of each compartment
n_1 = 400; #[mol] - Initial number of moles in compartment 1
V_1 = Vol_1/n_1; #[m**(3)/mol] - Molar volume at state 1
T_1 = 294; #[K]
Vol_2 = 0.2; #[m**(3)] - Final volume of the compartment after removing the partition.
n_2 = n_1; #[mol] - Number of moles remains the same
V_2 = Vol_2/n_2; #[m**(3)/mol] - Molar volume at state 2
a = 0.1362; #[Pa-m**(6)/mol**(2)]
b = 3.215*10**(-5); #[m**(3)/mol]
Cv_0 = 12.56; #[J/mol-K] - Heat capacity in ideal gas state
# Calculations
# For overall system q = 0, and no work is done, therefore delta_U = 0
# Therfore from the relation proved in part (1), we have
T_2 = T_1 + (a/Cv_0)*(1/V_2 - 1/V_1); #[K]
# Results
print " 2).The final temperatutre is %f K"%(T_2)
# Variables
P_1 = 1*10**(6); #[Pa] - Initial pressure
T_1 = 200 + 273.15; #[K] - Initial temperature
P_2 = 8*10**(6); #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
Y = 1.4; # Index of expansion
Cp_0 = 29.114; #[J/mol-K]
a = 0.55366; #[Pa-m**(6)/mol**(2)]
b = 3.049*10**(-5); #[m**(3)/mol]
# Calculations
V_1 = 3.816*10**(-3); #[m**(3)/mol]
Z_1 = (P_1*V_1)/(R*T_1);
T_2 = T_1*(P_2/P_1)**((Y-1)/Y); #[K]
# At 8 MPa and T_2,
# The molar volume of steam following van der Walls equation of state (as reported in the book) is
V_2 = 8.41*10**(-4); #[m**(3)/mol]
# And the compressibility factor is
Z_2 = (P_2*V_2)/(R*T_2);
# For van der Walls equation of state we know that
# delta_S_R/R = math.log(Z_2/Z_1) + math.log((V_2 - b)/V_2) - math.log((V_1 - b)/V_1)
delta_S_R = R*(math.log(Z_2/Z_1) + math.log((V_2 - b)/V_2) - math.log((V_1 - b)/V_1)); #[J/mol]
# delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1)
# The entropy change is therefore
# delta_S = delta_S_ig + delta_S_R
# Let us assume a temperature, say T = 870 K
# At 870 K the molar volume of steam following van der Walls equation of state (as reported in the book) is
# V_3 = 8.57*10**(-4); # [m**(3)/mol]
# Therefore
# Z_3 = (P_2*V_3)/(R*T_2);
# At this temperature,
# delta_S = Cp_0*math.log(T/T_1) - R*math.log(P_2/P_1) + R*(math.log(Z/Z_1) + R*math.log((V - b)/V) - R*math.log((V_1 - b)/V_1))
T = 800; #[K]
fault=10;
def f1(V):
return V**(3)-(b+(R*T)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2
def f2(T):
return Cp_0*math.log(T/T_1) - R*math.log(P_2/P_1) + R*(math.log(Z/Z_1) + R*(math.log((V - b)/V)) - R*(math.log((V_1 - b)/V_1)))
while(fault>0.3):
# At T and 8 MPa
V = fsolve(f1,1)
Z = (P_2*V)/(R*T);
T_exit = fsolve(f2,0.1)
fault=abs(T-T_exit);
T = T + 0.5;
Texit = T;
# W = - delta_H
# For van der Walls gas the enthalpy change is given by
delta_H_s = Cp_0*(T_exit - T_1) + (Z - 1)*R*T_exit - a/V - (Z_1-1)*R*T_1 + a/V_1; #[J/mol]
W = - delta_H_s; #[J/mol]
# Results
print " 1).The exit temperature is %f K"%(Texit);
print " The work required is given by W = %f J/mol"%(W);
#(2)
eff = 0.8; # Adiabatic efficiency
delta_H_a = eff*delta_H_s; #[J/mol] - Actual enthalpy change
W_2 = - delta_H_a;
# Let us assume a temperature, say
T_prime= 900; #[K]
fault1=10;
def f22(V):
return V**(3)-(b+(R*T_prime)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2
def f3(T_prime):
return Cp_0*(T_prime - T_1) + (Z_prime - 1)*R*T_prime - a/V_prime - 13230.49
while(fault1>0.3):
# At T_prime and 8 MPa
V_prime=fsolve(f22,1)
Z_prime = (P_2*V_prime)/(R*T_prime);
T_exit1 = fsolve(f3,100)
fault1=abs(T_prime-T_exit1);
T_prime = T_prime + 0.2;
Texit1 = T_prime;
print " 2).The exit temperature is %f K"%(Texit1);
print " The work required is given by W = %f J/mol"%(W_2);
# Variables
T = 100 + 273.15; #[K] - Temperature
Tc = 647.1; #[K] - Critical temperature of water
Pc = 220.55; #[bar] - Critical pressure of water
Pc = Pc*10**(5); #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations
a = (27*R**(2)*Tc**(2))/(64*Pc); #[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc); #[m**(3)/mol]
# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0
# For water vapour at 100 C under saturated conditions pressure is 1 atm, therefore
P = 1; #[atm]
P = P*101325; #[Pa]
# At 100 C and 1 atm
def f(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V_1 = fsolve(f,0.1)
V_1 = fsolve(f,10)
V_1 = fsolve(f,100)
V = V_1; #[m**(3)/mol]
f = P*(math.exp(math.log((R*T)/(P*(V-b))) + b/(V-b) - (2*a)/(R*T*V))); #[Pa]
f = f/101325; #[atm]
# Results
print " The molar volume is %f m**3)/mol"%(V);
print " The fugacity is %f atm"%(f);
P_1 = 6; #[bar] - Initial pressure
P_1 = P_1*10**(5); #[Pa]
T_1 = 100 + 273.15; #[T] - Initial temperature
P_2 = 12; #[bar] - Final pressure
P_2 = P_2*10**(5); #[Pa]
T_2 = 500 + 273.15; #[K] - Final temperature
R = 8.314; #[J/mol*K] - Universal gas constant
Y = 1.126; # Index of expansion
Cp_0 = (R*Y)/(Y-1); #[J/mol-K]
# For propane
Tc = 369.8; #[K]
Pc = 42.48; #[bar]
Pc = Pc*10**(5);
w = 0.152;
# Calculations and Results
#(1)
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc); #[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc); #[m**(3)/mol]
# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0
# At state 1 (100 C and 6 bar)
def f(V):
return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1
V_1_1 = fsolve(f,1)
V_1_2 = fsolve(f,10)
V_1_3 = fsolve(f,100)
V_1 = V_1_1; #[m**(3)/mol]
# Thus compressibility factor is
Z_1 = (P_1*V_1)/(R*T_1); #compressibility factor
H_R_1 = (Z_1 - 1)*R*T_1 - (a/V_1); # [J/mol]
S_R_1 = R*math.log((P_1*(V_1-b))/(R*T_1)); # [J/mol-K]
# At state 2 (500 C and 12 bar)
def f1(V):
return V**(3)-(b+(R*T_2)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2
V_2_1 = fsolve(f1,1)
V_2_2 = fsolve(f1,10)
V_2_3 = fsolve(f1,100)
# The largest root is considered because of molar volume of vapour phase is to determined
V_2 = V_2_1; #[m**(3)/mol]
# Thus compressibility factor is
Z_2 = (P_2*V_2)/(R*T_2); #compressibility factor
H_R_2 = (Z_2 - 1)*R*T_2 - (a/V_2); # [J/mol]
S_R_2 = R*math.log((P_2*(V_2-b))/(R*T_2)); # [J/mol-K]
# Ideal gas entropy change is given by
delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1); #[J/mol-K]
# Entropy change is given by
delta_S = delta_S_ig + (S_R_2 - S_R_1); #[J/mol-k]
# Ideal gas enthalpy change is given by
delta_H_ig = Cp_0*(T_2 - T_1); #[J/mol]
# Enthalpy change is given by
delta_H = delta_H_ig + (H_R_2 - H_R_1); #[J/mol]
print "1).The change in enthalpy is %f J/mol"%(delta_H);
print " The change in entropy is %f J/mol-K"%(delta_S);
#(2)
# Virial equation of state
# At state 1 (372.15 K, 6 bar) let us calculate B and dB/dT
Tr = T_1/Tc; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0+(w*B_1))*(R*Tc))/Pc; #[m**(3)/mol]
dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6); # (dB_0/dT) at state 1
dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2); # (dB_1/dT) at state 1
dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT)); # (dB/dT) at state 1
H_R_1_2 = B*P_1 - P_1*T_1*dB_dT; #[J/mol] - Residual enthalpy at state 1
S_R_1_2 = -P_1*(dB_dT); #[J/mol-K] - Residual entropy at state 1
# At state 2 (773.15 K, 12 bar)
Tr_2 = T_2/Tc; # Reduced temperature
B_0_2 = 0.083-(0.422/(Tr_2)**(1.6));
B_1_2 = 0.139-(0.172/(Tr_2)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B_2 = ((B_0_2+(w*B_1_2))*(R*Tc))/Pc; #[m**(3)/mol]
dB0_dT_2 = 0.422*1.6*Tc**(1.6)*T_2**(-2.6); # (dB_0/dT) at state 1
dB1_dT_2 = 0.172*4.2*Tc**(4.2)*T_2**(-5.2); # (dB_1/dT) at state 1
dB_dT_2 = ((R*Tc)/Pc)*((dB0_dT_2) + w*(dB1_dT_2)); # (dB/dT) at state 1
H_R_2_2 = B_2*P_2 - P_2*T_2*dB_dT_2; #[J/mol] - Residual enthalpy at state 1
S_R_2_2 = -P_2*(dB_dT_2); #[J/mol-K] - Residual entropy at state 1
delta_H_2 = delta_H_ig + (H_R_2_2 - H_R_1_2); #[J/mol]
delta_S_2 = delta_S_ig + (S_R_2_2 - S_R_1_2); #[J/mol]
print "2).The change in enthalpy is %f J/mol"%(delta_H_2);
print " The change in entropy is %f J/mol-K"%(delta_S_2);
# Variables
P = 2.76*10**(6); #[N/m**(2)] - Pressure
T = 310.93; #[K] - Temperature
R = 8.314; #[J/mol*K] - Universal gas constant
# For n-butane
Tc = 425.18; #[K] - Critical temperature
Pc = 37.97; #[bar] - Critical pressure
Pc = Pc*10**(5); #[Pa]
w = 0.193;
den = 0.61; #[g/cm**(3)]
mol_wt = 58; #[g/mol] - Molecular weight of butane
# Calculations and Results
# math.log(P_sat) = 15.7374 - 2151.63/(T-36.24)
P_sat = math.exp(15.7374 - 2151.63/(T-36.24)); #[mm Hg]
P_sat = (P_sat/760)*101325; #[N/m**(2)]
#(1)
# Let us determine the second virial coefficient at 310.93 K
Tr = T/Tc; # Reduced temperature
B_0 = 0.083-(0.422/(Tr)**(1.6));
B_1 = 0.139-(0.172/(Tr)**(4.2));
#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)
B = ((B_0+(w*B_1))*(R*Tc))/Pc; #[m**(3)/mol]
# Fugacity under saturated conditions is given by
# math.log(f_sat/P_sat) = (B*P_sat)/(R*T)
f_sat = P_sat*(math.exp((B*P_sat)/(R*T))); #[N/m**(2)]
# The molar volume is given by
V_liq = (1/(den*1000))*(mol_wt/1000); #[m**(3)/mol]
f = f_sat*math.exp(V_liq*(P-P_sat)/(R*T));
print " 1).The fugacity of n-butane is %e N/m**2)"%(f);
#(2)
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc); #[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc); #[m**(3)/mol]
# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0
# At 100 C and 1 atm
def f(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V_1 = fsolve(f,0.1)
V_1 = fsolve(f,10)
V_1 = fsolve(f,100)
# The above equation has only 1 real root, other two roots are imaginary
V = V_1; #[m**(3)/mol]
# math.log(f/P) = math.log((R*T)/(P*(V-b))) + b/(V-b) -(2*a)/(R*T*V)
f_2 = P*(math.exp(math.log((R*T)/(P*(V-b))) + b/(V-b) -(2*a)/(R*T*V)));
print " 2).The fugacity of n-butane is %e N/m**2)"%(f_2);
# Variables
T = 50+273.15; #[K] - Temperature
P = 25.*10**(3); #[Pa] - Pressure
y1 = 0.5; #[mol] - mole fraction of equimolar mixture
y2 = 0.5;
R = 8.314; #[J/mol*K] - Universal gas constant
#For component 1 (methyl ethyl ketone)
Tc_1 = 535.5; #[K] - Critical temperature
Pc_1 = 41.5*10**(5); #[N/m**(2)] - Critical pressure
Vc_1 = 267.; #[cm**(3)/mol] - Critical volume
Zc_1 = 0.249; # Critical compressibility factor
w_1 = 0.323; # acentric factor
#For component 2 (toluene)
Tc_2 = 591.8; #[K]
Pc_2 = 41.06*10**(5); #[N/m**(2)]
Vc_2 = 316.; #[cm**(3)/mol]
Zc_2 = 0.264;
w_2 = 0.262;
# Calculations and Results
# For equation of state Z = 1 + B/V
#For component 1, let us calculate B and dB/dT
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-K]
dB0_dT_1 = 0.422*1.6*Tc_1**(1.6)*T**(-2.6); # [m**(3)/mol-K] - (dB_0/dT)
dB1_dT_1 = 0.172*4.2*Tc_1**(4.2)*T**(-5.2); # [m**(3)/mol-K] - (dB_1/dT)
dB_dT_1 = ((R*Tc_1)/Pc_1)*((dB0_dT_1) + w_1*(dB1_dT_1)); #[m**(3)/mol-K] - (dB/dT)_
#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]
dB0_dT_2 = 0.422*1.6*Tc_2**(1.6)*T**(-2.6); # [m**(3)/mol-K] - (dB_0/dT)
dB1_dT_2 = 0.172*4.2*Tc_2**(4.2)*T**(-5.2); # [m**(3)/mol-K] - (dB_1/dT)
dB_dT_2 = ((R*Tc_2)/Pc_2)*((dB0_dT_2) + w_2*(dB1_dT_2)); #[m**(3)/mol-K] - (dB/dT)_
#For cross coeffcient, let us calculate B and dB/dT
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)]
#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]
dB0_dT_12 = 0.422*1.6*Tc_12**(1.6)*T**(-2.6); # [m**(3)/mol-K] - (dB_0/dT)
dB1_dT_12 = 0.172*4.2*Tc_12**(4.2)*T**(-5.2); # [m**(3)/mol-K] - (dB_1/dT)
dB_dT_12 = ((R*Tc_12)/Pc_12)*((dB0_dT_12) + w_12*(dB1_dT_12)); #[m**(3)/mol-K] - (dB/dT)_12
#For the mixture
B = y1**(2)*B_11 + 2*y1*y2*B_12 + y2**(2)*B_22; #[m**(3)/moL]
# The equation of state can be written as
# V**(2) - ((R*T)/P) - (B*R*T)/P = 0
# V**(2) - 0.1075*V + 1.737*10**(-4) = 0
def f(V):
return V**(2) - 0.1075*V + 1.737*10**(-4)
V1 = fsolve(f,0.1)
V2 = fsolve(f,1)
# We will consider the root which is near to R*T/P
V = V1;
# dB/dT = y_1**(2)*dB_11/dT + y_2**(2)*dB_22/dT + 2*y_1*y_2*dB_12/dT
dB_dT = y1**(2)*dB_dT_1 + y2**(2)*dB_dT_2 + 2*y1*y2*dB_dT_12; #[m**(3)/mol-K]
# For equation of state Z = 1 + B/V
H_R = (B*R*T)/V - ((R*T**(2))/V)*dB_dT; #[J/mol]
print " 1).The value of H_R for the mixture using virial equation of state is %f J/mol"%(H_R);
#(2)
# For van der Walls equation of state
a_11 = (27*R**(2)*Tc_1**(2))/(64*Pc_1); #[Pa-m**(6)/mol**(2)]
a_22 = (27*R**(2)*Tc_2**(2))/(64*Pc_2); #[Pa-m**(6)/mol**(2)]
a_12 = (a_11*a_22)**(1./2);
b_1 = (R*Tc_1)/(8*Pc_1); #[m**(3)/mol]
b_2 = (R*Tc_2)/(8*Pc_2); #[m**(3)/mol]
# 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]
# From the cubic form of van der Walls equation of state
def f1(V):
return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P
V2_1 = fsolve(f1,0.1)
V2_2 = fsolve(f1,10)
V2_3 = fsolve(f1,100)
# The largest root is considered
V_2 = V2_1;
# The residual enthalpy is given by
H_R_2 = P*V_2 - R*T -a/V_2; #[J/mol]
print " 2).The value of H_R for the mixture using van der Walls equation of state is %f J/mol"%(H_R_2);
from scipy.optimize import fsolve
import math
# Variables
T = 320 + 273.15; #[K]
R = 8.314; #[J/mol*K] - Universal gas constant
# For water
Tc = 647.1; #[K]
Pc = 220.55; #[bar]
Pc = Pc*10**(5); #[Pa]
# The cubic form of Redlich Kwong equation of state is given by,
# V**(3) - ((R*T)/P)*V**(2) - ((b_1**(2)) + ((b_1*R*T)/P) - (a/(T**(1/2)*P))*V - (a*b)/(T**(1/2)*P) = 0
# At 320 C and 70 bar pressure
P_1 = 70; #[bar]
P_1 = P_1*10**(5); #[Pa]
# Calculations and Results
a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc; #[Pa*m**(6)*K**(1/2)/mol]
b = (0.08664*R*Tc)/Pc; #[m**(3)/mol]
# Solving the cubic equation
def f1(V):
return V**(3)-((R*T)/P_1)*V**(2)-((b**(2))+((b*R*T)/P_1)-(a/(T**(1./2)*P_1)))*V-(a*b)/(T**(1./2)*P_1)
V1=fsolve(f1,1)
V2=fsolve(f1,10)
V3=fsolve(f1,100)
# The largest root is considered because at 320 C and 70 bar vapour phase exists.
V_1 = V1; #[m**(3)/mol]
# Thus compressibility factor is
Z_1 = (P_1*V_1)/(R*T);
# For Redlich-Kwong equation of state
# math.log(f/P) = Z - 1 - math.log(V_1/(V_1-b)) + (a/(b*R*(T**(3/2))))*math.log(V/(V+b))
f_1 = P_1*(math.exp(Z_1-1-math.log(Z_1)+math.log(V_1/(V_1-b))+(a/(b*R*(T**(3./2))))*math.log(V_1/(V_1+b)))); #[Pa]
f_1 = f_1*10**(-5); #[bar]
print " The fugacity of water vapour at 320 C and 70 bar pressure is %f bar"%(f_1);
# At 320 C and 170 bar pressure, we have
P_2 = 170; #[bar]
P_2 = P_2*10**(5); #[Pa]
# Solving the cubic equation
def f2(V):
return V**(3)-((R*T)/P_2)*V**(2)-((b**(2))+((b*R*T)/P_2)-(a/(T**(1./2)*P_2)))*V-(a*b)/(T**(1./2)*P_2)
V4 = fsolve(f2,1)
V5 = fsolve(f2,10)
V6 = fsolve(f2,100)
# The above equation has only 1 real root,other two roots are imaginary. Therefore,
V_2 = V6; #[m**(3)/mol]
# Thus compressibility factor is
Z_2 = (P_2*V_2)/(R*T);
# For Redlich-Kwong equation of state
# math.log(f/P) = Z - 1 - math.log(V_1/(V_1-b)) + (a/(b*R*(T**(3/2))))*math.log(V/(V+b))
f_2 = P_2*(math.exp(Z_2-1-math.log(Z_2)+math.log(V_2/(V_2-b))+(a/(b*R*(T**(3./2))))*math.log(V_2/(V_2+b)))); #[Pa]
f_2 = f_2*10**(-5); #[bar]
print " The fugacity of water vapour at 320 C and 170 bar pressure is %f bar"%(f_2);
from scipy.optimize import fsolve
import math
# Variables
Vol = 0.057; #[m**(3)] - Volume of car tyre
P_1 = 300.; #[kPa] - Initial pressure
P_1 = P_1*10**(3); #[Pa]
T_1 = 300.; #[K] - Initial temperature
P_2 = 330.; #[kPa] - Finnal pressure
P_2 = P_2*10**(3); #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
Cv_0 = 21.; #[J/mol-K] - Heat capacity for air
# For oxygen
Tc_O2 = 154.6; #[K] - Critical temperature
Pc_O2 = 50.43; #[bar] - Critical pressure
Pc_O2 = Pc_O2*10**(5); #[Pa]
y1 = 0.21; # - Mole fraction of oxygen
# For nitrogen
Tc_N2 = 126.2; #[K] - Critical temperature
Pc_N2 = 34.00; #[bar] - Critical pressure
Pc_N2 = Pc_N2*10**(5); #[Pa]
y2 = 0.79; # - Mole fraction of nitrogen
# Calculations and Results
# (1)
# Assuming ideal gas behaviour. The volume remains the same,therefore,we get
# P_1/T_1 = P_2/T_2
T_2 = P_2*(T_1/P_1); #[K]
n = (P_1*Vol)/(R*T_1); #[mol] - Number of moles
delta_U = n*Cv_0*(T_2-T_1); #[J]
print " 1).The change in internal energy for ideal gas behaviour) is %f J"%(delta_U);
#(2)
# For van der Walls equation of state
a_O2 = (27*R**(2)*Tc_O2**(2))/(64*Pc_O2); #[Pa-m**(6)/mol**(2)]
a_N2 = (27*R**(2)*Tc_N2**(2))/(64*Pc_N2); #[Pa-m**(6)/mol**(2)]
a_12 = (a_O2*a_N2)**(1./2);
b_O2 = (R*Tc_O2)/(8*Pc_O2); #[m**(3)/mol]
b_N2 = (R*Tc_N2)/(8*Pc_N2); #[m**(3)/mol]
# For the mixture
a = y1**(2)*a_O2 + y2**(2)*a_N2 + 2*y1*y2*a_12; #[Pa-m**(6)/mol**(2)]
b = y1*b_O2 + y2*b_N2; #[m**(3)/mol]
# From the cubic form of van der Walls equation of state
# At 300 K and 300 kPa,
def f1(V):
return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1
V_1 = fsolve(f1,0.1)
V_2 = fsolve(f1,10)
V_3 = fsolve(f1,100)
# The above equation has only 1 real root, other two roots are imaginary
V = V_1; #[m**(3)/mol]
# Now at P = 330 kPa and at molar volume V
# The van der Walls equation of state is
# (P + a/V**(2))*(V - b) = R*T
T_2_prime = ((P_2 + a/V**(2))*(V - b))/R; #[K]
n_prime = Vol/V; #[mol]
# Total change in internal energy is given by
# delta_U_prime = n_prime*delta_U_ig + n_prime*(U_R_2 - U_R_1)
# delta_U_prime = n_prime*Cv_0*(T_2_prime - T_1) + n_prime*a(1/V_2 - 1/V_1);
# Since V_1 = V_2 = V, therefore
delta_U_prime = n_prime*Cv_0*(T_2_prime - T_1);
print " 2).The change in internal energy for van der Walls equation of state) is %f J"%(delta_U_prime);
from scipy.optimize import fsolve
import math
# Variables
T_1 = 150 + 273.15; #[K] - Initial emperature
T_2 = T_1; # Isothermal process
P_1 = 100.*10**(3); #[Pa] - Initial pressure
P_2 = 450.*10**(3); #[Pa] - Final pressure
R = 8.314; #[J/mol*K] - Universal gas constant
# For water
Tc = 647.1; #[K] - Critical temperature
Pc = 220.55; #[bar] - Critical pressure
Pc = Pc*10.**(5);
w = 0.345;
Mol_wt = 18.015; #[g/mol] - Molecular weight of water
Cp_0 = 4.18; #[J/mol-K] - Smath.tan(math.radiansard heat capacity of water
# Calculations
# In Peng-Robinson equation of state
m = 0.37464 + 1.54226*w - 0.26992*w**(2);
# At T_1 and P_1, we have
Tr = T_1/Tc;
alpha = (1 + m*(1 - Tr**(1./2)))**(2);
a = ((0.45724*(R*Tc)**(2))/Pc)*alpha; #[Pa*m**(6)/mol**(2)]
b = (0.07780*R*Tc)/Pc; #[m**(3)/mol]
# Cubuc form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f(V):
return V**(3)+(b-(R*T_1)/P_1)*V**(2)-((3*b**(2))+((2*R*T_1*b)/P_1)-(a/P_1))*V+b**(3)+((R*T_1*(b**(2)))/P_1)-((a*b)/P_1)
V1 = fsolve(f,-1)
V2 = fsolve(f,0)
V3 = fsolve(f,1)
#The largest root is for vapour phase,
#The largest root is only considered as the systemis gas
V_1 = V3; #[m**(3)/mol]
# Thus compressibility factor is
Z_1 = (P_1*V_1)/(R*T_1);
# At T_2 and P_2, we have
# Cubuc form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f1(V):
return V**(3)+(b-(R*T_2)/P_2)*V**(2)-((3*b**(2))+((2*R*T_2*b)/P_2)-(a/P_2))*V+b**(3)+((R*T_2*(b**(2)))/P_2)-((a*b)/P_2)
V4 = fsolve(f1,-1)
V5 = fsolve(f1,0)
V6 = fsolve(f1,1)
V_2 = V6; #[m**(3)/mol]
# Thus compressibility factor is
Z_2 = (P_2*V_2)/(R*T_2);
# In the Peng-Robinson equation of stste
# da/dT = -(a*m)/((alpha*T*Tc)**(1/2))
# The residual enthalpy is given by
# H_R = R*T*(Z-1) + (((T*(da_dT))-a)/(2*2**(1/2)*b))*math.log((Z+(1+2**(1/2)*((P*b)/(R*T))))/(Z+(1-2**(1/2)*((P*b)/(R*T)))))
# At state 1
da_dT_1 = -(a*m)/((alpha*T_1*Tc)**(1./2)); #[Pa*m**(6)/mol**(2)]
H_R_1 = R*T_1*(Z_1-1) + (((T_1*(da_dT_1))-a)/(2*(2**(1./2))*b))* \
math.log((Z_1+(1+2**(1./2))*((P_1*b)/(R*T_1)))/(Z_1+(1-2**(1./2))*((P_1*b)/(R*T_1))));
# At state 2
da_dT_2 = -(a*m)/((alpha*T_2*Tc)**(1/2.)); #[Pa*m**(6)/mol**(2)]
H_R_2 = R*T_2*(Z_2-1) + (((T_2*(da_dT_2))-a)/(2*2**(1./2)* \
b))*math.log((Z_2+(1+2**(1./2))*((P_2*b)/(R*T_1)))/(Z_2+(1-2**(1./2))*((P_2*b)/(R*T_1))));
delta_H = H_R_2 - H_R_1; #[J/mol]
delta_H = delta_H/Mol_wt; #[kJ/kg]
# The residual entropy relation for a substance following Peng - Robinson equation of state ia
# S_R = R*math.log(Z - (P*b)/(R*T)) + (da_dT/(2*2**(1/2)*b))*math.log((Z+(1+2**(1/2))*((P*b)/(R*T)))/(Z+(1-2**(1/2))*((P*b)/(R*T))))
# The residual entropy at state 1 is
S_R_1 = R*math.log(Z_1 - (P_1*b)/(R*T_1)) + (da_dT_1/(2*2**(1./2)*b))* \
math.log((Z_1+(1+2**(1./2))*((P_1*b)/(R*T_1)))/(Z_1+(1-2**(1./2))*((P_1*b)/(R*T_1))));
# The residual entropy at state 2 is
S_R_2 = R*math.log(Z_2 - (P_2*b)/(R*T_2)) + (da_dT_2/(2*2**(1./2)*b)) \
*math.log((Z_2+(1+2**(1./2))*((P_2*b)/(R*T_2)))/(Z_2+(1-2**(1./2))*((P_2*b)/(R*T_2))));
delta_S_R = S_R_2 - S_R_1; #[J/mol-K]
# The ideal gas change in entropy is
delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1); #[J/mol-K]
# Therefore
delta_S = delta_S_R + delta_S_ig; #[J/mol-K]
# Results
print " The enthalpy change is given by delta_H = %f kJ/mol"%(delta_H);
print " The entropy change is given by delta_S = %f J/mol-K"%(delta_S);
from scipy.optimize import fsolve
import math
# Variables
Vol = 0.15; #[m**(3)]
T_1 = 170; #[K] - Initial emperature
P_1 = 100; #[bar] - Initial pressure
P_1 = P_1*10**(5); #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
# For nitrogen
Tc = 126.2; #[K] - Critical tempeature
Pc = 34; #[bar] - Critical pressure
Pc = Pc*10**(5); #[Pa]
w = 0.038;
# Cp_0 = 27.2+4.2*10**(-3)*T
# Calculations and Results
#(1)
# For van der Walls equation of state
a = (27*R**(2)*Tc**(2))/(64*Pc); #[Pa-m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc); #[m**(3)/mol]
# The cubic form of van der Walls equation of state is given by,
# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0
# On simplification the equation changes to
# V**(3) - 1.799*10**(4)*V**(2) + 1.366*10**(-8)*V - 5.269*10**(-13) = 0
# Solving the cubic equation
def f(V):
return V**(3)-1.799*10**(-4)*V**(2) + 1.366*10**(-8)*V - 5.269*10**(-13)
V1 = fsolve(f,1)
V2 = fsolve(f,10)
V3 = fsolve(f,100)
# The above equation has only 1 real root, other two roots are imagimnary
V_1 = V1; #[m**(3)/mol]
# Thus total number of moles is given by
n_1 = Vol/V_1; #[mol]
# After 500 mol are withdrawn, the final number of moles is given by
n_2 = n_1 - 500; #[mol]
# Thus molar volume at final state is
V_2 = Vol/n_2; #[m**(3)/mol]
# The ideal entropy change is guven by
def f24(T):
return 27.2+4.2*10**(-3)*T
# delta_S_ig = quad(f24,T_1,T_2) - R*math.log(P_2/P_1)[0]
# The residual entropy change is given by
# delta_S_R = R*math.log((P_2*(V_2-b))/(R*T_2)) - R*math.log((P_1*(V_1-b))/(R*T_1))
# delta_S = delta_S_ig = delta_S_R
def f25(T):
return 27.2+4.2*10**(-3)*T
# delta_S = quad(f25,T_1,T_2) + R*math.log((V_2-b)/(V_1-b))[0]
# During discharging delta_S = 0, thus on simplification we get
# 18.886*math.log(T_2) + 4.2*10**(-3)*T_2 - 92.937 = 0
# Solving the above equation we get
def f1(T_2):
return 18.886*math.log(T_2) + 4.2*10**(-3)*T_2 - 92.937
T_2 = fsolve(f1,1)
# Thus at T_2,
P_2 = (R*T_2)/(V_2-b) - a/V_2**(2); #[N/m**(2)]
P_2 = P_2*10**(-5); #[bar]
print " 1).The final temperature is %f K"%(T_2);
print " The final pressure is %f bar"%(P_2);
#(2)
# In Peng-Robinson equation of state
m = 0.37464 + 1.54226*w - 0.26992*w**(2);
# At T_1 and P_1, we have
Tr = T_1/Tc;
alpha = (1 + m*(1 - Tr**(1./2)))**(2);
a_2 = ((0.45724*(R*Tc)**(2))/Pc)*alpha; #[Pa*m**(6)/mol**(2)]
b_2 = (0.07780*R*Tc)/Pc; #[m**(3)/mol]
# Cubuc form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f2(V):
return V**(3)+(b_2-(R*T_1)/P_1)*V**(2)-((3*b_2**(2))+((2*R*T_1*b_2)/P_1)-(a_2/P_1))*V+b_2**(3)+((R*T_1*(b_2**(2)))/P_1)-((a_2*b_2)/P_1)
V4 = fsolve(f2,-1)
V5 = fsolve(f2,0)
V6 = fsolve(f2,0.006)
#The above equation has only 1 real root,the other two roots are imaginary
V_1_2 = V6; #[m**(3)/mol]
# The number of moles in the initial state is given by
n_1_2 = Vol/V_1_2; #[mol]
# After 500 mol are withdrawn, the final number of moles is given by
n_2_2 = n_1_2 - 500; #[mol]
# Thus molar volume at final state is
V_2_2 = Vol/n_2_2; #[m**(3)/mol]
# At the final state the relation between pressure and temperature is
# P_2_2 = (R*T_2_2)/(V_2_2-b_2) - a_2/V_2_2**(2)
# P_2_2 = 7.23*10**(4)*T_2 - 3.93*10**(7)*a_2
# Now let us calculate the residual entropy at initial state
Z_1 = (P_1*V_1_2)/(R*T_1);
da_dT_1 = -(a*m)/((alpha*T_1*Tc)**(1./2)); #[Pa*m**(6)/mol**(2)] - da/dT
# The residual entropy change for Peng-Robinson equatiob of state is given by
# S_R = R*math.log(Z-(P*b)/(R*T)) + (da_dT/(2*2**(1/2)*b))*math.log((V+(1+2**(1/2))*b))/((V+(1-2**(1/2)*b))));
S_R_1 = R*(math.log(Z_1-(P_1*b_2)/(R*T_1))) + (da_dT_1/(2*2**(1.2)*b_2))*(math.log((V_1_2+(1+2**(1./2))*b_2)/(V_1_2+(1-2**(1./2))*b_2)));
# The total entropy change is given by
# delta_S = delta_S_ig + delta_S_R
def f26(T):
return 27.2+4.2*10**(-3)*T
# where, delta_S_ig = quad(f26,T_1,T_2_2) - R*math.log(P_2_2/P_1)[0]
# and, P_2_2 = (R*T_2_2)/(V_2_2-b_2) - a_2/V_2_2**(2)
# On simplification we get
# delta_S = 27.2*math.log(T_2_2-T_1) + 4.2*10**(-3)*(T_2_2-T_1) - R*math.log(P_2_2/P_1) + R*math.log(Z_2-(P_2_2*b)/(R*T_2_2)) + 6226*(da_dT_2) + 9.22
# Now we have the determine the value of T_2_2 such that delta_S = 0
# Starting with a temperature of 150 K
T_prime = 100.; #[K]
error = 10.;
while(error>0.1):
Tr_prime = T_prime/Tc;
alpha_prime = (1 + m*(1 - Tr_prime**(1./2)))**(2);
a_prime = ((0.45724*(R*Tc)**(2))/Pc)*alpha_prime;
P_prime = 7.23*10**(4)*T_prime - 3.93*10**(7)*a_prime;
Z_prime = (P_prime*V_2_2)/(R*T_prime);
da_dT_prime = -(a_prime*m)/((alpha_prime*T_prime*Tc)**(1./2));
delta_S = 27.2*math.log(T_prime/T_1) + 4.2*10**(-3)*(T_prime-T_1) - R*math.log(P_prime/P_1) + R*math.log(Z_prime-((P_prime*b_2)/(R*T_prime))) + 6226*(da_dT_prime) + 9.22;
error=abs(delta_S);
T_prime = T_prime + 0.3;
T_2_2 = T_prime; #[K] - Final temperature
P_2_2 = P_prime*10**(-5); #[bar] - Final pressure
print " 2).The final temperature is %f K"%(T_2_2);
print " The final pressure is %f bar"%(P_2_2);
from scipy.optimize import fsolve
import math
# Variables
T = 373.15; #[K]
Tc = 562.16; #[K]
Pc = 48.98; #[bar]
Pc = Pc*10**(5); #[Pa]
R = 8.314; #[J/mol-K] - Universal gas constant
# The cubic form of Redlich Kwong equation of state is given by,
# V**(3) - ((R*T)/P)*V**(2) - ((b_1**(2)) + ((b_1*R*T)/P) - (a/(T**(1/2)*P))*V - (a*b)/(T**(1/2)*P) = 0
# Calculations
a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc; #[Pa*m**(6)*K**(1/2)/mol]
b = (0.08664*R*Tc)/Pc; #[m**(3)/mol]
# At 373.15 K, let us assume the pressure to be 2.5 bar and under these conditions
P_1 = 2.5; #[bar]
P_1 = P_1*10**(5); #[bar]
# Putting the values in Redlich Kwong equation of state, the equation becomes
# V**(3) - 0.0124*V**(2) + 8.326*10**(-6)*V - 7.74*10**(-10) = 0
# Solving the cubic equation
def f(V):
return V**(3) - 0.0124*V**(2) + 8.326*10**(-6)*V - 7.74*10**(-10)
V1=fsolve(f,-9)
V2=fsolve(f,10)
V3=fsolve(f,0.1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq = V1; #[m**(3)/mol] - Molar volume in liquid phase
V_vap = V3; #[m**(3)/mol] - Molar volume in vapour phase
# Let us calculate the fugacity of vapour phase
# math.log(f_vap/P) = b/(V-b) + math.log((R*T)/(P*(V-b))) - (a/(R*T**(1.5)))*(1/(V+b) - (1/b)*math.log(V/(V+b)))
f_vap = P_1*math.exp(b/(V_vap-b) + math.log((R*T)/(P_1*(V_vap-b))) - (a/(R*T**(1.5)))*(1/(V_vap+b) - (1/b)*math.log(V_vap/(V_vap+b)))); #[Pa]
# Let us calculate the fugacity of the liquid phase
f_liq = P_1*math.exp(b/(V_liq-b) + math.log((R*T)/(P_1*(V_liq-b))) - (a/(R*T**(1.5)))*(1/(V_liq+b) - (1/b)*math.log(V_liq/(V_liq+b))));
# The two fugacities are not same; therefore another pressure is to be assumed. The new pressure is
P_new = P_1*(f_liq/f_vap); #[Pa]
# At P_new
def f1(V):
return V**(3) - ((R*T)/P_new)*V**(2) - (b**(2) + ((b*R*T)/P_new) - a/(T**(1/2)*P_new))*V - (a*b)/(T**(1/2)*P_new)
V4=fsolve(f1,-9)
V5=fsolve(f1,10)
V6=fsolve(f1,0.1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq_2 = V4; #[m**(3)/mol] - Molar volume in liquid phase
V_vap_2 = V6; #[m**(3)/mol] - Molar volume in vapour phase
f_vap_prime = P_new*math.exp(b/(V_vap_2-b) + math.log((R*T)/(P_new*(V_vap_2-b))) - (a/(R*T**(1.5)))*(1/(V_vap_2+b) - (1/b)*math.log(V_vap_2/(V_vap_2+b)))); #[Pa]
f_liq_prime = P_new*math.exp(b/(V_liq_2-b) + math.log((R*T)/(P_new*(V_liq_2-b))) - (a/(R*T**(1.5)))*(1/(V_liq_2+b) - (1/b)*math.log(V_liq_2/(V_liq_2+b))));
# Since the fugacities of liquid and vapour phasesare almost same the assumed pressure may be taken as vapour pressure at 373.15 K
P_new = P_new*10**(-5); #[bar]
# Results
print " The vapour pressure of benzene using Redlich Kwong equation of state is %f bar"%(P_new);
from scipy.optimize import fsolve
import math
# Variables
T = 150 + 273.15; #[K]
Tc = 647.1; #[K]
Pc = 220.55; #[bar]
Pc = Pc*10**(5); #[Pa]
w = 0.345;
R = 8.314; #[J/mol-K] - Universal gas constant
# Let us assume a pressure of 100 kPa.
P_1 = 100.*10**(3); #[Pa]
# Calculations
# At 100 kPa and 423.15 K, from Peng-Robinson equation of stste
m = 0.37464 + 1.54226*w - 0.26992*w**(2);
Tr = T/Tc;
alpha = (1 + m*(1 - Tr**(1./2)))**(2);
a = ((0.45724*(R*Tc)**(2))/Pc)*alpha; #[Pa*m**(6)/mol**(2)]
b = (0.07780*R*Tc)/Pc; #[m**(3)/mol]
# Cubic form of Peng-Robinson equation of stste is given by
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f(V):
return V**(3)+(b-(R*T)/P_1)*V**(2)-((3*b**(2))+((2*R*T*b)/P_1)-(a/P_1))*V+b**(3)+((R*T*(b**(2)))/P_1)-((a*b)/P_1)
V1 = fsolve(f,-1)
V2 = fsolve(f,0)
V3 = fsolve(f,1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq = V1; #[m**(3)/mol] - Molar volume in liquid phase
V_vap = V3; #[m**(3)/mol] - Molar volume in vapour phase
# The compressibility factor is given by
Z_vap = (P_1*V_vap)/(R*T); # For liquid phase
Z_liq = (P_1*V_liq)/(R*T); # For vapour phase
# The math.expression for fugacity of Peng Robinson equation is
# math.log(f/P) = (Z-1) - math.log(Z-((P*b)/(R*T))) - (a/(2*2**(1/2)*b*R*T))*math.log((Z+(1+2**(1/2))*((P*b)/(R*T)))/((Z+(1-2**(1/2))*((P*b)/(R*T)))
# For vapour phase
f_P_vap = math.exp((Z_vap-1) - math.log(Z_vap-((P_1*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_vap+(1+2**(1./2))*((P_1*b)/(R*T)))/(Z_vap+(1-2**(1./2))*((P_1*b)/(R*T)))));
# For liquid phase
f_P_liq = math.exp((Z_liq-1) - math.log(Z_liq-((P_1*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_liq+(1+2**(1./2))*((P_1*b)/(R*T)))/(Z_liq+(1-2**(1./2))*((P_1*b)/(R*T)))));
# Therefore f_liq/f_vap can be calculated as
fL_fV = (f_P_liq/f_P_vap);
# The two values (f/P)_vap and (f/P)_vap are not same [ (f_P_liq/f_P_vap) >1 ]; therefore another pressure is to be assumed. The new pressure be
P_new = P_1*(f_P_liq/f_P_vap); #[Pa]
# At P_new and 423.15 K, from Peng-Robinson equation of stste
# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;
# Solving the cubic equation
def f(V):
return V**(3)+(b-(R*T)/P_new)*V**(2)-((3*b**(2))+((2*R*T*b)/P_new)-(a/P_new))*V+b**(3)+((R*T*(b**(2)))/P_new)-((a*b)/P_new)
V4 = fsolve(f,-1)
V5 = fsolve(f,0)
V6 = fsolve(f,1)
# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.
V_liq_2 = V4; #[m**(3)/mol] - Molar volume in liquid phase
V_vap_2 = V6; #[m**(3)/mol] - Molar volume in vapour phase
# The compressibility factor is given by
Z_vap_2 = (P_new*V_vap_2)/(R*T); # For liquid phase
Z_liq_2 = (P_new*V_liq_2)/(R*T); # For vapour phase
# For vapour phase
f_P_vap_2 = math.exp((Z_vap_2-1) - math.log(Z_vap_2-((P_new*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_vap_2+(1+2**(1./2))*((P_new*b)/(R*T)))/(Z_vap_2+(1-2**(1./2))*((P_new*b)/(R*T)))));
# For liquid phase
f_P_liq_2 = math.exp((Z_liq_2-1) - math.log(Z_liq_2-((P_new*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_liq_2+(1+2**(1./2))*((P_new*b)/(R*T)))/(Z_liq_2+(1-2**(1./2))*((P_new*b)/(R*T)))));
# Therefore f_liq/f_vap can be calculated as
fL_fV_2 = (f_P_liq_2/f_P_vap_2);
# And new pressure is given by
P_new_prime = P_new*(f_P_liq_2/f_P_vap_2); #[Pa]
P_new_prime = P_new_prime*10**(-5);
# Since the change in pressure is small, so we can take this to be the vapour pressure at 150 C
# Results
print " The vapour pressure of water using Peng-Robinson equation of stste is %f bar"%(P_new_prime);