Hydromagnetic Shock Structure in the Presence of Cosmic Rays [Drury1981]¶
Here, we reproduce some results shown in Drury1981. The main conclusion of this paper is that when cosmic ray pressure is comparable than the gas pressure, the shock structure is modified and a precursor and subshock region appear. Especially, when the CR pressure is large enough, the shock structure becomes a smoothed transition.
Starting from gas hydrodynamic equations coupled with CR pressure equation, we can write down as follows \begin{gathered} u{\partial _x}\rho + \rho {\partial _x}u = 0; \hfill \\ \rho u{\partial _x}u + {\partial _x}{p_g} + {\partial _x}{p_c} = 0; \hfill \\ u{\partial _x}{p_g} + {\gamma _g}{p_g}{\partial _x}u = 0; \hfill \\ u{\partial _x}{p_c} + {\gamma _c}{p_c}{\partial _x}u = {\partial _x}\left( {\kappa {\partial _x}{p_c}} \right), \hfill \\ \end{gathered} where $p_g$ and $p_c$ denote the gas and CR pressure, respectively.
To solve this second-order differential equations numerically, we introduce $y=\partial_xp_c$ to make the system to first-order differential equations. After doing this process the system becomes, \begin{gathered} u{\partial _x}\rho + \rho {\partial _x}u = 0; \hfill \\ \rho u{\partial _x}u + {\partial _x}{p_g} + y = 0; \hfill \\ u{\partial _x}{p_g} + {\gamma _g}{p_g}{\partial _x}u = 0; \hfill \\ uy + {\gamma _c}{p_c}{\partial _x}u - {\partial _x}y = 0; \hfill \\ {\partial _x}{p_c} = y, \hfill \\ \end{gathered} And writing this in a matrix form, we obtain, \begin{equation} \left( {\begin{array}{*{20}{c}} u&\rho &0&0&0 \\ 0&{\rho u}&1&0&0 \\ 0&{{\gamma _g}{p_g}}&u&0&0 \\ 0&{{\gamma _c}{p_c}}&0&0&{ - 1} \\ 0&0&0&1&0 \end{array}} \right){\partial _x}\left( {\begin{array}{*{20}{c}} \rho \\ u \\ {{p_g}} \\ {{p_c}} \\ y \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0 \\ { - y} \\ 0 \\ { - uy} \\ y \end{array}} \right). \end{equation} Here we normalized $\bar{x}=x/L$, $L\equiv\kappa/u_1$, $\bar{p}_{g,c}=p_{g,c}/(\rho_1u_1^2)$, where the subscript $1$ indicates the far upstream value.
We use the $4$th-order Runge-Kutta scheme by setting initial conditions at the left boundary (far upstream). The initial gas and CR pressures are determined by two dimensionless parameters, $M$ and $N$ which are defined as \begin{equation} M\equiv\sqrt{\frac{\rho u^2}{\gamma_gp_g+\gamma_cp_c}},\,\,\,N\equiv\frac{p_c}{p_g+p_c}. \end{equation}
Note that this method is only valid for a smooth solution.
%config InlineBackend.figure_format='retina'
from spacepy import pycdf
import math as mt
import matplotlib.pyplot as plt
import numpy as np
import scipy.fft as sf
import datetime
import bisect
import matplotlib.patches as patches
import matplotlib.dates as mdates
from IPython.display import display,Math
import matplotlib as mpl
from ai import cdas
import datetime
import pandas as pd
from matplotlib.gridspec import GridSpec
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset
from matplotlib.ticker import ScalarFormatter
color='#2c2c2c'
#color='#ffffff'
mpl.rcParams['font.family' ]='Arial'#'Gill Sans'
mpl.rcParams['mathtext.default' ]='regular'
mpl.rcParams['font.size' ]=12
mpl.rcParams['figure.figsize' ]=(4,2)
mpl.rcParams['figure.dpi' ]=120
mpl.rcParams['axes.linewidth' ]=1.2
mpl.rcParams['xtick.major.size' ]=5
mpl.rcParams['xtick.minor.size' ]=4
mpl.rcParams['ytick.major.size' ]=5
mpl.rcParams['ytick.minor.size' ]=3
mpl.rcParams['xtick.major.width']=1
mpl.rcParams['ytick.major.width']=1
mpl.rcParams['xtick.direction' ]='out'
mpl.rcParams['ytick.direction' ]='out'
mpl.rcParams['axes.facecolor' ]='None'
mpl.rcParams['figure.facecolor' ]='None'
mpl.rcParams['axes.edgecolor' ]=color
mpl.rcParams['xtick.color' ]=color
mpl.rcParams['ytick.color' ]=color
mpl.rcParams['axes.labelcolor' ]=color
mpl.rcParams['text.color' ]=color
mpl.rcParams['patch.edgecolor' ]=color
mpl.rcParams['savefig.bbox' ]='tight'
mpl.rcParams['animation.html' ]='jshtml'
###color scheme for lines: https://www.nature.com/articles/nmeth.1618 ###
from cycler import cycler
line_cycler = (cycler(color =[color, '#e69f00', '#56b4e9', '#009e73', '#f0e442', '#0072b2', '#d55e00', '#cc79a7']))
#line_cycler = (cycler(color =[color, color, color, color])+
# cycler(linestyle=["-" , "--" , "-." , ":" ]))
mpl.rcParams['axes.prop_cycle']=line_cycler
#########################################################################
def TWR(A,gamg,gamc,nx,dx,isav):
Ahst=np.zeros((nx,5))
for ix in range(nx):
gc=np.zeros((4,5))
#---4th-order Runge-Kutta---#
gc[0,:]=f(A )
gc[1,:]=f(A+0.5*dx*gc[0,:])
gc[2,:]=f(A+0.5*dx*gc[1,:])
gc[3,:]=f(A+ dx*gc[2,:])
A=A+dx*(gc[0,:]+2*gc[1,:]+2*gc[2,:]+gc[3,:])/6
if(ix%isav==0): Ahst[ix//isav,:]=A[:]
return locals()
def f(A):
ro=A[0]; u=A[1]; pg=A[2]; pc=A[3]; y=A[4]
B=np.array([[u, ro, 0, 0, 0],
[0, ro*u, 1, 0, 0],
[0, gamg*pg, u, 0, 0],
[0, gamc*pc, 0, 0, -1],
[0, 0, 0, 1, 0]])
b=np.array([0, -y, 0, -u*y, y])
Binv=np.linalg.inv(B)
x=Binv.dot(b)
#x=np.linalg.solve(B,b)
return x
Low Mach number case¶
nx=4000; dx=0.01; isav=1
gamg=5/3;gamc=5/3
M=2; N=0.8
pg=1.0/M**2/(gamg+gamc*N/(1-N))
pc=N*pg/(1-N)
rho=1; u=1
y=1e-10
A=[rho,u,pg,pc,y]
data=TWR(A,gamg,gamc,nx,dx,isav)
locals().update(data)
xx=dx*np.arange(nx)
fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True, figsize=(9, 5))
ax[0,0].plot(xx,Ahst[:,0],'-')
ax[0,1].plot(xx,Ahst[:,1],'-')
ax[0,1].plot(xx,np.sqrt(gamg*Ahst[:,2]/Ahst[:,0]),'--')
ax[1,0].plot(xx,Ahst[:,2],'-')
ax[1,1].plot(xx,Ahst[:,3],'-')
ax[0,0].set_ylabel(r'$\rho$')
ax[0,1].set_ylabel(r'$U$')
ax[1,0].set_ylabel(r'$p_g$')
ax[1,1].set_ylabel(r'$p_c$')
ax[1,0].set_xlabel(r'$X$')
ax[1,1].set_xlabel(r'$X$')
fig.suptitle('M=%.1f, N=%.1f' %(M,N),y=0.95)
plt.show()
High Mach number case¶
nx=2000; dx=0.01; isav=1
gamg=5/3;gamc=4/3
M=13; N=0.8
pg=1.0/M**2/(gamg+gamc*N/(1-N))
pc=N*pg/(1-N)
rho=1; u=1
y=1e-4
A=[rho,u,pg,pc,y]
data=TWR(A,gamg,gamc,nx,dx,isav)
locals().update(data)
xx=dx*np.arange(nx)
fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True, figsize=(9, 5))
ax[0,0].plot(xx,Ahst[:,0],'-')
ax[0,1].plot(xx,Ahst[:,1],'-')
ax[0,1].plot(xx,np.sqrt(gamg*Ahst[:,2]/Ahst[:,0]),'--')
ax[1,0].plot(xx,Ahst[:,2],'-')
ax[1,1].plot(xx,Ahst[:,3],'-')
ax[0,0].set_ylabel(r'$\rho$')
ax[0,1].set_ylabel(r'$U$')
ax[1,0].set_ylabel(r'$p_g$')
ax[1,1].set_ylabel(r'$p_c$')
ax[1,0].set_xlabel(r'$X$')
ax[1,1].set_xlabel(r'$X$')
fig.suptitle('M=%.1f, N=%.1f' %(M,N),y=0.95)
plt.show()
