In [1]:

```
import math
from scipy.misc import derivative
#Given:
r = 7. #Compression ratio
g = 1.4 #Specific heat ratio(gamma)
cv = 0.718 #(Assume)Specific heat at consmath.tant volume in kJ/kgK
dcv = 1*cv/100 #Change in specific heat in kJ/kgK
#Solution:
R = cv*(g-1) #Specific gas consmath.tant in kJ/kgK
eta = round(100*(1-1/r**(g-1)))/100 #Efficiency when there is no change in specific heat
def Otto(cv): #Defining efficiency as a function of specific heat
return 1-1/r**(R/cv)
detaBydcv = derivative(Otto,cv) #Derivative of efficiency wrt to specific heat at initial value of specific heat
detaByeta = detaBydcv*dcv/eta #Change in efficiency wrt to initial value of efficiency
#Results:
print " The percentage change in the efficiency of Otto cycle = %.3f percent"%(detaByeta*100)
if (detaByeta < 0):
print ("decrease")
```

In [12]:

```
from sympy import Symbol,solve
from numpy import roots
from math import floor
def horner(coeffs, x):
acc = 0
for c in reversed(coeffs):
acc = acc * x + c
return acc
#Given:
r = 6. #Compression ratio
CV = 44000. #Calorific value in kJ/kg of fuel
A_F = 15./1 #Air-fuel ratio
P1 = 1. #Pressure at 1 in bar
T1 = 60.+273 #Temperature at 1 in K
n = 1.32 #Index of compression
T = Symbol("T") #Defining temperature(T) as unknown in K
cv = 0.71+20e-5*T #Specific heat at consmath.tant volume as a function of temperature(T) in kJ/kgK
cv_c = 0.71 #Constant specific heat in kJ/kgK
cv1 = [0.71,20e-5]
#Solution:
#Refer fig 3.19
P2 = P1*r**n #Pressure at 2 in bar
T2 = floor(T1*r**(n-1)) #Temperature at 2 in K
T3 = Symbol("T3") #Defining temperature(T3) as unknown in K
T_av = (T3+T2)/2 #Average temperature during combustion of charge in K
cv_mean = horner(cv1,T_av) #Mean specific heat in kJ/kgK
#Assume cycle consumes 1 kg of air
m_a = 1.
m_f = m_a/A_F
m_c = m_f+m_a #Mass of air, fuel, and charge in kg
Q1 = CV*m_f #Heat added per kg of air in kJ/kg
T3_v = solve(Q1-cv_mean*m_c*(T3-T2))
T3_v = T3_v[1] #Temperature at 3 in K
P3_v = P2*T3_v/T2 #Pressure at 3 in bar
#For consmath.tant specific heat
T3_c = solve(Q1-cv_c*m_c*(T3-T2))[0] #Temperature at 3 for consmath.tant specific heat in K
P3_c = P2*T3_c/T2 #Pressure at 3 for consmath.tant specific heat in bar
#Results:
print " The maximum pressure in the cycle for variable specific heat, P3 = %.1f bar"%(P3_v)
print " The maximum pressure in the cycle for consmath.tant specific heat, P3 = %.1f bar"%(P3_c)
```

In [5]:

```
from scipy.optimize import fsolve
import math
from scipy.integrate import quad
#Given:
A_F = 28./1 #Air-fuel ratio
CV = 42000. #Calorific value in kJ/kg
R = 0.287 #Specific gas consmath.tant in kJ/kgK
r = 14./1 #Compression ratio
cv = 0.71+20-5*r #Specific heat at consmath.tant volume as a function of temperature(T) in kJ/kgK
T2 = 800. #Temperature at the end of the compression process (2) in K
#Solution:
#Refer fig 3.20
#Assume cycle consumes 1 kg of fuel
m_c = A_F*1+1 #Mass of charge in kg
cp = cv #addf(cv,R) #Specific heat at consmath.tant pressure as a function of temperature(T) in kJ/kgK
#Since, heat transfer at consmath.tant pressure, Q1 = integration(cp*dt) from T2 to T3
#Thus, Q1 is the function of T3. Defining the function Q1 of T3
def difference(T3):
def f0(T):
return cp
Q1 = quad(f0,T2,T3)[0]
return Q1-CV/m_c
#Since, heat transfer at consmath.tant pressure must be equal to calorific value per kg of charge
#Thus, their difference must be zero, function Q1toCV is solve for zero
V2 = 1.
T3 = 1939.5205
T3 = round(T3) #Temperature at the end of consmath.tant pressure proces (3) in K
rho = T3/T2 #Cut off ratio
V3 = rho #Volume at 3 in units
p = (V3-V2)*100/(r-V2) #percentage of stroke at which consmath.tant pressure process ends
#Results:
print " At %.2f percentage of stroke combustion is completed."%(p)
```

In [6]:

```
from scipy.optimize import fsolve
import math
from scipy.integrate import quad
from sympy import Symbol
#Given:
P1 = 1. #Pressure at 1 in bar
T1 = 90.+273 #Temperature at 1 in K
r = 13. #Compression ratio
Q1 = 1675. #Heat supplied per kg of air in kJ/kg
Q1_v = Q1/2
Q1_p = Q1/2 #Heat supplied at consmath.tant volume and pressure per kg of air in kJ/kg
g = 1.4 #Specific heat ratio(gamma)
R = 0.287 #Specific gas consmath.tant in kJ/kgK
T = Symbol("T")
cv = 0.71+20-5*T #Specific heat at consmath.tant volume as a function of temperature(T) in kJ/kgK
cp = 0.71+20-5*T+0.287
#Solution:
#Refer fig 3.21
P2 = P1*r**g #Pressure at 2 in bar
T2 = T1*r**(g-1) #Temperature at 2 in K
#Since, heat transfer at consmath.tant volume, Q1_v = integration(cv*dt) from T2 to T3
#Thus, Q1_v is the function of T3. Defining the function Q1_v of T3
def Volume(T3):
def f1(T):
return cv
Q1_v = quad(f1,T2,T3)[0]
return Q1_v-Q1/2
#Since, heat transfer at consmath.tant volume must be equal to half of total heat added
#Thus, their difference must be zero, function Q1_vtoQ1 is solve for zero
T3 = 1853.0823
P3 = P2*T3/T2 #Pressure at 3 in bar
#cp = addf(cv,R) #Specific heat at consmath.tant pressure as a function of temperature(T) in kJ/kgK
#Since, heat transfer at consmath.tant pressure, Q1_p = integration(cp*dt) from T3 to T4
#Thus, Q1_p is the function of T4. Defining the function Q1_p of T4
def Pressure(T4):
def f2(T):
return cp
Q1_p = quad(f2,T3,T4)[0]
return Q1_p-Q1/2
#Since, heat transfer at consmath.tant pressure must be equal to half of total heat added
#Thus, their difference must be zero, function Q1_ptoQ1 is solve for zero
T4 = 2440.2522
rho = T4/T3 #Cut off ratio
p = (rho-1)*100/(r-1) #Percentage of stroke at which cut off occurs
#Results:
print " The maximum pressure in the cycle, P3 = %.1f bar"%(P3)
print " The percentage of stroke at which cut off occurs = %.2f percent"%(p)
```

In [8]:

```
import math
#Given:
r = 7. #Compression ratio
CV = 44000. #Calorific value of the fuel in kJ/kg
A_F = 13.67 #Air fuel ratio of the mixture
cv = 0.718 #Specific heat at consmath.tant volume in kJ/kgK
n = 1.3 #Polytropic index
P1 = 1.
T1 = 67.+273 #Pressure and temperature at the beginning in bar and K
#Solution:
#Refer fig 3.22
C = 12. #Atomic mass of Carbon(C)
H = 1. #Atomic mass of Hydrogen(H)
O = 16. #Atomic mass of Oxygen(O)
p = 23. #Percentage of oxygen in air by mass
#Stoichiometric equation of combustion of fuel (C6H14)
# [C6H14] + x[O2] = y[CO2] + z[H2O]
#Equating coefficients
x = 9.5
y = 6.
z = 7. #Coefficients of stoichiometric equation
A_F_g = x*2*O/(6*C+14*H)*100/p #Gravimetric air fuel ratio
MS = A_F_g/A_F*100 #Actual mixture strength in percent
#Combustion is incomplete
#Stoichiometric equation of incomplete combustion of fuel (C6H14)
# MS/100[C6H14] + x[O2] = a[CO2] + b[CO] + c[H2O]
#Equating coefficients
a = 4.39
b = 2.36
c = 7.87 #Coefficients of stoichiometric equation
#Stoichiometric equation of combustion of fuel (C6H14) by adding Nitrogen
# MS/100[C6H14] + x[O2] + x*79/21[N2] = a[CO2] + b[CO] + c[H2O] + x*79/21[N2]
m1 = MS/100+x+x*79/21 #Moles before combustion
m2 = a+b+c+x*79/21 #Moles after combustion
Me = (m2-m1)/m1*100 #Molecular expansion in percent
T2 = T1*r**(n-1) #Temperature at 2 in K
m_c = A_F+1 #Mass of charge in kg
T3 = CV/(m_c*cv)+T2 #Temperature at 3 in K
T3 = round(T3)
P3 = P1*r*(T3/T1) #Pressure at 3 in bar (printing error)
#Temperature and pressure considering molecular expansion
T31 = T3 #Temperature remains same at 3 in K
P31 = P3*m2/m1 #Pressure at 3 in bar
#Results:
print "\t The molecular expansion = %.2f percent"%(Me)
print " a)Without considering the molecular contraction\t The maximum pressure, P3 = %.2f bar\
\n The maximum temperature, T3 = %.0f K"%(P3,T3)
print " b)Considering the molecular contraction\t The maximum pressure, P3 = %.2f bar \
\nThe maximum temperature, T3 = %.0f K"%(P31,T31)
#Answer in the book is wrong
```

In [9]:

```
import math
from numpy import interp
#Given:
p = 15. #Clearance volume in percentage of print lacement volume
V_s = 2.8 #Swept volume in litres
N = 2500. #Engine speed in rpm
Q1 = 1400. #Heat added in kJ/kg
T1 = 27.+273 #Temperature at inlet in K
P1 = 100. #Pressure at inlet in kPa
R = 0.287 #Specific gas consmath.tant in kJ/kgK
#Solution:
#Refer fig 3.23
#By umath.sing gas tables
#Refer Ideal-gas properties of air
V2 = (p/100)*(V_s/1000) #Volume at 2 (Clearance volume) in m**3
V3 = V2 #Volume at 3 in m**3
V1 = V_s/1000+V2
V4 = V1 #Volume at 1, 4 in m**3
# Process 1-2
vr1 = 621.2
pr1 = 1.3860
u1 = 214.09
phi1 = 5.7016 #Relative specific volume, relative pressure, specific internal energy(kJ/kg), specific entropy(kJ/kgK) at 1 (from air tables)
vr2 = vr1*(V2/V1) #Relative specific volume at 2
vr = [81.89, 78.61]
T = [660. ,670.]
pr = [23.13, 24.46]
u = [481.01, 488.81] #Relative specific volume, temperature(K), relative pressure, specific internal energy(kJ/kg) (extracted from air tables)
#Finding the corresponding temperature at vr2 by interpolation
#T2 = interp1d([vr,T],vr2) #Temperature at 2 in K
T2 = interp(vr2,vr,T)
#Finding the corresponding relative pressure at T2 by interpolation
pr2 = interp(T2,pr,T) #Relative pressure at 2
#Finding the corresponding specific internal energy at T2 by interpolation
u2 = interp(T2,u,T) #specific internal energy at 2 in kJ/kg
P2 = P1*(pr2/pr1) #Pressure at 2 in kPa
# Process 2-3
u3 = Q1+u2 #Specific internal energy at 3 in kJ/kg
vr = [2.356 ,2.175, 2.012]
T = [2100, 2150, 2200]
pr = [2559, 2837, 3138]
u = [1775.3 ,1823.8, 1872.8] #Relative specific volume, temperature(K), relative pressure, specific internal energy(kJ/kg) (extracted from air tables)
#Finding the corresponding relative specific volume at u3 by interpolation
vr3 = interp(u3,vr,u) #Relative specific volume at 3
#Finding the corresponding relative pressure at u3 by interpolation
pr3 = interp(u3,pr,u) #Relative pressure at 3
#Finding the corresponding temperature at u3 by interpolation
T3 = interp(u3,u,T) #Temperature at 3(maximum) in K (Round off error)
P3 = P2*(T3/T2) #Pressure at 3(maximum) in kPa
# Process 3-4
vr4 = vr3*(V4/V3) #Relative specific volume at 4
vr = [15.241 ,14.470]
T = [1180, 1200]
pr = [222.2 ,238.0]
u = [915.57, 933.33]
phi = [7.1586, 7.1684] #Relative specific volume, temperature(K), relative pressure, specific internal energy(kJ/kg), specific entropy(kJ/kgK) (extracted from air tables)
#Finding the corresponding temperature at vr4 by interpolation
T4 = interp(vr4,T,vr) #Temperature at 4 in K
#Finding the corresponding specific internal energy at T4 by interpolation
u4 = interp(T4,u,T) #Specific internal energy at 4 in kJ/kg
#Finding the corresponding relative pressure at T4 by interpolation
pr4 = interp(T4,pr,T) #Relative pressure at 4
P4 = P3*(pr4/pr3) #Pressure at 4 in kPa
#Finding the corresponding specific entropy at T4 by interpolation
phi4 = interp(T4,phi,T) #Specific entropy at 4 in kJ/kgK
# Process 4-1
Q2 = u1-u4 #Heat rejected in kJ/kg
W = Q1+Q2 #Work done in kJ/kg
eta = W/Q1 #Efficiency
m = P1*V1/(R*T1) #Mass of air in cycle in kg
W = m*W*N/60 #Rate of work in kW
Delta_s = phi1-phi4-R*math.log(P1/P4) #Change in specific entropy between 1 and 4 in kJ/kgK
AE = Q2-T1*(Delta_s) #Available portion of energy of Q2 in kJ/kg (Round off error)
p_AE = AE/Q2 #Available energy in percentage of Q2
# Without umath.sing gas tables
g = 1.4 #Specific heat ratio(gamma)
cv = 0.718 #Specific heat at consmath.tant volume in kJ/kgK
r = V1/V2 #Compression ratio
eta1 = 1-1/r**(g-1) #Efficiency
# Process 1-2
T2 = T1*(r)**(g-1) #Temperature at 2 in K
P2 = P1/100*(r)**g #Pressure at 2 in bar
# Process 2-3
T31 = Q1/cv+T2 #Temperature at 3(maximum) in K
P31 = P2*T31/T2 #Pressure at 3(maximum) in bar
# Process 3-4
T4 = T31*(1/r)**(g-1) #Temperature at 4 in K
Q2 = cv*(T1-T4) #Heat rejected in kJ/kg
W1 = Q1+Q2 #Work done in kJ/kg
et1a1 = W1/Q1 #Efficiency
power = m*W1*N/60 #Power in kW
Delta_s = cv*math.log(T1/T4) #Change in specific entropy between 1 and 4 in kJ/kgK
AE1 = Q2-T1*Delta_s #Available portion of energy of Q2 in kJ/kg (Round off error)
p_AE1 = AE1/Q2 #Available energy in percentage of Q2 (Round off error)
#Results:
print " Constant specific heat:\t Maximum temperature, Tmax = %.1f K\t Maximum pressure, \
\nPmax = %.1f bar\t Thermal efficiency, eta = %.2f percent\t \
\nPower = %.1f kW\t Available portion of heat rejected = %.1f kJ/kg %.1f percent"%(T31,P31,eta1*100,power,abs(AE1),p_AE1*100)
print " Variable specific heat:\t Maximum temperature, Tmax = %.0f K \
\n Maximum pressure, Pmax = %.1f bar\t Thermal efficiency eta = %.1f percent\t \
\nPower = %.1f kW\t Available portion of heat rejected = %.1f kJ/kg %.1f percent)"%(2210,56.5,49.8,108.6,384.2,54.6)
#Round off error in 'T3', 'AE', 'AE1', 'p_AE1'
```