In [1]:

```
from scipy.optimize import fsolve
import math
# Variables
V_vessel = 4.*10**(-3); #[m**(-3)] - Volume of vessel
T = 200+273.15; #[K] - Temperature
R = 8.314; #[J/mol*K] - Universal fas constant
P = 1.5*10**(6); #[Pa] - Pressure
Q = 40.*1000; #[J] - Heat input
# From steam table at 200 C,Psat=1.55549 MPa,therefore the steam is superheated.
# Calculations and Results
# (1)
# Using steam table,at 1.5 MPa and 200 C,
V_1 = 0.1325; #[m**(3)/mol] - Specific volume
U_1 = 2598.1; #[kJ/kg] - Specific internal energy
# From first law under constant pressure,
# Q - m*P*(V2 - V1) = m*(U2 - U1)
m = V_vessel/V_1; #[kg] - Mass of system
# Putting the values,the above equation becomes
# 45283*(V2 - 0.1325) + 30.1887*(U2 - 2598.1) = 40000
# From steam table at 700 C LHS is 33917.0 and at 800 C,it is 40925.3.
# Therefore the final temperature lies between 700 and 800 C
print " 1.From steam table the final temperature lies between 700 and 800 C";
# Alternate method
# Here we use first law at constant pressure,
# Q = m*(H_2 - H_1)
H_1 = 2796.8; #[kJ/kg]
# Substituting the values,
# 40 = 0.0301887*(H_2 - 2796.8)
H_2 = (40/0.0301887) + 2796.9; #[kJ/kg]
# Threfore,final enthalpy is (H2) 4121.8 [kJ/kg] and pressure is 1.5 [MPa].
# From steam table at 1.5 [MPa]and 700 C,enthalpy is 3920.3 [kj/kg] and at 1.5 [MPa]and 800 C,enthalpy is 4152.6 [kj/kg]
print "\tAlternate method";
print "\tBy linear interpolation we get the temperature at which enthlpy is 4121.8 kJ/kg to be 786.74 C";
# (2)
# Assuming ideal behaviour.
n = (P*V_vessel)/(R*T); #[mol] - No of moles
M = 18.015; # Molecular weight of water
m_2 = n*M; #[g] - mass of moles
Cp_1 = 7.7 + 0.04594*10**(-2)*T + 0.2521*10**(-5)*T**(2) - 0.8587*10**(-9)*T**(3); #[cal/mol*K] - Heat capacity at constant presure
R0 = 1.987; #[cal/mol*K] - universal gas constant
Cv_1 = Cp_1 - R0; #[cal/mol*K] - Heat capacity at constant volume
Cv_1 = Cv_1*4.184/M; #[J/g*K]
T1 = T;
# From 1st law energy balance for constant pressure, we have Q-W= m*(delta_U)
# Q = P*(V2 - V1)*m = m*Cv*(T2 - T1)
# Q = P*((T2/T1)-1)*V1*m = m*Cv*(T2-T1)
# But, (V1*Cv)=initial total volume of system = V_vessel
# Q-P((T2/T1)-1)*V_vessel = m_2*Cv_0*(T2-T1);
def f(T2):
return Q-P*((T2/T1)-1)*V_vessel-m_2*Cv_1*(T2-T1)
T2_1 = fsolve(f,1)
#The heat capacity should be evaluted at mean temperature
T_mean = (T1 + T2_1)/2;
Cp_2 = 7.7 + 0.04594*10**(-2)*T_mean+0.2521*10**(-5)*T_mean**(2) - 0.8587*10**(-9)*T_mean**(3); #[cal/mol*K] - Heat capacity at constant presure
Cv_2 = Cp_2-R0; #[cal/mol*K] - Heat capacity at constant volume
Cv_2 = Cv_2*4.184/M; #[J/g*K]
#Now again solving the equation Q=P*((T2/T1)-1)*V1*m = m*Cv*(T2-T1),for Cv=Cv_2
def f1(T2):
return Q-P*((T2/T1)-1)*V_vessel-m_2*Cv_2*(T2-T1)
T2_2 = fsolve(f1,1)
print " 2).The temperature assuming ideal behaviour is %f K"%(T2_2);
# Alternate method
# From 1st law at constant pressure, we have Q = m*Cp(T2-T1)
T2_3 = Q/(m_2*(Cp_1*4.184/M))+T1;
#We can calculate the mean temperature as done above
T_mean1 = (T1 + T2_3)/2; #[J/g*K]
#The heat capacity should be evaluted at mean temperature
Cp_3 = 7.7 + 0.04594*10**(-2)*T_mean1 + 0.2521*10**(-5)*T_mean1**(2)-0.8587*10**(-9)*T_mean1**(3); #[cal/mol*K] - Heat capacity at constant presure
Cp_3 = Cp_3*4.184/M; #[J/g*K]
# Again solving the equation Q = m*Cp(T2-T1), for Cp=Cp_3
T2_4 = Q/(m_2*Cp_3) + T1;
print "\tAlternate method";
print "\tThe temperature assuming ideal behaviour alternate method) is %f K"%(T2_4);
```

In [2]:

```
import math
# Variables
V_tank = 1; #[m**(3)] - Volume of the tank
V_liq = 0.05; #[m**(3)] - Volume of saturated water
V_vap = 0.95; #[m**(3)] - Volume of saturated vapour
P = 1; #[bar] - Pressure
V_liq_sat = 0.001043; #[m**(3)/kg] - Specific volume of saturated water
V_vap_sat = 1.6940; #[m**(3)/kg] - Specific volume of saturated vapour
U_liq_sat = 417.4; #[kJ/kg] - Saturated liquid internal energy
U_vap_sat = 2506.1; #[kJ/kg] - Saturated vapour internal energy
# Calculations
m = (V_liq/V_liq_sat) + (V_vap/V_vap_sat); #[kg] - Total mass of water
U_1 = (V_liq/V_liq_sat)*U_liq_sat + (V_vap/V_vap_sat)*U_vap_sat; #[kJ] - Total internal energy
# At final state,which is saturated vapour
V = V_tank/m; #[m**(3)/kg] - Molar volume
# From saturated steam table at 8 MPa,as reported in the book V_vap = 0.02352[m**(3)/kg] and U_vap = 2569.8[kJ/kg]
# At 9 MPa, Vv = 0.02048[m**(3)/kg] and Uv = 2557.8[kJ/kg]
# Therefore final state pressure of the system (from interpolation) is 8.954 [MPa] and internal energy of saturated vapour is 2558.35 [kJ/kg]
U_2 = m*2558.35; #[kJ] - Final total internal energy
del_Ut = U_2 - U_1; #[kJ]
#we have, del_U = Q - W
#Here work done is zero because volume is rigid.
Q = del_Ut; #[kJ]
Q = del_Ut*10**(-3); #[MJ]
# Results
print " The amount of heat to be added is %f MJ"%( Q);
```

In [3]:

```
# Variables
M_vap_sat = 0.22; #[kg] - mass of saturated vapour
M_liq_sat = 1.78; #[kg] - mass of saturated liquid
P = 700; #[kPa] - Pressure
#At P=700 kPa,the systen is saturated,from steam table as reported in the book
T_sat1 = 164.97; #[C]
V_liq_1 = 0.001108; #[m**(3)/kg]
V_vap_1 = 0.2729; #[m**(3)/kg]
# Calculations and Results
Vt_1 = V_liq_1*M_liq_sat + V_vap_1*M_vap_sat; #[m**(3)] - total volume
#At final state,P = 8 MPa
T_sat2 = 295.06; #[C]
V_liq_2 = 0.001384; #[m**(3)/kg]
V_vap_2=0.02352; #[m**(3)/kg]
Vt_2 = Vt_1; # Since the volume is rigid.
# Since the volume of 2 kg of vapour is 0.062 [m**(3)]
V = Vt_2/2; #[m**(3)/kg] - specific volume
# (a)
# From steam table at 8 [MPa]and 350 [C],V=0.02995[m**(3)/kg];
V_1 = 0.02995; #[m**(3)/kg]
# And at 8 [MPa]and 400 [C],
V_2 = 0.03432; #[m**(3)/kg]
# By interpolation,
T = ((V-V_1)/(V_2 - V_1))*(400-350)+350;
print " a).The final temperature is %f c"%(T);
# (b)
# From steam table
U_1 = 2747.7; #[kJ/kg]
H_1 = 2987.3; #[kJ/kg]
# And at 8 [MPa]and 400 C,
U_2 = 2863.8; #[kJ/kg]
H_2 = 3138.3; #[kJ/kg]
# Therefore at T = 362.01 C
U = U_1+((U_2 - U_1)/(400 - 350))*(T - 350);
print " b).The internal energy is %f kJ/kg"%(U);
#(c)
H = H_1+((H_2 - H_1)/(400 - 350))*(T - 350);
print " b).The enthalpy is %f kJ/kg"%(H);
```

In [4]:

```
from scipy.optimize import fsolve
import math
from scipy.integrate import quad
# Variables
T = 300; #[K] - Temperature
P1 = 1; #[bar] - Initial pressure
P1 = P1*10**(5); #[N/m**(2)]
P2 = 8; #[bar] - Final pressure
P2 = P2*10**(5); #[N/m**(2)]
R = 8.314; #[J/mol*K] - Universal gas constant
Tc = 126.2; #[K] - Critical temperature
Pc = 34; #[bar] - Critical pressure
Pc = Pc*10**(5); #[N/m**(2)]
w = 0.038; # Acentric factor
# Calculations and Results
# w = integral(Pdv)
# Z = 1 + (B/V)
# (P*V)/(R*T) = 1 + (B/V)
# P = (R*T)/V + (B*R*T)/V**(2)
def f29(V):
return (R*T/V) + (B*R*T)/V**(2)
# w = quad(f29,V1,V2)[0]
# Under isothermal conditions,
# w = R*T*math.log(V2/V1) - B*R*T*((1/V2) - (1/V1));
Tr = T/Tc;
B_0 = 0.083 - (0.422/(Tr)**(1.6));
B_1 = 0.139 - (0.172/(Tr)**(4.2));
B = ((B_0+(w*B_1))*(R*Tc))/Pc; #[m**(3)/mol]
# Now we have to calculate molar volume i.e V1 and V2 at given conditions
# At state 1,
def f(V):
return V**(2)-(R*T/P1)*V-(B*R*T)/P1
V_1 = fsolve(f,-1)
V_2 = fsolve(f,1)
# We will take root near to (R*T)/P, i.e V_2
V1 = V_2;
# At state 2,
def f1(V):
return V**(2)-(R*T/P2)*V-(B*R*T)/P2
V_3=fsolve(f1,-1)
V_4=fsolve(f1,1)
V2 = V_4;
# The work done is thus,
w = R*T*math.log(V2/V1) - B*R*T*((1/V2) - (1/V1)); #[J]
w = w*10**(-3); #[kJ]
print " The work done is %f kJ/mol"%(w);
print " Negative sign indicates that work is done on the gas";
```

In [5]:

```
# Variables
T = 300; #[K] - Temperature
P1 = 1; #[bar] - Initial pressure
P1 = P1*10**(5); #[N/m**(2)]
P2 = 8; #[bar] - Final pressure
P2 = P2*10**(5); #[N/m**(2)]
R = 8.314; #[J/mol*K] - Universal gas constant
y1 = 0.21; # Mole fraction of component 1 (oxygen)
y2 = 0.79; # Mole fraction of component 1 (nitroen)
# For component 1 (Oxygen)
Tc_1 = 154.6; #[K]
Pc_1 = 50.43*10**(5); #[N/m**(2)]
Vc_1 = 73.4; #[cm**(3)/mol]
Zc_1 = 0.288;
w_1 = 0.022;
#For component 2 (Nitrogen)
Tc_2 = 126.2; #[K]
Pc_2 = 34*10**(5); #[N/m**(2)]
Vc_2 = 89.2; #[cm**(3)/mol]
Zc_2 = 0.289;
w_2 = 0.038;
# 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)]
# 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 we have to calculate molar volume i.eV1 and V2 at given conditions
# At state 1,
def f(V):
return V**(2)-(R*T/P1)*V-(B*R*T)/P1
V_1=fsolve(f,-1)
V_2=fsolve(f,1)
# We will take root near to (R*T)/P, i.e V_2
V1 = V_2; #[m**(3)/mol]
# At state 2,
def f1(V):
return V**(2)-(R*T/P2)*V-(B*R*T)/P2
V_3=fsolve(f1,-1)
V_4=fsolve(f1,1)
V2 = V_4; #[m**(3)/mol]
# Work done per mole of air is given by, w=integral(Pdv)
# Z = 1 + (B/V)
# (P*V)/(R*T) = 1 +( B/V)
# P = (R*T)/V+(B*R*T)/V**(2)
def f43(V):
return (R*T/V)+(B*R*T)/V**(2)
# w = quad(f43,V1,V2)[0]
# Under isothermal conditions,
w = R*T*math.log(V2/V1)-B*R*T*((1/V2)-(1/V1));
w = w*10**(-3); #[kJ/mol]
# Results
print " The work done is %f kJ/mol"%(w);
```

In [6]:

```
# Variables
T = 125+273.15; #[K] - Temperature
P1 = 1; #[bar] - Initial pressure
P1 = P1*10**(5); #[N/m**(2)]
P2 = 60; #[bar] - Final pressure
P2 = P2*10**(5); #[N/m**(2)]
R = 8.314; #[J/mol*K] - Universal gas constant
Tc = 416.3; #[K] - Critical temperature
Pc = 66.80*10**(5); #[N/m**(2)] - Critical pressure
# (1)
# Virial equation of state, Z = 1 + (B/V)+(C/V**(2))
# (P*V)/(R*T) = 1 + (B/V)+(C/V**(2))
# P = (R*T)/V+(B*R*T)/V**(2)+(C*R*T)/V**(3)
# w = integral(PdV)=R*T*math.log(V2/V1)-(B*R*T)*(1/V2-1/V1)-(C*R*T/2)*(1/V2**(2)-1/V1**(2))
B = -207.5; #[cm**(3)/mol] - Second virial coefficient
B = -207.5*10**(-6); #[m**(3)/mol]
C = 18200; #[cm**(6)/mol**(2)] - Third virial coefficient
C = 18200*10**(-12); #[m**(6)/mol**(2)]
# Calculations and Results
# We need to calculate molar volume at state 1 and 2,
# At state 1,P = P1,
# V**(3)-(R*T/P)*V**(2)-(B*R*T/P)*V-(C*R*T/P)=0
# Solving the cubic equation
def f1(V):
return V**(3)-(R*T/P1)*V**(2)-(B*R*T/P1)*V-(C*R*T/P1)
V_1=fsolve(f1,-1)
V_2=fsolve(f1,0)
V_3=fsolve(f1,10)
# The cubic equation has only 1 real root,other two roots are imaginary.
V1 = V_3;
# Similarly at state 2,P=P2
# V**(3) - (R*T/P)*V**(2) - (B*R*T/P)*V - (C*R*T/P) = 0
# Solving the cubic equation
def f2(V):
return V**(3)-(R*T/P2)*V**(2)-(B*R*T/P2)*V-(C*R*T/P2)
V_4=fsolve(f2,-1)
V_5=fsolve(f2,0)
V_6=fsolve(f2,1)
V2 = V_6;
# Finally work done is given by,
w1 = R*T*math.log(V2/V1)-(B*R*T)*(1/V2-1/V1)-(C*R*T/2)*(1/V2**(2)-1/V1**(2)); #[J/mol]
w1 = w1*10**(-3); #[kJ/mol]
print " 1).The work done using given virial equation of state is %f kJ/mol"%(w1);
# (2)
# Virial equation of state, Z = 1+(B*P)/(R*T)+((C-B**(2))/(R*T)**(2))*P**(2)
# (P*V)/(R*T)= 1+(B*P)/(R*T)+((C-B**(2))/(R*T)**(2))*P**(2)
# V = (R*T)/P+B+((C-B**(2))/(R*T))*P
# Differentiating both sides by P and integrating we get,
# w = integral(PdV)=-(R*T)*math.log(P2/P1)+((C-B**(2))/(2*R*T))*(P2**(2)-P1**(2))
w2 = -(R*T)*math.log(P2/P1) + ((C-B**(2))/(2*R*T))*(P2**(2)-P1**(2)); #[J/mol]
w2 = w2*10**(-3); #[kJ/mol]
print " 2).The work done using given virial equation of state is %f kJ/mol"%(w2);
# (3)
# Van der Walls equation of state is given by,
a = (27*(R**(2))*(Tc**(2)))/(64*Pc); #[Pa*m**(6)/mol**(2)]
b = (R*Tc)/(8*Pc); #[m**(3)/mol]
# P = ((R*T)/(V-b))-a/(V**(2)); #[N/m**(2)]
# w = integral(PdV)=R*T*math.log((V2-b)/(V1-a))+a*(1/V2-1/V1)
# 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
# Solving the cubic equation for P=P1
def f3(V):
return V**(3)-(b+(R*T)/P1)*V**(2)+(a/P1)*V-(a*b)/P1
V2_1=fsolve(f3,1)
V2_2=fsolve(f3,10)
V2_3=fsolve(f3,100)
# The above equation has 1 real and 2 imaginary roots. We consider only real root (V2_3).
# Similarly at state 2,for P=P2,
def f4(V):
return V**(3)-(b+(R*T)/P2)*V**(2)+(a/P2)*V-(a*b)/P2
V2_4=fsolve(f4,1)
V2_5=fsolve(f4,10)
V2_6=fsolve(f4,100)
# The above equation has 1 real and 2 imaginary roots. We consider only real root (V2_6).
# Finally work done is given by
w3 = R*T*math.log((V2_6-b)/(V2_3-b))+a*(1/V2_6-1/V2_3); #[J/mol]
w3 = w3*10**(-3); #[kJ/mol]
print " 3).The work done using van der Walls equation of state is %f kJ/mol"%(w3);
# (4)
# Redlich Kwong equation of state,
a_1 = (0.42748*(R**(2))*(Tc**(2.5)))/Pc; #[Pa*m**(6)*K**(1/2)/mol]
b_1 = (0.08664*R*Tc)/Pc; #[m**(3)/mol]
# P = ((R*T)/(V-b_1))-(a_1/(T**(1/2)*V*(V+b_1))); #[N/m**(2)]
# Work done is given by
def f4(V):
return (V2-b)/(V1-b)-a/T**(1./2)*integrate(1/V*(V+b))
# w = R*T*math.log((V2-b)/(V1-b))-a/T**(1/2)* quad(f4,V1,V2)[0]
# Using the factorization 1/(V*(V+b))=(1/b)*((1/V)-(1/V+b)),we get
# w = R*T*math.log((V2-b)/(V1-b))-(a/(b*T**(1/2)))*(math.log(V2/V1)-math.log((V2+b)/(V1+b))
# Now we have calculate V1 and V2,
# 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
# Solving the cubic equation at state 1,
def f5(V):
return V**(3)-((R*T)/P1)*V**(2)-((b_1**(2))+((b_1*R*T)/P1)-(a_1/(T**(1./2)*P1)))*V-(a_1*b_1)/(T**(1./2)*P1)
V3_1=fsolve(f5,1)
V3_2=fsolve(f5,10)
V3_3=fsolve(f5,100)
# The above equation has 1 real and 2 imaginary roots. We consider only real root (V3_3).
# Similarly at state 2,for P = P2,
def f6(V):
return V**(3)-((R*T)/P2)*V**(2)-((b_1**(2))+((b_1*R*T)/P2)-(a_1/(T**(1./2)*P2)))*V-(a_1*b_1)/(T**(1./2)*P2)
V3_4=fsolve(f6,1)
V3_5=fsolve(f6,10)
V3_6=fsolve(f6,100)
# The above equation has 1 real and 2 imaginary roots. We consider only real root (V3_6).
# Finally work done is given by
w4 = R*T*math.log((V3_6-b_1)/(V3_3-b_1))-(a_1/(b_1*T**(1./2)))*(math.log(V3_6/V3_3)-math.log((V3_6+b_1)/(V3_3+b_1))); #[J/mol]
w4 = w4*10**(-3); #[kJ/mol]
print " 3).The work done using Redlich Kwong equation of state is %f kJ/mol"%(w4);
```

In [7]:

```
T = 20 + 273.15; #[K] - Temperature
P_1 = 140.; #[kPa] - Initial pressure
P_1 = P_1*10.**(3); #[Pa]
P_2 = 560.; #[kPa] - Final pressure
P_2 = P_2*10.**(3); #[Pa]
R = 1.987; #[cal/mol*K] - Universal gas constant
# Calculations
# Cp_0 = 1.648+4.124*10**(-2)*T - 1.53*10**(-5)*T**(2) + 1.74*10**(-9)*T**(3)
# Using adiabatic compression, P*V**(Y)=constant. For ideal gases
# P*(R*T/P)**(Y) = constant
# P**(1-Y)*T**(Y) = constant or,P1**(1-Y)*T1**(Y)=P2**(1-Y)*T2**(Y)
# Now,at state 1, i.e at T=20[C]
Cp_1 = 1.648+4.124*10**(-2)*T-1.53*10**(-5)*T**(2)+1.74*10**(-9)*T**(3); #[cal/mol*K] - Heat capacity at constant pressure
Cv_1 = Cp_1 - R; #[cal/mol*K] - Heat capacity at constant volume
Y1 = Cp_1/Cv_1; # Ratio of heat capacities
# Now calculatung the temperature at state 2 (T2)
# (T2/T1)=(P1/P2)**((1-Y1)/Y1)
T_1 = T;
T_2 = ((P_1/P_2)**((1-Y1)/Y1))*T_1; #[K]
# Now calculating the mean temperature
T_mean = (T_1 + T_2)/2; #[K]
# At mean temperature
Cp_2 = 1.648+4.124*10**(-2)*T_mean - 1.53*10**(-5)*T_mean**(2) + 1.74*10**(-9)*T_mean**(3); #[cal/mol*K] - Heat capacity at constant pressure
Cv_2 = Cp_2 - R; #[cal/mol*K] - Heat capacity at constant volume
Y2 = Cp_2/Cv_2;
# Calculating exit temperature
# Again using the realation,(T2/T1)=(P1/P2)**((1-Y1)/Y1)
T_exit = ((P_1/P_2)**((1-Y2)/Y2))*T_1; #[K]
# Since value of mean temperature has not changed much the molar heat capacity ratio can be assumed to be same.Therefore
# w = -delta(U)=Cv_0*(T2-T1)
w = Cv_2*(T_1 - T_exit); #[cal/mol]
w = w*4.184; #[J/mol]
# Results
print " The work done for adiabatic compression is %f J/mol"%(w);
```

In [8]:

```
m_ice = 1000; #[g] - Mass of ice
m_water = 1000; #[g] - Mass of water
T_ice = 273.15; #[K] - Temperature of ice
T_water = 373.15; #[K] - Temperature of water
L = 79.71; #[cal/g] - Latent heat of melting of ice.
# Calculations and Results
#(1)
Cp_1 = 1; #[cal/g-K] - Heat capacity at constant pressure
# Let the final temperature be T
# We assume that all of the ice melts.Energy taken up by ice is
# E1 = L*m_ice + m_ice*Cp_1*(T - T_ice)
# Energy given by hot water is,
# E2 = m_water*Cp_1*(T_water - T)
# No heat exchange with surrounding.Solving for T
T_1 = (m_ice*Cp_1*T_ice + m_water*Cp_1*T_water - L*m_ice)/(m_ice*Cp_1 + m_water*Cp_1); #[K]
T_1 = T_1 - 273.15; #[C]
print " 1).The final temperature taking Cp_water = 1 cal/g-K) is %f C"%(T_1);
#Since the final temperature is greater than 273.15 K,so our assumption that all of ice melts is correct
# (2)
# Cp_2 = 1.00874-0.7067*10**(-3)*T+15.93*10**(-6)*T**(2)-83.8*10**(-9)*T**(3);
def f15(T):
return Cp_2
# From energy balance,we get L*m_ice + m_ice* quad(f15,0,T) + m_water*integrate(Cp_2)[0]
# On putting the values and then simplifying we get
# 2.01748*T - 0.0007067*T**(2) + 1.062*10**(-5)*T**(3) - 4.19*10**(-8)*T**(4) - 20.8455 = 0
# Solving the above equation we get
def f1(T):
return 2.01748*T - 0.0007067*T**(2) + 1.062*10**(-5)*T**(3) - 4.19*10**(-8)*T**(4) - 20.8455
T_0 = fsolve(f1,1)
print " 2).The final temperature using specific heat capacity equation is %f C"%(T_0);
```

In [9]:

```
# Variables
n = 1.5; # - ratio of heat capacities
T_1 = 500.; #[K] - Initial temperature
T_2 = 1000.; #[K] - Final temperature
P_1 = 1.; #[bar] - Initial pressure
P_1 = P_1*10.**(5); #[Pa]
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations and Results
# The compression path is given by, P*V**(1.5) = constant
# P*(R*T/P)**(1.5) = constant
# P1**(-0.5)*T1**(1.5) = P2**(-0.5)*T2**(1.5)
P_2 = P_1*(T_1/T_2)**(-3); #[Pa]
P_2_final = P_2*10**(-5); #[bar] - Final pressure in bar
print " The final pressure is %f bar"%(P_2_final);
# From first law q - w = delta(U).
# First w and delta(U) are calculated and thereafter heat exchange is determined.
V_1 = R*T_1/P_1; #[m**(3)/mol] - Initial volume
V_2 = R*T_2/P_2; #[m**(3)/mol] - Final volume
w = ((P_1*V_1)/(n - 1))*(1 - (P_2/P_1)**(1 - 1/n)); #[J/mol] - work done
# Mean temperature is given by,
T_mean = (T_1 + T_2)/2; #[K]
#Now, heat capacity at T_mean is given by,
Cp_0 = R*(3.3 + 0.63*10**(-3)*T_mean); #[J/mol*K]
Cv_0 = Cp_0 - R; #[J/mol*K]
#Therefore delta(U) is given by
del_U = Cv_0*(T_2 - T_1); #[J/mol] - Change in internal energy
q = w + del_U; #[J/mol] - heat change
print " The amount of heat supplied to the system is %f J/mol"%(q);
```

In [10]:

```
P_1 = 150*10**(3); #[Pa] - Initial pressure
V_1 = 0.001; #[m**(3)] - Initial volume
P_2 = 1000*10**(3); #[Pa] - Final pressure
V_2 = 0.003; #[m**(3)] - Final volume
# Calculations and Results
# At x = 0, Vt(total volume) = 0.001 m**(3), therefore x = (V_t - V_1)/A; where A is area of cross section and x is length
# Force exerted b sprig is given by, F = Ps*A = k*x = k*(V_t - V_1)/A
# Ps = (k/A**(2))*(V_t - V_1)
# Total pressure = Initial pressre + Pressre due to spring
# P = P_1 + (k/A**(2))*(V_t - V_1)
# Let (k/A**(2)) = t (say)
# At state 2, i.e at P2 and V_t = V_2.
def f(t):
return P_2-P_1 - t*(V_2-V_1)
t = fsolve(f,1000)
# Therefore,pressure is related to total volume as P = P_1-t*(V_t - V_1)
# (a)
#slope = (k/A**(2))
print " a).The slope of the line on P-Vt diagram is %e N/m**5)"%(t);
# (b)
# Work done by the gas is given by w=integral(PdVt)
def f30(V_t):
return P_1+t*(V_t-V_1)
w = quad(f30,V_1,V_2)[0]
w = w*10**(-3); #[kJ]
print " b).The work done by gas is %f kJ"%(w);
```

In [11]:

```
V = 36; #[L] - Vol of gas on each side
P_1 = 1; #[atm] - pressure on left side of the piston
P_1 = P_1*101325; #[Pa]
T = 273.15; #[K]
P_2 = 3.375; #[atm] - Pressure on right side of the piston
P_2 = P_2*101325; #[Pa]
Y = 1.44; # Ratio of heat capacities
R = 8.314; #[J/mol*K] - Universal gas constnt
# Calculations and Results
# (a)
# For total system, del(U_total) = Q.
# T_2/T_1 = (P_2/P_1)**((Y - 1)/Y)
T_right = T*(P_2/P_1)**((Y - 1)/Y); #[K]
Cv_0 = R/(Y-1); #[J/mol*K] - Heat capacity at constant volume.
# Now work done on the gas on right hand side is given by
# W = (P_1*V_1 - P_2*V_2)/(Y - 1) = R*(T_2 - T_1)/(Y - 1) = Cv_0*(T_1 - T_2)
W_left = Cv_0*(T - T_right); #[J/mol]
# Negative sign for the work done on LHS gas implies work is done on the gas
# For right hand side of the gas
# P*Vt = n*R*T
n = P_1*(V*10**(-3))/(R*T); # number of moles
W_right = (-W_left)*n; #[J] - We used negative sign for 'W_left' because it is negative in magnitude.
W_right = W_right/1000; #[kJ]
print " a).Total work done on gas on the right hand side is %f kJ"%(W_right);
#(b)
print " b).The final temperature of the gas on right side is %f K"%(T_right);
# (c)
# We have (P_left*V_left)/T_left = (P_right*V_right)/T_right.
# Since P_left = P_right, (V_left/T_left) = (V_right/T-right) and also P*V**(Y) = constant.
V_right = V*(P_1/P_2)**(1/Y); #[L] - The total volume on right side
# The total volume on right side can also be calculated using P2*V2 = n*R*T2.
# Since total volume = 72 [L], therefore volume of left side is
V_left = 2*V - V_right; #[L]
T_left = T_right*(V_left/V_right);
print " c).Final temperature of the gas on the left side is %f K"%(T_left);
#(d)
#The first law applied to the total system (left side and right side) gives.
#Q - W = del(U_left) + del(U_right)
#There is no net work done by the total system as the cylinder is closed at both ends.
Q = n*Cv_0*(T_left-T) + n*Cv_0*(T_right-T); #[J]
Q = Q/1000; #[kJ]
print " d).Amount of heat added to the gas on the left side is %f kJ"%(Q);
```

In [12]:

```
P_2 = 0.2; #[bar]
P_2 = P_2*10**(5); #[Pa]
int_dia_2 = 2.4*10**(-2); #[m] - internal diameter at state 2.
Q = 5*10**(-3); #[cubic metre/s] - Flow rate at point 2.
den = 1000; #[kg/cubic metre] - density
delta_z = 1; #[m] - Difference in height
g = 9.81; #[m/s**(2)] - Acceleration due to gravity
# Calculations and Results
# (1)
# (delta(P)/den) = (P2-P1)/den = P2/den
Vel_2 = Q/(3.14*(int_dia_2/2)**(2)); #[m/s] - Velocity of water at state 2.
# Velocity at state 1 i negligible as compared to velocity at state 2,because the diameter of reservoir is very large as compared to diameter of pipe at state 2
# From bernaulli equation we get,
# -w = (delta(P)/den) + delta(v**(2))/2 + g*delta_z
w = -((P_2/den )+ (Vel_2**(2)/2) + (g*delta_z)); #[J/kg]
# w multiplied by m = (den*Q), will give the fluid power.
m = den*Q; #[kg/s]
W_net = m*w; #[Watt]
print " 1).The fluid power is %f Watt"%(W_net);
#(2)
# Total discharge head developed by the pump is given by
# h = (delta(P)/den*g) + (Vel_2**(2)/2*g) + delta_z
h = (P_2/(den*g)) + (Vel_2**(2)/(2*g)) + delta_z; #[m]
print " 2).Total discharge head developed by the pump is given by h = %f m"%(h);
```

In [1]:

```
T_1 = 1000.; #[K] - Temperature at entry
P_1 = 0.6; #[MPa] - Pressure at entry
P_2 = 0.2; #[MPa] - Exit pressure
Vel_1 = 50.; #[m/s] - Entry velocity
Y = 1.4; # Ratio of heat capacities
Mol_wt = 28.; #[g/mol] - Molecular weight of air
Cp = 29.099; #[J/mol-K] - Specific heat capacity at constant pressure
Cp = (Cp/Mol_wt)*1000; #[J/kg-K]
# Calculations
# We know that for a flow process
# delta_H + delta_V**(2)/2 + delta_(g*z) = q - w
# delta_H + delta_V**(2)/2 = 0
# For a reversible process P*V**(Y) = constant and thus (T_2/T_1) = (P_2/P_1)**((Y-1)/Y)
T_2 = T_1*(P_2/P_1)**((Y-1)/Y); #[K] - Exit temperature
# delta_H + delta_V**(2)/2 = 0
# Vel_2**(2)/2 - Vel_1**(2)/2 - (H_1 - H_2)= 0
# Vel_2**(2)/2 - Vel_1**(2)/2 - Cp*(T_1 - T_2) = 0
Vel_2_square = 2*(Vel_1**(2.)/2 - Cp*(T_2 - T_1)); #[m**(2)/s**(2)]
Vel_2 = (Vel_2_square)**(1./2); #[m/s]
# Results
print " The discharge velocity is %f m/s"%(Vel_2);
```

In [14]:

```
P_entry = 10; #[bar] - Pressure at entry
V_entry = 200; #[m/s] - Velocity at entry
P_exit = 1; #[bar] - Vressure at exit
V_exit = 800; #[m/s] - Velocity at exit
g = 9.81; #[m/s**(2)] - Acceleration due to gravity
# Calculations
#Heat balance gives
# delta_H + (delta_V**(2))/2 + g*delta_z = q - w
#delta_H = q - w - (delta_V**(2))/2
#From nozzle no work is extracted,therefore
delta_H = -(V_exit**(2)- V_entry**(2))/2; #[J/kg]
delta_H = delta_H*10**(-3); #[kJ/kg]
# Results
print " The change in enthalpy per kg of steam is %f kJ/kg"%(delta_H);
```

In [15]:

```
T_1 = 280; #[K] - Temperature at entry
P_1 = 100; #[kPa] - Pressure at entry
T_2 = 400; #[K] - Temperature at exit
P_2 = 600; #[kPa] - Pressure at exit
m = 0.02; #[kg/s] - Mass flow rate
m = m*10**(3); #[g/s]
heat_loss = 16; #[kJ/kg]
# Calculations and Results
#Cp_0 = 28.11 + 0.1967*10**(-2)*T + 0.4802*10**(-5)*T**(2) - 1.966*10**(-9)*T**(3)
#delta_H = q - w (neglecting kinetic and potential changes)
#delta_H = integral(Cp_0*dT)
def f22(T):
return 28.11 + 0.1967*10**(-2)*T + 0.4802*10**(-5)*T**(2) - 1.966*10**(-9)*T**(3)
delta_H = quad(f22,T_1,T_2)[0]
print " Change in enthalpy is %f J/mol"%(delta_H);
#Molecular weight of air(21 vol% O2 and 79 vol% N2)=(0.21*32)+(0.79*28)= 28.84 g/mol
Mol_wt = 28.84; #[g/mol]
q = - (heat_loss*Mol_wt); #[J/mol]
w = q - delta_H; #[J/mol]
print " The work done per mole of air is %f J/mol"%(w);
#the negative sign implies that work is done on the compressor.
n = m/Mol_wt; #[mol/s] - Mole flow rate
W_net = delta_H*n; #[W]
W_net = -W_net*10**(-3); #[kW]
print " And the necessary power input to the compressor is %f kW"%(W_net);
```

In [16]:

```
T_1 = 300; #[K] - Temperature at entry
P_1 = 100; #[kPa] - Pressure at entry
P_2 = 900; #[kPa] - Pressure at exit
R = 8.314; #[J/mol*K] - Universal gas constant
# Calculations and Results
# (a)
# Reversible adiabatic compression
Y = 1.4; #Ratio of specific heat capacities
# For ideal gas, P*V**(Y)= constant and it takes the form of (T_2/T_1) = (P_2/P_1)**((Y-1)/Y)
T_2 = T_1*(P_2/P_1)**((Y - 1)/Y); #[K]
# The work exchange for adiabatic process is given by
# W_adia = -delta_H = -Cp*(T2-T1) = Cp*(T1-T2) = ((Y*R)/(Y-1))*(T1-T2)
W_adia = ((Y*R)/(Y - 1))*(T_1 - T_2); #[J/mol] -work done
# Molecular weight of air(21 vol% O2 and 79 vol% N2)=(0.21*32)+(0.79*28)= 28.84 g/mol
Mol_wt = 28.84; #[g/mol]
W_adia = W_adia/Mol_wt; #[J/g]
print " a).The compressor work done for reversible adiabatic compession is %f J/g"%(W_adia);
#(b)
#Isothermal compression
#W_iso = -integral(V*dP) = -integral((R*T/P)*dP) = R*T*ln(P_2/P_1)
W_iso = -R*T_1*math.log(P_2/P_1); #[J/mol]
W_iso = W_iso/Mol_wt; #[J/g]
print " b).The compressor work done for isothermal compession is %f J/g"%(W_iso);
#Note that in isothermal compression between the same states work done is less as compared to reversible adiabatic compression.
#(c)
#Ideal two-stage compression
n = 1.3; #Polytropic exponent.
#Minimum work done in two stage compression is given by
#W_comp = ((2*n*R*T_1)/(n-1))*[1-(P_x/P_1)**(n-1)/n]
#where for minimum work, (P_x/P_1) = (P_x/P_2), and thus
P_x = (P_1*P_2)**(1./2); #[kPa]
#therefore, work done is given by,
W_comp = ((2*n*R*T_1)/(n-1))*(1-(P_x/P_1)**((n-1)/n)); #[J/mol]
W_comp = W_comp/Mol_wt; #[J/g]
print " c).The compressor work done for ideal two-stage compession is %f J/g"%(W_comp);
```

In [17]:

```
T_1 = 600; #[C] - Temperature at entry
P_1 = 15; #[MPa] - Pressure at entry
T_2 = 400; #[K] - Temperature at exit
P_2 = 100; #[kPa] - Pressure at exit
A_in = 0.045; #[metre square] - flow in area
A_out = 0.31; #[metre square] - flow out area
m = 30; #[kg/s] - mass flow rate.
# Calculations and Results
#At 15 MPa and 600 C,it has been reported in the book that the properties of steam are,
Vol_1 = 0.02491; #[m**(3)/kg] - Specific volume
H_1 = 3582.3; #[kJ/kg] - Enthalpy
# m = den*vel*A = (Vel*A)/Vol, substituting the values
vel_1 = (m*Vol_1)/A_in; #[m/s] - Velocity at point 1.
print " The inlet velocity is %f m/s"%(vel_1);
#At 100 MPa (saturated vapour),it has been reported in the book that the properties of steam are, T_sat = 99.63 C, and
Vol_vap_2 = 1.6940; #[m**(3)/kg] - specific volume of saturated vapour.
H_vap_2 = 2675.5; #[kJ/kg] - Enthalpy os saturated vapour.
vel_2 = (m*Vol_vap_2)/A_out; #[m/s] - Velocity at point 2.
print " The exit velocity is %f m/s"%(vel_2);
#From first law we get, q - w =delta_H + delta_V**(2)/2
#q = 0, therefore, -w = delta_H + delta_V**(2)/2
delta_H = H_vap_2 - H_1; #[kJ/kg] - change in enthalpy.
delta_V_square = (vel_2**(2) - vel_1**(2))/2; #[J/kg]
delta_V_square = delta_V_square*10**(-3); #[kJ/kg]
w = -(delta_H + delta_V_square); #[J/kg]
W_net = w*m; #[kW]
W_net = W_net*10**(-3); #[MW] - power produced.
print " The power that can be produced by the turbine is %f MW"%(W_net);
```

In [18]:

```
R = 8.314; #[J/mol-K] - Universal gas constant
Cp_0 = 2.5*R; #[J/mol-K] - Specific heat capacity at constant pressure
Cv_0 = 1.5*R; #[J/mol-K] - Specific heat capacity at constant volume
T_L = 300; #[K] - Temperature at which port properties are constant.
# Calculations
Y = Cp_0/Cv_0; # Ratio of heat capacities.
#From part(1) we obtained the relation,
# T_2 = 1/(((P_2-P_1)/(Y*T_L*P_2))+(P_1/(P_2*T_1)))
# Not that when P_2 >> P_1 ,T_2 approaches Y*T_L and thus
T_2 = Y*T_L; #[K]
# Results
print " b).The final temperature is %f K"%(T_2);
```

In [19]:

```
T_1 = 40 + 273.15; #[K] - Initial temperature.
P_1 = 1; #[bar] - Initial pressure.
P_1 = P_1*10**(5); #[Pa]
Vol_1 = 0.01; #[cubic metre] - Initial volume of the cylinder.
T_2 = 100 + 273.15; #[K] - Final temperature.
P_2 = 100; #[kPa] - Final pressure.
P_2 = P_2*10**(5); #[Pa]
Vol_2 = 0.02; #[cubic metre] - Final volume of the cylinder.
Cp = 1.005; #[J/g-K] - Specific heat capacity at constant pressure.
Cv = 0.718; #[J/g-K] - Specific heat capacity at constant volume.
Mol_wt = 28.84; #[g/mol] - Molecular weight of air.
R = 8.314; #[J/mol-K] - universal gas constant
# Calculations
delta_Vol = Vol_2 - Vol_1; # [cubic metre] - Change in volume.
# Assuming ideal gas P*V = R*T
V_1 = (R*T_1)/P_1; # [m**(3)/mol] - Initial specific volume.
# Therefore,the total number of moles initially in the system is,
n_1 = (Vol_1/V_1); # [mol]
m_1 = n_1*Mol_wt; # [g] - Initial mass of the system.
Y = Cp/Cv; #Ratio of heat capacities
# The energy balance equation is given by
# -P*delta_Vol + H_liq*(m_2 - m_1) = m_2*Cv*(P*V2)/R - m_1*Cv*T_1
# m_2*Cv*(P*V2)/R = (Cv*P_1*Vol_2)/R
# Cv/R = 1/(Y-1)
# Since pressure of the gas in system is assumed constant,therefore it remains at 1 bar and thus P = P_1,
H_liq = Cp*T_2; # [J/g] - Enthalpy of liquid
m_2 = (P_1*delta_Vol + ((P_1*Vol_2)/(Y-1)) + H_liq*m_1 - m_1*Cv*T_1)/H_liq; #[g]
#The mass entering the assembly during the filling process is given by
m = m_2 - m_1; #[g]
n_2 = m_2/Mol_wt; #[mol] - Number of moles in the final state.
V_2 = Vol_2/n_2; #[m**(3)/mol] - Final specific volume.
# Therfore,final temperature is given by,
T_2 = (P_1*V_2)/R; #[K] - Final temperature.
# Results
print " The final equilibrium temperature is %f K"%(T_2);
print " The mass entering through the valve is %f g"%(m);
```

In [20]:

```
V_total = 5; #[L] - Volume of pressure cooker.
V_total = V_total*10**(-3); #m**(3)
P_gauze = 15; #[psi] - Operating pressure (gauze)of pressure cooker.
P_gauze = (P_gauze/14.5)*10**(5); #[N/m**(2)]
P_atm = 0.966*10**(5); #[N/m**(2)] - Atmospheric pressure.
m_1 = 1; #[kg] - Initial mass.
t = 30*60; #[s] - Total time.
J = 500; #[W] - Rate of heat supply
# Calculations
P_abs = P_gauze + P_atm; #[N/m**(2)] - Absolute pressure.
#The energy balance equqtion gives,
# Q = m_e*H_e +(m_2*U_2 - m_1*U_1), where 'm_e' is the mass exit from the system and 'H_e' is enthalpy at exit conditions.
#It has been reported in the book that from steam table at P_abs,
T_sat = 120.23; #[K]- Saturated temperature
V_liq = 0.001061; #[m**(3)/kg] - specific volume of liquid.
V_vap = 0.8857; #[m**(3)/kg] - specific volume of vapour.
U_liq = 504.49; #[kJ/kg] - specific internal energy of liquid.
U_vap = 2529.5; #[kJ/kg] - specific internal energy of vapour.
H_liq = 504.70; #[kJ/kg] - specific enthalpy of liquid.
H_vap = 2706.7; #[kJ/kg] - specific internal energy of vapour.
#We know that total volume occupied by 1 kg of fluid is
#V_total = (1-x)*V_liq + x*V_vap
x1 = (V_liq - V_total)/(V_liq - V_vap); #[g]
#Internal energy at this state is
U_1 = (1-x1)*U_liq + x1*U_vap; #[kJ/kg] - specific internal energy
U_1_net = m_1*U_1; #[kJ] - Internal energy
#The amount of heat suplied is given by,
J_net = J*t; #[J] - Net heat supplied.
J_net = J_net*10**(-3); #[kJ]
#Let the dryness factor at the end of the process be x
#Let the total mass of H2O (liquid + vapour) at the end of the process be 'm' kg.
# V_total/m = (1-x)*(V_liq) + x*V_vap ......equqtion(1)
#The energy balance equqtion gives,
# Q = m_e*H_e +(m_2*U_2 - m_1*U_1), where 'm_e' is the mass exit from the system and 'H_e' is enthalpy at exit conditions.
# The second equation on simplification becomes
# x = ((0.005/m) - 0.001061)/0.884639
# Putting the expression of x in first equation and then simplifying, we get
# - 1293.2 = -2202.21*m + 11.445 - 2.429*m
m = (11.445+1293.2)/(2202.21+2.429); #[kg]
# Therefore x can be calculated as
x = ((0.005/m) - 0.001061)/0.884639;
# Therfore total water (liquid + vapour) present in the pressure cooker at the end of the process is m kg.
m_vapour = x*m; #[kg] - Mass of vapour
m_liquid = (1-x)*m; #[kg] - Mass of vapour
# Results
print " Total water liquid + vapour) present in the pressure cooker at the end of the process is %f kg"%(m);
print " The mass of vapour is %f kg"%(m_vapour);
print " The mass of liquid is %f kg"%(m_liquid);
```