Gas-Kinetic Scheme for Partially Ionized Plasma in Hydrodynamic Regime

For partially ionized plasma, the kinetic equation can be written as[35]:

where $f_{\alpha}=f_{\alpha}(t,\boldsymbol{x},\boldsymbol{u})$ is the distribution function for species $\alpha$ ( $\alpha=i$ for ion and $\alpha=e$ for electron, $\alpha=n$ for neutral) at space and time $(\boldsymbol{x},t)$ and microscopic translational velocity $\boldsymbol{u}$. $\boldsymbol{a}_{\alpha}$ is the Lorenz acceleration taking the form

For neutral species, $q_{n}=0$, thus $\boldsymbol{a}_{n}=0$. $Q_{\alpha}=\sum_{k=1}^{m}Q_{\alpha k}(f_{\alpha},f_{k})$ is the collision operator of species $\alpha$ between species $k$, where $m$ is total number of species in the system. In this work, $m=3$. In the Maxwell equations, $\boldsymbol{E}$ and $\boldsymbol{B}$ are the electric field strength and magnetic induction, $\boldsymbol{c}$ is the speed of light, and $\epsilon_{0}$ is the vacuum permittivity. $n_{\alpha}$ is number density of species $\alpha$.In this work, the charged species in the system are just protons and electrons, then the electric current is $\boldsymbol{J}=e\left(n_{i}\boldsymbol{U}_{i}-n_{e}\boldsymbol{U}_{e}\right)$ and the charge density is $q=e(n_{i}-n_{e})$, where $\boldsymbol{U}$ is macroscopic velocity and $e$ is the charge of a proton.

Collision between multiple species is modeled by the relaxation model by Andries, Aoki, and Perthanme[36], which is,

where $g_{\alpha}^{M}$ is a Maxwellian distribution,

and post-collision temperature and velocity are chosen as:

where $\mu_{\alpha k}=m_{\alpha}m_{k}/(m_{\alpha}+m_{k})$ is reduced mass, mean relaxation time $\tau_{\alpha}$ is determined by $1/\tau_{\alpha}=\sum_{k=1}^{m}\chi_{\alpha k}n_{k}$, and interaction coefficient $\chi_{\alpha k}$ for hard sphere model is [38]:

In this above formula, $d_{\alpha},d_{k}$are the diameters of the particles and can be approximated by

where Kn is Knudsen number, L is reference length

To satisfy the divergence constraint, the Perfect Hyperbolic Maxwell equations (PHM) are used to reformulate the Maxwell equations as

where $\phi,\psi$ are artificial correction potentials to accommodate divergence errors traveling at speed $\gamma c$ and $\chi c$[39, 40].

First, the kinetic governing equations are non-dimensionalized, resulting in a dimensionless kinetic system describing the behavior of electrons, protons, and neutrals. The Chapman-Enskog method can be applied to derive a three-fluid model from the kinetic system. Within this three-fluid model, the electron and ion fluids constitute a two-fluid subsystem. By analyzing an additional small parameter, the ideal Magnetohydrodynamic (MHD) limit is obtained from the two-fluid subsystem. In this limit, the electron and ion fluids are combined into a single conducting fluid. The ideal MHD equations describe the dynamics of this conducting fluid. The end result is a coupled system of equations for the ideal MHD conducting fluid and the separate neutral fluid. This coupled system describes the ideal partially ionized plasma in the MHD limit in Eq.(1).

The reference quantities are defined as follows:

where charateristic velocity $u_{0}$ is thermal velocity of ions, and $m_{0}=m_{i}$ is mass of ions. $\rho_{0}$ is reference density, $E_{0}$ and $B_{0}$ are reference electric and magnetic field strength, $a_{0}$ is reference acceleration, $f_{0}$ is reference distribution function, $c_{0}$ is reference sound speed and $\gamma$ is specific heat. Then variables are non-dimensionalized as:

Inserting the normalized variables into the BGK-Maxwell system, the following dimensionless BGK-Maxwell system is obtained:

where $\tilde{\lambda_{D}}$ is normalized Debye length, $\tilde{r}_{L}$ is normalized larmor radius. For simplicity, in the following of this paper, all the hats for nondimensionalized variables are omitted.

Zeroth-order asymptotic solution of $f$ based on Chapman-Enskog expansion is $f=g+O\left(\tau^{1}\right)$[33]. Substitute the solution into BGK equation and take moments, a three-fluid system composed of neutral species, ions, and electrons is obtained,

where $S_{\alpha}$ and $Q_{\alpha}$ are momentum and energy exchange between species $\alpha$ and other species in the system.

Except the first three equations for the neutral gas, the rest system in Eq.(11) forms an ion-electron two-fluid subsystem, based on which Hall-effect MHD and ideal MHD can be derived[35, 41].

Define center-of-mass velocity as

Denote mass ratio $\epsilon=m_{e}/m_{i}$, $(1+\epsilon)\boldsymbol{U}=\boldsymbol{U}_{i}+\epsilon\boldsymbol{U}_{e}$. For $\mathcal{O}(\epsilon^{0})$ balance,

Substituting the above approximation into an electron momentum equation in the two-fluid subsystem,

The right-hand side of the above equation contains four terms: the first term represents electric resistivity, the second term corresponds to the Hall effect, and the last two terms describe the effects of electron inertia and pressure. In the limit $\epsilon\rightarrow 0$, the electron’s inertia can be negected and the electron momentum equation gives the generalized Ohm’s law

For non-relativistic flows, the displacement current is negligible in view of,

where $U$ and $L$ is the characteristic macroscale velocity and spatial length of plasma, and $U^{2}/c^{2}\ll 1$. When $\lambda_{D}\sim c^{-1}\rightarrow 0$ (i.e., $\lambda c^{-1}=1$), the Ampère’s law then becomes

The above low-frequency Ampère’s law indicates that $\nabla\cdot\boldsymbol{J}=0$, and therefore in this regime, the plasma is quasi-neutral, namely $n_{i}\approx n_{e}$.

In such a regime, the two fluid equations reduce to one fluid Hall-MHD equation. The Hall term and the electron pressure term are on the order of Larmor radius. Then Hall-MHD can be written as,

In the limit $r_{L}\rightarrow 0$, the Hall current and electron pressure in Eq.(12) are negligible. And ignoring collisions between electrons and protons i.e. $\nu_{ie}=0$, electric resistivity thus can be dropped. Then the generalized Ohm’s law reduces to the ideal Ohm’s law,

So combining electron and ion momentum equations, ideal MHD equation can be written as,

In this limit, the three-fluid system in Eq.(11) turns to a neutral and ideal MHD two-fluid system.

In summary, as shown in Fig.1 in the limit of $r_{L}\rightarrow 0,\lambda_{D}\sim c^{-1}\rightarrow 0,m_{e}/m_{i}\rightarrow 0,\tau\rightarrow 0,\nu_{ie}\rightarrow 0$, the kinetic system in Eq.(10) can describe the ideal partially ionized plasma in Eq.(1).

In the framework of the FVM, the cell averaged conservative variables for species $\alpha$ is $(\boldsymbol{W}_{\alpha})_{i}=((\rho_{\alpha})_{i},(\rho_{\alpha}\boldsymbol{U}_{\alpha})_{i},(\rho_{\alpha}\mathscr{E}_{\alpha})_{i})$ on a physical cell $\Omega_{i}$ are defined as

where $|\Omega_{i}|$ is the volume of cell $\Omega_{i}$. For a discretized time step $\Delta t=t^{n+1}-t^{n}$, the evolution of $(\boldsymbol{W}_{\alpha})_{i}$ is

where $l_{s}\in\partial\Omega_{i}$ is the cell interface with center $\boldsymbol{x}_{s}$ and outer unit normal vector $\boldsymbol{n}_{s}$. $|l_{s}|$ is the area of the cell interface. $(\bar{\boldsymbol{W}}_{\alpha})_{i}=((\rho_{\alpha})_{i},(\rho_{\alpha}\bar{\boldsymbol{U}}_{\alpha})_{i},(\rho_{\alpha}\bar{\mathscr{E}}_{\alpha})_{i})$ where $(\bar{\boldsymbol{U}}_{\alpha})_{i}$ and $(\bar{\mathscr{E}}_{\alpha})_{i}$ are post-collision velocity and energy in AAP model as Eq.(5). $(\boldsymbol{S}_{\alpha})_{i}$ is source term due to electromagnetic force. The numerical flux across interface $(\mathscr{F}_{\boldsymbol{W}_{\alpha}})_{s}$ can be evaluated from distribution function at the interface,

where $\boldsymbol{\Psi}=(1,\boldsymbol{u},\frac{1}{2}(\boldsymbol{u}^{2}+\boldsymbol{\xi}^{2}))$ is the conservative moments of distribution functions with $\boldsymbol{\xi}=(\xi_{1},\xi_{2},\cdots,\xi_{n})$ the internal degree of freedom. $\mathrm{d}\boldsymbol{\Xi}=\mathrm{d}\boldsymbol{u}d\boldsymbol{\xi}$ is the volume element in the phase space. The evaluation of the distribution function will be presented in section 3.2.

To be specific, the elementwise equations in Eq.(13) are

In Eq.(15), the source term due to cross-species momentum and energy exchange will be evaluated by the operator splitting method. In the Euler limit,$\tau_{\alpha}\rightarrow 0$, $(\bar{\boldsymbol{W}_{\alpha}})\rightarrow\boldsymbol{W}_{\alpha})$ is obtained at every time step. Lorentz source term is split and coupled with source terms in PHM so as to get coupled evolution between fluid species and electromagnetic field. The interaction equations can be solve by Crank-Nicolson scheme introduced in Section 3.4.

The cell averaged quantities $\boldsymbol{Q}_{i}$ for electromagnetic variables $\boldsymbol{Q}=\left(E_{x},E_{y},E_{z},B_{x},B_{y},B_{z},\phi,\psi\right)$ in a cell are defined as

where $(\mathscr{F}_{\boldsymbol{Q}})_{s}$ is numerical flux across the cell interface $l_{s}$ which will be presented in section 3.3. The time evolution formula is

where $(\mathscr{F}_{\boldsymbol{Q}})_{s}$ is numerical flux across a cell interface, which will be presented in Section 3.3. $\boldsymbol{S}_{\boldsymbol{Q}}$ are sources terms in PHM equations. The componentwise equations are:

General numerical steps are listed as follows:

1. 1.

   Update conservative variable $\boldsymbol{W}_{\alpha}^{n}$ to $\boldsymbol{W}_{\alpha}^{*}$ considering net flux without force term across the cell interface.

2. 2.

   Update conservative variable $\boldsymbol{W}_{\alpha}^{*}$ to $\boldsymbol{W}_{\alpha}^{**}$ considering momentum and energy exchange between different species.

3. 3.

   Update electromagnetic field $\boldsymbol{E}^{n}\rightarrow\boldsymbol{E}^{*}$, $\boldsymbol{B}^{n}\rightarrow\boldsymbol{B}^{n+1}$, ${\phi}^{n}\rightarrow{\phi}^{*}$ and ${\psi}^{n}\rightarrow{\psi}^{n+1}$ by net flux across physical cell interface.

4. 4.

   Incorporate the interaction between electromagnetic field and charged species $\boldsymbol{W}_{\alpha}^{**}$ to $\boldsymbol{W}_{\alpha}^{n+1}$,$\boldsymbol{E}^{*}\rightarrow\boldsymbol{E}^{n+1}$,${\phi}^{*}\rightarrow{\phi}^{n+1}$.

After four steps, all variables are evolved from $t^{n}$ to $t^{n+1}$.

In this section, the solution distribution function at a cell interface is constructed so as to calculate numerical flux in Eq.(14) in the two-dimensional (2D) case. The 2D BGK equation without force term can be written as

where $\boldsymbol{u}=(u,v)$. For simplicity, species subscript $\alpha$ is omitted here. To get the distribution function at a cell interface in Eq.(14), the time evolution solution at the interface can be written as,

where $(x_{i+1/2,j}=0,y_{i+1/2,j})=(0,0)$ is point on the interface. $\boldsymbol{\xi}=(w,\xi_{1},\xi_{2},\cdots,\xi_{n})$ is the internal degree of freedom. $f_{0}$ is the initial gas distribution function at $t=0$, and $g$ is equilibrium distribution at $({x}^{{}^{\prime}},y^{{}^{\prime}},t^{{}^{\prime}})$. $(x^{{}^{\prime}},y^{{}^{\prime}})=(x_{i+1/2,j}-u(t-t^{{}^{\prime}}),y_{i+1/2,j}-v(t-t^{{}^{\prime}}))$ are the particle trajectory. $\tau$ is mean relaxation time between successive collisions. $\tau_{n}$ is the numerical collision time. For the inviscid flow, the collision time $\tau_{n}$ is

where $C_{1}=0.01$ and $C_{2}=5$. $p_{l}$ and $p_{r}$ denote the pressure on the left and right sides of the cell interface. The inclusion of the pressure jump term is to increase the non-equilibrium transport mechanism in the flux function to mimic the physical process in the shock layer[33].

For the inviscid flow computation, the initial distribution at the cell interface can be reconstructed from cell average value as,

$g^{l},g^{r}$ are the equilibrium states at left and right cells. $a^{l},a^{r},b^{l},b^{r}$ is the spatial slope of the initial distribution along the x and y directions at the left and right cells, i.e

The slope of the Maxwellian distribution can be evaluated as,

where the specific form of $a_{1}-a_{4}$ can be found in paper [33].

Equilibrium distribution can be approximated as

where $g_{0}$ is the equilibrium state at the interface,

where $H(x)$ is Heaviside function. $\bar{a}^{l},\bar{a}^{r}$ are microscopic spatial slopes of equilibrium distributions, $\bar{A}$ is the temporal slope, which can be obtained from macroscopic variables. $g^{l}$ and $g^{r}$ is the Maxwellian distribution at the left and right of the cell interface.

Finally, the gas distribution function $f$ at a cell interface can be written as,

Substituting the above formula into Eq.(14), the numerical flux can be calculated.

1D numerical flux is illustrated here, for 2D or 3D problems, simply rotating coordinate can be used to get flux in another direction. The general expression for 1D PHM system is:

where

The numerical flux across interface $(i-1/2,j)$ is [37]:

where $\mathcal{A}_{1}^{+}=R_{1}\Lambda^{+}R_{1}^{-1}$ and $\mathcal{A}_{1}^{-}=R_{1}\Lambda^{-}R_{1}^{-1}$. $R_{1}$ is the matrix composed of right eigenvectors of $A_{1}$, and $\Lambda^{+}=diag((\lambda^{1})^{+},(\lambda^{2})^{+},\cdots,(\lambda^{8})^{+})$ with $\lambda^{+}=max(\lambda,0)$ and $\Lambda^{-}=diag((\lambda^{1})^{-},(\lambda^{2})^{-},\cdots,(\lambda^{8})^{-})$ with $\lambda^{-}=min(\lambda,0)$. $\lambda^{p}$ is the $p$th eigenvalue of $A_{1}$. Besides,$\Delta\boldsymbol{Q}_{i-1/2}=\boldsymbol{Q}_{i+1}-\boldsymbol{Q}_{i}$. The flux slope in Eq.(19) is