from numpy import ones,array,sqrt,round
import math
# Variables
T = 573.15; #[K]
P = array([700,600,500,400,300,200]); #[KPa]
#values for H,V,S for various P from steam tables
H = array([3059.8,3020.4,2975.71,2923.5,2859.9,2777.35]); #[KJ/Kg]
V = array([371.39,418.25,481.26,571.23,711.93,970.04]); #[cm**3/g]
S = 7.29997*ones(6); #[KJ/Kg/K] Isentropic
u0 = 30.; #[m/s]
# Calculations
#u**2 = u1**2-2(H-H1)
u = round(sqrt((u0**2-2*(H-H[0])*10.**3)),1)
#Umath.sing Eq(2.27)
#A/A1 = u1*V/V1*u;
c = u[0]/V[0];
K = round((c*V/u),3); #K = A/A1 c = u1/V1
Ans = [P,V,u,K];
print ' P/[KPa] V/[cm**3/g] u/[m/s] A/A1 \n',Ans
import math
# Variables
T1 = 573.15; #[K]
R = 8314;
P1 = 700; #[KPa]
M = 18.015;
Gamma = 1.3;
u0 = 30; #[m/s]
# Calculations and Results
#(a)
#Umath.sing Eqn(7.12)
#K = P2/P1 = (2/(Gamma+1))**(Gamma/(Gamma-1))
K = round((2/(Gamma+1))**(Gamma/(Gamma-1)),2); #rounding to 2 decimal places
P1V1 = round(R*T1/M); #m**2/s**2
#Umath.sing Eqn(7.11)
#u_throat**2 = u**2+2(Gamma)(P1V1)/(Gamma-1)[1-(P2/P1)**((Gamma-1)/Gamma))]
u_throat = round(math.sqrt(u0**2+((2*Gamma*P1V1/(Gamma-1))*(1-(K**((Gamma-1)/Gamma))))),2);
print '(a)Critical Pressure ratio(P2/P1)',K
print ' Velocity at the throat',u_throat,' m/s'
#(b)Mach No 2.0
u = 2*u_throat;
K = (1-((u**2-u0**2)*(Gamma-1)/(2*Gamma*P1V1)))**(Gamma/(Gamma-1)); #K = P2/P1
P2 = round(K*P1);
print '(b)Discharge Pressure for Mach Number of 2.0',P2,'KPa'
import math
def HRB(Tr,Pr,omega):
B0 = 0.083-(0.422/(Tr**1.6));
diffr_B0 = 0.675/(Tr**2.6); #dB0/dTr
B1 = 0.139-(0.172/(Tr**4.2));
diffr_B1 = 0.722/(Tr**5.2); #dB0/dTr
return Pr*(B0-(Tr*diffr_B0)+(omega*(B1-(Tr*diffr_B1))));
def SRB(Tr,Pr,omega):
B0 = 0.083-(0.422/(Tr**1.6));
diffr_B0 = 0.675/(Tr**2.6); #dB0/dTr
B1 = 0.139-(0.172/(Tr**4.2));
diffr_B1 = 0.722/(Tr**5.2); #dB0/dTr
return -Pr*(diffr_B0+(omega*diffr_B1));
# Variables
P1 = 20.; #[bar]
T = 400.; #[K]
P2 = 1.; #[bar]
R = 8.314;
# Calculations
#using Eq(6.84)
#del_H = Cp(T2-T1)+Hr2-Hr1 = 0 but Hr2 = 0
#T2 = Hr1/Cp + T1
Tc = 369.8; #[K]
Pc = 42.48; #[bar]
omega = 0.152;
a = T; #Initial
for i in range(2):
Tr = a/Tc
Pr = P1/Pc;
Hr1 = R*Tc*HRB(Tr,Pr,omega); #[J/mol]
Cp = R*(1.213+(28.785*10**-3*a)-(8.824*10**-6*a*a)); #[J/mol/K]
T2 = (Hr1/Cp)+a;
Tm = (a+T2)/2;
i = i+1;
a = Tm;
Tm = a;
T2 = round(Tm) #[K]
Tr = T/Tc;
Sr = R*SRB(Tr,Pr,omega);
del_S = round((Cp*math.log(T2/T))-(R*math.log(P2/P1))-Sr,2);
# Results
print 'Entropy',del_S,'J/mol/K'
print ('Positive Entropy represents the irreversibility of Throttling Process')
print 'Final Temperature',T2,'K'
# Variables
P1 = 8600.; #[KPa]
T1 = 773.15; #[K]
#values of Enthalpy and Entropy from Steam tables
H1 = 3391.6; #[KJ/Kg]
S1 = 6.6858; #[KJ/Kg/K]
eta = 0.75;
P2 = 10000.; #[KPa]
rW = 56400.; #[KW] or [KJ/s]
S2i = S1; #Isentropic
S2_liquid = 0.6493;
S2_vapor = 8.1511;
H2_liquid = 191.8;
H2_vapor = 2584.8;
# Calculations
x2 = (S2i-S2_liquid)/(S2_vapor-S2_liquid);
H2i = H2_liquid+(x2*(H2_vapor-H2_liquid));
del_Hs = H2i-H1; #[KJ/Kg]
del_H = eta*del_Hs;
H2 = round(H1+del_H,0); #[KJ/Kg]
x2 = (H2-H2_liquid)/(H2_vapor-H2_liquid);
S2 = round(S2_liquid+(x2*(S2_vapor-S2_liquid)),4);
rm = round(-rW/(H2-H1),2); #[Kg/s]
# Results
print 'Enthalpy',H2,'KJ/Kg'
print 'Entropy',S2,'KJ/Kg/K'
print 'Rate of mass change',rm,'Kg/s'
import math
#To find Approx Value
def approx(V,n):
return (V*10**n)/10**n; #V-Value n-To what place
def MCPS(T0,T,A,B,C,D):
t = T/T0;
if t == 1:
t = 0.9997383
return (A)+(((B*T0)+(((C*T0*T0)+(D/(t*t*T0*T0)))*(t+1)/2.))*((t-1)/math.log(t)))
def MCPH(T0,T,A,B,C,D):
t = T/T0;
return (A+((B/2)*T0*(t+1))+((C/3)*T0*T0*((t**2)+t+1))+(D/(t*T0*T0)))
def HRB(Tr,Pr,omega):
B0 = 0.083-(0.422/(Tr**1.6));
diffr_B0 = 0.675/(Tr**2.6); #dB0/dTr
B1 = 0.139-(0.172/(Tr**4.2));
diffr_B1 = 0.722/(Tr**5.2); #dB0/dTr
return Pr*(B0-(Tr*diffr_B0)+(omega*(B1-(Tr*diffr_B1))));
def SRB(Tr,Pr,omega):
B0 = 0.083-(0.422/(Tr**1.6));
diffr_B0 = 0.675/(Tr**2.6); #dB0/dTr
B1 = 0.139-(0.172/(Tr**4.2));
diffr_B1 = 0.722/(Tr**5.2); #dB0/dTr
return -Pr*(diffr_B0+(omega*diffr_B1));
# Variables
T1 = 573.15; #[K]
P1 = 45. #[bar]
P2 = 2. #[bar]
Tc = 282.3 #[K]
Pc = 50.4; #[bar]
omega = 0.087;
A = 1.424;
B = 14.394*10**-3;
C = -4.392*10**-6;
D = 0;
R = 8.314;
# Calculations and Results
#Using Eqn(6.84)
#del_H = <Cp>h (T2-T1)+Hr2-Hr1
#Umath.sing Eqn(6.85))
#del_S = <Cp>s ln(T2/T1) - R*ln(P2/P1)+Sr2-Sr1
#(a) equations for Ideal gas
#No residuals terms, whence
#del_H = <Cp>h(T2-T1)
#del_S = <Cp>s ln(T2/T1) - R*ln(P2/P1)
del_S = 0 #isentropic
#Whence K = <Cp>s/R ln(T2/T1) = ln(P2/P1)
K = math.log(P2/P1);
#let c = <Cp>s/R
#T2 = exp(K/c+ln(T1))
i = -1;
a = (T1); #Initial
while (i == -1):
b = MCPS(T1,a,A,B,C,D);
temp = math.exp((K/b)+math.log(T1));
flag = a-temp;
if(flag<= 0.1):
T2 = a;
i = 1;
else:
a = temp-0.1;
i = -1;
print ('(a)by Equations for an Ideal gas')
print ('K',approx(T2,1),'Temp = ')
Cp_h = R*MCPH(T1,T2,A,B,C,D);
del_Hs = Cp_h*(T2-T1);
Ws_a = approx(del_Hs,0);
print ('J/mol',Ws_a,'Work')
#(b)-Appropriate Generalized correlations
Tr1 = T1/Tc;
Pr1 = P1/Pc;
Hr1 = R*Tc*HRB(Tr1,Pr1,omega); #[J/mol]
Sr1 = R*SRB(Tr1,Pr1,omega); #[J/mol/K]
Tr2 = T2/Tc;
Pr2 = P2/Pc;
Sr2 = R*SRB(Tr2,Pr2,omega);
#Using Eqn(6.85))
#del_S = <Cp>s ln(T2/T1) - R*ln(P2/P1)+Sr2-Sr1
#del_S = 0 isentropic
#K = <Cp>s ln(T2/T1) = Rln(P2/P1)-Sr2+Sr1
K = R*math.log(P2/P1)-Sr2+Sr1;
#T2 = exp((K/<Cp>s)+ln T1)
i = -1;
a = (T1); #Initial
while (i == -1):
b = R*MCPS(T1,a,A,B,C,D);
temp = math.exp((K/b)+math.log(T1));
flag = a-temp;
if(flag <= 0.1):
T2 = a;
i = 1;
else:
a = temp-0.1;
i = -1;
print ('(b)by Appropriate generalized correlations')
print 'Temp = ',approx(T2,1),'K'
Tr2 = T2/Tc;
Sr2 = R*SRB(Tr2,Pr2,omega); #[J/mol/K]
Hr2 = R*Tc*HRB(Tr2,Pr2,omega); #[J/mol]
Cp_h = R*MCPH(T1,T2,A,B,C,D);
del_Hs = Cp_h*(T2-T1)+Hr2-Hr1;
Ws_b = approx(del_Hs,-1);
print 'Work',Ws_b,'J/mol'
# Program to Find the Work Required and Properties of Discharge Steam
import math
#To find Approx Value
def approx(V,n):
return (V*10**n)/10**n; #V-Value n-To what place
# Variables
P1 = 100.; #[KPa] (Tsat/tsat) = 327.78K/99.63`C)
#From Steam Tables @ 100KPa
S1 = 7.3598; #[KJ/Kg/K]
H1 = 2675.4; #[KJ/Kg]
P2 = 300.; #[KPa]
#From Steam Tables @ 300KPa
S2 = S1; #Isentropic
H2i = 2888.8; #[KJ/Kg]
eta = 0.75; #Efficiency
# Calculations
del_H = H2i-H1;
del_H = del_H/eta;
H2 = approx(H1+del_H,1); #[KJ/Kg]
#From Steam Tables w.r.t H2
T2 = 519.25; #[K]
S2 = 7.5019; #[KJ/Kg/K]
Ws = approx(del_H,1); #[KJ/Kg] Work Reqd
# Results
print 'Enthalpy',H2,'KJ/Kg'
print 'Entropy',S2,'KJ/Kg/K'
print 'Temperature',T2,'K'
print 'Work Done',Ws,'KJ/Kg'
# Program to Find Work Reqiured and Discharge Temperature of Methane
import math
#To find Approx Value
def approx(V,n):
return (V*10**n)/10**n; #V-Value n-To what place
def MCPH(T0,T,A,B,C,D):
t = T/T0;
return (A+((B/2)*T0*(t+1))+((C/3)*T0*T0*((t**2)+t+1))+(D/(t*T0*T0)))
def MCPS(T0,T,A,B,C,D):
t = T/T0;
if t==1:
t = 0.9997383
return (A)+(((B*T0)+(((C*T0*T0)+(D/(t*t*T0*T0)))*(t+1)/2))*((t-1)/math.log(t)))
# Variables
R = 8.314;
T1 = 293.15; #[K]
P1 = 140.; #[KPa]
P2 = 560.; #[KPa]
eta = 0.75; #[Efficiency]
A = 1.702;
B = 9.081*10**-3;
C = -2.164*10**-6;
D = 0;
# Calculations
i = -1;
a = (T1); #Initial
while (i == -1):
b = MCPS(T1,a,A,B,C,D);
b = b**-1;
c = T1*((P2/P1)**b);
flag = c-a;
if(flag <= 0.0001):
T2i = a;
i = 1;
else:
a = a+0.01;
i = -1;
Cps = R*MCPS(T1,T2i,A,B,C,D);
Cph = approx(R*MCPH(T1,T2i,A,B,C,D),3);
#from Eqn(7.19)
Ws = approx(Cph*(T2i-T1),0) #[J/mol]
Ws = approx(Ws/eta,0) #Actual work
del_H = Ws;
#From eqn(7.21) Actual discharge Temperature
#T2 = T1+(del_H/Cph)
i = -1;
a = (T2i); #Initial
chk = 1;
while (i == -1):
b = R*MCPH(T2i,a,A,B,C,D);
c = del_H/(a-T1);
flag = c-b;
if(flag<= 0.001):
T2 = a;
i = 1;
else:
a = a+0.001;
i = -1;
Cph_T2 = approx(R*MCPH(T2i,T2,A,B,C,D),2);
# Results
print 'Temperature',T2,'K'
print 'Enthalpy',Cph_T2,'J/mol/K'
print 'Actual Work',Ws,'J/mol'
print ('Note: The answer in the Book varies with that of this code because the Calculation\
in the Book does not leads to the answer given')
# Program to Find Work,Temperature Change and Entropy Change in Pump
import math
#To find Approx Value
def approx(V,n):
return (V*10**n)/10**n; #V-Value n-To what place
# Variables
T1 = 318.15; #[K]
P1 = 10.; #[KPa]
P2 = 8600.; #[KPa]
eta = 0.75; #Efficiency
#Properties of saturated liquid water @ 318.15K
V = 1010.; #[cm**3/Kg]
V = 1010.*10**-6; #[m**3/Kg]
Beta = 425*10**-6; #[K**-1]
Cp = 4.178; #[KJ/Kg/K]
# Calculations
#From Eqn(7.24)
Ws = V*(P2-P1); #[KPa m**3/Kg]
del_H = Ws;
#From Eqn(7.17)
del_H = del_H/eta;
Ws = approx(del_H,2);
#From Eqn(7.25)
del_T = approx((del_H-(V*(1-(Beta*T1))*(P2-P1)))/Cp,2);
#From Eqn(7.26)
T2 = T1+del_T;
del_S = approx(Cp*math.log(T2/T1)-(Beta*V*(P2-P1)),3);
# Results
print 'Work Done',Ws,'KJ/Kg'
print 'Change in Temperature',del_T,'K'
print 'Change in Entropy',del_S,'KJ/Kg/K'