# Example 3.1.py
# At a point in the flow over an F-15 high performance fighter airplane,
# the pressure, temperature and Mach number are 1890 lb/ft^2, 450 deg R,
# and 1.5, respectively. At this point calculate To, po, T*, p* and the
# flow velocity.
# Variable declaration
p = 1890.0 # pressure (lb/ft^2)
T = 450.0 # temperature (R)
M = 1.5 # Mach number
R = 1716.0 # gas constant (ft lbf/slug/R)
gamma = 1.4 # ratio of specific heat capacity for air
# Calculations
# from isentropic relations
def temperature_ratio(M):
return 1 + (gamma-1)*M*M/2.0
def pressure_ratio(M):
return pow(temperature_ratio(M), gamma/(gamma-1))
po_by_p = pressure_ratio(M) # ratio of stagnation pressure to actual pressure at mach 1.5
To_by_T = temperature_ratio(M) # ratio of stagnation temperature to actual temperature at mach 1.5
po = po_by_p * p # stagnation pressure po = po / p * p (lb/ft^2)
To = To_by_T * T # stagnation temperature To = To / T * T (R)
po_by_pstar = pressure_ratio(1.0) # ratio of stagnation pressure to actual pressure(star) at mach 1.0
To_by_Tstar = temperature_ratio(1.0) # ratio of stagnation temperature to actual temperature(star) at mach 1.0
pstar = po_by_p / po_by_pstar * p # pstar(lb/ft^2), pstar/po * po/p * p
Tstar = To_by_T / To_by_Tstar * T # Tstar(R), Tstar/To * To/T * T
a = pow(gamma*R*T, 0.5) # speed of sound(ft/s)
V = M*a # velocity(ft/s)
# Result
print "To = %.2f R, po = %.2f lb/ft^2, Tstar = %.2f R, pstar = %.2f lb/ft^2" %(To,po,Tstar,pstar)
print "Velocity at the required point is %.2f ft/s" %(V)
# Example 3.2.py
# Return to Example 1.6, Calculate the Mach Number and velocity at the exit of the rocket
# nozzle.
# Variable declaration from example 1.6
pc = 15.0 # pressure combustion chamber (atm)
Tc = 2500.0 # temperature combustion chamber (K)
mol_wt = 12.0 # molecular weight (gm)
cp = 4157.0 # specific heat at constant pressure (J/Kg/K)
Tn = 1350.0 # temperature at nozzle exit (K)
# Calculations
R = 8314.0/mol_wt # gas constant = R_prime/mo_wt, R_prime = 8314 J/K
cv = cp - R # specific heat at constant volume (J/Kg k)
gamma = cp/cv # ratio of specific heat
pn_by_pc = pow(Tn/Tc, gamma/(gamma-1)) # ratio of pressure for isentropic process, pn/pc
Mn = pow(2/(gamma-1)*(pow(1/pn_by_pc,(gamma-1)/gamma) - 1), 0.5) # Mach number at exit, from isentropic flow relation
an = pow(gamma*R*Tn, 0.5) # Speed of sound at exit (m/s)
Vn = Mn*an # Velocity at exit (m/s)
# Result
print "Mach number at the exit of the rocket nozzle is %.3f"%(Mn)
print "Velocity at the exit of the rocket nozzle is %.1f m/s"%(Vn)
# Example 3.3.py
# Return to Example 1.1, calculate the percentage density change between the given point
# on the wing and the free-stream, assuming compressible flow.
# Variable declaration from example 1.1
rho_1 = 0.002377 # density at sea level (slug/ft^3)
T_1 = 519.0 # temperature at sea level (R)
v_1 = 100.0 # velocity far ahead of the wing (mi/h)
v_2 = 150.0 # velocity at some point on the wing (mi/h)
gamma = 1.4 # ratio of specific heat capacity for air
R = 1716.0 # gas constant (ft lbf/slug/R)
# Calculations
cp = gamma*R/(gamma-1) # specific heat capacity at constant pressure (ft lb/ slug / R)
u_1 = v_1 * 88.0/60.0 # converting v_1 in ft/s
u_2 = v_2 * 88.0/60.0 # converting v_2 in ft/s
T_2 = T_1 + (u_1*u_1 - u_2*u_2)/cp/2.0 # temperature at a point from energy equation (R)
rho_2_by_rho_1 = pow((T_2/T_1), 1/(gamma-1))# density ratio from isentropic flow relation
rho_2 = rho_2_by_rho_1 * rho_1 # density at the point (slug/ ft^3)
delrho = rho_1 - rho_2 # change in density (slug/ ft^3)
fracrho = delrho/rho_1 # fractional change in density
# Result
print "Percentage change in density is %.1f"%(fracrho*100)
# Example 3.4.py
# A normal shock wave is standing in the test section of a supersonic wind tunnel.
# Upstream of the wave, M1 = 3, p1 = 0.5 atm, and T1 = 200 K. Find M2, p2, T2 and
# u2 downstream of the wave
# Variable declaration from example 1.1
M1 = 3.0 # upstream mach number
p1 = 0.5 # upstream pressure (atm)
T1 = 200.0 # upstream temperature (K)
R = 287.0 # gas constant (J/Kg/K)
gamma = 1.4 # ratio of specific heats for air
# Calculations
# from shock relation (Table A2) for M1 = 3.0
# subscript 2 means downstream of the shock
p2_by_p1 = 10.33 # p2/p1
T2_by_T1 = 2.679 # T2/T1
M2 = 0.4752 # M2
p2 = p2_by_p1 * p1 # downstream pressure (atm)
T2 = T2_by_T1 * T1 # downstream temperature (K)
a2 = pow(gamma*R*T2, 0.5) # speed of sound downstream of the shock (m/s)
u2 = M2*a2 # downstream velocity (m/s)
# Result
print "M2 = %.4f"%(M2)
print "p2 = %.3f atm"%(p2)
print "T2 = %.1f K"%(T2)
print "u2 = %.1f m/s"%(u2)
# Example 3.5.py
# A blunt-nosed missile is flying at Mach 2 at standard sea level. Calculate the
# temperature and pressure at the nose of the missile.
# Variable declaration
M1 = 2.0 # upstream mach number
p1 = 2116.0 # upstream pressure(sea level) (lb/ft^2)
T1 = 519.0 # upstream temperature(sea level) (R)
gamma = 1.4 # ratio of specific heat
# Calculations
# subscript 2 means downstream of the shock
# from isentropic relations
def temperature_ratio(M):
return 1 + (gamma-1)*M*M/2.0
def pressure_ratio(M):
return pow(temperature_ratio(M), gamma/(gamma-1))
To1_by_T1 = temperature_ratio(M1) # ratio of stagnation temperature to actual temperature at mach M1
po1_by_p1 = pressure_ratio(M1) # ratio of stagnation pressure to actual pressure at mach M1
To2 = To1 = To1_by_T1 * T1 # To2 and To1 = To1/T1 * T1 (in R)
# from shock relation (Table A2) for M1 = 2.0
po2_by_po1 = 0.7209 # po2/po1
po2 = po2_by_po1 * po1_by_p1 * p1 # po2 = po2/po1 * po1/p1 * p1 (in lb/ft^2)
# Result
print "Pressure at the nose of the missile is %.1f lb/ft^2" %(po2)
print "Temperature at the nose of the missle is %.1f deg R"%(To2)
# Example 3.6.py
# Consider a point in a supersonic flow where the static pressure is 0.4 atm. When
# a pitot tube is inserted in the float at this point, the pressure measured by the
# pitot tube is 3 atm. Calculate the mach number at this point.
# Variable declaration
p1 = 0.4 # static pressure (in atm)
po2 = 3.0 # pressure measured by the pitot tube (in atm)
#Calculations
# from table A2 for po2/p1 = 7.5
M1 = 2.35
# Results
print "Mach number is %.2f"%(M1)
# Example 3.7.py
# For the normal shock that occurs in front of the pitot tube in Example 3.6,
# calculate the entropy change across the shock.
# Variable declaration
M1 = 2.34 # mach number from example 3.6
R = 1716.0 # gas constant (ft lbf/slug/R)
# Calculations
# subscript 2 means downstream of the shock
po2_by_po1 = 0.5615 # from shock table A2 for mach M1
from numpy import log # importing log
dels = -R*log(po2_by_po1) # s2 - s1 (in lb/slug R)
# Result
print "Change is entropy is %.1f lb/slug R" %(dels)
# Example 3.8.py
# Air enters a constant area duct at M1 = 0.2, p1 = 1 atm, and T1 = 273 K.
# Inside the duct the heat is added per unit mass is q = 1.0 * 10^6 J/kg.
# Calculate the flow properties M2, p2, T2, rho2, To2, po2 at the exit of the
# duct.
# Variable declaration
M1 = 0.2 # upstream mach number
p1 = 1.0 # upstream pressure (in atm)
T1 = 273.0 # upstream temperature (in K)
q = 1e6 # heat added (in J/Kg)
R = 287.0 # gas constant (in J/Kg K)
gamma = 1.4 # ratio of specific heats
# from Table A3, for M1 = 0.2
T1_by_Tstar = 0.2066
p1_by_pstar = 2.273
po1_by_postar = 1.235
To1_by_Tostar = 0.1736
# Calculations
# from isentropic relations
def temperature_ratio(M):
return 1 + (gamma-1)*M*M/2.0
def pressure_ratio(M):
return pow(temperature_ratio(M), gamma/(gamma-1))
po1_by_p1 = pressure_ratio(M1) # ratio of stagnation pressure to actual pressure at mach M1
To1_by_T1 = temperature_ratio(M1) # ratio of stagnation temperature to actual temperature at mach M1
po1 = po1_by_p1 * p1 # po1 = po1/p1 * p1 (in atm)
To1 = To1_by_T1 * T1 # To1 = To1/T1 * T1 (in K)
cp = gamma * R / (gamma-1) # specific heat at constant pressure (J/ Kg K)
To2 = q / cp + To1 # total temperature after heat addition (K)
# hence
To2_by_Tostar = To2/To1 * To1_by_Tostar # To2/Tostar = To2/To1 * To1/Tostar
M2 = 0.58 # from table A3
# from Table A3, for M2 = 0.58
T2_by_Tstar = 0.8955
p2_by_pstar = 1.632
po2_by_postar = 1.083
T2 = T2_by_Tstar / T1_by_Tstar * T1 # T2 = T2/Tstar * Tstar/T1 * T1 (in K)
p2 = p2_by_pstar / p1_by_pstar * p1 # p2 = p2/pstar * pstar/p1 * p1 (in atm)
po2 = po2_by_postar / po1_by_postar * po1 # po2 = po2/postar * postar/po1 * po1 (in atm)
rho2 = p2*101325.0 / R / T2 # rho (in Kg/m^3) = P/(RT), where p is in N/m^2
# Result
# M2, p2, T2, rho2, To2, po2 at the exit
print "M2 = %.2f "%(M2)
print "p2 = %.3f atm"%(p2)
print "T2 = %.2f K"%(T2)
print "rho2 = %.3f Kg/m^3"%(rho2)
print "To2 = %.1f K"%(To2)
print "po2 = %.3f atm"%(po2)
# Example 3.9.py
# Air enters a constant area duct at M1 = 3.0, p1 = 1 atm, and T1 = 300 K.
# Inside the duct the heat is added per unit mass is q = 3.0 * 10^5 J/kg.
# Calculate the flow properties M2, p2, T2, rho2, To2, po2 at the exit of the
# duct.
# Variable declaration
M1 = 3.0 # upstream mach number
p1 = 1.0 # upstream pressure (in atm)
T1 = 300.0 # upstream temperature (in K)
q = 3e5 # heat added (in J/Kg)
R = 287.0 # gas constant (in J/Kg K)
gamma = 1.4 # ratio of specific heats
# from Table A3, for M1 = 3.0
T1_by_Tstar = 0.2803
p1_by_pstar = 0.1765
To1_by_Tostar = 0.6540
po1_by_postar = 3.424
# Calculations
# from isentropic relations
def temperature_ratio(M):
return 1 + (gamma-1)*M*M/2.0
def pressure_ratio(M):
return pow(temperature_ratio(M), gamma/(gamma-1))
po1_by_p1 = pressure_ratio(M1) # ratio of stagnation pressure to actual pressure at mach M1
To1_by_T1 = temperature_ratio(M1) # ratio of stagnation temperature to actual temperature at mach M1
po1 = po1_by_p1 * p1 # po1 = po1/p1 * p1 (in atm)
To1 = To1_by_T1 * T1 # To1 = To1/T1 * T1 (in K)
cp = gamma * R / (gamma-1) # specific heat at constant pressure (J/ Kg K)
To2 = q / cp + To1 # total temperature after heat addition (K)
# hence
To2_by_Tostar = To2/To1 * To1_by_Tostar # To2/Tostar = To2/To1 * To1/Tostar
M2 = 1.58 # from table A3
# from Table A3, for M2 = 1.58
T2_by_Tstar = 0.7117
p2_by_pstar = 0.5339
po2_by_postar = 1.164
T2 = T2_by_Tstar / T1_by_Tstar * T1 # T2 = T2/Tstar * Tstar/T1 * T1 (in K)
p2 = p2_by_pstar / p1_by_pstar * p1 # p2 = p2/pstar * pstar/p1 * p1 (in atm)
po2 = po2_by_postar / po1_by_postar * po1 # po2 = po2/postar * postar/po1 * po1 (in atm)
rho2 = p2*101325.0 / R / T2 # rho (in Kg/m^3) = P/(RT), where p is in N/m^2
# Result
# M2, p2, T2, rho2, To2, po2 at the exit
print "M2 = %.2f "%(M2)
print "p2 = %.3f atm"%(p2)
print "T2 = %.2f K"%(T2)
print "rho2 = %.3f Kg/m^3"%(rho2)
print "To2 = %.1f K"%(To2)
print "po2 = %.3f atm"%(po2)
# Example 3.10.py
# In example 3.9, how much heat per unit mass must be added to choke the flow?
# Variable declaration from example 3.9
To1 = 840 # upstream total temperature (in K)
M1 = 3.0 # upstream mach number
To1_by_Tostar = 0.6540 # To1/Tostar from Table A3
cp = 1004.5 # specific heat at constant pressure for air (in J/Kg K)
# Calculations
Tostar = To1 / To1_by_Tostar # Tostar = To1 * Tostar/To1 (in K)
M2 = 1.0 # for choked flow
To2 = Tostar # since M2 = 1.0
q = cp * (To2 - To1) # required heat = cp(To2 - To1) (in J/kg)
# Result
print "Heat require to choke the flow is %.2e J/kg" %(q)
# Example 3.11.py
# Consider the flow of air through a pipe of diameter = 0.15 m and length = 30 m.
# The inlet flow conditions are M1 = 0.3, p1 = 1.0 atm and T1 = 273 K. Assuming
# f = constant = 0.005, calculate the flow conditions at the exit, M2, p2, T2 and
# po2.
# Variable declaration
M1 = 0.3 # upstream mach number
p1 = 1.0 # upstream pressure (in atm)
T1 = 273.0 # upstream temperature (in K)
R = 287.0 # gas constant (in J/Kg K)
gamma = 1.4 # ratio of specific heats
D = 0.15 # diameter of pipe (in m)
L = 30.0 # length of pipe (in m)
f = 0.005 # friction coefficient
# Calculations
# from isentropic relations
def temperature_ratio(M):
return 1 + (gamma-1)*M*M/2.0
def pressure_ratio(M):
return pow(temperature_ratio(M), gamma/(gamma-1))
po1_by_p1 = pressure_ratio(M1) # ratio of stagnation pressure to actual pressure at mach M1
To1_by_T1 = temperature_ratio(M1) # ratio of stagnation temperature to actual temperature at mach M1
po1 = po1_by_p1 * p1 # po1 = po1/p1 * p1 (in atm)
To1 = To1_by_T1 * T1 # To1 = To1/T1 * T1 (in K)
# from Table A4 for M1 = 0.3
# C = 4*f*Lstar/D
C1 = 5.299
p1_by_pstar = 3.619 # p1 / pstar
T1_by_Tstar = 1.179 # T1 / Tstar
po1_by_pstar = 2.035 # po1 / pstar
C2 = C1 - 4*f*L/D
# from Table A4 for 4*f*Lstar/D = 1.2993
M2 = 0.475
T2_by_Tstar = 1.148 # T2 / Tstar
p2_by_pstar = 2.258 # p2 / pstar
po2_by_postar = 1.392 # po2 / postar
T2 = T2_by_Tstar / T1_by_Tstar * T1 # T2 = T2/Tstar * Tstar/T1 * T1 (in K)
p2 = p2_by_pstar / p1_by_pstar * p1 # p2 = p2/pstar * pstar/p1 * p1 (in atm)
po2 = po2_by_postar / po1_by_pstar * po1 # po2 = po2/postar * postar/po1 * po1 (in atm)
# Results
print "M2 = %.3f "%(M2)
print "p2 = %.3f atm"%(p2)
print "T2 = %.2f K"%(T2)
print "po2 = %.3f atm"%(po2)
# Example 3.12.py
# Consider the flow of air through a pipe of diameter = 0.4 ft and length = 5 ft.
# The inlet flow conditions are M1 = 3.0, p1 = 1.0 atm and T1 = 300 K. Assuming
# f = constant = 0.005, calculate the flow conditions at the exit, M2, p2, T2 and
# po2.
# Variable declaration
M1 = 3.0 # upstream mach number
p1 = 1.0 # upstream pressure (in atm)
T1 = 300.0 # upstream temperature (in K)
gamma = 1.4 # ratio of specific heats
D = 0.4 # diameter of pipe (in ft)
L = 5.0 # length of pipe (in ft)
f = 0.005 # friction coefficient
# Calculations
# from isentropic relations
def temperature_ratio(M):
return 1 + (gamma-1)*M*M/2.0
def pressure_ratio(M):
return pow(temperature_ratio(M), gamma/(gamma-1))
po1_by_p1 = pressure_ratio(M1) # ratio of stagnation pressure to actual pressure at mach M1
To1_by_T1 = temperature_ratio(M1) # ratio of stagnation temperature to actual temperature at mach M1
po1 = po1_by_p1 * p1 # po1 = po1/p1 * p1 (in atm)
To1 = To1_by_T1 * T1 # To1 = To1/T1 * T1 (in K)
# from Table A4 for M1 = 3.0
# C = 4*f*Lstar/D
C1 = 0.5222
p1_by_pstar = 0.2182 # p1 / pstar
T1_by_Tstar = 0.4286 # T1 / Tstar
po1_by_pstar = 4.235 # po1 / pstar
C2 = C1 - 4*f*L/D
# from Table A4 for 4*f*Lstar/D = 0.2722
M2 = 1.9
T2_by_Tstar = 0.6969 # T2 / Tstar
p2_by_pstar = 0.4394 # p2 / pstar
po2_by_postar = 1.555 # po2 / postar
T2 = T2_by_Tstar / T1_by_Tstar * T1 # T2 = T2/Tstar * Tstar/T1 * T1 (in K)
p2 = p2_by_pstar / p1_by_pstar * p1 # p2 = p2/pstar * pstar/p1 * p1 (in atm)
po2 = po2_by_postar / po1_by_pstar * po1 # po2 = po2/postar * postar/po1 * po1 (in atm)
# Results
print "M2 = %.1f "%(M2)
print "p2 = %.3f atm"%(p2)
print "T2 = %.1f K"%(T2)
print "po2 = %.2f atm"%(po2)
# Example 3.13.py
# In example 3.12, what is the length of the duct required to choke the flow?
# Variable declaration from example 3.12
M1 = 3.0 # mach number
C1 = 0.5222 # C1 = 4*f*L1star/D
f = 0.005 # friction coefficient
D = 0.4 # diameter of pipe (in ft)
# Calculations
L1star = 0.5222 * D/4.0/f
# Result
print "Length required to choke the flow is %.2f ft" %(L1star)