# Chapter 3 : The First Law and Its Applications¶

### Example - 3.1 Page number - 80¶

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);

 1.From steam table the final temperature lies between 700 and 800 C
Alternate method
By linear interpolation we get the temperature at which enthlpy is 4121.8 kJ/kg to be 786.74 C
2).The temperature assuming ideal behaviour is 1141.732355 K
Alternate method
The temperature assuming ideal behaviour alternate method) is 1141.738180 K

/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:227: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last ten iterations.
warnings.warn(msg, RuntimeWarning)


### Example - 3.2 Page number - 83¶

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);

 The amount of heat to be added is 102.663530 MJ


### Example 3.3 Page number - 83¶

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);

 a).The final temperature is 362.072311 c
b).The internal energy is 2775.731907 kJ/kg
b).The enthalpy is 3023.758380 kJ/kg


### Example - 3.4 Page number - 85¶

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";

 The work done is -5.186547 kJ/mol
Negative sign indicates that work is done on the gas


### Example - 3.5 Page number - 86¶

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);

 The work done is -5.186545 kJ/mol


### Example - 3.6 Page number - 88¶

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);

 1).The work done using given virial equation of state is -13.848149 kJ/mol
2).The work done using given virial equation of state is -13.688301 kJ/mol
3).The work done using  van der Walls equation of state is -14.802420 kJ/mol
3).The work done using Redlich Kwong equation of state is -14.704965 kJ/mol

/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:227: RuntimeWarning: The iteration is not making good progress, as measured by the
improvement from the last five Jacobian evaluations.
warnings.warn(msg, RuntimeWarning)


### Example - 3.8 Page number - 93¶

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);

 The work done for adiabatic compression is -3198.427072 J/mol


### Example - 3.9 Page number - 93¶

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);

 1).The final temperature taking Cp_water = 1 cal/g-K) is 10.145000 C
2).The final temperature using specific heat capacity equation is 10.364452 C


### Example - 3.11 Page number - 97¶

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);

 The final pressure is 8.000000 bar
The amount of heat supplied to the system is 3211.282500 J/mol


### Example - 3.12 Page number - 99¶

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 = w*10**(-3);			#[kJ]
print " b).The work done by gas is %f kJ"%(w);

 a).The slope of the line on P-Vt diagram is 4.250000e+08 N/m**5)
b).The work done by gas is 1.150000 kJ


### Example - 3.13 Page number - 99¶

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);

 a).Total work done on gas on the right hand side is 3.731958 kJ
b).The final temperature of the gas on right side is 396.112176 K
c).Final temperature of the gas on the left side is 1447.650324 K
d).Amount of heat added to the gas on the left side is 39.378580 kJ


### Example - 3.14 Page number - 105¶

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);

 1).The fluid power is -454.750210 Watt
2).Total discharge head developed by the pump is given by h = 9.271156 m


### Example - 3.15 Page number - 106¶

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);

 The discharge velocity is 749.965327 m/s


### Example - 3.16 Page number - 107¶

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);

 The change in enthalpy per kg of steam is -300.000000 kJ/kg


### Example - 3.17 Page number - 111¶

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)

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);

 Change in enthalpy is 3511.197066 J/mol
The work done per mole of air is -3972.637066 J/mol
And the necessary power input to the compressor is -2.434949 kW


### Example - 3.18 Page number - 112¶

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]
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);

 a).The compressor work done for reversible adiabatic compession is -264.386412 J/g
b).The compressor work done for isothermal compession is -190.024880 J/g
c).The compressor work done for ideal two-stage compession is -216.284501 J/g


### Example - 3.19 Page number - 113¶

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);

 The inlet velocity is 16.606667 m/s
The exit velocity is 163.935484 m/s
The power that can be produced by the turbine is 26.805014 MW


### Example - 3.20 Page number - 117¶

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);

 b).The final temperature is 500.000000 K


### Example - 3.21 Page number - 119¶

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);

 The final equilibrium temperature is 339.343160 K
The mass entering through the valve is 9.367211 g


### Example - 3.22 Page number - 122¶

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);

 Total water liquid + vapour) present in the pressure cooker at the end of the process is 0.591773 kg
The mass of vapour is 0.004942 kg
The mass of liquid is 0.586830 kg