Physics-Based-Adaptive Plasma Model for High-Fidelity Numerical Simulations

A physics-based-adaptive plasma model and an appropriate computational algorithm are developed to numerically simulate plasma phenomena in high fidelity. The physics-based-adaptive plasma model can be dynamically refined based on the local plasma conditions to increase model fidelity uniformity throughout the domain at all times of the simulation. The adaptive plasma model uses continuum representations of the plasma, which include a kinetic Vlasov model for the highest fidelity, multi-fluid 5N-moment plasma model, and single-fluid MHD model for the lowest fidelity. The models include evolution equations for the electromagnetic fields, electron species, ion species, and neutral species. A nodal discontinuous Galerkin finite element method is implemented and is coupled with various implicit and explicit Runge-Kutta methods. Various model coupling techniques are investigated for a 5N-moment multi-fluid models with a Vlasov-Maxwell model, and a 5N-moment two-fluid model with an MHD model. Continuum plasma models using consistent normalizations and identical spatial representations provide straightforward and accurate coupling between the models. The solution approach offers high-order accuracy and computational efficiency. Target compute platforms are heterogeneous computer architectures using a compute model that minimizes data movement.

1. Introduction

There are a wide variety of computational plasma models available which balance physical accuracy with simplifying approximations. When the assumptions for the reduced models are justified, the predictions from the reduced models approach those from higher fidelity models while requiring a fraction of the computational resources. However, there are a wide range of scientific and engineering applications which may violate the assumptions of a given reduced model only in a localized subset of the domain. The domain boundaries between the valid physical models may evolve as the dynamics progress.

Some notable examples include the current sheet behavior during geomagnetic reconnection [1, 2], plasma photonic crystals [3, 4], field-reversed configurations (FRC's) [5–7], and Z-pinches [8–10]. Reddell [2] showed that the distribution plasma function within the thin sheath boundary was not near Maxwellian, but regions inside the reconnected region or far away from the sheath were still near Maxwellian. Tummel et al. [10] demonstrated that kinetic effects can have a stabilizing effect on a flow-shear stabilized Z-pinch, while the bulk plasma may still be mostly well described by a fluid model. However, some physical phenomena are not adequately captured with reduced models and require models with higher fidelity. Drift turbulence [11] is accurately represented by multi-fluid plasma models [12] but not by single-fluid plasma models.

Spatial coupling of multiple plasma models offers the potential to minimize the computational cost for the required degree of physical fidelity throughout the domain. The investigation is performed using a high-order discontinuous Galerkin finite element method [13–15], specifically involving the continuum kinetic multi-species plasma model, 5N-moment multi-fluid plasma model, and magnetohydrodynamics (MHD) models, implemented in the WARPXM (Washington Approximate Riemann Plasma) codes [12, 16, 17], which provides a general framework for performing parallel computational plasma physics simulations. It is demonstrated that transition mixing layers do not necessarily guarantee conservation of important physical quantities such as mass, momentum, or energy. Nevertheless, coupling between the models can be performed consistently. Two such numerical fluxes are derived, of which one flux is consistent with the MHD and two-fluid plasma models, and a second flux guarantees conservation properties but without a consistent definition of the underlying variables. Investigations of spatial coupling of the continuum kinetic and 5N-moment multi-fluid plasma models are also presented. The continuum kinetic model is implemented on a mixed structured/unstructured phase space finite element mesh to be able to handle complex physical space geometry while maintaining the computational efficiency of a rectilinear mesh in velocity space.

Related research switches between fluid and kinetic equations solved with a finite volume method using the Knudsen number and gradients of primitive variables as switching criteria between models within the domain [18, 19] or have a static region solved with a kinetic particle method embedded in domain modeled with the MHD model [20]. Other work considers a hybrid approach where charged particles may be simulated using a fluid model while neutrals are simulated kinetically [21] or the ions are modeled as kinetic particles and the electrons as a fluid [22]. There has also been work to couple different numerical methods [23] investigating a blended finite element method to solve the five moment multi-fluid plasma model [12].

In this work, coupling procedures between plasma models are presented, including between the MHD and 5N-moment multi-fluid models as well as between the 5N-moment multi-fluid model and continuum kinetics. An implementation of the kinetic model on a mixed structured/unstructured mesh is also described that facilitates coupling to the 5N-moment multi-fluid plasma model. Section 2 presents the governing equations of the plasma models in normalized form. The numerical method using discontinuous Galerkin for the spatial representation and Runge-Kutta time integration are described in section 3. Validation of the continuum kinetic plasma model is presented in section 5. Coupling the plasma models is described in section 4 and results from numerical testing are presented in sections 6 and 7.

2. Plasma Models

The plasma model with the highest physical fidelity among continuum models statistically represents the plasma constituents with continuous probability distribution functions f_α in phase space (r,v) for each species α. The evolution of f_α is described by the Boltzmann transport equation,

$\begin{array}{cc}\frac{\partial f_{\alpha }}{\partial t}+\nabla ·(\mathit{v}{f}_{\alpha })+\left[\frac{q_{\alpha }}{m_{\alpha }}(\mathit{E}+\mathit{v}\times \mathit{B})f_{\alpha }\right]={\left(\frac{\partial f_{\alpha }}{\partial t}\right)}_{coll},& (1)\end{array}$

where q_α and m_α are the charge and mass of species α. Charged species are accelerated by the Lorentz force generated by the electric and magnetic fields, E and B. The right hand side of Equation (1) accounts for collisional effects, which are often negligible in high temperature plasmas. The resulting Vlasov equation is expressed in normalized form as

$\begin{array}{cc}\frac{\partial f_{\alpha }}{\partial t}+\nabla ·(\mathit{v}{f}_{\alpha })+\left(\frac{L}{{\delta }_p}\right){\nabla }_{\mathit{v}}.\left[\frac{Z_{\alpha }}{A_{\alpha }}(\mathit{E}+\mathit{v}\times \mathit{B})f_{\alpha }\right]=0,& (2)\end{array}$

where L is the characteristic length, δ_p is the proton skin depth, Z_α is the ionization state of species α, and A_α is the atomic mass of species α. Each species is evolved with an instance of Equation (2). The normalization is described in section 2.1. The zeroth and first moments of the distribution functions f_α provides the charge density and current density

$\begin{array}{cc}\begin{array}{rc}{\rho }_c=& {\displaystyle \sum _{\alpha }}{Z}_{\alpha }\int f_{\alpha }(\mathit{r},\mathit{v})d\mathit{v},\end{array}& (3)\end{array}$

$\begin{array}{cc}\begin{array}{rc}\mathit{j}=& {\displaystyle \sum _{\alpha }}{Z}_{\alpha }\int \mathit{v}{f}_{\alpha }(\mathit{r},\mathit{v})d\mathit{v},\end{array}& (4)\end{array}$

which couples the Vlasov equations to Maxwell's equations, given in normalized form as

$\begin{array}{cc}\begin{array}{r}\frac{\partial \mathit{B}}{\partial t}+\nabla \times \mathit{E}=0,\end{array}& (5)\end{array}$

$\begin{array}{cc}\begin{array}{r}-\frac{1}{{({\omega }_p\tau )}^2}{\left(\frac{L}{{\delta }_p}\right)}^2\frac{\partial \mathit{E}}{\partial t}+\nabla \; \times \mathit{B}=\left(\frac{L}{{\delta }_p}\right)\mathit{j},\end{array}& (6)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{1}{{({\omega }_p\tau )}^2}\frac{L}{{\delta }_p}\nabla ·\mathit{E}={\rho }_c,\end{array}& (7)\end{array}$

$\begin{array}{cc}\begin{array}{rc}\nabla ·\mathit{B}=& 0,\end{array}& (8)\end{array}$

where ω_p is the proton plasma frequency and τ is the characteristic time. Pressure can also be related by the second moment of the distribution function

$\begin{array}{cc}\begin{array}{r}{P}_{\alpha }=n_{\alpha }{T}_{\alpha }=\frac{m_{\alpha }}{3}\int {(\mathit{v}-{\mathit{u}}_{\alpha })}^2{f}_{\alpha }(\mathit{r},\mathit{v})d\mathit{v},\end{array}& (9)\end{array}$

where n_α and u_α are the fluid density and velocity defined by the zeroth and first moments of the distribution function.

Moments of Equation (1) combined with appropriate closure relations result in the fluid plasma models, which have reduced dimensionality from 6D phase space to 3D physical space. The 5N-moment plasma multi-fluid model is common [12], but higher moment models, such as the 13N-moment multi-fluid plasma model [24], are also possible.

The 5N-moment multi-fluid plasma model for a species α is described by the normalized equation system

$\begin{array}{cc}\begin{array}{r}\frac{\partial {\rho }_{\alpha }}{\partial t}+\nabla ·{\mathit{p}}_{\alpha }=0,\end{array}& (10)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial {\mathit{p}}_{\alpha }}{\partial t}+\nabla ·\left(\frac{{\mathit{p}}_{\alpha }\otimes {\mathit{p}}_{\alpha }}{{\rho }_{\alpha }}+P_{\alpha }\overleftrightarrow{I}\right)=\left(\frac{L}{{\delta }_p}\right)\left(\frac{Z_{\alpha }}{A_{\alpha }}\right)\left({\rho }_{\alpha }\mathit{E}+{\mathit{p}}_{\alpha }\times \mathit{B}\right),\end{array}& (11)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial e_{\alpha }}{\partial t}+\nabla ·\left((e_{\alpha }+P_{\alpha })\frac{{\mathit{p}}_{\alpha }}{{\rho }_{\alpha }}\right)=\left(\frac{L}{{\delta }_p}\right)\left(\frac{Z_{\alpha }}{A_{\alpha }}\right){\mathit{p}}_{\alpha }·\mathit{E},\end{array}& (12)\end{array}$

where ρ_α is the mass density, P_α is the momentum, P_α is the scalar pressure, and e_α is the total fluid energy. These are coupled with the electric and magnetic fields E and B, which are evolved with Maxwell's equations, Equations (5), (6). Note this form of the multi-fluid plasma model neglects collisions and atomic reactions, which can be included [25]. The multi-fluid plasma model is a significant reduction in the physical fidelity of the kinetic plasma model and has a corresponding lower computational cost.

Further reductions in physical fidelity are possible by restriction to a two-fluid system of ions and electrons [26], i.e., α = i, e, and application of the asymptotic approximations of single-fluid MHD [17]—infinite ion/electron mass ratio and infinite speed of light. The resulting reduced plasma model is the two-temperature Hall-MHD model, which is given in normalized form as

$\begin{array}{cc}\begin{array}{r}\frac{\partial \rho }{\partial t}+\nabla ·\mathit{p}=0,\end{array}& (13)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial \mathit{p}}{\partial t}+\nabla ·\left(\frac{\mathit{p}\otimes \mathit{p}}{\rho }-\mathit{B}\otimes \mathit{B}+\left(P+\frac{\mathit{B}·\mathit{B}}{2}\right)\overleftrightarrow{\text{I}}\right)=0,\end{array}& (14)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial e_{\alpha }}{\partial t}+\nabla ·\left(\left(e_{\alpha }+P_{\alpha }\right)\frac{{\mathit{p}}_{\alpha }}{{\rho }_{\alpha }}\right)=\left(\frac{L}{{\delta }_p}\right)\left(\frac{Z_{\alpha }}{A_{\alpha }}\right){\mathit{p}}_{\alpha }·\mathit{E},\end{array}& (15)\end{array}$

$\begin{array}{cc}\begin{array}{r}\rho \mathit{E}+\mathit{p}\times \mathit{B}=\left(\frac{A_i}{Z_i}\right)\left(\left(\frac{{\delta }_p}{L}\right)\nabla P_e+\mathit{j}\times \mathit{B}\right)\mathrm{\text{.}}\end{array}& (16)\end{array}$

The last expression is the generalized Ohm's law, which states the relationship between the electric field and other MHD variables. The electric field is eliminated by combining the generalized Ohm's law with Faraday's law, Equation (5).

The MHD variables are related to the two-fluid variables through a center-of-mass formulation to maintain consistency between the two plasma models, with the definitions

$\begin{array}{cc}\begin{array}{r}\rho ={\rho }_i+{\rho }_e,\end{array}& (17)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{p}={\mathit{p}}_i+{\mathit{p}}_e,\end{array}& (18)\end{array}$

$\begin{array}{cc}\begin{array}{r}P=P_i+P_e,\end{array}& (19)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{j}=\frac{Z_i}{A_i}{\mathit{p}}_i+\frac{Z_e}{A_e}{\mathit{p}}_e=\left(\frac{{\delta }_p}{L}\right)\nabla \; \times \mathit{B},\end{array}& (20)\end{array}$

where the last expression is the low-frequency form of Ampère's law that results from the infinite speed of light approximation.

The plasma model can be reduced further by assuming a negligible skin depth, resulting in the ideal MHD equations. The species temperatures are assumed equal (T_i = T_e). This is justified since the terms which could lead to different temperatures depend either on non-thermal distribution functions and collisional effects, or from an appreciable skin depth. The governing equations for the ideal MHD plasma model can be written in conservative form as

$\begin{array}{cc}\begin{array}{r}\frac{\partial \rho }{\partial t}+\nabla ·\mathit{p}=0,\end{array}& (21)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial \mathit{p}}{\partial t}+\nabla ·\left(\frac{\mathit{p}\otimes \mathit{p}}{\rho }-\mathit{B}\otimes \mathit{B}+\left(P+\frac{\mathit{B}·\mathit{B}}{2}\right)\overleftrightarrow{\text{I}}\right)=0,\end{array}& (22)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial e}{\partial t}+\nabla ·\left(\left(e+P+\frac{\mathit{B}·\mathit{B}}{2}\right)\frac{\mathit{p}}{\rho }-\left(\frac{\mathit{p}·\mathit{B}}{\rho }\right)\mathit{B}\right)=0,\end{array}& (23)\end{array}$

$\begin{array}{cc}\begin{array}{r}\rho \mathit{E}+\mathit{p}\times \mathit{B}=0,\end{array}& (24)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{\partial \mathit{B}}{\partial t}+\nabla ·\left(\frac{\mathit{p}\otimes \mathit{B}-\mathit{B}\otimes \mathit{p}}{\rho }\right)=0,\end{array}& (25)\end{array}$

where the total energy is defined as

$e=\frac{P}{\gamma -1}+\frac{\mathit{p}·\mathit{p}}{2\rho }+\frac{\mathit{B}·\mathit{B}}{2},$

which contains the internal, kinetic, and magnetic energies.

2.1. Normalized Equation System

The normalization used in the model equations described above has been chosen such that the transition of Ampere's Law to its pre-Maxwell form (loss of displacement current) is well-behaved under the approximations used to derive the MHD model. Consider the non-normalized Ampere's Law

$\begin{array}{cc}\begin{array}{r}-{ϵ}_0\frac{\partial \mathit{E}}{\partial t}+\frac{1}{{\mu }_0}\nabla \; \times \mathit{B}=\mathit{j}\mathrm{\text{.}}\end{array}& (26)\end{array}$

It can be seen that the displacement term disappears as the speed of light approaches ∞, which is consistent in the chosen normalization. Vacuum permittivity and permeability are eliminated by introducing the reference proton plasma frequency and reference Alfvén speed through

$\begin{array}{cc}\begin{array}{r}{ϵ}_0=\frac{e^2{n}_0}{m_0{\omega }_p^2}\end{array}& (27)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mu }_0=\frac{B_0^2}{m_0{n}_0{V}_A^2},\end{array}& (28)\end{array}$

which leads to

$\begin{array}{cc}\begin{array}{r}-\frac{e^2{n}_0}{m_p{\omega }_p^2}\frac{E_0}{\tau }\frac{\partial \tilde{\mathit{E}}}{\partial \tilde{t}}+\frac{m_p{n}_0{V}_A^2}{B_0^2}\frac{B_0}{L}\tilde{\nabla }\times \tilde{\mathit{B}}=e{n}_0{v}_0\tilde{\mathit{j}},\end{array}& (29)\end{array}$

where all variables with a tilde are normalized by reference values, e.g., $\tilde{E}$ = E/E₀. Relating the reference values defines a set of normalizations of electromagnetic fields and scales,

$\begin{array}{cc}\begin{array}{r}{E}_0=v_0{B}_0\end{array}& (30)\end{array}$

$\begin{array}{cc}\begin{array}{r}{v}_0=\frac{L}{\tau }=V_A\end{array}& (31)\end{array}$

and introducing proton cyclotron frequency,

$\begin{array}{cc}\begin{array}{r}{\omega }_c=\frac{e{B}_0}{m_p},\end{array}& (32)\end{array}$

leads to a normalized equation

$\begin{array}{cc}\begin{array}{r}-\frac{{({\omega }_c\tau )}^2}{{({\omega }_p\tau )}^2}\frac{\partial \mathit{E}}{\partial t}+\nabla \times \mathit{B}=({\omega }_c\tau )\mathit{j}\mathrm{\text{,}}\end{array}& (33)\end{array}$

where tildes have been dropped for clarity. This form of Ampere's Law remains consistent since μ₀ is related to the Alfvén speed, ϵ₀ → 0 as c → ∞ ($c^2=\frac{1}{{\mu }_0{ϵ}_0}$), implying ω_p → ∞ eliminating the displacement current in Equation (33). While maintaining the desired consistency, this choice of normalization also implies a magnetized plasma through the cyclotron frequency, which is not always appropriate. Thus, the normalized cyclotron frequency is replaced with a normalized proton skin depth by noticing that

$\begin{array}{cc}\begin{array}{r}{\omega }_c\tau =\frac{e{B}_0}{m_p}\frac{L}{V_A}=\frac{e{B}_0}{m_p}\frac{L\sqrt{{\mu }_0{m}_p{n}_0}}{B_0}=\sqrt{\frac{e^2{n}_0}{{ϵ}_0{m}_p}}\frac{L}{c}=\frac{{\omega }_{p}L}{c}=\frac{L}{{\delta }_p}.\end{array}& (34)\end{array}$

Applying this to Equation (33) yields Equation (6). All other equations are then normalized using the same reference values and non-dimensional terms.

3. Numerical Discretization Using Discontinuous Galerkin

The equations of the various plasma models described in section 2 can be expressed as a general set of coupled partial differential equations as

$\begin{array}{cc}\begin{array}{r}\mathit{R}(\mathit{q},\mathit{r},t)=\mathit{S}(\mathit{q},\mathit{r},t)+\nabla ·(\overleftrightarrow{\text{F}}(\mathit{q},\mathit{r},t))+\overleftrightarrow{\text{A}}(\mathit{q},\mathit{r},t)·\nabla \mathit{q}=0,\end{array}& (35)\end{array}$

where q is the vector of dependent variables. Any of the terms $\overleftrightarrow{\text{F}}$, $\overleftrightarrow{\text{A}}$, and S may be a function of temporal derivatives of q, but not spatial derivatives. A typical PDE of interest might have

$\begin{array}{cc}\begin{array}{r}\mathit{S}(\mathit{q},\mathit{r},t)=\frac{\partial \mathit{q}}{\partial t}+\mathit{G}(\mathit{q},\mathit{r},t)\end{array}& (36)\end{array}$

The general numerical procedure is to discretize all spatial derivatives to obtain a semi-discrete form containing only temporal derivatives (if any). Any temporal discretizations are handled using a method of lines approach and solved using standard ODE/DAE solvers. For the purposes of discussing the spatial discretization method these two terms are grouped together since they are handled in the same fashion by the spatial discretization. Equation (35) is spatially discretized using the discontinuous Galerkin (DG) method [13–15]. This type of PDE contains three primary categories of terms to be handled:

$\begin{array}{cc}\begin{array}{r}\nabla ·(\overleftrightarrow{\text{F}}(\mathit{q},\mathit{r},t))=\frac{\partial F_{ij}}{\partial r_j},\end{array}& (37)\end{array}$

$\begin{array}{cc}\begin{array}{r}\overleftrightarrow{\text{A}}(\mathit{q},\mathit{r},t)·\nabla \mathit{q}=A_{ijk}\frac{\partial q_k}{\partial r_j},\end{array}& (38)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{S}(\mathit{q},\mathit{r},t)=S_i,\end{array}& (39)\end{array}$

which correspond to conservation form, non-conservation form, and sources, respectively. Note the introduction of index notation which is used for compact representations of vectors and higher rank tensors. While it is theoretically possible to construct a system that mixes conservation form and non-conservation form fluxes, in practice it is difficult to produce a consistent numerical method for this case using the DG method. It is preferable to write the equation set in conservation form; however, it may not be possible or practical to do so for all systems of interest. As an alternative, the system can be written in non-conservation form. The conservation form flux tensor can be converted into a flux Jacobian and combined with ${\overleftrightarrow{\text{A}}}^{\prime }$ such that

$\begin{array}{cc}\begin{array}{r}\frac{\partial F_{ij}}{\partial r_j}+A_{ijk}^{\prime }\frac{\partial q_k}{\partial r_j}=\left(A_{ijk}^{\prime }+\frac{\partial F_{ij}}{\partial q_k}\right)\frac{\partial q_k}{\partial r_j}=A_{ijk}\frac{\partial q_k}{\partial r_j}.\end{array}& (40)\end{array}$

Thus, the only PDE's that need to be treated for this research are of the forms

$\begin{array}{cc}\begin{array}{r}{S}_i+\frac{\partial F_{ij}}{\partial r_j}=0,\end{array}& (41)\end{array}$

$\begin{array}{cc}\begin{array}{r}{S}_i+A_{ijk}\frac{\partial q_k}{\partial r_j}=0.\end{array}& (42)\end{array}$

The specific details of the spatial discretization of these two classes of PDE's are described in sections 3.1 and 3.2, respectively.

The computation domain Ω is discretized into a finite element mesh ${T}$, as shown in Figure 1. The elements are grouped into subdomains corresponding to the plasma model that will be solved. For example, elements where the MHD model will be used are denoted as the subdomain ${{T}}_{MHD}$, and the elements where the two-fluid model will be used are denoted as ${{T}}_{TF}$. Note that these subdomains do not necessarily have to be contiguous. Additionally, while normally there is no reason why these two subdomains cannot overlap with a transition region, for the present work it is assumed that no elements overlap.

FIGURE 1

Figure 1. Example 2D discretized domain ${T}$ with blue elements $K\in {{T}}_{MHD}$ and red elements $K\in {{T}}_{TF}$.

For an element K, denote its boundary as ∂K, and define the shared interface between two elements K₁ and K₂ to be ϵ = K₁ ∩ K₂. The set of all shared faces is denoted as ${{E}}^0$, the set of all exterior domain boundary faces is denoted as ∂Ω, and the set of all faces is denoted ${E}={{E}}^0\cup \partial \Omega$. Define the finite element spaces

$\begin{array}{cc}\begin{array}{r}{W}_{\text{Γ}}=\left\{w\in L^2(\text{Γ}):w{|}_K\in {{P}}_m(K),\forall K\in {{T}}_{\text{Γ}}\right\},\end{array}& (43)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{V}}_{\text{Γ}}=\left\{\mathit{v}\in {(L^2(\text{Γ}))}^d:\mathit{v}{|}_K\in {({{P}}_m(K))}^d,\forall K\in {{T}}_{\text{Γ}}\right\},\end{array}& (44)\end{array}$

where ${{P}}_m(K)$ is the set of all polynomials up to degree m restricted to the element K. A scalar variable of the original PDE is approximated by a field in W_Γ, and vector variables are approximated by a field in V_Γ. For example, mass density ρ ≈ ρ_h ∈ W_Γ, while momentum ρ ≈ ρ_h ∈ V_Γ. The complete solution for q is then composed of a tensor product of variables found in these various spaces, for example q = [ρ p e] ≈ q_h ∈ W_Γ ⊗ V_Γ ⊗ W_Γ. Define this tensor product finite element space to be Z_Γ.

A general DG method for Equation (35) consisting of a linear combination of terms given by forms of Equations (37)–(39) then seeks to find q_h ∈ Z_Γ such that

$\begin{array}{cc}\begin{array}{c}{\displaystyle {\int }_{\text{Γ}}}{R}_i({\mathit{q}}_h,\mathit{r},t)z_{i}dV=0,\forall \mathit{z}\in {\mathit{Z}}_{\text{Γ}},\end{array}& (45)\end{array}$

which is a general formulation. However, the research presented here focuses on a global Cartesian space, which is equivalent to stating that each component i must satisfy

$\begin{array}{cc}\begin{array}{rc}{\displaystyle {\int }_{\text{Γ}}}{R}_i({\mathit{q}}_h,\mathit{r},t)wdV& =0,\forall w\in W_{\text{Γ}}.\end{array}& (46)\end{array}$

A practical implementation would then choose a basis Φ_Γ which spans the entire space W_Γ such that ${\varphi }^{(\alpha )}\in {\Phi }_{\Gamma }$ and ${\varphi }^{(\beta )}\in {\Phi }_{\Gamma }$ are linearly independent for α ≠ β. The basis used in this research are the interpolation basis, producing a nodal DG scheme. In 1D nodes are co-located with the Gauss-Lobatto quadrature nodes. In higher dimensions, cube-like elements can be formed by taking a tensor product of 1D node locations. This is used to form the 1D-1V basis for Vlasov's equation. Nodes for triangular elements are located using the warping function method described in Hesthaven and Warburton [27].

The treatment for each of the different types of terms represented by Equations (37)–(39) is discussed in the following sections. The complete algorithm is a combination of the individual treatments. After the spatial discretization has been applied, the system becomes a coupled system of ordinary differential equations, which are discretized temporally using explicit or implicit Runge-Kutta methods. These methods are presented in section 3.5.

3.1. Conservation Form Fluxes

Consider the PDE

$\begin{array}{cc}\begin{array}{r}{S}_i+\frac{\partial F_{ij}}{\partial r_j}=0.\end{array}& (47)\end{array}$

Integrating over a test function ϕ^(α), each component i must satisfy

$\begin{array}{cc}\begin{array}{r}{\displaystyle \int }{\varphi }^{(\alpha )}\left(S_i+\frac{\partial F_{ij}}{\partial r_j}\right)dV=0,\forall {\varphi }^{(\alpha )}\in \mathbf{\Phi }\end{array}& (48)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\displaystyle \int }{\varphi }^{(\alpha )}{S}_{i}dV+\oint {\varphi }^{(\alpha )}{(F_{ij}{n}_j)}^{*}dS-{\displaystyle \int }{F}_{ij}\frac{\partial {\varphi }^{(\alpha )}}{\partial r_j}dV=0,\forall {\varphi }^{(\alpha )}\in \mathbf{\Phi }\end{array}& (49)\end{array}$

where ${(F_{ij}{n}_j)}^{*}$ is an appropriately chosen numerical flux. Expanding terms over the same basis,

$\begin{array}{cc}\begin{array}{r}{S}_i\approx S_i^{(\beta )}{\varphi }^{(\beta )},\end{array}& (50)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(F_{ij}{n}_j)}^{*}\approx {(F_{ij}{n}_j)}_{(\beta )}^{*}{\varphi }^{(\beta )},\end{array}& (51)\end{array}$

$\begin{array}{cc}\begin{array}{r}{F}_{ij}\approx F_{ij}^{(\beta )}{\varphi }^{(\beta )},\end{array}& (52)\end{array}$

where $S_i^{(\beta )}$, ${(F_{ij}{n}_j)}_{(\beta )}^{*}$, and $F_{ij}^{(\beta )}$ are not spatially varying. This allows the PDE to be written in the semi-discrete form

$\begin{array}{cc}\begin{array}{r}{S}_i^{(\beta )}{\displaystyle \int }{\varphi }^{(\alpha )}{\varphi }^{(\beta )}dV+{(F_{ij}{n}_j)}_{(\beta )}^{*}\oint {\varphi }^{(\alpha )}{\varphi }^{(\beta )}dS\\ -F_{ij}^{(\beta )}{\displaystyle \int }{\varphi }^{(\beta )}\frac{\partial {\varphi }^{(\alpha )}}{\partial r_j}dV=0,\forall {\varphi }^{(\alpha )}\in \mathbf{\Phi }.\end{array}& (53)\end{array}$

where all spatial variation of the unknowns has been removed, and only temporal variation remains.

3.2. Non-Conservation Form Fluxes

Consider the PDE

$\begin{array}{cc}\begin{array}{r}{S}_i+A_{ijk}\frac{\partial q_k}{\partial r_j}=0,\end{array}& (54)\end{array}$

and introduce ${\overleftrightarrow{\text{A}}}^{-1}\overleftrightarrow{\text{A}}=\overleftrightarrow{\text{I}}$, which is purely notational, since ${\overleftrightarrow{A}}^{-1}$ may not necessarily exist. Integrating over a test function ϕ^(α), each component i must satisfy

$\begin{array}{cc}\begin{array}{r}{\displaystyle \int }{\varphi }^{(\alpha )}\left(A_{ijk}^{-1}{S}_k+I_{ik}\frac{\partial q_k}{\partial r_j}\right)dV=0,\forall {\varphi }^{(\alpha )}\in \mathbf{\Phi },\end{array}& (55)\end{array}$

$\begin{array}{cc}\begin{array}{l}{\displaystyle \int {\varphi }^{(\alpha )}}{A}_{ijk}^{-1}{S}_{k}dV+I_{ik}({\displaystyle \oint {(q_k{n}_j)}^{*}}{\varphi }^{(\alpha )}dS\\ +{\displaystyle \int q_k}\frac{\partial {\varphi }^{(\alpha )}}{\partial r_j}dV)=0,\forall {\varphi }^{(\alpha )}\in \mathbf{\Phi },\end{array}& (56)\end{array}$

A similar expansion as performed in section 3.1 for spatial varying unknowns can be performed such that

$\begin{array}{cc}\begin{array}{r}{S}_i{(A_{ijk})}^{-1}\approx S_i^{(\beta )}{(A_{ijk})}_{(\alpha )}^{-1}{\varphi }^{(\beta )},\end{array}& (57)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(q_k{n}_j)}^{*}\approx {(q_k{n}_j)}_{(\beta )}^{*}{\varphi }_{(\beta )},\end{array}& (58)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}_k\approx q_k^{(\beta )}{\varphi }^{(\beta )}.\end{array}& (59)\end{array}$

Again this allows the PDE to be written in a semi-discrete form where all unknowns have no spatial variation. Finally, multiplication by $\overleftrightarrow{\text{A}}$ gives

$\begin{array}{cc}\begin{array}{l}{S}_i^{(\beta )}{\displaystyle \int {\varphi }^{(\alpha )}}{\varphi }^{(\beta )}dV+\Delta A_{ijk}^{(\gamma )}({(q_k{n}_j)}_{(\xi )}^{*}{\displaystyle \oint {\varphi }^{(\alpha )}}{\varphi }^{(\beta )}dS\\ +q_k^{(\xi )}{\displaystyle \int {\varphi }^{(\beta )}}\frac{\partial {\varphi }^{(\alpha )}}{\partial r_j}dV)=0,\forall {\varphi }^{(\alpha )}\in \mathbf{\Phi },\end{array}& (60)\end{array}$

$\begin{array}{cc}\begin{array}{r}\Delta ={\delta }^{(\alpha \gamma )}{\delta }^{(\beta \xi )}.\end{array}& (61)\end{array}$

A key difference between Equations (60) and (53) is the fact that A_ijk is evaluated at different locations than the locations where the internal and numerical fluxes are evaluated. Additionally, the characteristics of A_ijk at an interface affect the numerical flux ${(q_k{n}_j)}^{*}$ at that interface. Derivation of an appropriate numerical flux which is stable and convergent can be done using techniques demonstrated in Hesthaven and Warburton [27] and Rhebergen [28].

3.3. Domains for Continuum Kinetic Plasma Models

Kinetic models are useful for plasmas in regions where fluid models break down. However, continuum kinetic models can be computationally expensive, as they solve the kinetic equation for the distribution function, which exists on up to a six-dimensional phase space while fluid and field variables exist on up to only three dimensions. Continuum kinetic models also require high resolution to resolve small phase-space features that can be present in collisionless plasmas, though previous work using higher order elements has shown success in this regard [2, 29]. A Particle-In-Cell (PIC) approach can readily resolve these features and is thus a common choice for applications for collisionless problems. However, PIC models are susceptible to noise and may not always resolve high energy tails of a distribution function as well as a continuum model due to sparsity of particles in low density regions of phase space. The coupling of a continuum kinetic model to a fluid model is also simpler and more accurate than between a PIC and fluid model since the fluid variables are moments which can be directly converted to a distribution function. The conversion to discrete particles is more complex and may require a Monte Carlo procedure.

With this in mind, the continuum kinetic model was chosen in this work. For this model, the phase space domain is discretized into a mixed structured/unstructured mesh, where velocity space is discretized in a structured manner using rectilinear hypercubes (lines, rectangles, and rectangular prisms), and physical space is represented using simplexes (lines, triangles, and tetrahedrons) in an unstructured manner in order to handle potentially complex physical geometries. The resulting phase space element is a tensor product of these two constituent elements. Only distribution functions f_α are discretized over the entire phase space mesh; electric and magnetic fields are expressed only over the co-located physical space mesh.

Numerical quadrature using Gauss-Lobatto rules in the interval [0, 1] for each dimension is used to compute moments of the distribution function to couple the plasma species with Maxwell's equations, as expressed in Equations (3) and (4). The integration is performed over all of velocity space at a physical space location,

$\begin{array}{cc}\begin{array}{ll}{\displaystyle \iiint F}(v_x,v_y,v_z)d{v}_{x}d{v}_{y}d{v}_z\approx \hfill & {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{j=1}^N{\displaystyle \sum _{k=1}^N{w}_i}}}{w}_j{w}_{k}F\left(v_x,v_y,v_z\right)\hfill \\ \hfill & \left|\overleftrightarrow{J}({\xi }_i,{\xi }_j,{\xi }_k)\right|,\hfill \end{array}& (62)\end{array}$

where F is any functional to be integrated, N is number of quadrature locations in each dimension, each with weight w. Though not a required limitation, the hypercubes are implemented as subparametric elements with linearly varying position coordinates and M vertices. Thus, the velocities at quadrature nodes are determined by the mapping between the global coordinate system into a canonical local element space

$\begin{array}{cc}\begin{array}{r}{v}_x=P({\xi }_i,{\xi }_j,{\xi }_k)={\displaystyle \sum _{m=1}^M}{v}_{x_m}{\psi }_m({\xi }_1,{\xi }_2,{\xi }_3),\end{array}& (63)\end{array}$

$\begin{array}{cc}\begin{array}{r}{v}_y=Q({\xi }_i,{\xi }_j,{\xi }_k)={\displaystyle \sum _{m=1}^M}{v}_{y_m}{\psi }_m({\xi }_1,{\xi }_2,{\xi }_3),\end{array}& (64)\end{array}$

$\begin{array}{cc}\begin{array}{r}{v}_z=R({\xi }_i,{\xi }_j,{\xi }_k)={\displaystyle \sum _{m=1}^M}{v}_{z_m}{\psi }_m({\xi }_1,{\xi }_2,{\xi }_3).\end{array}& (65)\end{array}$

The Jacobian defines this mapping

$\begin{array}{cc}\begin{array}{r}{J}_{ij}=\frac{\partial v_i}{\partial {\xi }_j}.\end{array}& (66)\end{array}$

Note that since velocity space is rectilinear, $\overleftrightarrow{\text{J}}$ is diagonal.

3.4. Numerical Fluxes

The numerical fluxes are necessary for computing the terms identified with asterisks in Equations (53) and (60). A further numerical simplification can be found by noting that the lifting matrix ∮ϕ^(α)ϕ^(β)dS is only non-zero if nodes α and β are on the surface ∂K of an element K. Thus, the numerical fluxes ${(F_{ij}{n}_j)}_{(\beta )}^{*}$ and ${(q_k{n}_j)}_{(\beta )}^{*}$ only need to be evaluated at surface nodes β ⊇ λ. The numerical flux is a function of q both on the interior and exterior of ∂K, denoted using the notation q^(λ) and q^(μ) respectively.

An appropriate upwinding numerical flux [30] for the Vlasov equation is given by

$\begin{array}{cc}\begin{array}{r}{(F_{ij}{n}_j)}_{(\lambda )}^{*}=\left[\frac{1}{2}({\widehat{F}}_{ij}^{(\lambda )}+{\widehat{F}}_{ij}^{(\mu )})n_j+\frac{\left|c_j{n}_j\right|}{2}(f_i^{(\lambda )}-f_i^{(\mu )})\right],\end{array}& (67)\end{array}$

where on velocity flux faces ${\widehat{F}}_{ij}=f_i{v}_j$ and c_j = v_j while on acceleration flux faces ${\widehat{F}}_{ij}=f_i{(\text{E}+\text{v}\times \text{B})}_j$ and c_j = (E + v × B)_j.

The numerical fluxes for Maxwell's equations are given [31, 27] by

$\begin{array}{cc}\begin{array}{l}{\left(F_{ij}{n}_j\right)}_{{(\lambda )}_E}^{*}=[\frac{1}{2}\left(-{ϵ}_{ijk}{c}^2\left(B_k^{(\lambda )}+B_k^{(\mu )}\right)\right)n_j\\ +\frac{c{n}_j}{2}{ϵ}_{ijk}{ϵ}_{klm}\left(E_l^{(\lambda )}-E_l^{(\mu )}\right)n_m],\end{array}& (68)\end{array}$

$\begin{array}{cc}\begin{array}{l}{\left(F_{ij}{n}_j\right)}_{{(\lambda )}_B}^{*}=[\frac{1}{2}\left({ϵ}_{ijk}\left(E_k^{(\lambda )}+E_k^{(\mu )}\right)\right)n_j\\ +\frac{c{n}_j}{2}{ϵ}_{ijk}{ϵ}_{klm}\left(B_l^{(\lambda )}-B_l^{(\mu )}\right)n_m],\end{array}& (69)\end{array}$

where $c=({\omega }_p\tau )(\frac{{\delta }_p}{L})$ is the speed of light. ${(F_{ij}{n}_j)}_{{(\lambda )}_{\text{E}}}^{*}$ refers to the numerical flux on the evolution of E (Ampère's Law) where the analytical flux is $F_{ij}=-{ϵ}_{ijk}{c}^2{B}_k$. ${(F_{ij}{n}_j)}_{{(\lambda )}_{\text{B}}}^{*}$ refers to the numerical flux on the evolution of B (Faraday's Law) where the analytical flux is F_ij = ϵ_ijkE_k.

Appropriate forms of numerical fluxes for plasma fluid models [15, 26] have been described previously and are not repeated here. Rusanov [32], HLL [33], or Roe [34] numerical fluxes are used for the MHD and five-moment single species Euler fluid systems.

3.5. Temporal Advance

The time integration of the plasma models is accomplished using implicit and explicit Runge-Kutta methods, including strong stability-preserving Runge Kutta methods [35, 36]. These methods are implemented in a form in which each stage is updated successively over a timestep. The methods solve ordinary differential equations of the form

$\begin{array}{cc}\begin{array}{r}\frac{dq}{dt}={L}(q,t),\end{array}& (70)\end{array}$

where ${L}$ is the right-hand-side spatial discretization of the equations as described in sections 3.1 and 3.2. Multiple methods have been employed, including second and third order strong stability preserving Runge Kutta methods [37]. The second-order total variation bounded Runge Kutta method (Huen's method) [38] is written as

$\begin{array}{cc}\begin{array}{r}{q}^{*}=q^n+\Delta t·{L}(q^n,t_n),\end{array}& (71)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}^{n+1}=\frac{1}{2}{q}^{*}+\frac{1}{2}{q}^n+\frac{1}{2}\Delta t·{L}(q^{*},t_n+\Delta t).\end{array}& (72)\end{array}$

The third-order total variation diminishing Runge Kutta method is written as

$\begin{array}{cc}\begin{array}{r}{q}^{*}=q^n+\Delta t·{L}(q^n,t_n),\end{array}& (73)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}^{**}=\frac{3}{4}{q}^n+\frac{1}{4}{q}^{*}+\frac{1}{4}\Delta t·{L}(q^{*},t_n+\Delta t),\end{array}& (74)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}^{n+1}=\frac{1}{3}{q}^n+\frac{2}{3}{q}^{**}+\frac{2}{3}\Delta t·{L}(q^{**},t_n+\frac{1}{2}\Delta t).\end{array}& (75)\end{array}$

A fourth order method is also used, written as

$\begin{array}{cc}\begin{array}{r}{q}^{*}=q^n+\frac{1}{2}\Delta t·{L}(q^n,t_n),\end{array}& (76)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}^{**}=q^n+\frac{1}{2}\Delta t·{L}(q^{*},t_n+\frac{1}{2}\Delta t),\end{array}& (77)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}^{***}=q^n+\Delta t·{L}(q^{**},t_n+\frac{1}{2}\Delta t),\end{array}& (78)\end{array}$

$\begin{array}{cc}\begin{array}{r}{q}^{n+1}=\frac{1}{3}(-q^n+q^{*}+2q^{**}+q^{***})+\frac{1}{6}\Delta t·{L}(q^{***},t_n+\Delta t).\end{array}& (79)\end{array}$

4. Model Compatibility and Conservation Properties

The MHD and two-fluid plasma models are coupled by defining numerical fluxes at the shared interface. There are several possible ways to define this numerical flux. The first approach is to derive coupling variables based on the two-fluid plasma model using as few assumptions as possible. These coupling variables can then be used to construct a numerical flux for the two-fluid plasma model, which can then be combined to give a conservative numerical flux for the MHD model as well. This approach is examined in section 4.1.

Another approach is to define coupling variables based on all the local assumptions so that these coupling variables are consistent with the underlying local plasma model. This approach complicates the derivation of a numerical flux that guarantees conservation. However, provided the model interface is located in a region where the reduced plasma model, e.g., MHD, is valid, the lack of conservation is insignificant. This approach is examined in section 4.2.

A third approach is to define an intermediate overlapping transition domain where both models influence the plasma properties. Any coupling impedance mismatches are smoothed over several elements. However, this overlap region requires additional computational work since the transition domain must reside entirely within a region where both plasma models are valid. This approach is analyzed in section 4.3.

4.1. Approach 1: Minimize Assumptions

The derivation of the MHD equations provides a clear approach by which the underlying variables are constructed: the continuity equation for each species are summed together, then the momentum equation, etc. The MHD variables q are then defined, Equations (17)–(20), such that the governing equations for the MHD model can be written as

$\begin{array}{cc}\begin{array}{r}\frac{\partial \mathit{q}}{\partial t}+\mathit{R}=0,\end{array}& (80)\end{array}$

where R contains no temporal derivatives. Simplifying assumptions are applied to the operator R. However, if these assumptions are not applied to q, then the definitions of q can be used as the coupling variables definitions. Rewriting Equations (17), (18), (20) gives

$\begin{array}{cc}\begin{array}{r}{\rho }_{\alpha }=\frac{\rho }{1-\frac{Z_{\alpha }{A}_{\beta }}{Z_{\beta }{A}_{\alpha }}},\end{array}& (81)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{p}}_{\alpha }=\frac{A_{\alpha }(Z_{\beta }\mathit{p}-A_{\beta }\mathit{j})}{A_{\alpha }{Z}_{\beta }-A_{\beta }{Z}_{\alpha }},\end{array}& (82)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{u}}_{\alpha }=\mathit{u}-\frac{A_{\beta }\mathit{j}}{Z_{\beta }\rho },\end{array}& (83)\end{array}$

where α, β = i, e. In the case of the single temperature MHD equations, an additional relation is needed to recover the two-fluid energy. This is given by

$\begin{array}{cc}\begin{array}{r}{e}_{\alpha }=\frac{1}{1-\frac{Z_{\alpha }}{Z_{\beta }}}\left(e-\frac{1}{2}{\mathit{p}}_{\beta }·{\mathit{u}}_{\beta }-\frac{1}{2}\frac{Z_{\alpha }}{Z_{\beta }}{\mathit{p}}_{\alpha }·{\mathit{u}}_{\alpha }-\frac{\mathit{B}·\mathit{B}}{2}\right).\end{array}& (84)\end{array}$

It is possible to derive a numerical flux which is guaranteed to be conservative by only computing the numerical flux of the higher fidelity model (in this case, the two-fluid plasma model) by using Equations (81)–(83) to apply “boundary” conditions for the two-fluid plasma model, then the two-temperature MHD numerical flux can be computed using

$\begin{array}{cc}\begin{array}{r}{(\mathit{p}·\mathit{n})}^{*}={({\mathit{p}}_i·\mathit{n})}^{*}+{({\mathit{p}}_e·\mathit{n})}^{*},\end{array}& (85)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(\mathit{H}·\mathit{n})}^{*}={({\mathit{H}}_i·\mathit{n})}^{*}+{({\mathit{H}}_e·\mathit{n})}^{*},\end{array}& (86)\end{array}$

$\begin{array}{cc}\begin{array}{r}{({\mathit{q}}_{\alpha }·\mathit{n})}_{\text{MHD}}^{*}={({\mathit{q}}_{\alpha }·\mathit{n})}_{\text{two}-\text{fluid}}^{*},\end{array}& (87)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(\mathit{E}·\mathit{n})}^{*}={(\mathit{E}·\mathit{n})}^{*},\end{array}& (88)\end{array}$

and for the single temperature MHD numerical fluxes, this is modified slightly to be

$\begin{array}{cc}\begin{array}{r}{(\mathit{p}·\mathit{n})}^{*}={({\mathit{p}}_i·\mathit{n})}^{*}+{({\mathit{p}}_e·\mathit{n})}^{*},\end{array}& (89)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(\overleftrightarrow{\text{H}}·\mathit{n})}^{*}={({\overleftrightarrow{\text{H}}}_i·\mathit{n})}^{*}+{({\overleftrightarrow{\text{H}}}_e·\mathit{n})}^{*},\end{array}& (90)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(\mathit{q}·\mathit{n})}_{\text{MHD}}^{*}=\mathit{B}·{(\overleftrightarrow{\text{E}}·\mathit{n})}_{\text{two}-\text{fluid}}^{*}+{\displaystyle \sum _{\alpha }}{({\mathit{q}}_{\alpha }·\mathit{n})}_{\text{two}-\text{fluid}}^{*},\end{array}& (91)\end{array}$

$\begin{array}{cc}\begin{array}{r}{(\overleftrightarrow{\text{E}}·\mathit{n})}_{\text{MHD}}^{*}={(\overleftrightarrow{\text{E}}·\mathit{n})}_{\text{two}-\text{fluid}}^{*},\end{array}& (92)\end{array}$

where $\overleftrightarrow{\text{H}}$, q, and $\overleftrightarrow{\text{E}}$ are the terms inside the divergence operators. Namely,

$\begin{array}{cc}\begin{array}{r}{\overleftrightarrow{\text{H}}}_{\alpha }={\mathit{p}}_{\alpha }\otimes {\mathit{u}}_{\alpha }+P_{\alpha }\overleftrightarrow{\text{I}},\end{array}& (93)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{q}}_{\alpha }=(e_{\alpha }+P_{\alpha }){\mathit{u}}_{\alpha },\end{array}& (94)\end{array}$

$\begin{array}{cc}\begin{array}{r}\nabla ·\overleftrightarrow{\text{E}}=\nabla \; \times \mathit{E}.\end{array}& (95)\end{array}$

Note that the flux of just magnetic field energy cannot be written in conservation form; it must be implemented using the approach described in section 3.2.

Extending this concept to coupling the continuum kinetic plasma model with the 5N-moment multi-fluid plasma model assumes that the distribution function is well approximated by a Maxwellian distribution function

$\begin{array}{cc}\begin{array}{r}{f}_{\alpha ,\text{two}-\text{fluid}}=n_{\alpha }{\left(\frac{A_{\alpha }}{2\pi T_{\alpha }}\right)}^{\frac{3}{2}}exp\left(\frac{-A_{\alpha }(\mathit{v}-{\mathit{u}}_{\alpha })·(\mathit{v}-{\mathit{u}}_{\alpha })}{2T_{\alpha }}\right).\end{array}& (96)\end{array}$

Once the numerical fluxes for the continuum kinetic plasma model are computed, these can be integrated over the faces in velocity space to construct conservative fluxes for the 5N-moment multi-fluid plasma model.

4.2. Approach 2: Consistent Assumptions

Another approach for coupling the numerical fluxes is to define the coupling variables using consistent assumptions in each region. For the MHD model, the asymptotic approximations are applied to the definitions of the variables.

$\begin{array}{cc}\begin{array}{r}\rho ={\rho }_i,\end{array}& (97)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{p}={\mathit{p}}_i,\end{array}& (98)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{u}=\frac{\mathit{p}}{\rho },\end{array}& (99)\end{array}$

$\begin{array}{cc}\begin{array}{r}e=\frac{P_i+P_e}{\gamma -1}+\frac{\mathit{p}·\mathit{u}}{2}+\frac{\mathit{B}·\mathit{B}}{2},\end{array}& (100)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\rho }_e=-\frac{A_e{Z}_i}{A_i{Z}_e}\rho ,\end{array}& (101)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{p}}_e=\frac{A_e}{A_i{Z}_e}(A_i\mathit{j}-Z_i\mathit{p}),\end{array}& (102)\end{array}$

$\begin{array}{cc}\begin{array}{r}{T}_i=T_e=\frac{e-\frac{\mathit{p}·\mathit{u}}{2}-\frac{\mathit{B}·\mathit{B}}{2}}{n_i+n_e}\end{array}& (103)\end{array}$

While this form is consistent with the definitions of the MHD variables, it complicates formulating a conservative numerical flux due to the definition of e that omits the electron kinetic energy, which must then be absorbed into the ion kinetic energy.

The extension of this procedure for coupling the continuum plasma model the 5N-moment multi-fluid model is again to assume a Maxwellian distribution at the model interface but to apply Equation (96) using the 5N-moment variables on the fluid side to the construct the distribution function while taking moments of the distribution function on the kinetic side, yielding number density, velocity, and temperature for the fluid side. As with MHD to 5N-moment coupling, this provides consistency but does not necessarily provide a conservative flux. This method is referred to as a direct variable translation.

4.3. Approach 3: Transition Region Conservation Properties

For a sharp model interface, the conservation properties are tied directly to the jump of the numerical flux across the shared interface. An interesting question is how does this extend to an overlapping transition region. Suppose the coupling variable set associated with Equations (81)–(83) are used to couple the two systems. The plasma properties can then be defined as a mixture of the two constituent model properties as

$\begin{array}{cc}\begin{array}{r}{\rho }_i=\zeta {\tilde{\rho }}_i+(1-\zeta )\frac{\tilde{\rho }}{1-\frac{Z_i{A}_e}{Z_e{A}_i}},\end{array}& (104)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\rho }_e=\zeta {\tilde{\rho }}_e+(1-\zeta )\frac{\tilde{\rho }}{1-\frac{Z_e{A}_i}{Z_i{A}_e}},\end{array}& (105)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{p}}_i=\zeta {\tilde{\mathit{p}}}_i+(1-\zeta )\left(\frac{A_i(Z_e\tilde{\mathit{p}}-A_e\tilde{\mathit{j}})}{A_i{Z}_e-A_e{Z}_i}\right),\end{array}& (106)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\mathit{p}}_e=\zeta {\tilde{\mathit{p}}}_e+(1-\zeta )\left(\frac{A_e(Z_i\tilde{\mathit{p}}-A_i\tilde{\mathit{j}})}{A_e{Z}_i-A_i{Z}_e}\right),\end{array}& (107)\end{array}$

$\begin{array}{cc}\begin{array}{r}{e}_i=\zeta {ẽ}_i+(1-\zeta ){\widetilde{ϵ}}_i,\end{array}& (108)\end{array}$

$\begin{array}{cc}\begin{array}{r}{e}_e=\zeta {ẽ}_e+(1-\zeta ){\widetilde{ϵ}}_e,\end{array}& (109)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{B}=\zeta \tilde{\mathit{B}}+(1-\zeta )\tilde{{B}},\end{array}& (110)\end{array}$

$\begin{array}{cc}\begin{array}{r}\mathit{E}=\zeta \tilde{\mathit{E}}+(1-\zeta )\tilde{{E}},\end{array}& (111)\end{array}$

where ζ ∈ [0, 1] is a mixing parameter, and the tilde quantities are the plasma parameters predicted by the respective models (note: $\tilde{{E}}$ and $\tilde{{B}}$ are the electric and magnetic fields predicted by the MHD model). The relations in Equations (104)–(111) are guaranteed to be locally instantaneously compatible and conservative. Additionally, it is assumed that evolving each individual model is conservative. However, the resulting scheme is not necessarily conservative. Consider a scenario where no mass enters or leaves the transition domain. Then

$\begin{array}{cc}\begin{array}{r}{\displaystyle \int }{\tilde{\rho }}_i^{(n)}dV={\displaystyle \int }{\rho }_i^{(n)}dV,\end{array}& (112)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\displaystyle \int }\frac{{\tilde{\rho }}^{(n+1)}}{1-\frac{Z_i{A}_e}{Z_e{A}_i}}dV={\displaystyle \int }{\rho }_i^{(n)}dV,\end{array}& (113)\end{array}$

$\begin{array}{cc}\begin{array}{l}{{\displaystyle \int }}^{\text{}}{\rho }_i^{(n+1)}dV={{\displaystyle \int }}^{\text{}}\zeta {\tilde{\rho }}_i^{(n)}+(1-\zeta )\frac{{\tilde{\rho }}^{(n+1)}}{1-\frac{Z_i{A}_e}{Z_e{A}_i}}dV,\\ ={{\displaystyle \int }}^{\text{}}\zeta \left({\tilde{\rho }}_i^{(n)}-\frac{{\tilde{\rho }}^{(n+1)}}{1-\frac{Z_i{A}_e}{Z_e{A}_i}}\right)+\frac{{\tilde{\rho }}^{(n+1)}}{1-\frac{Z_i{A}_e}{Z_e{A}_i}}dV,\\ ={{\displaystyle \int }}^{\text{}}\zeta \left({\tilde{\rho }}_i^{(n)}-\frac{{\tilde{\rho }}^{(n+1)}}{1-\frac{Z_i{A}_e}{Z_e{A}_i}}\right)dV+{{\displaystyle \int }}^{\text{}}{\rho }_i^{(n)}dV.\end{array}& (114)\end{array}$

However, if no mass enters or leaves the transition region then

$\begin{array}{cc}\begin{array}{r}0={\displaystyle \int }\zeta \left({\tilde{\rho }}_i^{(n)}-\frac{{\tilde{\rho }}^{(n+1)}}{1-\frac{Z_i{A}_e}{Z_e{A}_i}}\right)dV.\end{array}& (115)\end{array}$

Consider ζ = x for a transition domain x ∈ [0, 1] and the initial/ending density distributions are

$\begin{array}{cc}\begin{array}{r}{\rho }_i^{(n)}=1,\end{array}& (116)\end{array}$

$\begin{array}{cc}\begin{array}{l}{\tilde{\rho }}_i^{(n+1)}=1+\frac{\sin (2\pi x)}{4},\\ {\tilde{\rho }}^{(n+1)}=\left(1+\frac{\cos (2\pi x)}{4}\right)\left(1-\frac{Z_i{A}_e}{Z_e{A}_i}\right).\end{array}& (117)\end{array}$

Clearly these quantities satisfy the single model conservation properties. However, the integral of Equation (115) is

$\begin{array}{cc}{\displaystyle {{\displaystyle \int }}_0^{1}x}\left(1+\frac{\sin (2\pi x)}{4}-\left(1+\frac{\cos (2\pi x)}{4}\right)\right)dV=\frac{3}{8}-\frac{1}{8\pi }\ne 0& (118)\end{array}$

Thus, the hypothetical scenario does not conserve ion mass. By similar analyses, it can be shown that the other quantities (momentum, energy, etc.) are also not necessarily conserved. This can also be extended to apply different functions ζ. There is only one obvious general function ζ which guarantees conservation inside the transition domain: ζ is spatially constant in the entire transition domain.

If the two-fluid model and MHD model predict the exact same flow behavior then presumably the MHD model is equally valid to use as the two-fluid plasma model for the entire transition region. The deviation from this ideal point, as described by Equation (115), governs the degree to which the resulting model will conserve mass. Similar deviation equations can be derived for the other important conserved quantities momentum and energy. This lack of guaranteed conservation holds regardless of the definition of the coupling variables chosen. However, this property does not appear to be significant since the integral in Equation (115) may be negligibly small in practice. Furthermore, since the lack of conservation can be locally determined, the coupling can be modified to account for any deficit or surplus and correct any lack of conservation.

5. Continuum Kinetic Plasma Model Validation

As described in section 3.3, implementation of the continuum kinetic multi-species plasma model into the DG method and the computational framework in a mixed structured/unstructured formulation requires special considerations. Previous work by Reddell [2] implemented a continuum kinetic plasma model in a fully structured phase-space mesh with the DG implementation without coupling to any fluid model. This work incorporates the mixed phase-space formulation that allows coupling to the unstructured finite element spatial mesh used for the fluid plasma models. The implementation is validated with standard kinetic plasma problems of Landau damping and the two-stream instability.

5.1. Landau Damping

Landau damping in the weak (linear) and strong (nonlinear) regime is validated for the 1D1V case by setting a Maxwellian distribution of f in velocity space with sinusoidal variation in position space, given by the form

$\begin{array}{cc}\begin{array}{r}f(x,v_x)=\frac{n_0}{{(2\pi )}^{\frac{1}{2}}{v}_{\mathrm{\text{th}}}}exp\left(-\frac{v_x^2}{2v_{\mathrm{\text{th}}}^2}\right)(1+\alpha cos(kx)),\end{array}& (119)\end{array}$

where $v_{\mathrm{\text{th}}}=\sqrt{T/A}$. The simulation is performed on a domain of x ∈ [−2π, 2π] and v ∈ [−5, 5] subject to periodic boundary conditions in x and a no flux condition on the velocity boundaries. A consistent electric field is initialized by solving the Poisson equation subject to the periodic boundary condition

$\begin{array}{cc}\begin{array}{l}-\frac{1}{{({\omega }_p\tau )}^2}\left(\frac{L}{{\delta }_p}\right)\frac{{\partial }^2\varphi }{\partial x^2}=Z_{i}n-Z_i{n}_0\\ =Z_i{n}_0(1+\alpha \cos (kx))-Z_i{n}_0,\end{array}& (120)\end{array}$

where a uniformly distributed oppositely charged distribution is set to enforce net charge neutrality. Solving for the electrostatic potential and using the electric field definition gives

$\begin{array}{cc}\begin{array}{r}{E}_x(x)=+Z_i{n}_0{({\omega }_p\tau )}^2\left(\frac{{\delta }_p}{L}\right)\frac{\alpha }{k}sin(kx)\end{array}& (121)\end{array}$

With these initial conditions the Vlasov-Maxwell system is solved using the described DG method. The parameters for this simulation are n₀ = 1, v_th = 1 α = 0.01, k = 0.5, ω_pτ = 1, δ_p/L = 1, which is run for 80 plasma periods using 40 × 40 elements with second order polynomials for the spatial and velocity domains, using the 4th order explicit Runge-Kutta timestepping scheme as written in Equations (76) – (79). The damping rate is measured by plotting an integral of the electric field energy over time until recurrence occurs, as shown in Figure 2. The calculated rate of −0.3106 is compared to the theoretical result of −0.3066 [39, 2]. Increasing the resolution of the simulation leads to the calculated damping rate converging to the theoretical value.

FIGURE 2

Figure 2. Evaluation of electric field energy for the weak (Left) and strong (Right) Landau damping problem. For the weak damping problem, a damping rate of −0.3106 is measured, compared with a theoretical rate of −0.3066. For the strong damping problem, a damping and growth rate of −0.5634 and 0.1742 are measured, compared with −0.562 and 0.168, respectively from published results.

For strong Landau damping, the perturbation is increased to α = 0.5. The damping rate of −0.5634 and growth rate of 0.1742 shown in Figure 2 match well with values from the published literature of −0.562 and 0.168, respectively [39].

5.2. Two-Stream Instability

The two-stream instability is another electrostatic 1D1V test problem with a theoretical growth rate that can be used to validate the implementation of the continuum kinetic plasma model. Two counter-streaming plasma beams are initialized as

$\begin{array}{cc}\begin{array}{l}f(x,v_x)=\frac{1}{2}\frac{n_0}{{\left(2\pi \right)}^{\frac{1}{2}}{v}_{\text{th}}}[\exp \left(-\frac{1}{2}\frac{{\left(v_x-v^{\prime }\right)}^2}{v_{\text{th}}^2}\right)\\ +\exp \left(-\frac{1}{2}\frac{{\left(v_x+v^{\prime }\right)}^2}{v_{\text{th}}^2}\right)]\left(1+\alpha \cos \left(kx\right)\right).\end{array}& (122)\end{array}$

With this form of f, the same initial electric field is used, Equation (121). The separation of each beam from the axis is v′ = π/2 with a perturbation α = 0.01. The beams have v_th = 0.1. The simulation is run using 80 × 160 elements with second order polynomials and 4th order explicit Runge-Kutta timestepping using Equations (76)–(79), on x ∈ [−2π, 2π/k] and v ∈ [−10, 10] with n₀ = 1, α = 0.01, k = 0.5, ω_pτ = 1, δ_p/L = 1. With these parameters, $k{v}^{\prime }/{\omega }_p=0.7854$, yielding a theoretical growth rate of 0.4952. The calculated growth rate of this simulation, as shown in Figure 3 is 0.4879. Figure 4 shows the evolution of the distribution function in phase space. Growth of the instability occurs until a characteristic time of 20 after which there is the expected saturation.

FIGURE 3

Figure 3. Evaluation of electric field energy for the two-stream instability problem. A growth rate of 0.4879 is measured, compared with the theoretical rate of 0.4952.

FIGURE 4

Figure 4. Distribution function f∈[0, 2] for the two-stream instability problem for v_th with separation v′ = π/2 at different times. x ∈ [0, 4π] is the bottom axis while v_x ∈ [−5, 5] is the left axis, though the simulation was performed for v_x ∈ [−10, 10]. Growth of the instability occurs at the expected rate until saturation at $({\omega }_p\tau )\tilde{\tau }=20$. See Figure 3.

6. Numerical Tests With Direct Variable Translation

6.1. 5N-Moment to MHD

A 2D planar plasma opening switch is implemented to test the direct variable translation method for coupling two-fluid and MHD. Five different domain subdivisions are tested:

(a) Uniform MHD: A single MHD plasma model is used in the entire domain.

(b) Uniform two-fluid: A two-fluid plasma model is used in the entire domain.

(c) Mixed.25: A two-fluid model is used in the region x ∈ [0, 0.25], and an MHD model is used in the region x ∈ [0.25, 1].

(d) Mixed.5: A two-fluid model plasma model is used in the region x ∈ [0, 0.5], and an MHD model is used in the region x ∈ [0.5, 1].

(e) Adaptive: Uses the subdivision Mixed.25 for t ∈ [0, 0.8], then switches to Mixed.5 for t ≥ 0.8.

The temporal domain is discretized using Heun's method [38]. A high density slab of plasma is initialized with a notch perturbation in the middle. This configuration is shown in Figure 5. The plasma is accelerated to the right by imposing B_z = 1 on the left boundary. This drives a current in the plasma, producing a Lorentz force. The testing parameters used are:

$\begin{array}{cc}\begin{array}{r}{\omega }_p\tau =2000\end{array}& (123)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{{\delta }_p}{L}=10^{-2}\end{array}& (124)\end{array}$

$\begin{array}{cc}\begin{array}{r}{Z}_i=-Z_e=1\end{array}& (125)\end{array}$

$\begin{array}{cc}\begin{array}{r}{A}_i=100A_e=1\end{array}& (126)\end{array}$

$\begin{array}{cc}\begin{array}{r}\gamma =\frac{5}{3}\end{array}& (127)\end{array}$

$\begin{array}{cc}\begin{array}{r}{P}_i=P_e=\frac{1}{2}\end{array}& (128)\end{array}$

FIGURE 5

Figure 5. Planar plasma switch domain and initial plasma density. Dashed lines denote the various decomposition points used. Mesh consists of 400 × 100 rectangles subdivided into two triangle elements.

Figure 6 shows the plasma density at t = 1, shortly after the adaptive method has remapped the domain. An instability begins developing around where the initial perturbation was in all cases except for the uniform MHD case. At t = 2, Figure 7 shows that the adaptive case is capable of matching the fine-scale instability structure of the Uniform two-fluid and Mixed.5 cases. However, the instabilities which developed on the two-fluid side of the Mixed.25 case balloons on the MHD side and loses the fine-scale structure. This is especially prevalent at t = 2.5 in Figure 8. At this point the adaptive case is still capable of maintaining fine-scale instability structures on the trailing edge of the bulk plasma; however, the leading edge does not have similar instability structures as the Uniform two-fluid and Mixed.5 cases. It is suspected that this is due to the leading edge being on the MHD side initially before the adaptive method moves the model interface ahead of the leading edge.

FIGURE 6

Figure 6. Planar plasma switch density at t = 1. The adaptive method has re-partitioned the domain variables. An instability at the density perturbation starts developing in all models except for the Uniform MHD case.

FIGURE 7

Figure 7. Planar plasma switch density at t = 2. The bulk plasma has transitioned past x = 0.25 for all models. Some instabilities which have developed in the two-fluid region of the Mixed.25 case balloon on the MHD side, unlike the Uniform two-fluid or Mixed.5 cases. The adaptive method maintains similar instability structures as the Uniform two-fluid and Mixed.5 cases.

FIGURE 8

Figure 8. Planar plasma switch density at t = 2.5. Ballooning of the instability in the Mixed.25 case is significant. The adaptive case maintains most of the fine-scale structure found in the Uniform two-fluid and Mixed.5 cases on the trailing edge, however the leading edge has noticeable differences.

The test cases have been run with a few different fixed timesteps using identical initial conditions. No new numerical instabilities which result from coupling the models when the coupling assumptions are well matched was noticed. Stability appears to be limited by the limits of the constituent models used. Additional stabilization mechanisms such as artificial viscosity or limiters were not needed for this particular problem.

6.2. Continuum Kinetic to 5N-Moment

6.2.1. Double Rarefaction Waves

Exploration of the direct variable translation method for coupling described in section 4.2 between a neutral fluid with no heat flux or viscosity in 1D and the kinetic model in 1D1V is performed. A Riemann problem producing double rarefaction waves is initialized where

$\begin{array}{cc}\begin{array}{l}\left(\rho ,v_x,p\right)=\{\begin{array}{l}\left(1.0,-0.2,0.4\right)\text{for }x<0.5\hfill \\ \left(1.0,+0.2,0.4\right)\text{for }x\ge 0.5\hfill \end{array}\hfill \end{array}& (129)\end{array}$

on a domain for x ∈ [0, 1] up to t = 0.15 with the ratio of specific heats, γ = 3. A problem is run in which the domain is split with the 5N-moment model on the left and the kinetic model on the right, with v ∈ [−10, 10]. No limiter is used in this coupling problem, but may be necessary for larger jumps and the resulting stronger waves. On the kinetic side, these fluid variables are converted to a Maxwellian distribution. Dirichlet conditions are imposed on fluid variables while copy-out is imposed on the distribution function in phase space, acceptable up to this time since waves do not reach the boundary. Two cases are run, one where the model interface is at the jump (x = 0.5) and another where the interface is at x = 2/3, allowing for the waves to propagate in the fluid region for some time before coming into contact with the kinetic region. In each case, a mesh of 128 × 320 elements (64 physical elements on either side in both cases) are used with first order polynomials on the physical and phase space with second order Runge-Kutta timestepping using Equations (71) – (72). Also n₀ = 1, ω_pτ = 1, and δ_p/L = 1. Rusanov numerical fluxes are used for the fluid side of the domain.

To drive the distribution function toward a Maxwellian to be consistent with the fluid model , a BGK collision operator is introduced for the collision term in the kinetic equation,

$\begin{array}{cc}\begin{array}{r}{\left(\frac{\partial f_{\alpha }}{\partial t}\right)}_{coll}=-({\nu }_p\tau )(f-f_M)\end{array}& (130)\end{array}$

where f_M is the Maxwellian distribution calculated from the zeroth, first, and second moments of f. The normalized collisional relaxation time is set to ν_pτ = 1000. Results are given in Figure 9. The solution is shown on the physical-space domain with fluid variables (ρ,v_x,p) shown directly on the left. The same variables are shown in the right, which are the result of moments taken of the distribution function as given in Equations (3), (4), and (9).

FIGURE 9

Figure 9. Double rarefaction waves problem with 5N-moment model on left of model interface and continuum kinetic model with BGK operator on right of model interface at t = 0.15 on x ∈ [0, 1] subject to direct variable translation matching. The top case has the model interface at x = 0.5 while the bottom case has the model interface at x = 0.667. Black dashed lines are the analytic solution to the fluid Riemann problem [40].

6.2.2. 1D1V Sheath

The direct variable translation method for coupling is tested on a non-neutral plasma using a simulation of a sheath. The Vlasov-Maxwell system is used to solve this problem as described in Cagas et al. [41] with a domain of x ∈ [−128λ_D, 128λ_D] with 256 elements. The outer 32 elements from both left and right are solved using the Vlasov-Maxwell system for ion and electron distributions while the interior elements are solved with the 5N-moment model for fluid ions and electrons. A realistic mass ratio is used, such that A_i = 1, A_e = 1/1836, with charges Z_i = +1, Z_e = −1. Phase space is spans v_i ∈ [−6v_thi, +6v_thi] and v_i ∈ [−6v_the, +6v_the] for ions and electrons, respectively, using 48 square elements for each physical element and as with the double rarefaction waves problem, first order polynomials and second order timestepping is used.

The normalization is set by defining a reference length (L), reference number density (n₀), and reference temperature (T₀), from which δ_p/L and ω_pτ are calculated. The choice is of a typical number density ($n_0=10^{19}{\mathrm{\text{m}}}^{-3}$), T₀ = 1eV, and $L=L_D={(\frac{{ϵ}_0{T}_0}{n_0{e}^2})}^{\frac{1}{2}}=2.35\mu \mathrm{\text{m}}$. From these parameters, using the pressure normalization, $p_0=n_0{T}_0=\frac{B_0^2}{{\mu }_0}$ leads to $v_0=v_A=\frac{B_0}{\sqrt{{\mu }_0{m}_p{n}_0}}=n\sqrt{\frac{T_0}{m_p}}=9.79\times 10^3\mathrm{\text{m/s}}$ where reference $m_p=1.67\times 10^{-27}\mathrm{\text{kg}}$ (proton mass). Then $\tau =\frac{L}{v_0}=2.4\times 10^{-10}\mathrm{\text{s}}$ leading to the normalization constants

$\begin{array}{cc}\begin{array}{r}{\omega }_p\tau =\sqrt{\frac{e^2{n}_0}{{ϵ}_0{m}_p}}\tau =1\end{array}& (131)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{{\delta }_p}{L}=\frac{1}{{\omega }_c\tau }=\frac{m_p}{e{B}_0}\frac{1}{\tau }=3.265\times 10^{-5}\end{array}& (132)\end{array}$

using fundamental constants ${ϵ}_0=8.85\times 10^{-12}\mathrm{\text{F/m}}$, ${\mu }_0=4\pi \times 10^{-7}\mathrm{\text{H/m}}$, and e = 1.602 × 10⁻¹⁹C.

The initial conditions for ions and electrons are for normalized n_i = n_e = 1, v_xi = v_xe = 0, T_i = T_e = 10. In the fluid domains, these variables are set directly, while in the kinetic domain they are converted to equivalent Maxwellian distributions, using γ = 3, appropriate for one spatial degree of freedom. The simulations are run to t = 20/ω_pe. The boundary condition of the sheath in physical space is that of a conducting wall on Maxwell's equations and are set as outflow conditions on the distribution function in phase space. Additionally, locally-dependent BGK operators are applied on both ion and electron species to replenish electron tails to achieve steady-state, as shown in Cagas et al. [41], where the form of the BGK operator is

$\begin{array}{cc}\begin{array}{r}{\left(\frac{\partial f_{\alpha }}{\partial t}\right)}_{coll}=-({\nu }_p\tau ){\nu }_{\alpha }(f-f_M)\end{array}& (133)\end{array}$

with

$\begin{array}{cc}\begin{array}{r}{\nu }_{\alpha }=\frac{n_{\alpha }}{A_{\alpha }^{\frac{1}{2}}{T}_{\alpha }^{\frac{3}{2}}}ln\text{Λ}.\end{array}& (134)\end{array}$

For this problem, ν_pτ was chosen to be 2, with ln Λ = 10.

Results comparing fluid moments are shown in Figure 10 (ions) and Figure 11 (electrons). For comparison, a full kinetic simulation is also performed (dashed black lines). The phase space for both ions and electrons, respectively, are shown in Figures 12, 13. The solutions for the domain-hybridized plasma model and for the full kinetic model match closely until the electrons deviate significantly from a Maxwellian distribution at model boundaries. However, as long as the 5N-moment plasma model is physically valid at the model interface, i.e., the distribution functions remain Maxwellian, the domain-hybridized plasma model produces accurate results with reduced computational effort. To simulate this particular problem longer in time, one might choose to move the model boundary farther into the middle portion of the domain where the electron distribution remains closer to a Maxwellian, or to solve the electrons with a full kinetic model while solving the ions with the hybrid approach, which would still reduce the computational effort. The boundary may also be moved dynamically during a simulation using a metric determining the departure of the distribution from a Maxwellian. An example is given by Reddell [2], by defining a spatially local parameter

$\begin{array}{cc}\begin{array}{ll}{\chi }_{\alpha }=\hfill & \frac{{\displaystyle {\int }_v\left|f_{\alpha }-f_{M_{\alpha }}\right|}dv}{n_{\alpha }}\hfill \end{array}& (135)\end{array}$

FIGURE 10

Figure 10. Ion sheath density, velocity, and pressure moments at t = 20/ω_pe with direct variable translation. The middle domain is solved using the 5N-moment two-fluid plasma model while the left and right domains are solved using the continuum kinetic two-species plasma model. Simulation results from applying the continuum kinetic model for the entire domain are plotted in black dashed lines for comparison.

FIGURE 11

Figure 11. Electron sheath density, velocity, and pressure moments at t = 20/ω_pe with direct variable translation. The middle domain is solved using the 5N-moment two-fluid plasma model while the left and right domains are solved using the continuum kinetic two-species plasma model. Simulation results from applying the continuum kinetic model for the entire domain are plotted in black dashed lines for comparison. Some solution mismatch is seen due to departure of the electron distribution function away from a Maxwellian where the model interface assumes validity of the fluid approximation.

FIGURE 12

Figure 12. Contours of the ion distribution function at t = 20/ω_pe solved using the domain-hybridized plasma model. The left and right subdomains x ∈ [−128, −96] and x ∈ [96, 128] are modeled with the continuum kinetic two-species plasma model, and the middle subdomain x ∈ [−96, 96] is modeled with the 5N-moment two-fluid plasma model.

FIGURE 13

Figure 13. Contours of the electron distribution function at t = 20/ω_pe solved using the domain-hybridized plasma model. The left and right subdomains x ∈ [−128, −96] and x ∈ [96, 128] are modeled with the continuum kinetic two-species plasma model, and the middle subdomain x ∈ [−96, 96] is modeled with the 5N-moment two-fluid plasma model.

where f_Mα is the Maxwellian determined by moments of f_α.

7. Comparison of Conservative Numerical Fluxes and Direct Variable Translation

One question of interest is how different are the numerical fluxes of different models assuming both models are individually valid locally. Computing independent numerical fluxes at an interface is the simplest method, and if it is sufficient then there may not be any reason to use coupling numerical fluxes which guarantee conservation.

To investigate this issue, a 1% jump discontinuity centered about x = 0 in density is initialized in an isothermal plasma. The subdomain is subdivided such that $-\frac{1}{5}<x<\frac{1}{5}$ is simulated using two-fluid, while the remainder of the domain is simulated using MHD.

Two different test regimes are initialized: one where the MHD assumptions are reasonable, and the second where the MHD assumption of quasi-neutrality breaks down. In the latter case, Langmuir waves propagate out from the jump discontinuity, and comparisons are made on the plasma behavior as this wave interacts with the transition boundaries. Heun's method is used for the temporal discretization.

The testing parameters for “good” matching are

$\begin{array}{cc}\begin{array}{r}{A}_i=100\end{array}& (136)\end{array}$

$\begin{array}{cc}\begin{array}{r}{A}_e=1\end{array}& (137)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{{\delta }_p}{L}=\frac{1}{10}\end{array}& (138)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\omega }_p\tau =1000\end{array}& (139)\end{array}$

The testing parameters for “poor” matching are

$\begin{array}{cc}\begin{array}{r}{A}_i=100\end{array}& (140)\end{array}$

$\begin{array}{cc}\begin{array}{r}{A}_e=1\end{array}& (141)\end{array}$

$\begin{array}{cc}\begin{array}{r}\frac{{\delta }_p}{L}=1\end{array}& (142)\end{array}$

$\begin{array}{cc}\begin{array}{r}{\omega }_p\tau =100\end{array}& (143)\end{array}$

Figure 14 demonstrate that when the MHD assumptions are valid, the direct variable translation method conserves mass and energy. The conservative numerical flux is able to preserve mass and energy several orders of magnitude better.

FIGURE 14

Figure 14. Conservation of mass and energy testing. (Top Left) Mass conservation when MHD assumptions somewhat valid. (Top Right) Mass conservation when MHD assumptions invalid. (Bottom Left) Energy conservation when MHD assumptions somewhat valid. (Bottom Right) Energy conservation when MHD assumptions invalid. Validity of the MHD assumption produces an improvement in the conservation properties for direct variable translation.

For the small magnitude electrostatic shock test problem, Figure 15 shows that the parameter mismatch results in oscillations at the model interface. Using the direct variable translation approach eliminates the oscillations. Both methods couple across the domain boundary between the models without generating deleterious effects.

FIGURE 15

Figure 15. Total plasma mass density of the small electrostatic shock test at t = 1.5 with invalid MHD assumptions. Both the conservative numerical flux coupling and direct variable translation methods are able to maintain reasonable coupling behavior, and the impedance mismatch shows up as oscillations at the model interface. The direct variable translation method seems to average over some of these oscillations to produce a visually smoother coupling.

As the parameter regime better satisfies the MHD assumptions, both approaches produce similar results. In the context of a practical adaptive method, this needs to be satisfied anyways in order to accurately capture the correct physics. This implies that the direct variable translation method should work better in practice on the account that it is significantly simpler to implement.

8. Conclusion

Spatial coupling of different plasma models is facilitated by derivations of the governing equations that use a consistent formulation and normalization, which allows direct translation between models of higher and lower physical fidelity. The DG method provides a uniform structure for spatial discretization to represent the continuum plasma models used for this development. Using the same spatial discretization simplifies coupling the plasma models. Care must be exercised when treating the fluxes in conservative and non-conservative forms that appear in the governing equations.

The continuum kinetic plasma model is also implemented with the DG method by exploiting the rectilinear hypercubes used to represent velocity space. Validation test cases for the continuum kinetic plasma model in 1D1V agree with theoretical and published numerical results for weak and strong Landau damping and for the two-stream instability. Moments of the distribution function provide a means to couple the continuum kinetic plasma model to fluid plasma models, such as the multi-fluid plasma model.

The spatial coupling of different plasma fluid models has been implemented and successfully tested on 1D and 2D test problems. Model coupling accuracies are investigated for different approaches. Devising a consistent normalization for all of the plasma models used facilitates more accurate coupling and does not require modification of the plasma models. The coupling method does not appear to introduce any stricter stability requirements than are already required for the constituent models. It is demonstrated that direct coupling of model variables does not guarantee conservation of key physical properties, but the lack of exact conservation produces no significant detrimental effect. Conservation improves as the models approach their regions of validity, which indicates that an adaptive plasma model that smoothly transitions between levels of physical fidelity can produce accurate simulations at lower computational costs.

Initial work on adaptively changing the domain model partitioning is demonstrated on a 2D planar plasma switch, and demonstrates that the variable translation method can be used to remap the domain, and the resulting simulations can maintain better physical fidelity. Future work on conditions for automatically determining where re-mapping should take place is being investigated.

Data Availability Statement

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

AH and ID implemented and tested the various plasma models. AH investigated the plasma model coupling and its conservation properties. US designed the approach and methodology for a physics-based-adaptive kinetic-fluid plasma model and developed consistent representations for the plasma models and normalizations. AH, ID, and US are primary architects of the WARPXM plasma simulation framework, analyzed and interpreted the numerical and physical results, and wrote the manuscript.

Funding

This material is based upon work supported by the Air Force Office of Scientific Research under award numbers FA9550-15-1-0271 and FA9550-14-1-0317.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors acknowledge helpful discussion from other members of the University of Washington's Computational Plasma Dynamics Laboratory: D. W. Crews, A. R. Johansen, E. T. Meier, S. Taheri, Y. Takagaki, and W.R. Thomas.

References

1. Birn J, Drake JF, Shay MA, Rogers BN, Denton RE, Hesse M, et al. Geospace Environmental Modeling (GEM) magnetic reconnection challenge. J Geophys Res. (2001) 106:3715. doi: 10.1029/1999JA001005

CrossRef Full Text | Google Scholar

2. Reddell N. A Kinetic Vlasov Model for Plasma Simulation Using Discontinuous Galerkin Method on Many-Core Architectures. University of Washington. Seattle, WA (2016).

Google Scholar

3. Sakai O, Tachibana K. Plasmas as metamaterials: a review. Plasma Sour Sci Technol. (2012) 21:1–18. doi: 10.1088/0963-0252/21/1/013001

CrossRef Full Text | Google Scholar

4. Wang B, Capelli MA. A plasma photonic crystal bandgap device. Appl Phys Lett. (2016) 108:161101. doi: 10.1063/1.4946805

CrossRef Full Text | Google Scholar

5. Wang MY, Miley GH. Particle orbits in field-reversed mirrors. Nuclear Fusion (1979) 19:39–49.

Google Scholar

6. Tuszewski M. Field reversed configurations. Nuclear Fusion (1988) 28:2033.

Google Scholar

7. Steinhauer LC. Review of field-reversed configurations. Phys Plasmas (2011) 18:1–38. doi: 10.1063/1.3613680

CrossRef Full Text | Google Scholar

8. Haines MG, Lebedev SV, Chittenden JP, Beg FN, Bland SN, Dangor AE. The past, present, and future of Z pinches. Phys Plasmas (2000) 7:1672–80. doi: 10.1063/1.874047

CrossRef Full Text | Google Scholar

9. Shumlak U, Nelson BA, Claveau EL, Forbes EG, Golingo RP, Hughes MC, et al. Increasing plasma parameters using sheared flow stabilization of a Z-pinch. Phys Plasmas (2017) 24:1–10. doi: 10.1063/1.4977468

CrossRef Full Text | Google Scholar

10. Tummel K, Higginson D, Link A, Schmidt A, McLean H, Shumlak U, et al. 2D kinetic particle in cell simulations of a flow-shear stabilized Z-pinch. Bull Am Physical Soc. (2017) 62:UP11.109.

Google Scholar

11. Loverich J, Shumlak U. Nonlinear full two-fluid study of m = 0 sausage instabilities in an axisymmetric Z pinch. Phys Plasmas (2006) 13:082310. doi: 10.1063/1.2220009

CrossRef Full Text | Google Scholar

12. Shumlak U, Lilly R, Reddell N, Sousa E, Srinivasan B. Advanced physics calculations using a multi-fluid plasma model. Comput Phys Commun. (2011) 182:1767–70. doi: 10.1016/j.cpc.2010.12.048

CrossRef Full Text | Google Scholar

13. Cockburn B, Shu CW. The Runge-Kutta discontinuous Galerkin method for conservation laws. J Comput Phys. (1998) 141:199–224.

Google Scholar

14. Loverich J, Shumlak U. A discontinuous Galerkin method for the full two-fluid plasma model. Comput Phys Commun. (2005) 169:251–5. doi: 10.1016/j.cpc.2005.03.058

CrossRef Full Text | Google Scholar

15. Loverich J, Hakim A, Shumlak U. A discontinuous Galerkin method for ideal two-fluid plasma equations. Commun Comput Phys. (2011) 9:240–68. doi: 10.4208/cicp.250509.210610a

CrossRef Full Text | Google Scholar

16. Srinivasan B, Hakim A, Shumlak U. Numerical methods for two-fluid dispersive fast MHD phenomena. Commun Comput Phys. (2011) 10:183–215. doi: 10.4208/cicp.230909.020910a

CrossRef Full Text | Google Scholar

17. Srinivasan B, Shumlak U. Analytical and computational study of the ideal full two-fluid plasma model and asymptotic approximations for Hall-magnetohydrodynamics. Phys Plasmas (2011) 18:092113. doi: 10.1063/1.3640811

CrossRef Full Text | Google Scholar

18. Kolobov VI, Arslanbekov RR, Aristov VV, Frolova AA, Zabelok SA. Unified solver for rarefied and continuum flows with adaptive mesh and algorithm refinement. J Comput Phys. (2007) 223:589–608. doi: 10.1016/j.jcp.2006.09.021

CrossRef Full Text | Google Scholar

19. Kolobov VI, Arslanbekov RR. Towards adaptive kinetic-fluid simulations of weakly ionized plasmas. J Comput Phys. (2012) 231:839–69. doi: 10.1016/j.jcp.2011.05.036

CrossRef Full Text | Google Scholar

20. Tóth G, Jia X, Markidis S, Peng IB, Chen Y, Daldorff LKS, et al. Extended magnetohydrodynamics with embedded particle-in-cell simulation of Ganymede's magnetosphere. J Geophys Res. Space Phys. (2016) 121:1273–93. doi: 10.1002/2015JA021997

CrossRef Full Text | Google Scholar

21. Alouani-Bibi F, Opher M, Alexashov D, Izmodenov V, Toth G. Kinetic versus multi-fluid approach for interstellar neutrals in the heliosphere: exploration of the interstellar magnetic field effects. Astrophys J. (2011) 734:1–13. doi: 10.1088/0004-637X/734/1/45

CrossRef Full Text | Google Scholar

22. Winske D, Yin L, Omidi N, Karimabadi H, Quest K. Hybrid Simulation Codes: Past, Present and Future—A Tutorial. Berlin; Heidelberg: Springer (2003). p. 136–65.

Google Scholar

23. Sousa EM, Shumlak U. A blended continuous-discontinuous finite element method for solving the multi-fluid plasma model. J Comput Phys. (2016) 326:56–75. doi: 10.1016/j.jcp.2016.08.044

CrossRef Full Text | Google Scholar

24. Miller ST, Shumlak U. A multi-species 13-moment model for moderately collisional plasmas. Phys Plasmas (2016) 23:082303. doi: 10.1063/1.4960041

CrossRef Full Text | Google Scholar

25. Meier ET, Shumlak U. A general nonlinear fluid model for reacting plasma-neutral mixtures. Phys Plasmas (2012) 19:1–11. doi: 10.1063/1.4736975

CrossRef Full Text | Google Scholar

26. Shumlak U, Loverich J. Approximate Riemann solver for the two-fluid plasma model. J Comput Phys. (2003) 187:620–38. doi: 10.1016/S0021-9991(03)00151-7

CrossRef Full Text | Google Scholar

27. Hesthaven JS, Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. New York, NY: Springer Science+Business Media, LLC (2008).

Google Scholar

28. Rhebergen S. Discontinuous Galerkin Finite Element Methods for (Non)Conservative Partial Differential Equations. University of Twente (2010).

29. Vogman G. Fourth-Order Conservative Vlasov-Maxwell Solver for Cartesian and Cylindrical Phase Space Coordinates. University of California, Berkeley, CA (2016).

Google Scholar

30. Cheng Y, Christlieb AJ, Zhong X. Energy-conserving discontinuous Galerkin methods for the Vlasov-Maxwell System. J Comput Phys. (2014) 279:145–73. doi: 10.1016/j.jcp.2014.08.041

CrossRef Full Text | Google Scholar

31. Cheng Y, Gamba IM, Li F, Morrison PJ. Discontinuous galerkin methods for the Vlasov-Maxwell equations. SIAM J Numer Anal. (2013) 52:1–37. doi: 10.1137/130915091

CrossRef Full Text | Google Scholar

32. Rusanov VV. The calculation of the interaction of non-stationary shock waves with barriers. USSR Comput Math Math Phys. (1962) 1:304–20.

Google Scholar

33. Harten A, Lax PD, Leer BV. On upstream differencing and godunov-type schemes for hyperbolic conservation laws. SIAM Rev. (1983) 25:35–61.

Google Scholar

34. Roe PL. Approximate riemann solvers, parameter vectors, and difference schemes. J Comput Phys. (1981) 43:357–72.

Google Scholar

35. Gottlieb S, Shu CW. Total variation diminishing runge-kutta schemes. Math Comput. (1998) 267:73–85.

Google Scholar

36. Gottlieb S, Shu CW, Tadmor E. Strong stability-preserving high-order time discretization methods. SIAM Rev. (2001) 43:89–112. doi: 10.1137/S003614450036757X

CrossRef Full Text | Google Scholar

37. Shu CW, Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys. (1988) 77:439–71.

Google Scholar

38. Hirsch C. Numerical Computatoin of Internal and External Flows. vol. 1. 2nd ed. Oxford: Butterworth-Heinemann (2007).

Google Scholar

39. Cheng CZ, Knorr G. The integration of the vlasov equation in configuration space. J Comput Phys. (1976) 22:330–51.

Google Scholar

40. Leveque RJ. Finite Volume Methods for Hyperbolic Problems. Cambridge, UK: Cambridge University Press (2002).

Google Scholar

41. Cagas P, Hakim A, Juno J, Srinivasan B. Continuum Kinetic and multi-fluid simulations of classical sheaths. Phys Plasmas (2017) 24:1–11. doi: 10.1063/1.4976544

CrossRef Full Text | Google Scholar