Rankine-Hugoniot, Shock Waves and Wall Reflexion¶
The current notebook aims at computing attached shock wave and reflexions on walls.
import numpy as np
import matplotlib.pyplot as plt
from aerokit.common import defaultgas
from aerokit.aero import degree as deg # import trigo functions with degree unit support
from aerokit.aero import ShockWave as sw # import functions for shockwave computation
%matplotlib inline
First, fluid definition ($\gamma$) and upstream conditions are defined. A function is also defined to plot the geometry.
# definition of problem parameters
gam = 1.4 ; defaultgas.set_gamma(gam)
M0 = 2.8
wdev = 18.
print("wall deviation (deg): {} with upwstream Mach number {}".format(wdev, M0))
# function to plot the geometry
#
def plot_geom(xneg=-.5, length=2., dev=wdev, zoom=1):
fig = plt.figure(figsize=(14*zoom,8*zoom))
ax = fig.add_subplot(111)
#plt.axis([xneg, length])
ax.set(aspect="equal", xlim=[xneg, length], ylim=[-.1, 1.1])
plt.plot([xneg, length], [1, 1], color="black", linewidth=2)
plt.plot([xneg, 0, length], [0, 0, length*deg.tan(dev)], color="black", linewidth=2)
#
# test de la fonction de tracé
plot_geom(dev=wdev, zoom=.5)
wall deviation (deg): 18.0 with upwstream Mach number 2.8
First attached shock wave¶
Using aerokit.aero.ShockWave
module, one can check the maximum allowed deviation for this upstream Mach number (either sonic or max deviation limit). Then, the shock angle $\sigma$ is computed, followed by the normal Mach number and all ratios of quantities.
devmax = sw.dev_Max(M0)
devsonic = sw.dev_Sonic(M0)
print(("For upstream Mach number M0= {:1.4},\n* maximum deviation is {:1.4}°\n"+
"* limit for downstream supersonic flow is {:1.4}°").format(M0, devmax, devsonic))
For upstream Mach number M0= 2.8, * maximum deviation is 32.59° * limit for downstream supersonic flow is 32.5°
sig1 = sw.sigma_Mach_deflection(M0, wdev)
Mn0 = M0*deg.sin(sig1)
p1p0 = sw.Ps_ratio(Mn0)
Mn1 = sw.downstream_Mn(Mn0)
M1 = Mn1/deg.sin(sig1-wdev)
print("shock with {:1.4}° deviation and angle {:1.4}\ndownstream Mach number is M1= {:1.4}\nCompression ratio is p1/p0= {:1.4}".format(wdev, sig1, M1, p1p0))
shock with 18.0° deviation and angle 37.14 downstream Mach number is M1= 1.961 Compression ratio is p1/p0= 3.168
Shock reflexion¶
Again, existence of reflected shock is checked before computing $\sigma_2$. Then, the reflected shock wave is computed and downstream state 3.
devmax = sw.dev_Max(M1)
if wdev > devmax:
print("no attached solution, the reflexion will be irregular")
# if irregular reflexion, assume that reflected shock wave is at downstream sonic limit
sig2 = sw.sigma_Sonic(M1)
else:
sig2 = sw.sigma_Mach_deflection(M1, wdev)
Mn1 = M1*deg.sin(sig2)
p2p1 = sw.Ps_ratio(Mn1)
Mn2 = sw.downstream_Mn(Mn1)
M2 = Mn2/deg.sin(sig2-wdev)
print("shock reflexion with {:1.4}° deviation and angle {:1.4}° (max deviation is {:1.4}°)".format(wdev, sig2, devmax))
print("downstream Mach number is M2= {:1.4}\nCompression ratio is p2/p1= {:1.4}".format(M2, p2p1))
p2p0 = p2p1*p1p0
shock reflexion with 18.0° deviation and angle 50.96° (max deviation is 22.29°) downstream Mach number is M2= 1.274 Compression ratio is p2/p1= 2.54
All geometric angles are available to draw the shocks solution (the reflexion is assumed regular).
yend=.75 # parameter to draw reflected shock
plot_geom(dev=wdev, zoom=.8)
xup = 1./deg.tan(sig1) # impact of first shock on top wall
xbot = xup + (1.-yend)/deg.tan(sig2-wdev) # abscissa of reflected shock at yend
plt.plot([0, xup, xbot],
[0, 1, yend], 'red', linewidth=2)
[<matplotlib.lines.Line2D at 0x7fa8bde21330>]
Representation in the shock polar $\theta/p$¶
The black line is the polar curve which depends on $M_0$. The solution of the shock wave is at $\theta=\mbox{wdev}$. From this point, the new polar curve for $M_1$ (in red) is drawn ; the deviation is negative up to $\theta=0$. If this intersection exists, the shock reflexion is regular and the solution is state 0. The last curve in blue is the polar curve for $M_2$ which may not reach the wall deviation.
import aerokit.aero.plot.shockpolar as swplt
fig=swplt.figure_theta_pressure(figsize=(14,8))
fig.suptitle('Polar of Shock-Waves, $\gamma = %.1f$'%gam, fontsize=18, y=0.93)
plt.xlabel('flow angle', fontsize=14)
plt.ylabel('normalized static pressure', fontsize=14)
if p2p0 < 20.: plt.yscale('linear') # default is logarithmic
#
# plot polar curves
swplt.plot_theta_pressure(M0, devmax=True, sonic=True)
swplt.plot_theta_pressure(M1, thet_init=wdev, p_init=p1p0, color='red')
swplt.plot_theta_pressure(M2, thet_init=0., p_init=p2p0, color='blue')
# plot symbols for flow regions
plt.plot(0, 1., 'bo')
plt.plot(wdev, p1p0, 'wo')
plt.plot(0., p2p0, 'go')
[<matplotlib.lines.Line2D at 0x7fa8bdebebf0>]