In this notebook we will study the effect of a shock on a choked converging diverging nozzle. After the flow reaches the throat the flow will transfer to supersonic, meaning the flow is travelling faster than the speed of sound and there can be no more communication between the inlet and exit as it can only move as fast as the speed of sound. This means that fluids that normally react to changes in pressure and temperature can no longer do so. For sound waves at supersonic speeds the source is travelling faster than the waves can propogate (at the speed of sound) and the waves produced by the source will constructvely interfere and form areas of high pressure that lead to a shock. As the pressure outside the nozzle (back pressure) is decreased the shock will move closer to the exit, until the back pressure is lower than the pressure at the throat, at which point a shock will no longer occur within the nozzle and the flow will either be underexpanded, exit pressure is higher than back pressure, or overexpanded, exit pressure is lower than back pressure. The diagram below shows the different cases that can occur in the nozzle, with case C being what we will study
In this problem the nozzle is frictionless this means that both before and after the shock the flow is isentropic – no change in entropy / no heat transfer – therefore up until the shock the flow is isentropic and afterwards the flow is also isentropic, across the shock however the flow is not isentropic and the entropy increases.
Sources: Nozzle Diagram
#import relevant libraries
import numpy as np
import pandas as pd
from scipy.optimize import root_scalar
from scipy import interpolate
import plotly.graph_objects as go
import plotly.figure_factory as ff
import plotly.io as pio
pio.renderers.default = 'iframe'Let us create functions for the duct curvature.
# create functions for duct curvature and an array of x values to use in graphing
x = np.linspace(0, 1, 200)
alpha = (x/2 - 0.3)**2 + 0.05
beta = -alpha
# define given variables and constants
p_0i = 140 # kPa
T_0i = 290 # K
R = 287 # J/kgK
c_p = 1006 # J/kgK
g = 1.4
shock_x = 0.8
throat_x = 0.6
# assume that flow is choked so the throat Mach number can be defined as 1
M_t = 1Information to keep in mind about the flow: - The flow is given to be in a frictionless duct and is a perfect gas - The flow before and after the shock will be isentropic and therefore isentropic flow relations can be used - Across a normal shock the stagnation pressure will drop and stagnation temperature will remain the same - Due to conservation of mass the mass flow rate is constant through out the duct
The curvature of the top of the nozzle is given by $ (x) $ and the curvature of the bottom is the negative of the top: \[ \alpha(x) = \left(\frac{x}{2}-0.3\right)^{2} + 0.05 = -\beta(x) \]
\(x\) is given to be \(x \in [0, 1]\)
The area distribution above is captured below
# Define a function that solves for an area per unit width (can also solve for different areas such as circular)
def area(x):
area = 2 * ((x/2 - 0.3)**2 + 0.05)
return areaTo find the Mach number in the nozzle before the shock we can use the isentropic area ratio Mach number relation[1]:
\[ \large \require{color}\frac{A}{A^{*}} = \frac{\left(1 + \frac{\gamma+1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}}}{{\color[rgb]{0.041893,0.355669,0.727621}M}} \left(\frac{\gamma +1}{2}\right)^{-\frac{\gamma +1}{2(\gamma -1)}} \]
There is no analytical solution to find the Mach number from the area ratio so another method must be used, in this notebook optimization is used
# Define the area at the throat, which is needed in order to solve for Mach number
A_t = area(throat_x)
# Define a mach function to calculate the roots of the area ratio mach relation set to zero
def mach_function(M, x_pos):
A = area(x_pos)
return np.sqrt((1/M**2)*((2/(g+1))*(1+((g-1)/2)*M**2))**((g+1)/(g-1))) - A/A_tIn order to define the Mach number after the shock based on the Mach number before the shock we need to use the relationship [2]:
\[ \large \require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_{as}^2 = \frac{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_s^2 + 2}{2\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_s^2 - (\gamma -1)} \]
Where \({\color[rgb]{0.041893,0.355669,0.727621}M}_{as}\) is the Mach number after the shock and \({\color[rgb]{0.041893,0.355669,0.727621}M}_{s}\) is the Mach number directly before the shock.
# Define the Mach number at the shock position, min Mach value at the entrance, and Mach number after shock
M_s_root = root_scalar(mach_function, args = (shock_x), x0 = 1, x1 = 2.5)
M_s = M_s_root.root
M_min_root = root_scalar(mach_function, args = (0), x0 = 0.2, x1 = 1)
M_min = M_min_root.root
M_as = np.sqrt(((g-1)*M_s**2 + 2)/(2*g*M_s**(2)-(g-1)))To find the Mach number after the shock we need the mass flow rate[3]:
\[ \large \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}} = \frac{A{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}_{0}}} \sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}} {\color[rgb]{0.041893,0.355669,0.727621}M} \left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-\frac{\gamma+1}{2(\gamma -1)}} \]
This equation also cannot be solved analytically so the same process as before can be used. Since the mass flow rate must remain constant to obey conservation of mass and the stagnation temperature also remains constant, the area and the stagnation pressure are the only changing variables
The stagnation pressure will remain constant from the inlet until the flow reaches the shock. As the stagnation pressure after the shock only depends on the Mach number at the shock and the inlet stagnation pressure we can go ahead and define a function to calculate stagnation pressure[4] based on the normal shock relation:
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_{0s} = {\color[rgb]{0.315209,0.728565,0.037706}p}_{0i} \left(\frac{(\gamma+1){\color[rgb]{0.041893,0.355669,0.727621}M}_{s}^{2}}{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_{s}^{2}-2}\right)^{\frac{\gamma}{\gamma-1}} \left(\frac{\gamma+1}{2 \gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_{s}^{2} - (\gamma -1)}\right)^{\frac{1}{\gamma-1}}, \]
where \({\color[rgb]{0.041893,0.355669,0.727621}M}_{s}\) is the value of the Mach number at the shock.
# Define a function to calculate the stagnation pressure before and after the shock -- based on shock relations
def stag_pressure(x_pos):
if x_pos < shock_x:
stag_pressure = p_0i
else:
stag_pressure = p_0i * ((((g + 1) * M_s**2) / ((g - 1) * M_s**2 +2))**(g/(g - 1))) * ((g + 1) / (2*g* M_s**2 - (g -1)))**(1/(g-1))
return stag_pressure
# Define the mass flow rate at the throat where M=1 for ease, this number will remain constant
m_dot = ((A_t)*(p_0i*1000)/(np.sqrt(T_0i)))*np.sqrt(g/R)*((g+1)/2)**(-(g+1)/(2*(g-1)))
#p0 is multipled by 1000 because it is given in kPa and we need units of Pa
# Define a function that finds the roots of the mass flow rate equation set to zero
def mach_function_as(M, x_pos):
A = area(x_pos)
p_0 = stag_pressure(x_pos)*1000
return (A*p_0/np.sqrt(T_0i))*(np.sqrt(g/R))*M*(1+((g-1)/2)*M**2)**(-(g+1)/(2*(g-1))) - m_dot# Define a function to calculate the Mach value at an x value, x0 and x1 are lower and upper limits on the Mach number that
# are based on location in the nozzle
def mach_number(x_pos):
if x_pos <= 0.5:
M_root = root_scalar(mach_function, args = (x_pos), x0 = M_min, x1 = 0.8)
M = M_root.root
elif 0.5 < x_pos <= throat_x:
M_root = root_scalar(mach_function, args = (x_pos), x0 = M_min, x1 = 1)
M = M_root.root
elif x_pos < shock_x:
M_root = root_scalar(mach_function, args = (x_pos), x0 = 1, x1 = M_s)
M = M_root.root
else:
M_root = root_scalar(mach_function_as, args = (x_pos), x0 = M_as, x1 = 0.2)
M = M_root.root
return MAfter you have a function that defines the stagnation pressure and Mach number the static pressure [5] can be calculated using isentropic flow relations, as flow is only momentarily nonisentropic right at the shock so there is a discontinuity:
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p} = {\color[rgb]{0.315209,0.728565,0.037706}p}_{0} \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-\frac{\gamma}{\gamma-1}} \]
# Define a function to calculate the static pressure -- based on isentropic flow relations
def static_pressure(x_pos):
M = mach_number(x_pos)
p_0 = stag_pressure(x_pos)
static_pressure = p_0 * (1 + ((g - 1)/2)*M**2)**(-g/(g-1))
return static_pressureTo calculate the static temperature[6] isentropic flow relations can also be used, keeping in mind that the stagnation temperature remains constant over a shock:
\[ \large \require{color}{\color[rgb]{0.121820,0.954406,0.966585}T} = {\color[rgb]{0.121820,0.954406,0.966585}T}_{0} \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-1} \]
Where \({\color[rgb]{0.121820,0.954406,0.966585}T}_{0}\) is the value of the inlet stagnation temperature
# Define a function to calculate the static temperature -- based on isentropic flow relations
def static_temperature(x_pos):
M = mach_number(x_pos)
static_temperature = T_0i * (1 + ((g - 1)/2)*M**2)**(-1)
return static_temperatureNow that the Mach number and static temperature can be calculated at every point, the speed of sound and thus the flow speed at each point can be calculated. To calculate the speed of sound [7] we can use:
\[ \large \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a} = \sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R} {\color[rgb]{0.121820,0.954406,0.966585}T}} \]
To calculate the speed we can use the relation between Mach number and speed of sound, and solve for speed:
\[ \large \require{color}{\color[rgb]{0.041893,0.355669,0.727621}M} = \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.989013,0.435749,0.811750}a}} \]
# Define a function to calculate the speed of the particle at each position
def speed(x_pos):
M = mach_number(x_pos)
T = static_temperature(x_pos)
a = np.sqrt(g*R*T)
v = M * a
return v# now using the five defined functions from above the values of the five fields of interest can be found at x values
M_values = []
p0_values = []
p_values = []
T_values = []
v_values = []
for x_val in x:
M_values.append(mach_number(x_val))
p0_values.append(stag_pressure(x_val))
p_values.append(static_pressure(x_val))
T_values.append(static_temperature(x_val))
v_values.append(speed(x_val))# using the dataframe we can also create a plot with a dropdown menu between the fields
fig = go.Figure()
fig.add_scatter(x=x, y=M_values, line_color = 'plum')
fig.add_scatter(x=x, y=p0_values, line_color = 'mediumaquamarine', visible=False)
fig.add_scatter(x=x, y=p_values, line_color = 'cornflowerblue', visible=False)
fig.add_scatter(x=x, y=T_values, line_color = 'tomato', visible=False)
fig.add_scatter(x=x, y=v_values, line_color = 'darkmagenta', visible=False)
fig.update_layout(updatemenus = [dict(active=0, buttons = [dict(label='Mach Number', method='update', args=[{'visible':
[True, False, False, False, False]}, {'yaxis':{'title':'Mach Number'}}]), dict(label='Stagnation Pressure',
method='update', args=[{'visible': [False, True, False, False, False]}, {'yaxis':{'title':
'Stagnation Pressure (kPa)'}}]), dict(label='Static Pressure', method='update', args=[{'visible': [False,
False, True, False, False]}, {'yaxis':{'title':'Static Pressure (kPa)'}}]), dict(label='Static Temperature'
, method='update', args=[{'visible': [False, False, False, True, False]}, {'yaxis':{'title':
'Static Temperature (K)'}}]), dict(label='Velocity', method='update', args=[{'visible':[False, False,
False, False, True]}, {'yaxis':{'title':'Velocity (m/s)'}}])])])
fig.update_yaxes(title_text='Mach Number')
fig.update_xaxes(title_text="x position")
fig.show()Next we can make contour plots of the fields within the nozzle. As there is no \(y\) variation these plots will be uniform in the \(y\)-direction, however, when we get to topics such as boundary layers this will change and will no longer be uniform. In reality, isentropic flow is an idealized process as is a frictionless duct. Even in flows with a very low viscosity the viscous forces will never be negligible near a wall and will cause a boundary layer to form thus the field contours would differ along \(y\) and can lead to the potential for flow separation. Friction also causes a loss of heat meaning the flow would no longer be adiabatic and causes a drop in enthalpy, therefore reducing efficency as well as the exit Mach number.
# create contour maps of the nozzle
fig = go.Figure()
fig.add_contour(z=[M_values, M_values], x=x, contours_coloring = 'heatmap', line_width = 0, colorscale = 'deep',
reversescale = True, colorbar=dict(title='Mach Number', titleside='right'))
fig.add_contour(z=[p0_values, p0_values], x=x, contours_coloring ='heatmap', line_width = 0, colorscale = 'Teal',
visible = False, colorbar=dict(title='Stagnation Pressure (kPa)', titleside='right'))
fig.add_contour(z=[p_values, p_values], x=x, contours_coloring = 'heatmap', line_width = 0, colorscale = 'viridis',
visible = False, colorbar=dict(title='Static Pressure (kPa)', titleside='right'))
fig.add_contour(z=[T_values, T_values], x=x, contours_coloring = 'heatmap', line_width = 0, colorscale = 'blackbody',
reversescale = True, visible = False, colorbar=dict(title='Static Temperature (K)', titleside='right'))
fig.add_scatter(x=x, y=alpha+0.15, line_color = 'white')
fig.add_scatter(x=x, y=x-x+0.3, fill='tonexty', fillcolor='white', line_color='white')
fig.add_scatter(x=x, y=beta+0.15, fill='tozeroy', fillcolor='white', line_color='white')
fig.update_layout(showlegend=False, updatemenus=[dict(active=0, buttons = [dict(label='Mach Number', method='update', args=
[{'visible':[True, False, False, False, True, True, True]}]), dict(label='Stagnation Pressure',
method='update', args=[{'visible':[False, True, False, False, True, True,True]}]), dict(label=
'Static Pressure', method='update', args=[{'visible':[False, False, True, False, True, True, True]}]),
dict(label='Static Temperature', method='update', args=[{'visible':[False, False, False, True, True, True,
True]}])], y=1.2, x=1.1)])
fig.update_yaxes(range=[0, 0.3], showticklabels = False)
fig.update_xaxes(title='x position')
fig.show()From the velocity function we can visualize how fluid would move through the nozzle. By utilizing the equation \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v} = \frac{dx}{dt}\) we can find the approximate amount of time between each equally spaced \(x\) value within the defined array. From time we can find \(x(t)\) and linerally interpolate for equally spaced values of \(t\) to find approximately how far a particle would move within that time interval.
# find the interval of time between each equally spaced x position
dx = np.diff(x)
vav = []
for i in range(0, len(dx)):
av = sum(v_values[i: i+1]) / 2
vav.append(av)
dt = dx / vav
# find the time at each equally spaced x position
time = []
for j in range(0, len(x)):
time.append(sum(dt[0:j]))
# interpolate using the x and time values then use the interpolation function to find x positions at equal intervals of time
f = interpolate.interp1d(x=time, y=x)
tnew = np.linspace(0, max(time), 500)
xnew = f(tnew)# using the interpolated data we can animate how a particle would move through the nozzle
fig = go.Figure()
fig.add_scatter(x=x, y=np.zeros(500))
fig.add_scatter(x=x, y=np.zeros(500), line_dash = 'dot', line_color = 'lightsteelblue')
fig.add_scatter(x=x, y=alpha, line_color = 'black')
fig.add_scatter(x=x, y=beta, line_color = 'black')
fig.add_vline(x = shock_x, line_dash = 'dash', line_color = 'red', annotation_text = 'shock')
fig.add_vline(x = throat_x, line_dash = 'dash', line_color = 'blue', annotation_text = 'throat')
fig.update_layout(xaxis=dict(range=[0,1], autorange=False), yaxis=dict(range=[-0.15, 0.15], autorange=False), updatemenus=
[dict(type='buttons', buttons=[dict(label='play', method='animate', args=[None, {'frame':{'duration':10}}])])])
fig.update(frames = [go.Frame(data=[go.Scatter(x=[xnew[k]], y=[np.zeros(500)[k]], mode='markers', marker=dict(color='teal',
size=10))]) for k in range(500)], layout_showlegend=False)
fig.show()# define x and v values for 18 different positions to create arrows
x_18 = np.linspace(0,1,18)
v_18 = []
for x_val in x_18:
v_18.append(speed(x_val) / 1000)
alpha_2 = (x_18/3 -0.2)**2 + 0.025
beta_2 = -alpha_2
slope_2 = 2*(x_18/3 - 0.2)/3 * abs(np.append(np.sin(np.arctan(np.diff(alpha_2) / np.diff(x_18))), 0))Now let’s plot the velocity vector field.
# create velocity vector field
fig = ff.create_quiver(x=np.linspace(0,1,18), y=np.zeros(18), u=v_18, v=np.zeros(18), line_color = 'blue')
fig2 = ff.create_quiver(x=x_18, y=alpha_2, u=v_18, v=slope_2, line_color = 'blue')
fig3 = ff.create_quiver(x=x_18, y=beta_2, u=v_18, v=-slope_2, line_color = 'blue')
fig.add_scatter(x=x, y=alpha, line_color = 'black')
fig.add_scatter(x=x, y=beta, line_color='black')
fig.update_layout(xaxis=dict(range=[0,1]), yaxis=dict(range=[-0.15, 0.15]), showlegend = False)
fig.add_traces(data = fig2.data)
fig.add_traces(data = fig3.data)
fig.show()Using conservation of mass between the flow at a point to the throat (designated by *) we have that:
\[ \large \require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} A{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.814433,0.253157,0.091125}\rho}^* A^* {\color[rgb]{0.059472,0.501943,0.998465}v}^* \longrightarrow \frac{A}{A^*} = \frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}^* {\color[rgb]{0.059472,0.501943,0.998465}v}^*}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}} \]
Next we multiply and divide the right side by the initial stagnation denisty \({\color[rgb]{0.814433,0.253157,0.091125}\rho}_0\) to get
\[ \large \require{color}\frac{A}{A^*}=\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}^*}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^*}{{\color[rgb]{0.059472,0.501943,0.998465}v}} \]
Recognizing that at the throat \({\color[rgb]{0.041893,0.355669,0.727621}M}=1\) so \({\color[rgb]{0.989013,0.435749,0.811750}a}^*={\color[rgb]{0.059472,0.501943,0.998465}v}^*\) and multiplying and dividing the right side by the speed of sound a at the initial point we have:
\[ \large \require{color}\frac{A}{A^*}=\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}^*}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} \]
Then using the isentropic flow relation between initial stagnation density and denisty [8] \(\require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}} = \left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{\frac{1}{\gamma -1}}\) using this for both the initial condition and at the throat where \(\require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}^*}=\left(1+\frac{\gamma-1}{2}\right)^{\frac{1}{\gamma -1}}\) we have:
\[ \large \require{color}\frac{A}{A^*}=\left(1+\frac{\gamma-1}{2}\right)^{-\frac{1}{\gamma -1}}\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{\frac{1}{\gamma -1}}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}}{{\color[rgb]{0.059472,0.501943,0.998465}v}} \]
Then recognizing that \(\require{color}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}}{{\color[rgb]{0.059472,0.501943,0.998465}v}}=\frac{1}{{\color[rgb]{0.041893,0.355669,0.727621}M}}\) we have:
\[ \large \require{color}\frac{A}{A^*}=\left(1+\frac{\gamma-1}{2}\right)^{-\frac{1}{\gamma -1}}\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{\frac{1}{\gamma -1}}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}\frac{1}{{\color[rgb]{0.041893,0.355669,0.727621}M}} \] \[ \large = \left(\frac{\gamma+1}{2}\right)^{-\frac{1}{\gamma -1}}\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{\frac{1}{\gamma -1}}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}\frac{1}{{\color[rgb]{0.041893,0.355669,0.727621}M}} \]
Next using that \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) [7] we can get that \(\require{color}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}=\frac{\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}^*}}{\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}}=\sqrt{\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}^*}{{\color[rgb]{0.121820,0.954406,0.966585}T}}}\) then multiplying and dividing the right side by \(\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}\) we get that \(\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}=\sqrt{\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_0{\color[rgb]{0.121820,0.954406,0.966585}T}^*}{{\color[rgb]{0.121820,0.954406,0.966585}T}_0{\color[rgb]{0.121820,0.954406,0.966585}T}}}\) and by using the isentropic flow relation between stagnation and static temperature[6] \(\require{color}\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}{{\color[rgb]{0.121820,0.954406,0.966585}T}}=1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\) which at the throat would be \(\require{color}\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}{{\color[rgb]{0.121820,0.954406,0.966585}T}^*}=1+\frac{\gamma -1}{2}\) we get that:
\[ \large \require{color}\frac{{\color[rgb]{0.989013,0.435749,0.811750}a}^*}{{\color[rgb]{0.989013,0.435749,0.811750}a}}=\sqrt{\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)\left(1+\frac{\gamma -1}{2}\right)^{-1}} \]
Plugging this into the area ratio equation we get:
\[ \large \require{color}\frac{A}{A^*}=\frac{1}{{\color[rgb]{0.041893,0.355669,0.727621}M}}\left(\frac{\gamma+1}{2}\right)^{-\frac{1}{\gamma -1}}\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{\frac{1}{\gamma -1}}\sqrt{\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)\left(1+\frac{\gamma -1}{2}\right)^{-1}} \]
Simplifying this we get the area ratio mach number relationship:
\[ \large \require{color}\frac{A}{A^{*}} = \frac{\left(1 + \frac{\gamma+1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{\frac{\gamma+1}{2(\gamma-1)}}}{{\color[rgb]{0.041893,0.355669,0.727621}M}} \left(\frac{\gamma +1}{2}\right)^{-\frac{\gamma +1}{2(\gamma -1)}} \]
First we compare the speeds before and after the shock using the relationships \(\require{color}{\color[rgb]{0.041893,0.355669,0.727621}M}=\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.989013,0.435749,0.811750}a}}\) and \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) [7] :
\[ \large \require{color}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}{{\color[rgb]{0.059472,0.501943,0.998465}v}_2}=\frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_1{\color[rgb]{0.989013,0.435749,0.811750}a}_1}{{\color[rgb]{0.041893,0.355669,0.727621}M}_2{\color[rgb]{0.989013,0.435749,0.811750}a}_2}=\frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_1\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_1}}{{\color[rgb]{0.041893,0.355669,0.727621}M}_2\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_2}}=\frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_1}{{\color[rgb]{0.041893,0.355669,0.727621}M}_2}\sqrt{\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_1}{{\color[rgb]{0.121820,0.954406,0.966585}T}_2}} \]
From conservation of mass \(\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1 A_1 {\color[rgb]{0.059472,0.501943,0.998465}v}_1 = {\color[rgb]{0.814433,0.253157,0.091125}\rho}_2 A_2 {\color[rgb]{0.059472,0.501943,0.998465}v}_2\) and since the change in Mach number at a shock happens instantaneously \(A_1=A_2\) therefore \(\require{color}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}{{\color[rgb]{0.059472,0.501943,0.998465}v}_2}=\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\)
Using the ideal gas law \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}\) we have that \(\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) therefore:
\[ \large \require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2 {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_1}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1 {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_2} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2 {\color[rgb]{0.121820,0.954406,0.966585}T}_1}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1 {\color[rgb]{0.121820,0.954406,0.966585}T}_2} \]
We can find the pressure ratio by using Euler’s equation \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_1 +{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1 {\color[rgb]{0.059472,0.501943,0.998465}v}_1^2 = {\color[rgb]{0.315209,0.728565,0.037706}p}_2 +{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2 {\color[rgb]{0.059472,0.501943,0.998465}v}_2^2\) and substituting \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}\) and using that \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}=\sqrt{\gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}}\) [7]:
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_1 + {\color[rgb]{0.814433,0.253157,0.091125}\rho}_1 \left({\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 \gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right) = {\color[rgb]{0.315209,0.728565,0.037706}p}_2 + {\color[rgb]{0.814433,0.253157,0.091125}\rho}_2 \left({\color[rgb]{0.041893,0.355669,0.727621}M}_2^2 \gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}\right) \] \[ \large \Rightarrow {\color[rgb]{0.315209,0.728565,0.037706}p}_1(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2) = {\color[rgb]{0.315209,0.728565,0.037706}p}_2(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2) \] \[ \large \Rightarrow \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1} = \frac{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2)}{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2)} \]
We can find the temperature ratio using isentropic flow relations [6] even though the flow is not isentropic over the shock it is both directly before and after so we can find the temperatures directly before and after and compare them, remembering that across a shock the stagnation pressure remains constant
\[ \large \require{color}{\color[rgb]{0.121820,0.954406,0.966585}T}_1 = {\color[rgb]{0.121820,0.954406,0.966585}T}_0 \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2\right)^{-1}, \] \[ \large {\color[rgb]{0.121820,0.954406,0.966585}T}_2 = {\color[rgb]{0.121820,0.954406,0.966585}T}_0 \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\right)^{-1} \] \[ \large \Rightarrow \frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_1}{{\color[rgb]{0.121820,0.954406,0.966585}T}_2} = \frac{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2\right)^{-1}}{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\right)^{-1}} \]
Then we can solve for \(\require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\) by making these subsitutions:
\[ \large \require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1} = \frac{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2)}{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2)}\frac{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\right)}{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2\right)} \]
We can also use the temperature ratio to solve for \(\require{color}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}{{\color[rgb]{0.059472,0.501943,0.998465}v}_2}\):
\[ \large \require{color}\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}_1}{{\color[rgb]{0.059472,0.501943,0.998465}v}_2} = \frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_1}{{\color[rgb]{0.041893,0.355669,0.727621}M}_2}\sqrt{\frac{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\right)}{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2\right)}} \]
Setting these two equivalent we have:
\[ \large \require{color}\frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_1}{{\color[rgb]{0.041893,0.355669,0.727621}M}_2}\sqrt{\frac{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\right)}{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2\right)}} = \frac{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2)}{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2)}\frac{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\right)}{\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2\right)} \]
Moving terms of \({\color[rgb]{0.041893,0.355669,0.727621}M}_1\) to one side and \({\color[rgb]{0.041893,0.355669,0.727621}M}_2\) to the other we have:
\[ \large \require{color}\frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_1}{1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2}\sqrt{1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^2} = \frac{{\color[rgb]{0.041893,0.355669,0.727621}M}_2}{1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2}\sqrt{1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2} \]
Next with some algebra we can solve for \({\color[rgb]{0.041893,0.355669,0.727621}M}_2^2\)
\[ \large \require{color} {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2 = \frac{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 + 2}{2\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 - (\gamma -1)} \]
We begin again with continuity of mass \(\require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}}={\color[rgb]{0.814433,0.253157,0.091125}\rho} A{\color[rgb]{0.059472,0.501943,0.998465}v}\) and use the relationship between Mach number and velocity \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}\) to get that \(\require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}}={\color[rgb]{0.814433,0.253157,0.091125}\rho} A{\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}\) then we use the forumla for speed of sound \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) [7] to get that:
\[ \large \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}}={\color[rgb]{0.814433,0.253157,0.091125}\rho} A{\color[rgb]{0.041893,0.355669,0.727621}M}\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}} \]
Then we use the ideal gas law \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} \longrightarrow {\color[rgb]{0.814433,0.253157,0.091125}\rho} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) to substitute for \({\color[rgb]{0.814433,0.253157,0.091125}\rho}\):
\[ \large \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}}=\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}} A{\color[rgb]{0.041893,0.355669,0.727621}M}\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}= A{\color[rgb]{0.041893,0.355669,0.727621}M}\sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}}} \]
Then by using the relation that \(\require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_0}=\left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}\right)^{\frac{\gamma}{\gamma -1}}\) which is derived by dividing the two isentropic flow relations for \(\require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_0}\) and \(\require{color}\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}\) we can substitute for \({\color[rgb]{0.315209,0.728565,0.037706}p}\):
\[ \large \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}}= A{\color[rgb]{0.041893,0.355669,0.727621}M}\sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_0\left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}\right)^{\frac{\gamma}{\gamma -1}}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}}}=A{\color[rgb]{0.041893,0.355669,0.727621}M}\sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_0\left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}\right)^{\frac{\gamma+1}{2(\gamma -1)}}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}} \]
Then using the isentropic flow relation \(\require{color}\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}{{\color[rgb]{0.121820,0.954406,0.966585}T}}=1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\) [6] we get that:
\[ \large \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}}=A{\color[rgb]{0.041893,0.355669,0.727621}M}\sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_0\left(1+\frac{\gamma -1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{-\frac{\gamma+1}{2(\gamma -1)}}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}} \]
Rearraging this we get the mass flow equation:
\[ \large \require{color}\dot{{\color[rgb]{0.501961,0.250953,0.010028}m}} = \frac{A{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}}{\sqrt{{\color[rgb]{0.121820,0.954406,0.966585}T}_{0}}} \sqrt{\frac{\gamma}{{\color[rgb]{0.986252,0.007236,0.027423}R}}} {\color[rgb]{0.041893,0.355669,0.727621}M} \left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-\frac{\gamma+1}{2(\gamma -1)}} \]
First we find the pressure ratio between the static pressure before and after the shock by using Euler’s equation \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_1 +{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1 {\color[rgb]{0.059472,0.501943,0.998465}v}_1^2 = {\color[rgb]{0.315209,0.728565,0.037706}p}_2 +{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2 {\color[rgb]{0.059472,0.501943,0.998465}v}_2^2\) and substituting \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}\) and using that \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}=\sqrt{\gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}}\) [7]:
\[ \large \require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}_1 + {\color[rgb]{0.814433,0.253157,0.091125}\rho}_1 \left({\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 \gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right) = {\color[rgb]{0.315209,0.728565,0.037706}p}_2 + {\color[rgb]{0.814433,0.253157,0.091125}\rho}_2 \left({\color[rgb]{0.041893,0.355669,0.727621}M}_2^2 \gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}\right) \] \[ \large \Rightarrow {\color[rgb]{0.315209,0.728565,0.037706}p}_1(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2) = {\color[rgb]{0.315209,0.728565,0.037706}p}_2(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2) \] \[ \large \Rightarrow \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1} = \frac{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2)}{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2)} \]
Next to find the pressure ratio \(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}}\) we can multiply the static pressure ratio by \(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}\) and \(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}}\) to get that \(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\), then we can use isentropic flow relations between stagnation and static pressure [5] as both directly before and directly after the shock the flow is isentropic so the stagnation and static pressure before can be related and the stagnation and static pressure after can also be related:
\[ \large \require{color}\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}} = \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^{2}\right)^{-\frac{\gamma}{\gamma-1}} \] \[ \large \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}} = \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^{2}\right)^{-\frac{\gamma}{\gamma-1}} \]
Plugging this into the relation of stagnation pressures and using the previously derived ratio between static pressures before and after we have:
\[ \large \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}}= \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^{2}\right)^{-\frac{\gamma}{\gamma-1}}\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^{2}\right)^{\frac{\gamma}{\gamma-1}}\frac{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2)}{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_2^2)} \]
Next using the relation between \({\color[rgb]{0.041893,0.355669,0.727621}M}_1\) and \({\color[rgb]{0.041893,0.355669,0.727621}M}_2 \longrightarrow \require{color}{\color[rgb]{0.041893,0.355669,0.727621}M}_2^2 = \frac{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 + 2}{2\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 - (\gamma -1)}\) [2] we can find the stagnation pressure ratio in terms of only \({\color[rgb]{0.041893,0.355669,0.727621}M}_1\) the Mach number before the shock:
\[ \large \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}}= \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}_1^{2}\right)^{-\frac{\gamma}{\gamma-1}}\left(1+ \frac{\gamma-1}{2}\left(\frac{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 + 2}{2\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 - (\gamma -1)}\right)\right)^{\frac{\gamma}{\gamma-1}} \times \] \[ \large \frac{(1+\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2)}{(1+\gamma \left(\frac{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 + 2}{2\gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_1^2 - (\gamma -1)}\right))} \]
By doing lots of algebra and simplification we arrive at:
\[ \large \require{color} \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_{02}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{01}} = \left(\frac{(\gamma+1){\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}}{(\gamma-1){\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2}-2}\right)^{\frac{\gamma}{\gamma-1}} \left(\frac{\gamma+1}{2 \gamma {\color[rgb]{0.041893,0.355669,0.727621}M}_{1}^{2} - (\gamma -1)}\right)^{\frac{1}{\gamma-1}} \]
By starting with the equations \(\require{color}d{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501945,0.999984,0.031006}V}d{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.501958,0.501942,0.014744}c_p}d{\color[rgb]{0.121820,0.954406,0.966585}T} \longrightarrow d{\color[rgb]{0.121820,0.954406,0.966585}T}=\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}d{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}\) and \(\require{color}de=-{\color[rgb]{0.315209,0.728565,0.037706}p}d{\color[rgb]{0.501945,0.999984,0.031006}V}=c_vd{\color[rgb]{0.121820,0.954406,0.966585}T}\) and substituting for \(d{\color[rgb]{0.121820,0.954406,0.966585}T}\) as well as using the relation \(\require{color}\gamma = \frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{c_v}\) we get:
\[ \large \require{color}-{\color[rgb]{0.315209,0.728565,0.037706}p}d{\color[rgb]{0.501945,0.999984,0.031006}V} = \frac{c_{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{\color[rgb]{0.501945,0.999984,0.031006}V}d{\color[rgb]{0.315209,0.728565,0.037706}p} \longrightarrow -\gamma \frac{1}{{\color[rgb]{0.501945,0.999984,0.031006}V}} d{\color[rgb]{0.501945,0.999984,0.031006}V} = \frac{1}{{\color[rgb]{0.315209,0.728565,0.037706}p}}d{\color[rgb]{0.315209,0.728565,0.037706}p} \]
Then by integrating from initial to final conditions:
\[ \large \require{color}-\gamma \int_{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}^{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}\frac{1}{{\color[rgb]{0.501945,0.999984,0.031006}V}}d{\color[rgb]{0.501945,0.999984,0.031006}V} = \int_{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}^{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}\frac{1}{{\color[rgb]{0.315209,0.728565,0.037706}p}}d{\color[rgb]{0.315209,0.728565,0.037706}p} \longrightarrow -\gamma ln\left(\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}\right) = ln\left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \] \[ \large \Rightarrow ln\left(\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}\right)^{-\gamma} = ln\left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}\right)^{-\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \]
Then by using the fact that denisty and volue are inversely proportional we can get that:
\[ \large \require{color}\left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}\right)^{-\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \longrightarrow \left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right)^{\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \]
Next we can use the ideal gas law \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} \longrightarrow {\color[rgb]{0.814433,0.253157,0.091125}\rho} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) to substitute for \({\color[rgb]{0.814433,0.253157,0.091125}\rho}\):
\[ \large \require{color}\left(\frac{\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_2}}{\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}{{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_1}}\right)^{\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \longrightarrow \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2{\color[rgb]{0.121820,0.954406,0.966585}T}_1}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1{\color[rgb]{0.121820,0.954406,0.966585}T}_2}\right)^{\gamma} \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \longrightarrow \left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_1}{{\color[rgb]{0.121820,0.954406,0.966585}T}_2}\right) = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right)^{\frac{1}{\gamma}}\left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1} \right)^{-1} =\left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right)^{\frac{1}{1-\gamma}} \]
Solving for the pressure ratio we have:
\[ \large \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right)=\left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_2}{{\color[rgb]{0.121820,0.954406,0.966585}T}_1} \right)^{\frac{\gamma}{\gamma-1}} \]
Then we can use the isentropic flow relation between static and stagnation temperature [6] and define the initial condition as being stagnation (we can make this assumption as we can imagine a point in isentropic flow where v=0 and therefore static=stagnation):
\[ \large \left(\left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-1}\right)^{\frac{\gamma}{\gamma-1}}=\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}} \]
Simplifying this we get the isentropic flow relation between static and stagnation pressure:
\[ \large \require{color} \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.315209,0.728565,0.037706}p}_{0}} = \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-\frac{\gamma}{\gamma-1}} \]
In isentropic flow stagnation enthalpy is conserved and is defined as \(\require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.064095,0.501831,0.501977}h}+\frac{1}{2}{\color[rgb]{0.059472,0.501943,0.998465}v}^2\) thand enthalpy can be defined as \(\require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}\) plugging this in we get:
\[ \large \require{color}{\color[rgb]{0.064095,0.501831,0.501977}h}_0 = {\color[rgb]{0.501958,0.501942,0.014744}c_p}{\color[rgb]{0.121820,0.954406,0.966585}T}+ \frac{1}{{\color[rgb]{0.059472,0.501943,0.998465}v}^2} \] \[ \large \Rightarrow \frac{{\color[rgb]{0.064095,0.501831,0.501977}h}_0}{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}={\color[rgb]{0.121820,0.954406,0.966585}T} + \frac{1}{2{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{\color[rgb]{0.059472,0.501943,0.998465}v}^2 \] \[ \large \Rightarrow {\color[rgb]{0.121820,0.954406,0.966585}T}_0 = {\color[rgb]{0.121820,0.954406,0.966585}T} + \frac{1}{2{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{\color[rgb]{0.059472,0.501943,0.998465}v}^2 \]
Recognizing that \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}={\color[rgb]{0.041893,0.355669,0.727621}M}{\color[rgb]{0.989013,0.435749,0.811750}a}\) and that \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}=\sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}\) [7] we can replace \({\color[rgb]{0.059472,0.501943,0.998465}v}^2\) with \(\require{color}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}\) to get:
\[ \large \require{color}{\color[rgb]{0.121820,0.954406,0.966585}T}_0 = {\color[rgb]{0.121820,0.954406,0.966585}T} + \frac{1}{2{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} \longrightarrow \frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}{{\color[rgb]{0.121820,0.954406,0.966585}T}} = \left(1+ \frac{{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}}{2{\color[rgb]{0.501958,0.501942,0.014744}c_p}}\right) \]
Next we can substitute \(\require{color}{\color[rgb]{0.986252,0.007236,0.027423}R}={\color[rgb]{0.501958,0.501942,0.014744}c_p}-c_v\) and use the fact that \(\require{color}\gamma = \frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{c_v}\) to get
\[ \large \require{color}\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}{{\color[rgb]{0.121820,0.954406,0.966585}T}} = \left(1+ \frac{{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma ({\color[rgb]{0.501958,0.501942,0.014744}c_p} - c_v)}{2{\color[rgb]{0.501958,0.501942,0.014744}c_p}}\right) \] \[ \large = \left(1+ \frac{{\color[rgb]{0.041893,0.355669,0.727621}M}^2 \gamma}{2}\left(1 - \frac{1}{\gamma}\right)\right)= \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right) \]
Rearranging to get \(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.121820,0.954406,0.966585}T}_0}\) we get that:
\[ \large \require{color} \frac{{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.121820,0.954406,0.966585}T}_{0}} = \left(1+ \frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^{2}\right)^{-1} \]
Starting with continuity of mass of a system with decreasing speed: \(\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho} A{\color[rgb]{0.059472,0.501943,0.998465}v}=({\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho})A({\color[rgb]{0.059472,0.501943,0.998465}v}-d{\color[rgb]{0.059472,0.501943,0.998465}v})\) we can solve for the change in velocity:
\[ \large \require{color}d{\color[rgb]{0.059472,0.501943,0.998465}v} = \frac{{\color[rgb]{0.059472,0.501943,0.998465}v}d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho}} \]
Then using conservation of momentum we have: \(\require{color}Ad{\color[rgb]{0.315209,0.728565,0.037706}p}=A{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 -A({\color[rgb]{0.814433,0.253157,0.091125}\rho} +d{\color[rgb]{0.814433,0.253157,0.091125}\rho})({\color[rgb]{0.059472,0.501943,0.998465}v}-d{\color[rgb]{0.059472,0.501943,0.998465}v})^2\) then using the equation for \(d{\color[rgb]{0.059472,0.501943,0.998465}v}\) we get:
\[ \large \require{color}Ad{\color[rgb]{0.315209,0.728565,0.037706}p}=A{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 -A({\color[rgb]{0.814433,0.253157,0.091125}\rho} +d{\color[rgb]{0.814433,0.253157,0.091125}\rho})\left({\color[rgb]{0.059472,0.501943,0.998465}v}-\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\right)^2 \] \[ \large A{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 - A({\color[rgb]{0.814433,0.253157,0.091125}\rho} +d{\color[rgb]{0.814433,0.253157,0.091125}\rho})\left(\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}+d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\right)^2 = A{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2 - A\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2 {\color[rgb]{0.814433,0.253157,0.091125}\rho}^2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho}} \] \[ \large A\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.059472,0.501943,0.998465}v}^2({\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho})}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho}} -A\frac{{\color[rgb]{0.059472,0.501943,0.998465}v}^2 {\color[rgb]{0.814433,0.253157,0.091125}\rho}^2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}=A{\color[rgb]{0.059472,0.501943,0.998465}v}^2 \frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho} d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}+d{\color[rgb]{0.814433,0.253157,0.091125}\rho}} \]
Solving for \({\color[rgb]{0.059472,0.501943,0.998465}v}^2\) we get:
\[ \large \require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}^2 = d{\color[rgb]{0.315209,0.728565,0.037706}p}\left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho} + d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho} d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\right) \] \[ \large =\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\left(1+\frac{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\right) \]
The speed of sound is defined as the speed of the wave for small perturbations meaning that \(\require{color}d{\color[rgb]{0.315209,0.728565,0.037706}p} \rightarrow 0\), \(\require{color}d{\color[rgb]{0.814433,0.253157,0.091125}\rho} \rightarrow 0\) and \(\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v} \rightarrow {\color[rgb]{0.989013,0.435749,0.811750}a}\) thus we get:
\[ \large \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}^2 = \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}} \]
For an adiabatic process we have that \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.501945,0.999984,0.031006}V}^{\gamma} = const\) or equivalently \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}{\color[rgb]{0.814433,0.253157,0.091125}\rho}^{-\gamma}=(const) \longrightarrow {\color[rgb]{0.315209,0.728565,0.037706}p}=const{\color[rgb]{0.814433,0.253157,0.091125}\rho}^{-\gamma}\) if we take the derivative we get:
\[ \large \require{color}d{\color[rgb]{0.315209,0.728565,0.037706}p} = \gamma {\color[rgb]{0.814433,0.253157,0.091125}\rho}^{\gamma-1} d{\color[rgb]{0.814433,0.253157,0.091125}\rho} (const) \] \[ \large \Rightarrow \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}=\gamma \frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}^{\gamma}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}(const) \]
Then using that \(\require{color}{\color[rgb]{0.814433,0.253157,0.091125}\rho}^{\gamma} = \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{const}\) we have:
\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}=\gamma \frac{\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{const}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}(const)=\gamma \frac{{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}} \]
Next we use the ideal gas law \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}\) and substitute for \({\color[rgb]{0.315209,0.728565,0.037706}p}\):
\[ \large \require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}=\gamma \frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}} = \gamma {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T} \]
Then using the equation for \(a\) that we derived earlier \(\require{color}{\color[rgb]{0.989013,0.435749,0.811750}a}^2 = \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{d{\color[rgb]{0.814433,0.253157,0.091125}\rho}}\) we get that:
\[ \large \require{color}{\color[rgb]{0.989013,0.435749,0.811750}a} = \sqrt{\gamma {\color[rgb]{0.986252,0.007236,0.027423}R} {\color[rgb]{0.121820,0.954406,0.966585}T}} \]
By starting with the equations \(\require{color}d{\color[rgb]{0.064095,0.501831,0.501977}h}={\color[rgb]{0.501945,0.999984,0.031006}V}d{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.501958,0.501942,0.014744}c_p}d{\color[rgb]{0.121820,0.954406,0.966585}T} \longrightarrow d{\color[rgb]{0.121820,0.954406,0.966585}T}=\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}d{\color[rgb]{0.315209,0.728565,0.037706}p}}{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}\) and \(\require{color}de=-{\color[rgb]{0.315209,0.728565,0.037706}p}d{\color[rgb]{0.501945,0.999984,0.031006}V}=c_vd{\color[rgb]{0.121820,0.954406,0.966585}T}\) and substituting for \(d{\color[rgb]{0.121820,0.954406,0.966585}T}\) as well as using the relation \(\require{color}\gamma = \frac{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{c_v}\) we get:
\[ \large \require{color}-{\color[rgb]{0.315209,0.728565,0.037706}p}d{\color[rgb]{0.501945,0.999984,0.031006}V} = \frac{c_{\color[rgb]{0.059472,0.501943,0.998465}v}}{{\color[rgb]{0.501958,0.501942,0.014744}c_p}}{\color[rgb]{0.501945,0.999984,0.031006}V}d{\color[rgb]{0.315209,0.728565,0.037706}p} \]
\[ \large \Rightarrow -\gamma \frac{1}{{\color[rgb]{0.501945,0.999984,0.031006}V}} d{\color[rgb]{0.501945,0.999984,0.031006}V} = \frac{1}{{\color[rgb]{0.315209,0.728565,0.037706}p}}d{\color[rgb]{0.315209,0.728565,0.037706}p} \]
Then by integrating from initial to final conditions:
\[ \large \require{color}-\gamma \int_{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}^{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}\frac{1}{{\color[rgb]{0.501945,0.999984,0.031006}V}}d{\color[rgb]{0.501945,0.999984,0.031006}V} = \int_{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}^{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}\frac{1}{{\color[rgb]{0.315209,0.728565,0.037706}p}}d{\color[rgb]{0.315209,0.728565,0.037706}p} \]
\[ \large \Rightarrow -\gamma ln\left(\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}\right) = ln\left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \] \[ \large \Rightarrow ln\left(\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}\right)^{-\gamma} = ln\left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.501945,0.999984,0.031006}V}_2}{{\color[rgb]{0.501945,0.999984,0.031006}V}_1}\right)^{-\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \]
Then by using the fact that denisty and volue are inversely proportional we can get that:
\[ \large \require{color}\left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}\right)^{-\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right)^{\gamma} = \left(\frac{{\color[rgb]{0.315209,0.728565,0.037706}p}_2}{{\color[rgb]{0.315209,0.728565,0.037706}p}_1}\right) \]
Next we can use the ideal gas law \(\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}={\color[rgb]{0.814433,0.253157,0.091125}\rho} {\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}\):
\[ \large \require{color}\left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right)^{\gamma} = \left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1{\color[rgb]{0.986252,0.007236,0.027423}R}{\color[rgb]{0.121820,0.954406,0.966585}T}_1}\right) \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right)^{\gamma}\left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right)^{-1} = \left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_2}{{\color[rgb]{0.121820,0.954406,0.966585}T}_1}\right) \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right)^{\gamma-1} = \left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_2}{{\color[rgb]{0.121820,0.954406,0.966585}T}_1}\right) \] \[ \large \Rightarrow \left(\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_2}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_1}\right) = \left(\frac{{\color[rgb]{0.121820,0.954406,0.966585}T}_2}{{\color[rgb]{0.121820,0.954406,0.966585}T}_1}\right)^{\frac{1}{\gamma-1}} \]
Then we can use the isentropic flow relation between static and stagnation temperature [6] and define the initial condition as being stagnation (we can make this assumption as we can imagine a point in isentropic flow where v=0 and therefore static=stagnation):
\[ \large \require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0} = \left(\left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{-1}\right)^{\frac{1}{\gamma -1}} \]
Simplifying this we get the relationship between static and stagnation density
\[ \large \require{color}\frac{{\color[rgb]{0.814433,0.253157,0.091125}\rho}}{{\color[rgb]{0.814433,0.253157,0.091125}\rho}_0} = \left(1+\frac{\gamma-1}{2}{\color[rgb]{0.041893,0.355669,0.727621}M}^2\right)^{-\frac{1}{\gamma -1}}. \]
Sources: Mass Flow Rate, Isentropic Flow Relations, Normal Shock Relations, AE 3450 Slides, Class Notes