Quasi-Neutral Two Fluid Code (QNTF): 1D Brio-Wu Shock Tube

Reference: Amano2015CPC
Equations which we solve are written in a conservative form as

\begin{gathered} \frac{\partial }{{\partial t}}{\mathbf{U}} + \nabla \cdot {\mathbf{F}} = 0 \hfill \\ {\mathbf{U}} = \left( {\begin{array}{*{20}{c}} {{\rho _i} + {\rho _e}} \\ {{\rho _i}{{\mathbf{v}}_i}{\rho _e}{{\mathbf{v}}_e}} \\ {\frac{1}{2}{\rho _i}{\mathbf{v}}_i^2 + \frac{1}{2}{\rho _e}{\mathbf{v}}_e^2 + \frac{1}{{\gamma - 1}}\left( {{p_i} + {p_e}} \right) + \frac{{{{\mathbf{B}}^2}}}{{8\pi }}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} D \\ {\mathbf{M}} \\ K \end{array}} \right) \hfill \\ {\mathbf{F}} = \left( {\begin{array}{*{20}{c}} {{\rho _i}{{\mathbf{v}}_i} + {\rho _e}{{\mathbf{v}}_e}} \\ {{\rho _i}{{\mathbf{v}}_i}{{\mathbf{v}}_i} + {\rho _e}{{\mathbf{v}}_e}{{\mathbf{v}}_e} + \left( {{p_i} + {p_e} + \frac{{{{\mathbf{B}}^2}}}{{8\pi }}} \right){\mathbf{I}} - \frac{{{\mathbf{BB}}}}{{4\pi }}} \\ {\left( {\frac{1}{2}{\rho _i}{\mathbf{v}}_i^2 + \frac{\gamma }{{\gamma - 1}}{p_i}} \right){{\mathbf{v}}_i} + \left( {\frac{1}{2}{\rho _e}{\mathbf{v}}_e^2 + \frac{\gamma }{{\gamma - 1}}{p_e}} \right){{\mathbf{v}}_e} + \frac{c}{{4\pi }}{\mathbf{E}} \times {\mathbf{B}}} \end{array}} \right). \hfill \\ \end{gathered}

The subscripts $i$ and $e$ indicate ion and electron, respectively.

Mgnetic fields are obtained from \begin{gathered} \frac{1}{c}\frac{\partial }{{\partial t}}{\mathbf{B}} = - \nabla \times {\mathbf{E}} \hfill \\ \nabla \times {\mathbf{B}} = \frac{{4\pi }}{c}{\mathbf{J}}, \hfill \\ \end{gathered} where we assume the displacement current is ignored.

Electric fields are determined from the generalized Ohm's law, \begin{gathered} \left( {\Lambda + {c^2}\nabla \times \nabla \times } \right){\mathbf{E}} = - \frac{{\mathbf{\Gamma }}}{c} \times {\mathbf{B}} + \nabla \cdot {\mathbf{\Pi }} + \eta \Lambda {\mathbf{J}} \hfill \\ \Lambda = 4\pi {\rho _i}\frac{{q_i^2}}{{m_i^2}} + 4\pi {\rho _e}\frac{{q_e^2}}{{m_e^2}} = \omega _{pi}^2 + \omega _{pe}^2 \hfill \\ {\mathbf{\Gamma }} = 4\pi {\rho _i}\frac{{q_i^2}}{{m_i^2}}{{\mathbf{v}}_i} + 4\pi {\rho _e}\frac{{q_e^2}}{{m_e^2}}{{\mathbf{v}}_e} = \omega _{pi}^2{{\mathbf{v}}_i} + \omega _{pe}^2{{\mathbf{v}}_e} \hfill \\ {\mathbf{\Pi }} = \frac{{4\pi {q_i}}}{{{m_i}}}\left( {{\rho _i}{{\mathbf{v}}_i}{{\mathbf{v}}_i} + {p_i}{\mathbf{I}}} \right) + \frac{{4\pi {q_e}}}{{{m_i}}}\left( {{\rho _e}{{\mathbf{v}}_e}{{\mathbf{v}}_e} + {p_e}{\mathbf{I}}} \right) \hfill \\ \end{gathered}

Fluid values between ion and electron are related with the following equations, \begin{gathered} {\rho _i} = \frac{{{q_e}}}{{{m_e}}}D/\left( {\frac{{{q_e}}}{{{m_e}}} - \frac{{{q_i}}}{{{m_i}}}} \right),\,\,\,{\rho _e} = - {\rho _i}\frac{{{q_i}/{m_i}}}{{{q_e}/{m_e}}} \hfill \\ {{\mathbf{v}}_i} = \left( {\frac{{{q_e}}}{{{m_e}}}{\mathbf{M}} - \frac{c}{{4\pi }}\nabla \times {\mathbf{B}}} \right)/\frac{{{q_e}}}{{{m_e}}}D,\,\,{{\mathbf{v}}_e} = {v_i} - \frac{{{m_i}}}{{{q_i}}}\frac{c}{{4\pi {\rho _i}}}\nabla \times {\mathbf{B}} \hfill \\ {p_i} = \frac{{\gamma - 1}}{{1 + 1/\tau }}\left( {K - \frac{1}{2}{\rho _i}{\mathbf{v}}_i^2 - \frac{1}{2}{\rho _e}{\mathbf{v}}_e^2 - \frac{{{{\mathbf{B}}^2}}}{{8\pi }}} \right),\,\,{p_e} = {p_i}/\tau. \hfill \\ \end{gathered}

In the reference paper, HLL-UCT scheme is adopted.

Test Problem: Brio-Wu Shock Tube

The initial condition is \begin{gathered} {\left( {\begin{array}{*{20}{c}} \rho \\ p \\ {{B_x}/\sqrt {4\pi } } \\ {{B_y}/\sqrt {4\pi } } \end{array}} \right)_{{\text{left}}}} = \left( {\begin{array}{*{20}{c}} {1.0} \\ {1.0} \\ {0.75} \\ {1.0} \end{array}} \right),\,\,\,\,\,\,\,\,{\left( {\begin{array}{*{20}{c}} \rho \\ p \\ {{B_x}/\sqrt {4\pi } } \\ {{B_y}/\sqrt {4\pi } } \end{array}} \right)_{{\text{right}}}} = \left( {\begin{array}{*{20}{c}} {0.1} \\ {0.125} \\ {0.75} \\ { - 1.0} \end{array}} \right), \end{gathered} where $\rho\equiv\rho_e+\rho_i$ and $p\equiv p_e+p_i$.

The normalized ion inertial length ($\rho=1$) is $\lambda_i/L=10^{-3}$ where $L$ is a unit length. The ion to electron mass ratio is $m_i/m_e=100$.

In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 14
plt.rcParams['axes.linewidth'] = 1.5

Case 1: $\Delta x=2.5\times10^{-3}$

In [8]:
nx=800
f1 =open('data_Nx=800/rhoihst.d'  ,'rb')
f2 =open('data_Nx=800/vxihst.d'   ,'rb')
f3 =open('data_Nx=800/vyihst.d'   ,'rb')
f4 =open('data_Nx=800/rhoehst.d'  ,'rb')
f5 =open('data_Nx=800/vxehst.d'   ,'rb')
f6 =open('data_Nx=800/vyehst.d'   ,'rb')
f7 =open('data_Nx=800/bxhst.d'    ,'rb')
f8 =open('data_Nx=800/byhst.d'    ,'rb')
f9 =open('data_Nx=800/bzhst.d'    ,'rb')

rhoi =np.fromfile(f1, dtype='float64').reshape(-1,nx)
vix  =np.fromfile(f2, dtype='float64').reshape(-1,nx)
viy  =np.fromfile(f3, dtype='float64').reshape(-1,nx)
rhoe =np.fromfile(f4, dtype='float64').reshape(-1,nx)
vex  =np.fromfile(f5, dtype='float64').reshape(-1,nx)
vey  =np.fromfile(f6, dtype='float64').reshape(-1,nx)
bx   =np.fromfile(f7, dtype='float64').reshape(-1,nx)
by   =np.fromfile(f8, dtype='float64').reshape(-1,nx)
bz   =np.fromfile(f9, dtype='float64').reshape(-1,nx)

xx=np.linspace(0,2,nx)
plt.figure(figsize=(12,8))
plt.plot(xx,(rhoi+rhoe)[92,:])
plt.xlabel('X/L');plt.ylabel(r'$\rho$')
plt.xlim(0.5,1.5)
plt.show()

plt.figure(figsize=(12,8))
plt.subplot(311);plt.plot(xx,vix[19,:]);plt.ylabel(r'$v_{ix}$');plt.xlim(0.9,1.3);plt.ylim(-0.4,1.2)
plt.subplot(312);plt.plot(xx,viy[19,:]);plt.ylabel(r'$v_{iy}$');plt.xlim(0.9,1.3);plt.ylim(-2.5,0.5)
plt.subplot(313);plt.plot(xx,vey[19,:]);plt.ylabel(r'$v_{ey}$');plt.xlim(0.9,1.3);plt.ylim(-2.5,0.5)
plt.xlabel('X/L')
plt.show()

Case 2: $\Delta x=6.25\times10^{-3}$

In [15]:
nx=3200
f1 =open('data_Nx=3200/rhoihst.d'  ,'rb')
f2 =open('data_Nx=3200/vxihst.d'   ,'rb')
f3 =open('data_Nx=3200/vyihst.d'   ,'rb')
f4 =open('data_Nx=3200/rhoehst.d'  ,'rb')
f5 =open('data_Nx=3200/vxehst.d'   ,'rb')
f6 =open('data_Nx=3200/vyehst.d'   ,'rb')
f7 =open('data_Nx=3200/bxhst.d'    ,'rb')
f8 =open('data_Nx=3200/byhst.d'    ,'rb')
f9 =open('data_Nx=3200/bzhst.d'    ,'rb')

rhoi =np.fromfile(f1, dtype='float64').reshape(-1,nx)
vix  =np.fromfile(f2, dtype='float64').reshape(-1,nx)
viy  =np.fromfile(f3, dtype='float64').reshape(-1,nx)
rhoe =np.fromfile(f4, dtype='float64').reshape(-1,nx)
vex  =np.fromfile(f5, dtype='float64').reshape(-1,nx)
vey  =np.fromfile(f6, dtype='float64').reshape(-1,nx)
bx   =np.fromfile(f7, dtype='float64').reshape(-1,nx)
by   =np.fromfile(f8, dtype='float64').reshape(-1,nx)
bz   =np.fromfile(f9, dtype='float64').reshape(-1,nx)

xx=np.linspace(0,2,nx)
plt.figure(figsize=(12,8))
plt.plot(xx,(rhoi+rhoe)[216,:])
plt.xlabel('X/L');plt.ylabel(r'$\rho$')
plt.xlim(0.5,1.5)
plt.show()

plt.figure(figsize=(12,8))
plt.subplot(311);plt.plot(xx,vix[42,:]);plt.ylabel(r'$v_{ix}$');plt.xlim(0.9,1.3);plt.ylim(-0.4,1.2)
plt.subplot(312);plt.plot(xx,viy[42,:]);plt.ylabel(r'$v_{iy}$');plt.xlim(0.9,1.3);plt.ylim(-2.5,0.5)
plt.subplot(313);plt.plot(xx,vey[42,:]);plt.ylabel(r'$v_{ey}$');plt.xlim(0.9,1.3);plt.ylim(-2.5,0.5)
plt.xlabel('X/L')
plt.show()

Case 3: $\Delta x=1.5625\times10^{-4}$

In [13]:
nx=12800
f1 =open('data_Nx=12800/rhoihst.d'  ,'rb')
f2 =open('data_Nx=12800/vxihst.d'   ,'rb')
f3 =open('data_Nx=12800/vyihst.d'   ,'rb')
f4 =open('data_Nx=12800/rhoehst.d'  ,'rb')
f5 =open('data_Nx=12800/vxehst.d'   ,'rb')
f6 =open('data_Nx=12800/vyehst.d'   ,'rb')
f7 =open('data_Nx=12800/bxhst.d'    ,'rb')
f8 =open('data_Nx=12800/byhst.d'    ,'rb')
f9 =open('data_Nx=12800/bzhst.d'    ,'rb')

rhoi =np.fromfile(f1, dtype='float64').reshape(-1,nx)
vix  =np.fromfile(f2, dtype='float64').reshape(-1,nx)
viy  =np.fromfile(f3, dtype='float64').reshape(-1,nx)
rhoe =np.fromfile(f4, dtype='float64').reshape(-1,nx)
vex  =np.fromfile(f5, dtype='float64').reshape(-1,nx)
vey  =np.fromfile(f6, dtype='float64').reshape(-1,nx)
bx   =np.fromfile(f7, dtype='float64').reshape(-1,nx)
by   =np.fromfile(f8, dtype='float64').reshape(-1,nx)
bz   =np.fromfile(f9, dtype='float64').reshape(-1,nx)

xx=np.linspace(0,2,nx)
plt.figure(figsize=(12,8))
plt.plot(xx,(rhoi+rhoe)[449,:])
plt.xlabel('X/L');plt.ylabel(r'$\rho$')
plt.xlim(0.5,1.5)
plt.show()

plt.figure(figsize=(12,8))
plt.subplot(311);plt.plot(xx,vix[89,:]);plt.ylabel(r'$v_{ix}$');plt.xlim(0.9,1.3);plt.ylim(-0.4,1.2)
plt.subplot(312);plt.plot(xx,viy[89,:]);plt.ylabel(r'$v_{iy}$');plt.xlim(0.9,1.3);plt.ylim(-2.5,0.5)
plt.subplot(313);plt.plot(xx,vey[89,:]);plt.ylabel(r'$v_{ey}$');plt.xlim(0.9,1.3);plt.ylim(-2.5,0.5)
plt.xlabel('X/L')
plt.show()
In [18]:
!jupyter nbconvert --to html --template tmp.tpl QNTF_1D.ipynb

[NbConvertApp] Converting notebook QNTF_1D.ipynb to html
[NbConvertApp] Writing 816979 bytes to QNTF_1D.html