Rankine-Hugoniot Relation¶
The jump conditions in the de Hoffmann-Teller frame can be written $$\begin{array}{l} {\rho _{m1}}{U_{n1}} = {\rho _{m2}}{U_{n2}},\\ {\rho _{m1}}U_{n1}^2 + {\rho _{m1}}\frac{{V_{s1}^2}}{\gamma } + \frac{{B_{t1}^2}}{{2{\mu _0}}} = {\rho _{m2}}U_{n2}^2 + {\rho _{m2}}\frac{{V_{s2}^2}}{\gamma } + \frac{{B_{t2}^2}}{{2{\mu _0}}},\\ {\rho _{m1}}{U_{n1}}{U_{t1}} - {B_{n1}}\frac{{{B_{t1}}}}{{{\mu _0}}} = {\rho _{m2}}{U_{n2}}{U_{t2}} - {B_{n2}}\frac{{{B_{t2}}}}{{{\mu _0}}},\\ {\rho _{m1}}{U_{n1}}\left( {\frac{{V_{s1}^2}}{{\gamma - 1}} + \frac{{U_{t1}^2 + U_{n1}^2}}{2}} \right) + \frac{{{B_{t1}}}}{{{\mu _0}}}\left( {{B_{t1}}{U_{n1}} - {B_n}{U_{t1}}} \right)\\ \,\,\,\,\,\, = {\rho _{m2}}{U_{n2}}\left( {\frac{{V_{s2}^2}}{{\gamma - 1}} + \frac{{U_{t2}^2 + U_{n2}^2}}{2}} \right) + \frac{{{B_{t2}}}}{{{\mu _0}}}\left( {{B_{t2}}{U_{n2}} - {B_n}{U_{t2}}} \right),\\ {U_{n1}}{B_{t1}} - {U_{t1}}{B_{n1}} = {U_{n2}}{B_{t2}} - {U_{t2}}{B_{n2}},\\ {U_{n2}}{B_{t2}} - {U_{t2}}{B_{n2}} = 0\,\,\,\,(i.e.,\,\,the\,{E_z} = 0\,condition),\\ {B_{n1}} = {B_{n2}} \end{array}$$ Here $V_s^2=\frac{\gamma P}{\rho_m}$.
The Alfven Mach number is defined as $M_A=\frac{U_n}{V_{A}}=\frac{U_n\sqrt{\mu_0\rho_m}}{B}$.
The sound Mach number is defined as $M_s=\frac{U_n}{V_s}=U_n\sqrt{\frac{\rho_m}{\gamma P}}$.
Once five upstream parameters ($\rho_{m1}, U_{n1}, V_{s1}, B_{t1}$, and $B_{n1}$) are given, the all downstream parameters can be determined.
$$\begin{array}{l}
\frac{{{B_{t2}}}}{{{B_{t1}}}} = \frac{{M_{A1}^2 - {{\cos }^2}{\theta _{Bn1}}}}{{\left( {M_{A1}^2/r} \right) - {{\cos }^2}{\theta _{Bn1}}}},\\
{{{U_{t2}}}} = U_{t1}+\frac{(r-1)\tan\theta_1}{M_{I1}^2-r}U_{n1},\\
\frac{{{P_2}}}{{{P_1}}} = 1 + 2 \frac{{M_{A1}^2}}{\beta }\frac{{r - 1}}{r}\left[ {1 - \frac{{r[(r + 1)M_{A1}^2 - 2r{{\cos }^2}{\theta _{Bn1}}]}}{{{{2\left( {M_{A1}^2 - r{{\cos }^2}{\theta _{Bn1}}} \right)}^2}}}{{\sin }^2}{\theta _{Bn1}}} \right].
\end{array}$$
Here $U_{n1}$ corresponds to the shock speed.
The shock compression rate $r=\rho_{m2}/\rho_{m1}$ is determined by the shock adiabatic,
$$
{\left( {M_{A1}^2 - r{{\cos }^2}{\theta _{Bn1}}} \right)^2}\left( {M_{A1}^2 - \frac{{2r{\beta _1}}}{{(1 - \gamma )r + (1 + \gamma )}}} \right) - r{\sin ^2}{\theta _{Bn1}}M_{A1}^2\left( {\frac{{(2 - \gamma )r + \gamma }}{{(1 - \gamma )r + (1 + \gamma )}}M_{A1}^2 - r{{\cos }^2}{\theta _{Bn1}}} \right) = 0.
$$
$M_{I1}=M_{A1}/\cos\theta_{Bn1}$
note: The de Hoffmann-Teller frame does not mean that it is in the shock rest frame.
It does not exist at $\theta_{Bn}=90^{\circ}$.
from "Introduction to plasma physics"
%config InlineBackend.figure_format='retina'
import numpy as np
import matplotlib.pyplot as plt
from IPython import display
import math as mt
def shockadiabatic(MA,theta,beta,gam):
theta=theta*mt.pi/180.0
MI=MA/mt.cos(theta)
coeff=np.zeros(4)
if theta==90:
coeff=np.zeros(4)
coeff[0]=0.0#gam*beta+MI**2*(gam-1)
coeff[1]=-(gam-2.0)#-MI**2*(2*gam*beta+gam+1+(gam-2+gam*mt.cos(theta)**2)*MI**2)
coeff[2]=gam*beta+gam+(gam-1.0)*MA**2#MI**4*(gam*beta+gam+(gam+2)*mt.cos(theta)**2+(gam-1)*MA**2)
coeff[3]=-(gam+1.0)*MA**2#-(gam+1)*MA**2*MI**4
else:
coeff[0]=gam*beta+MI**2*(gam-1.0)
coeff[1]=-MI**2*(2.0*gam*beta+gam+1.0+(gam-2.0+gam*mt.cos(theta)**2)*MI**2)
coeff[2]=MI**4*(gam*beta+gam+(gam+2.0)*mt.cos(theta)**2+(gam-1.0)*MA**2)
coeff[3]=-(gam+1.0)*MA**2*MI**4
A=MA**4
B=-2*MA**2*mt.cos(theta)**2
C=mt.cos(theta)**4
D=MA**2*(1-gam)-2*beta
E=MA**2*(1+gam)
F=-MA**2*mt.sin(theta)**2
G=-(1-gam)*mt.cos(theta)**2
H=MA**2*(2-gam)-mt.cos(theta)**2*(1+gam)
I=gam*MA**2
#coeff[0]=C*D+F*G
#coeff[1]=C*E+D*B+F*H
#coeff[2]=E*B+A*D+F*I
#coeff[3]=A*E
root=np.roots(coeff)
rcmp=np.zeros(3,dtype=complex)
for ir in range(3):
if np.real(root[ir])>0 and np.isreal(root[ir])==True:
rcmp[ir]=root[ir]
else:
rcmp[ir]=0
return rcmp
def dwnstream(MA,theta,beta,gam,r):
theta=theta*mt.pi/180.0
MI=MA/mt.cos(theta)
bt2=(MI**2-1.0)/(MI**2/r-1.0)
ut2=(r-1.0)*mt.tan(theta)/(MI**2-r)*MA
p2 =1+2*MA**2/beta*(r-1.)/r*(1.-(r*(r+1)*MA**2-2*r**2*mt.cos(theta)**2)/2/(MA**2-r*mt.cos(theta)**2)**2*mt.sin(theta)**2)
#ds =np.log(p2*r**gam)
return bt2,ut2,p2
def MA_r(MAmax,nma,theta,beta,gam):
MA=np.linspace(0,MAmax,nma)
r =np.zeros((3,nma))
bt2 =np.zeros((3,nma))
ut2 =np.zeros((3,nma))
p2 =np.zeros((3,nma))
for ima in range(nma):
r[:,ima]=np.real(shockadiabatic(MA[ima],theta,beta,gam))
#for ir in range(3):
# for ima in range(nma):
# if r[ir,ima]>1:
# #print(r[ir,ima])
# bt2[ir,ima],ut2[ir,ima],p2[ir,ima]=dwnstream(MA[ima],theta,beta,gam,r[ir,ima])
plt.figure()
plt.xlabel(r'$M_A$');plt.ylabel('r')
plt.plot(MA,r[0,:],'.')
plt.plot(MA,r[1,:],'.')
plt.plot(MA,r[2,:],'.')
plt.ylim(1,4)
#plt.figure()
#plt.xlabel(r'$M_A$');plt.ylabel('bt2/bt1')
#plt.plot(MA,bt2[0,:],'.')
#plt.plot(MA,bt2[1,:],'.')
#plt.plot(MA,bt2[2,:],'.')
plt.show()
def th_r(thmax,nth,MA,beta,gam):
th=np.linspace(0,thmax,nth)*mt.pi/180
r =np.zeros((3,nth))
for ith in range(nth):
r[:,ith]=np.real(shockadiabatic(MA,th[ith],beta,gam))
plt.figure()
plt.xlabel(r'$\theta_{Bn1}$');plt.ylabel('r')
plt.plot(th*180/mt.pi,r[0,:],'.')
plt.plot(th*180/mt.pi,r[1,:],'.')
plt.plot(th*180/mt.pi,r[2,:],'.')
plt.ylim(1,4)
plt.show()
def beta_r(betamax,nbt,MA,theta,gam):
beta=np.linspace(0,betamax,nbt)
r =np.zeros((3,nbt))
for ibt in range(nbt):
r[:,ibt]=np.real(shockadiabatic(MA,theta,beta[ibt],gam))
plt.figure()
plt.xlabel(r'$\beta_1$');plt.ylabel('r')
plt.plot(beta,r[0,:],'.')
plt.plot(beta,r[1,:],'.')
plt.plot(beta,r[2,:],'.')
plt.ylim(1,5)
plt.show()
MA=6.5;theta=30;beta=1.;gam=5/3
r=shockadiabatic(MA,theta,beta,gam)
print(r)
#bt2, ut2, p2 note: ut2 is normalized by VA
dwnstream(MA,theta,beta,gam,np.real(r[2]))
MAmax=10;nma=2000
theta=60.
gam=5/3
beta=1
MA_r(MAmax,nma,theta,beta,gam)
thmax=90;nth=2000
MA=1.2
gam=5/3
beta=0.1
th_r(thmax,nth,MA,beta,gam)
betamax=10;nbt=2000
MA=1.2
theta=10
gam=5/3
beta_r(betamax,nbt,MA,theta,gam)
