import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Welcome to Aerospace Propulsion. The full Github repository is available at https://github.com/wlamont/aerospace_propulsion. Please feel free to contribute
With sufficient thrust, pigs fly just fine.
-R. Callon, The Twelve Networking Truths
The background required for the subsequent chapters falls into two categories: gas dynamics and combustion. If the reader has experience with these subjects this chapter can be safely skipped. Note both gas dynamics and combustion are vast fields with many textbooks and full courses dedicated to each. For a more in-depth treatment of gas dynamics I would recommend Prof. Zucrow's Gas Dynamics and for a good detailed book on combustion I would recommend Prof. Glassman's Combustion.
Gas dynamics is the overlapping area of fluid mechanics and thermodynamics and is concerned with the motion of gases and its effects on physical structures. Before we begin, we must define our approach (statistical or continuum) and perspective (Eulerian or Lagrangian) for defining the governing equations and our framework for solving problems. For this text, we will use the continuum approach -- we will ignore individual molecules and treat the substance as being continuously distributed in space. The continuum approach is valid when the smallest macroscopic length scale is much larger than the microscopic length scale (usually the mean free path of a molecule). Conversely, the statistical approach considers individual molecules and the solution requires probability and statistics (this is the approach taken in Statistical Mechanics). In terms of perspective, a moving fluid can be studied by either looking at a control volume (Eulerian point of view) or a control mass (Lagrangian point of view). In the control volume (Eulerian) perspective, a stationary control volume is created in space and we observe the change in fluid passing through the control volume. In the control mass (Lagrangian) perspective, we consider a specific mass as it moves and changes with time. Note: the control mass is allowed to deform and change with time but once defined the control volume is not. For this text, we will use the control volume (Eulerian) perspective.
The most important parameter in aerospace propulsion (and arguably gas dynamics) is the Mach number (M), which is a dimensionless number representing the ratio of the flow velocity (V) to the local speed of sound (a):
As we will see, the Mach number can be used to recast the governing equations into forms that are better suited for aerospace propulsion. The Mach number can also be used to classify flows:
Before we define the speed of sound we first need to introduce the equation of state which describes the relation between various thermodynamic properties. For this text we will use the ideal gas law:
where P is the pressure, V is the volume, n is the number of moles, Ru is the universal gas constant(Ru=8314 kJ /(K⋅mol)), and T is the temperature. This equation can be recast into a more useful form by: (1) replacing the number of moles (n) with the mass (m) divided by the molar mass (MW), (2) dividing the mass (m) by the volume (V) and replacing by the density (ρ=m/V), and (3) replacing the universal gas constant by the specific gas constant (R=Ru/MW):
The conservation of mass (also sometimes called continuity) states that matter can not be created nor destroyed. For control volume analysis, the conservation of mass states that mass change within the control volume must equal the net rate at which mass enters and exits the control volume:
For a steady state process, as depicted in figure, the time dependence drops out and the conservation of mass reduces to the simple:
The mass flow rate by definition is:
where ρ is the density, u is the velocity and A is the area normal to the flow.
The above two equations can be programmed into Python to create the D Function.
def d_function(gamma, MW, M):
'''
Calculates D(M) and D(M)/dM
where
gamma is the ratio of specific heat
MW is the molecular weight
M is the Mach Number
'''
D = M / (1+(gamma-1)/2*M**2) ** ((gamma+1)/(2*gamma-2))
dD = D / M*(1-M**2) / (1+(gamma-1)/2*M**2)
return D, dD
gamma_air = 1.4
MW_air = 28.97
M_air = []
D_air = []
for M in np.arange(0.1, 10, 0.1):
d = d_function(gamma_air, MW_air, M)
D_air.append(d[0])
M_air.append(M)
plt.plot(M_air, D_air, 'k')
plt.xlabel('Mach number [-]')
plt.xticks(range(1, 11, 1))
plt.ylabel('D(M) [-]')
plt.grid()
where ∑Fx is the total force on the control volume and comprised of pressure and viscous forces:
def g_function(gamma, MW, M):
'''
Calculates G(M)
where
gamma is the ratio of specific heat
MW is the molecular weight
M is the Mach Number
'''
G = (1 + gamma * M **2) / (1 + (gamma - 1) / 2 * M **2) ** (gamma/(gamma-1))
return G
M_air = []
G_air = []
for M in np.arange(0, 10, 0.1):
g = g_function (gamma_air, MW_air, M)
G_air.append(g)
M_air.append(M)
plt.plot(M_air, G_air, 'k')
plt.xlabel('Mach number [-]')
plt.xticks(range(1, 11, 1))
plt.ylabel('G(M) [-]')
plt.grid()
Heat addition in constant area duct: A1=A2
Continuity:
Momentum Equation:
Divide continuity by momentum equation:
where: N=D/G
Use Newton's Method:
def n_function(gamma, MW, M):
'''
Calculates N(M)
where
gamma is the ratio of specific heat
MW is the molecular weight
M is the Mach Number
'''
N = M*(1+(gamma-1)/2*M**2)**0.5/(1+gamma*M**2)
dN = N*(1/M+((gamma-1)/2)*M/(1+(gamma-1)/2*M**2)-2*gamma*M/(1+gamma*M**2))
return N, dN
M_air = []
N_air = []
for M in np.arange(0.1, 10, 0.1):
n, dn = n_function (gamma_air, MW_air, M)
N_air.append(n)
M_air.append(M)
plt.plot(M_air, N_air, 'k')
plt.xlabel('Mach number [-]')
plt.xticks(range(1, 11, 1))
plt.ylabel('N(M) [-]')
plt.grid()
T1 = 300 #inlet temperature
T2 = 2000 #outlet temperature
gamma = 1.4;
M1 = 0.1
T1o = T1*(1+(gamma-1)/2*M1**2)
N1, dN1 = n_function (gamma_air, MW_air, M1)
if M1 > 1:
scale = 0.9
else:
scale = 1.05
M2 = M1 * scale # guess M2
T2o = T2*(1+(gamma-1)/2*M2**2)
N2, dN2 = n_function (gamma_air, MW_air, M2)
ERROR = []
M = []
M.append(M2)
ERROR.append(np.abs(N1/T1o**0.5 - N2/T2o**0.5))
while ERROR[-1] > 1e-8:
M2new = M2 + (N1*((T2o)/T1o)**0.5-N2)/dN2
if M2new > 1.1*M2:
M2new = 1.1*M2
elif M2new < 0.9 * M2:
M2new = 0.9 * M2
M2 = M2new
M.append(M2)
T2o = T2*(1+(gamma-1)/2*M2**2)
N2, dN2 = n_function (gamma_air, MW_air, M2)
ERROR.append(np.abs(N1/T1o**0.5 - N2/T2o**0.5))
plt.plot(M, ERROR)
print 'M2 = ', M2
print 'N1 / T1o**0.5=', N1 / T1o**0.5
print 'N2 / T2o**0.5=', N2 / T2o**0.5
M2 = 0.283231351005 N1 / T1o**0.5= 0.00569378963698 N2 / T2o**0.5= 0.00569378765277
from IPython.core.display import HTML
def css_styling():
styles = open("/home/wlamont/Documents/Textbooks/aerospace_prop_python/styles/custom.css", "r").read()
return HTML(styles)
css_styling()