import aerokit.aero.unsteady1D as uq
import aerokit.instance.riemann as riem
Shock tube and interface interaction¶
This test case is a classical application of a shock tube facility to study shock/interface interaction. Most interactions between locally uniform states are theoretically available to computation: waves speeds and intermediate states are computed. The $x/t$ diagram is then provided at the end. The limitation of these computation is given by expansion/discontinuities interaction where non-uniform flows cannot be simply described.
Shock generation¶
The first Riemann problem is given the initial interaction between 2 zones with a pressure ratio. The driver gas q20
is air. The driven gas q10
(which will propagate the shock wave) is air too. Gas states will be denoted qxn
where x
is the initial zone number and 'n' the successive states in time.
Pinit = 1.e5 # reference pressure
Pratio = 2.38 # pressure ratio (left over right)
Tisoth = 293.15 # same temperature for both gas
r_air = 287.1
r_He = 2078.
gam_air = 1.4
gam_He = 5./3.
#
# driver gas (q20)
q20 = uq.unsteady_state(rho = Pinit*Pratio/(r_air*Tisoth),
u = 0.,
p = Pinit*Pratio, gamma=gam_air)
# driven gas (q10)
q10 = uq.unsteady_state(rho = Pinit/(r_air*Tisoth),
u = 0.,
p = Pinit, gamma=gam_air)
# Riemann problem
pb1 = riem.riemann_pb(q20, q10)
q21 = pb1.qstarL()
q11 = pb1.qstarR()
print("shock compression ratio: {:.3f}".format(q11.p/q10.p))
print("shock Mach number : {:.3f}".format(pb1.right_fastest()/q10.asound()))
shock compression ratio: 1.521 shock Mach number : 1.203
Shock / Interface interaction¶
The previous computed shock wave will reach an interface between q10
(driven gas) and q10
. The Riemann problem when the shock interacts is defined by q11
(shock downstream flow) and q10
.
# Helium
q00 = uq.unsteady_state(rho = Pinit/(r_He*Tisoth),
u = 0.,
p = Pinit, gamma=gam_He)
# Riemann problem
pb2 = riem.riemann_pb(q11, q00)
q12 = pb2.qstarL()
q01 = pb2.qstarR()
print("transmitted shock compression ratio: {:.3f}".format(q01.p/q00.p))
print("interface speed : {:.3f}".format(q12.u))
transmitted shock compression ratio: 1.277 interface speed : 151.399
Shock reflexion¶
The transmitted shock wave will bounce or reflect on the wall. Even if the reflected shock wave could have been computed directly, the trick here is to defined a symmetric state and a new Riemann problem. Helium is stopped and compressed to q02
state.
# define symmetric (ghost state)
q01gh = q01.copysymmetric()
# Wall interaction
pbw = riem.riemann_pb(q01, q01gh)
q02 = pbw.qstarL()
print("shock compression ratio: {:.3f}".format(q02.p/q01.p))
print("shock Mach number : {:.3f}".format(pbw.left_fastest()/q01.asound()))
shock compression ratio: 1.262 shock Mach number : -0.957
Reshock of interface¶
The reflected shock will reach the interface with a new complex Riemann problem: the transmitted shock wave will change the air from q12
to q13
. A reflected right wave will be produced: it can be an expansion or a shock depending on the initial condition and correlated to the new interface velocity.
# Riemann problem
pb3 = riem.riemann_pb(q12, q02)
q13 = pb3.qstarL()
q03 = pb3.qstarR()
#print("interface speed : {:.3f}".format(q13.u))
Position of interactions¶
Given all waves velocities, all interaction positions can be computed and drawn in $x/t$ diagram.
# positions
x_int = 2.5
x_Rwall = 3.
x_Lwall = -1.
#
def xt_intersect(x1, t1, u1, x2, t2, u2):
ti = (-x2+x1 + t2*u2-t1*u1)/(u2-u1)
return x1+u1*(ti-t1), ti
#
# shock / interface
t_shock_int = x_int / pb1.right_fastest()
#print("shock / interface interaction at (x,t)=(%.3f,%.5f)"%(x_int, t_shock_int))
# shock / wall
t_shock_wall = (x_Rwall-x_int) / pb2.right_fastest() + t_shock_int
#print("shock / wall interaction at (x,t)=(%.3f,%.5f)"%(x_Rwall, t_shock_wall))
# interface / reshock
x_reshock, t_reshock = xt_intersect(x_int, t_shock_int, q12.u,
x_Rwall, t_shock_wall, pbw.left_fastest())
#print("interface / reshock interaction at (x,t)=(%.3f,%.5f)"%(x_reshock, t_reshock))
t_max = 1.2*t_reshock
x/t diagram¶
import matplotlib.pyplot as plt
%matplotlib inline
#
plt.figure(figsize=(15,8))
plt.axis([-1, 3, 0, t_max])
plt.xlabel(r'$x$', fontsize=18)
plt.ylabel(r"$t$", fontsize=18)
plt.title(u"shocktube in x/t diagram", fontsize=18)
#
# plot shocks
plt.plot([0, x_int, x_Rwall, x_reshock, x_reshock+pb3.left_fastest()*(t_max-t_reshock)],
[0, t_shock_int, t_shock_wall, t_reshock, t_max],
color="red", linewidth=2)
# plot interface
plt.plot([x_int, x_int, x_reshock, x_reshock+q13.u*(t_max-t_reshock)],
[0, t_shock_int, t_reshock, t_max],
color="green", linewidth=2, linestyle="dashed")
plt.plot([0, q11.u*t_max],
[0, t_max],
color="green", linewidth=2, linestyle="dashed")
# plot expansion
t_exp_wall = x_Lwall/q20.left_acoustic()
plt.plot([0, x_Lwall, x_Lwall+q21.right_acoustic()*(t_max-t_exp_wall)],
[0, t_exp_wall, t_max],
color="orange", linewidth=1)
plt.plot([0, x_Lwall],
[0, x_Lwall/q21.left_acoustic()],
color="orange", linewidth=1)
[<matplotlib.lines.Line2D at 0x7f9b97eaa860>]