Chapter 11:Compressible flow¶
Example 11.1 page no 617.¶
In [16]:
D=4.0 #in
T1=540.0 #degree R
p1=100.0 #psia
T2=453.0 #degree R
p2=18.4 #psia
k=1.4
R=1716/32.174 #ft*lb/(lbm*(degree R))
cv=R/(k-1) #ft*lb/(lbm*(degree R))
udiff=cv*(T2-T1) #ft*lb/lbm change in internal energy
#Result
print "a)The change in internal energy between (1) and (2)=",round(udiff,2),"ft*lb/lbm"
cp=k*round(cv,0) #ft*lb/(lbm*(degree R))
hdiff=cp*(T2-T1) #ft*lb/lbm change in enthalpy
print "b)The change in enthalpy between (1) and (2)=",round(hdiff,0),"ft*lb/lbm"
ddiff=(1/R)*((p2*144/T2)-(p1*144/T1)) #lbm/(ft**3) change in density
#Result
print "The change in density betwenn (1) and (2)=",round(ddiff,3),"ft*lb/lbm"
Example 11.2 page no.619¶
In [15]:
D=4.0 #in
T1=540.0 #degree R
p1=100.0 #psia
T2=453.0 #degree R
p2=18.4 #psia
#Calculation
import math
dratio=(p1/T1)*(T2/p2)
sdif=(cv*(math.log(T2/T1)))+(R*(math.log(dratio)))#ft*lb/lbm*(degree R) change in entropy
#Result
print "The change in entropy between (1) and (2)=",round(sdif,1),"ft*lb/lbm*(degree R)"
Example 11.3 page no.623¶
In [16]:
T=0 #degree C
R=286.9 #J/(kg*K)
k=1.401
c=(R*(T+273.15)*k)**0.5 #m/s
#Result
print "The speed of sound for air at 0 degree C =",round(c,1),"m/s"
Example 11.4 page no.628¶
In [1]:
%matplotlib inline
In [4]:
z=1000 #m
Ma=1.5
T=20 #degree C
#alpha=atan(z/x), x=V*t,and Ma=(1/sin(alpha)) where alpha is the angle of the Mach cone
#V=Ma*c
#calculation
import math
c=343.3 #m/s found from the value of temperature
V=Ma*c #m/sec
t=z/(Ma*c*math.tan(math.asin(1/Ma))) #sec
print"The number of seconds to wait after the plane passes over-head before it is heard=",round(t,2),"s"
#Plot
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111)
Ma=[1,1.1,1.3,1.5,2,4]
t=[0,1.5,2,2.17,2.3,2.75]
xlabel("Ma")
ylabel("t (s)")
plt.xlim((0,4))
plt.ylim((0,3))
ax.plot([1.5], [2.17], 'o')
ax.annotate('(1.5,2.17s)', xy=(1.7,2.1))
a=plot(Ma,t)
Example 11.5 page no.635¶
In [20]:
A=1*(10**(-4)) #m**2
p1=80.0 #kPa(abs)
p2=40.0 #kPa(abs)
p0=101.0 #kPa(abs)
pcritical=0.528*p0 #kPa(abs)
k=1.4
#for (a) pth=p1>pcritical
Math1=((((p0/p1)**((k-1)/k))-1)/((k-1)/2))**(0.5) #Math=Mach number at throat
#dth/d0=p1/p0 dth=density at throat
dth1=(1.23)*(1/(1+(((k-1)/2)*(Math1**2))))**(1/(k-1)) #kg/(m**3) density at throat
Tth1=(288)*(1/(1+(((k-1)/2)*(Math1**2)))) #K temperature at throat
Vth1=Math1*(286.9*269*k)**(0.5) #m/sec
m1=dth1*A*Vth1 #kg/sec
#Result
print "a) The mass flowrate through the duct=",round(m1,4),"kg/s"
#for (b) pth=p2<pcritical, hence
Math2=1
dth2=1.23*(1/(1+(((k-1)/2)*(Math2**2))))**(1/(k-1)) #kg/(m**3) density at throat
Tth2=(288)*(1/(1+(((k-1)/2)*(Math2**2)))) #K temperature at throat
Vth2=Math2*(286.9*Tth2*k)**0.5 #m/sec
m2=dth2*A*Vth2 #kg/sec
print "b) The mass flowrate through the duct=",round(m2,3),"kg/s"
Example 11.6 page no.637¶
In [15]:
A=1*(10**(-4)) #m**2
p1=80 #kPa(abs)
p2=40 #kPa(abs)
p0=101 #kPa(abs)
k=1.4
#for (a)
pratio1=p1/p0
#for this value of p1/p0,
Math1=0.59
Tratio1=0.94 #=Tth/T0
dratio1=0.85 #=dth/d0
Tth1=Tratio1*(288) #K
dth1=dratio1*(1.23) #kg/(m**3)
Vth1=Math1*(286.9*Tth1*k)**0.5 #m/sec
m1=(dth1*A*Vth1) #kg/sec
print "a)The mass flowrate=",round(m1,3),"kg/s"
#for (b)
Math2=1
Tratio2=0.83 #=Tth/T0
dratio2=0.64 #=dth/d0
Tth2=Tratio2*(288) #K
dth2=dratio2*(1.23) #kg/(m**3)
Vth2=Math2*(286.9*Tth2*k)**0.5 #m/sec
m2=(dth2*A*Vth2) #kg/sec
#Result
print "b)The mass flowrate=",round(m2,3),"kg/s"
Example 11.7 page no.638¶
In [18]:
#Given
pratio=0.82 #ratio of static to stagnation pressure
T=68 #degree F
#for (a)
#for the value of pratio given Ma is calculated as
Ma1=0.54
k1=1.4
Tratio1=0.94 #T/T0
T1=Tratio1*(T+460) # degree R
V1=(Ma1*(53.3*T1*k1)**0.5)*(32.2**0.5) #ft/sec
#for (b)
k2=1.66
Ma2=((((1/pratio)**((k2-1)/k2))-1)/((k2-1)/2))**0.5
Tratio2=1/(1+(((k2-1)/2)*(Ma2**2))) #T/T0
T2=Tratio2*(T+460) #degree R
V2=(Ma2*(386*T2*k2)**0.5)*(32.2**0.5) #ft/sec
#Result
print "The flow velocity if fluid is air=",round(V1,0),"ft/s"
print "The flow velocity if fluid is helium=",round(V2,0),"ft/s"
Example 11.11 page no.648¶
In [14]:
#Given
k=1.4
T0=518.67 #degree R
T1=514.55 #degree R
p1=14.3 #psia
R=53.3 #(ft*lb)/(lbm* degree R)
#calculation
from scipy.optimize import fsolve
import math
cp=R*k/(k-1) #(ft*lb/(lbm* degree R))
Tratio=T1/T0
Ma=(((1/Tratio)-1)/((k-1)/2))**0.5
x=(R*T1*k*32.2)**0.5 #ft/sec where x=(R*T1*k)**0.5
y=p1*144/(R*T1)*(Ma*x) # lbm/((ft**2)*sec) where y=d*V
#for p=7 psia
p=7 #psia
def f(T):
return(6.5*10**-5*T**2+T-518.67)
T=fsolve(f,500)
sdif=(cp*math.log(T/T1))-(R*math.log(p/p1)) #(ft*lb)/(lbm* degree R)
#Result
print "The corrosponding value of temperature for Fanno for downstream pressure of 7psia=",round(T,1),"K"
print "The corrosponding value of entropy change for Fanno for downstream pressure of 7psia=",round(sdif,1),"(ft*lb)/(lbm* degree R)"
#Plot
s=[33.6,39.8,46.3,52.6,57.3,57.9,55.5,53.0,47.0,44.2]
T=[502.3,496.8,488.3,474.0,447.7,432.1,394.7,378.1,347.6,335.6]
a=plot(s,T)
xlabel("s-s1 (ft lb/lbm R)")
ylabel("T (R)")
plt.xlim((30, 65))
plt.ylim((300,550))
show(a)
Example 11.12 page no.655¶
In [1]:
#Given
T0=288 #K
p0=101 #kPa(abs)
l=2 #m
D=0.1 #m
f=0.02
k=1.4
x=f*l/D
Tratio=2/(k+1) #where Tratio is Tcritical/T0
Tcritical=Tratio*T0 #K = T2
Vcritical=(286.9*Tcritical*k)**0.5 #m/sec = V2
#from value of x, the following are found
Ma=0.63
Trat=1.1 #where Trat=T1/Tcritical
Vrat=0.66 #where Vrat=V1/Vcritical
prat=1.7 #where prat=p1/pcritical
pratio=1.16 #where pratio=p0,1/p0critical
#from value of Ma, the following are found
Tfraction=0.93 #whereTfraction=T1/T0
pfraction=0.76 #where pfraction=p1/p0,1
dfraction=0.83 #where dfraction=d1/d0,1
#hence,
#calculation
import math
V1=Vrat*Vcritical #m/sec
d1=dfraction*(1.23) #kg/(m**3)
m=d1*math.pi*(D**2)*V1/4 #kg/sec
T1=Tfraction*T0 #K
p1=pfraction*p0 #kPa(abs)
T01=T0 #K and T01=T02
p01=p0 #kPa(abs)
p2=(1/prat)*(pfraction)*p01 #kpa(abs)
p02=(1/pratio)*p01 #kPa(abs)
#Result
print "Critical temperature=",round(Tcritical,2),"K"
print "Critical velocity=",round(Vcritical,0),"m/s"
print "Velocity at inlet=",round(V1,0),"m/s"
print "Maximum mass flowrate=",round(m,2),"Kg/s"
print "Temperature at inlet=",round(T1,2),"K"
print "Pressure at inlet=",round(p1,0),"kPa(abs)"
print "stagnation temperature at inlet and exit=",round(T01,2),"K"
print "The stagnation pressure at inlet=",round(p01,2),"kPa(abs)"
print "Pressure at exit=",round(p2,2),"kPa(abs)"
print "The stagnation pressure at exit=",round(p02,2),"kPa(abs)"
#Plot
#import numpy as np
#import matplotlib.pyplot as plt
#fig = plt.figure()
#ax = fig.add_subplot(111)
s=[-5,50]
T=[288,288]
s1=[0,0]
T1=[288,268]
s2=[40,40]
T2=[288,240]
s3=[0,10,20,30,40]
T3=[268,265,261,252,240]
s4=[30,50]
T4=[240,240]
s5=[-3,50]
T5=[288,288]
s6=[-2,8]
T6=[284,294]
s7=[35,45]
T7=[281,295]
s8=[-2,8]
T8=[265,275]
s9=[35,45]
T9=[235,245]
a=plot(s,T)
a1=plot(s1,T1)
a2=plot(s2,T2)
a3=plot(s3,T3,linestyle='--')
a4=plot(s4,T4)
a5=plot(s5,T5)
a6=plot(s6,T6)
a7=plot(s7,T7)
a8=plot(s8,T8)
a9=plot(s9,T9)
xlabel("s-s1 (J/kg K)")
ylabel("T (K)")
plt.xlim((-5, 70))
plt.ylim((230,300))
plt.text(50,240,'T2=240K')
plt.text(50,288,'T0=288K')
plt.text(15,255,'fanno line')
plt.text(10,280,'p1=70kpa(abs)')
plt.text(15,268,'T1=268K')
plt.text(10,295,'p01=101kpa(abs)')
plt.text(50,295,'p02=84kpa(abs)')
plt.text(45,245,'p2=45kpa(abs)')
show(a)
show(a1)
show(a2)
show(a3)
show(a4)
show(a5)
show(a6)
show(a7)
show(a8)
show(a9)
Example 11.13 page no.657¶
In [10]:
#Given
T0=288 #K Temprature
p0=101 #kPa(abs) pressure
l=2 #m length of duct
D=0.1 #m diameter
f=0.02
pd=45 #kPa(abs)
f=0.02
k=1.4
lnew=(50/100)*l #m new length of duct
x=lnew*f/D
#Calculation
#from this value of x, following are found
Ma=0.7
prat=1.5 #where prat=p1/pcritical
#from this value of Ma, following are found
pratio=0.72 #where pratio=p1/p0
dratio=0.79 #where dratio=d1/d0,1
Vratio=0.73 #where Vratio=V1/Vcritical
#hence,
p2=(1/prat)*pratio*p0 #kPa(abs)
pcritical=p2
#we find that pd<pcritical
d1=dratio*(1.23) #kg/(m**3)
Vcritical=310 #m/sec = V2, from example 11.12
V1=Vratio*Vcritical #m/sec
m=d1*math.pi*(D**2)*V1/4 #kg/sec
#Result
print "The flowrate for the smaller tube is =",round(m,2),"kg/s"
Example 11.14 page no.658¶
In [14]:
#given
T0=288 #K
p0=101 #kPa(abs)
l=2 #m
D=0.1 #m
f=0.02
pd=45 #kPa(abs)
f=0.02
m=1.65 #kg/sec
lnew=l/2 #m
x=f*l/D
#from this value of x, Ma at exit is found as
Ma=0.7
#and p2/pcritical is found as
pratio=1.5
#and, from example 11.12,
prat=1.7 #where prat=p1/pcritical
pfraction=0.76 #where pfraction=p1/p0,1
#Hence,
p2=pratio*(1/prat)*pfraction*p0 #kPa(abs)
#Result
print "The Mach number at the exit=",round(Ma,3)
print "The back pressure required=",round(p2,0),"kPa(abs)"
Example 11.15 page no.659¶
In [9]:
k=1.4
T0=518.67 #degree R
T1=514.55 #degree R
p1=14.3 #psia
R=53.3 #(ft*lb)/(lbm*degree R)
cp=R*k/(k-1) #(ft*lb/(lbm* degree R))
Tratio=T1/T0
Ma=(((1/Tratio)-1)/((k-1)/2))**0.5
x=(R*T1*k*32.2)**0.5 #ft/sec where x=(R*T1*k)**0.5
y=p1*144/(R*T1)*(Ma*x) #lbm/((ft**2)*sec) where y=d*V
z=R*T1/(p1*144) #(ft**3)/lbm
c=(p1)+(y*y*z/(32.2*144)) #psia =constant
#when downstream pressure p=13.5 psia
p=13.5 #psia
a=(y**2)*R/(p*144*32.2*144) #(lb/(in**2))/degree R
T=(c-p)/a
sdif=(cp*math.log(T/T1))-(R*math.log(p/p1)) #ft*lb/(lbm*degree R)
#Result
print "The corrosponding value of temperature for the downstream pressure of 13.5 psia=",round(T,0),"degree R"
print "The corrosponding value of change in entropy for the downstream pressure of 13.5 psia=",round(sdif,0), "ft*lb/(lbm*degree R)"
s=[93.2,202,251,285,317,330,333,334,336,338,338,336,333,328,321,259,181]
T=[969,1459,1859,2168,2464,2549,2558,2558,2544,2488,2450,2369,2260,2140,1992,1175,633]
a=plot(s,T)
xlabel("s-s1 (ft lb/lbm R)")
ylabel("T (R)")
plt.xlim((90,350))
plt.ylim((500,3000))
show(a)
Example 11.18 page no.670¶
In [24]:
#Given
p=60 #psia total pressure
T=1000 #degree R Temprature
px=12 #psia
k=1.4
R=53.3 #ft*lb/(lbm*degree R)
pratio=p/px
#for this value of pratio, Max is calculated as(from fig D.4)
Max=1.9
#using this value of Max, Tx/T0,x is found as
Tratio=0.59
#T=T0,x=T0,y
#Calculation
Tx=Tratio*T #degree R
cx=(R*Tx*k)**0.5 #ft/sec
Vx=1.87*cx*(32.2**0.5) #ft/sec
#Result
print "The Mach number for the flow=",round(Max,1)
print "The velocity of the flow=",round(Vx,0),"ft/sec"
Example 11.19 page no.671¶
In [10]:
x1=0.5 #m exit of duct
x2=0.3 #m entring of duct
Acritical=0.1 #m**2 area
#at x1, Max1 is found as
Max1=2.8
#and px/p0,x is found as
pratio1=0.04
#For this value of Max, py/px is found as
prat1=9
pfraction1=prat1*pratio1 #where pfraction=py/p0,x = pIII/p0,x
#at x2, Max2 is found as
Max2=2.14
#for this value of Max2, the following are found
prat2=5.2
prat22=0.66 #where prat22=p0,y/p0,x
May=0.56
#for this valur of May, Ay/Acritical is found as
Aratio=1.24
Arat=(Acritical+(x1**2))/(Acritical+(x2**2)) #where Aratio=A2/Ay
Afraction=Aratio*Arat #where Afraction=A2/Acritical
A2=Acritical+(x1**2) #m**2
Acritical1=A2/Afraction #where Acritical1 critical area for the isentropic flow downstream of the shock
#with the value of Afraction, the following are found
Ma2=0.26
pfraction=0.95 #where pfraction=p2/p0,y
#hence,
pfrac=pfraction*prat22 #where pfrac=p2/p0,x
#Result
print "The ratio of back pressure to inlet stagnation pressure \n that will result in a normal shock at the exit of the duct=",round(pfraction1,2)
print "The value of back pressure to inlet stagnation pressure required to position the shock at (x=0.3 m)=",round(pfrac,2)
#Plot
s=[0,0]
T=[288,112]
s1=[0,360]
T1=[288,288]
s2=[0,280]
T2=[112,275]
s3=[240,320]
T3=[275,275]
s4=[-10,80]
T4=[112,112]
s5=[280,280]
T5=[288,275]
a=plot(s,T)
a1=plot(s1,T1)
a2=plot(s2,T2,linestyle='--')
a3=plot(s3,T3)
a4=plot(s4,T4)
a5=plot(s5,T5)
xlabel("s-s1 (J/Kg K)")
ylabel("T (K)")
plt.xlim((-40,500))
plt.ylim((60,340))
plt.text(82,112,'Tx=112K')
plt.text(60,310,'p0x=101kpa(abs)')
plt.text(300,315,'p0y=38kpa(abs)')
plt.text(340,270,'py=36kpa(abs)')
plt.text(10,275,'0,x')
plt.text(10,110,'x')
plt.text(300,290,'0,y')
plt.text(280,270,'y')
plt.text(280,375,'Ty=275K')
plt.text(350,285,'Tox=Toy=288K')
plt.text(40,125,'px=4kpa(abs)')
plt.text(80,90,'shock at nozzle exit plane x=(0.5m)')
show(a)
show(a1)
show(a2)
show(a3)
show(a4)
show(a5)
# For 2nd flow
s=[0,0]
T=[288,150]
s1=[0,360]
T1=[288,288]
s2=[80,200]
T2=[284,284]
s3=[80,160]
T3=[271,271]
s4=[0,120]
T4=(150,271)
s5=[-10,80]
T5=[150,150]
s6=[120,120]
T6=[271,284]
s7=[120,250]
T7=[288,320]
s8=[120,250]
T8=[284,300]
a=plot(s,T)
a1=plot(s1,T1)
a2=plot(s2,T2)
a3=plot(s3,T3)
a4=plot(s4,T4,linestyle='--')
a5=plot(s5,T5)
a6=plot(s6,T6)
a7=plot(s7,T7)
a8=plot(s8,T8)
plt.text(100,140,'Tx=150K')
plt.text(40,310,'p0x=101 kpa(abs)')
plt.text(250,320,'p0y=67kpa(abs)')
plt.text(120,250,'py=52kpa(abs)')
plt.text(10,275,'0,x')
plt.text(10,140,'x')
plt.text(120,290,'0,y')
plt.text(80,270,'y')
plt.text(160,265,'Ty=271K')
plt.text(200,280,'T2=284K')
plt.text(380,288,'Tox=Toy=288K')
plt.text(40,160,'px=10kpa(abs)')
plt.text(240,300,'p2=64kpa(abs)')
plt.text(80,100,'shock with nozzle x=(0.3m)')
plt.text(20,200,'normal shock')
xlabel("s-s1 (J/Kg K)")
ylabel("T (K)")
plt.xlim((-40,500))
plt.ylim((60,340))
show(a)
show(a1)
show(a2)
show(a3)
show(a4)
show(a5)
show(a6)
show(a7)
show(a8)
