Ten-moment fluid model for low-temperature magnetized plasmas

Kuldinow, Derek Amur; Yamashita, Yusuke; Hara, Kentaro

doi:10.1063/5.0240993

In this paper, a one-dimensional 10-moment multi-fluid plasma model is developed and applied to low-temperature magnetized plasmas. The 10-moment model solves for six anisotropic pressure terms, in addition to density and three components of fluid momentum, which allows the model to capture finite kinetic effects. The results are benchmarked against a 5-moment model, which assumes that the gas constituents follow a Maxwellian velocity distribution function (VDF), and a particle-in-cell Monte Carlo collision model that allows for arbitrary non-Maxwellian VDFs. The models are compared in a one-dimensional, low-temperature, partially magnetized plasma test case. The 10-moment results accurately reproduce the anisotropic temperature profile in low-temperature magnetized plasmas, where shear gradients exist due to the $E \times B$ drift. We discuss the mechanisms by which the anisotropic pressure can be generated in low-temperature magnetized plasmas. In addition, the importance of a self-consistent heat flux closure to the 10-moment model is studied, showing consistency with other models only when the assumptions of the underlying model are met. The 10-moment model allows for study of electron inertia effects and non-Maxwellian VDFs without the need for kinetic methods that are more computationally expensive.

I. INTRODUCTION

Low-temperature plasmas (LTPs) play an important role in natural phenomena, such as space weather,^1–3 and engineering applications including space propulsion^4–6 and semiconductor manufacturing.⁷ LTPs are a rich area of study because the particle dynamics include a wide range of spatial and temporal scales arising from the mass difference between different species, e.g., ions, electrons, and neutral gases. The addition of a magnetic field can greatly affect the dynamics within LTPs due to particle gyromotion and by introducing a direction of anisotropy. The atomic-scale interactions between particles, including short-range collisions and long-range field interactions, directly cascade up to device-scale behavior and performance. Therefore, the basis for any macroscopic modeling of LTPs must be faithful to the underlying microscopic dynamics.

A Hall-effect thruster (HET) is an electric propulsion (EP) device that employs a plasma in an annular channel with applied crossed electric and magnetic fields. The crossed fields trap the electrons and generate an $E \times B$ drift in the azimuthal direction, and the high-energy trapped electrons can then collide with and ionize a neutral gas. High electron resistivity increases the electric field which then accelerates the ions, generating thrust. The electron trapping is key to electron heating, the ionization of neutral atoms which sustain the plasma, and the instantiation of the strong accelerating field. An HET discharge chamber is a prudent test problem for the 10-moment model because it contains crossed electric and magnetic fields, multiple relevant species, a number of reactions and contains physics that occur over a wide range of spatial and temporal scales. In particular, anisotropy can be important in such a plasma as the magnetic field modifies the electron transport in various directions.⁸ Electron transport is one of the most important phenomena in this type of crossed-field discharges and other plasma phenomena. Through shear-induced transport, and a generalized diamagnetic drift, anisotropic pressure components can greatly affect the electron transport.

Fluid moment models are an attractive option for modeling plasmas because macroscopic system-scale behavior can be described without needing to track individual particle trajectories and may thus have a significant computational cost advantage over kinetic models. However, the key challenge of fluid models is the closure problem. That is, the evolution of the calculated moments (e.g., density, bulk velocity, and temperature) depends on higher-order moments; consequently, some ansatz must be made to determine the higher-order moments from the lower-order moments. One common assumption for fluid models of plasma and gas dynamics is that the gaseous particles follow a near-Maxwellian velocity distribution function (VDF). All gaseous systems relax to a Maxwellian distribution in the limit of collisional equilibrium. Thus, the Maxwellian assumption is valid when the timescale of intermolecular collisions is fast compared to any other process. Such a system would be fully described by the spatial profiles of density, n, bulk velocity, $u$ ⁠, and isotropic pressure, p, i.e., a 5-moment model. Another common approach in LTPs is to use the two-term spherical harmonic expansion of the Boltzmann equation, leading to the local field and local mean energy approximations.⁹ However, particularly in the presence of a magnetic field, non-local effects and collisional processes may lead to anisotropic transport and deviation from local equilibrium; thus, the 5-moment assumption may be invalid.

Experimental measurements have found evidence of anisotropic electron VDFs, particularly near phase boundaries (e.g., walls, channel exit) where shear terms are large.^10,11 Anisotropic pressures (and temperatures) can affect the dynamics of a system in two main ways, which are not captured in equilibrium models. First, anisotropy modifies particle transport properties; in general, in addition to changing pressure gradients, anisotropic temperatures can lead to anisotropic resistivities and mobilities.^8,12 Furthermore, a non-Maxwellian VDF can modify reaction rates. For instance, consider a VDF which has a higher temperature in the x direction than the other two directions. If the thermal energy is near the activation energy of a reaction, an equilibrium model may underpredict the reaction rate caused by the particles with higher x velocity.

The fluid equations are derived by taking velocity moments of the kinetic transport equation (KTE), such as the Boltzmann equation, and fluid models are categorized by which moments they determine and the type of closure. The typical 5-moment fluid model uses the zeroth, first, and contracted second moments of the KTE. Higher-order moment models (HOMMs) are an attractive approach for modeling fluids in the transition between equilibrium fluid and kinetic regimes by accounting for non-Maxwellian VDFs, while still remaining a manageable set of first-order partial differential equations (PDEs). Previous works have used 10-moment models,^3,13,14 13-moment models,^15,16 and 14-moment models,^17,18 which consider a maximum-entropy ansatz and Pearson-IV distribution for closure. Intuitively, having more moments, i.e., more information, should only increase the fidelity of the simulation; however, as described in previous work, there are practical advantages to maintaining a hyperbolic set of PDEs. In Ref. 14, a closure based on the H-theorem was developed for the 10-moment model, which retains a hyperbolic wave structure and performed favorably against 5-moment models and previously developed 10-moment closures when compared to kinetic simulations for gas dynamics problems.

In this study, the 10-moment model is extended to charged species and applied to the partially ionized, partially magnetized plasmas that are observed in an HET discharge channel, and benchmarked against a 5-moment fluid model and a kinetic model. In Sec. II, we discuss the governing equations for our fluid and particle simulations. In Sec. III, we introduce the numerical methods, the boundary conditions, and the solver for the electric field. Section IV describes the details of the simplified HET model. The results of our simulations are presented and analyzed in Sec. V.

II. GOVERNING EQUATIONS

The governing equation to describe a distribution of particles of species s is the kinetic transport equation, given by

\frac{\partial f_{s}}{\partial t} + v_{i} \frac{\partial f_{s}}{\partial x_{i}} + a_{i} \frac{\partial f_{s}}{\partial v_{i}} = C_{s},

(1)

where

f_{s} (x_{i}, v_{i}, t)

is the velocity distribution function of species s as a function of position

x_{i}

⁠, velocity

v_{i}

⁠, and time t,

a_{i}

is the acceleration acting on the particles,

C_{s} (f) = \sum_{s^{'}} C_{s s^{'}} (f_{s}, f_{s^{'}})

is the collision operator of species s under collisions with species

s^{'}

⁠, and subscripts

i = {x, y, z}

indicate Einstein summation notation.

A. The fluid approximation

Fluid models are derived by solving for certain moments of Eq. (1). In the limit of retaining infinite moments, the function f is resolved exactly (i.e., moment problem); however, in practice, usually only a limited set of moments are considered to resolve f to sufficient accuracy.

1. The 10-moment equations

The 10-moment equations are derived by taking zeroth, first, and second velocity moments of Eq. , which yields the following transport laws:

\frac{\partial}{\partial t} ρ + \frac{\partial}{\partial x_{i}} Γ_{i} = C^{(0)},

(2)

\frac{\partial}{\partial t} Γ_{i} + \frac{\partial}{\partial x_{j}} P_{i j} = \frac{q}{m} (ρ E_{i} + ϵ_{ijk} Γ_{j} B_{k}) + C_{i}^{(1)},

(3)

\begin{matrix} \frac{\partial}{\partial t} P_{i j} + \frac{\partial}{\partial x_{k}} Q_{ijk} = \frac{q}{m} [Γ_{i} E_{j} + Γ_{j} E_{i} + (ϵ_{i ℓ k} P_{j ℓ} + ϵ_{j ℓ k} P_{i ℓ}) B_{k}] + C_{i j}^{(2)}, \end{matrix}

(4)

where we have dropped the species subscript s and taken only the acceleration due to electromagnetic forces:

a_{i} = \frac{q}{m} (E_{i} + ϵ_{ijk} v_{j} B_{k})

⁠, where q is the species charge, m is the species mass,

ϵ_{ijk}

is the Levi-Civita symbol,

E_{i}

is the electric field,

B_{i}

is the magnetic field, and

C^{(0)}

⁠,

C_{i}^{(1)}

⁠, and

C_{i j}^{(2)}

are the zeroth, first, and second moments of the collision operator, respectively. The mass density, momentum density, total stress, and energy flux are defined as

ρ = m \int f d^{3} v, Γ_{i} = m \int v_{i} f d^{3} v,

(5)

P_{i j} = m \int v_{i} v_{j} f d^{3} v, Q_{ijk} = m \int v_{i} v_{j} v_{k} f d^{3} v .

(6)

The set of variables

U = {ρ, Γ_{i}, P_{i j}, Q_{ijk}}

can be called conservative variables because, without collisions and external forces, they are conserved within the domain. Under elastic collisions only,

{ρ, Γ_{i}, P_{i i}}

are conserved quantities, corresponding to conservation of mass, momentum, and energy. The conservative variables are in contrast to the so-called primitive variables

Π = {ρ, u_{i}, p_{i j}, q_{ijk}}

⁠, where the bulk velocity, pressure tensor, and heat flux tensor can be written, respectively, as

u_{i} = Γ_{i} / ρ, p_{i j} = m \int w_{i} w_{j} f d^{3} v = P_{i j} - ρ u_{i} u_{j},

(7)

\begin{matrix} q_{ijk} = m \int w_{i} w_{j} w_{k} f d^{3} v \\ = Q_{ijk} - p_{i j} u_{k} - p_{i k} u_{j} - p_{j k} u_{i} - ρ u_{i} u_{j} u_{k} . \end{matrix}

(8)

Here,

w_{i} = v_{i} - u_{i}

is the peculiar velocity. Then, through a change of variables, the governing equations [Eqs. ] can be written using the primitive variables as

\frac{D ρ}{D t} + ρ \frac{\partial u_{i}}{\partial x_{i}} = c^{(0)},

(9)

\frac{D u_{i}}{D t} + \frac{1}{ρ} \frac{\partial p_{i j}}{\partial x_{j}} = \frac{q}{m} (E_{i} + ϵ_{ijk} u_{j} B_{k}) + c_{i}^{(1)},

(10)

\begin{matrix} \frac{D p_{i j}}{D t} + p_{i j} \frac{\partial u_{k}}{\partial x_{k}} + p_{i k} \frac{\partial u_{j}}{\partial x_{k}} + p_{j k} \frac{\partial u_{i}}{\partial x_{k}} + \frac{\partial q_{ijk}}{\partial x_{k}} \\ = \frac{q}{m} (ϵ_{i ℓ k} p_{j ℓ} + ϵ_{j ℓ k} p_{i ℓ}) B_{k} + c_{i j}^{(2)}, \end{matrix}

(11)

where

\frac{D}{D t} = \frac{\partial}{\partial t} + u_{i} \frac{\partial}{\partial x_{i}}

denotes the material derivative and the primitive sources can be derived to be

c^{(0)} = C^{(0)},

(12)

c_{i}^{(1)} = \frac{1}{ρ} (C_{i}^{(1)} - u_{i} C^{(0)}),

(13)

c_{i j}^{(2)} = C_{i j}^{(2)} - u_{i} C_{j}^{(1)} - u_{j} C_{i}^{(1)} + u_{i} u_{j} C^{(0)} .

(14)

It is to be noted that the conductive heat flux

q_{ijk}

is not solved for by the model and some ansatz must be made to approximate it from the solved quantities. Our closure model for heat flux will be described in Sec. II D.

2. The 5-moment equations

The 5-moment equations can be derived from the 10-moment equations under the assumption that

p_{i j} = p δ_{i j}

⁠, where

p = p_{i i} / 3

⁠, i.e., that the pressure tensor can be reduced to a scalar, isotropic pressure in which case the primitive equations can be written

\frac{D ρ}{D t} + ρ \frac{\partial u_{i}}{\partial x_{i}} = c^{(0)},

(15)

\frac{D u_{i}}{D t} + \frac{1}{ρ} \frac{\partial p}{\partial x_{i}} = \frac{q}{m} (E_{i} + ϵ_{ijk} u_{j} B_{k}) + c_{i}^{(1)},

(16)

\frac{3}{2} \frac{D p}{D t} + \frac{5}{2} p \frac{\partial u_{k}}{\partial x_{k}} + \frac{\partial q_{k}}{\partial x_{k}} = c^{(2)},

(17)

where the heat flux vector

q_{k} = q_{iik} / 2

and the pressure source term

c^{(2)} = c_{i i}^{(2)} / 2

⁠. The conservative equations are equivalent in form to those in Eqs. , but only one second-moment equation, i.e., the contracted second-moment equation, must be solved for.

B. The particle-in-cell method

In particle-based kinetic methods, like particle-in-cell (PIC), instead of evolving the VDF, f, directly or taking its moments, the VDF is approximated as a sum of discrete computational particles, each represented as a delta function

f ≃ \sum_{p} f_{p} = \sum_{p} w_{p} δ (x_{p, i}, v_{p, i}, t),

(18)

where

f_{p}

is a kernel function for each particle, assumed in our case to be a Dirac delta function in phase space,

δ (x_{i}, v_{i}, t)

⁠,

x_{p, i}

and

v_{p, i}

are the position and velocity of the

p^{th}

particle, and

w_{p}

is the weight of each particle.

The time evolution of each computational particle can be determined by splitting the left-hand side of Eq. and integrating to yield

\frac{d x_{p, i}}{d t} - v_{p, i} (t) = 0,

(19)

\frac{d v_{p, i}}{d t} - a_{i} (x_{p, i}, v_{p, i}, t) = 0.

(20)

For coupling with the field equations, it is often necessary to determine the number density and current density of the system of particles on a grid for Gauss's law and Ampère's law, respectively. If the time rate change of the induced magnetic field due to current density is small, the problem can be considered electrostatic, which is further explained in Sec. II E.

C. Collisions

For a complete description of a system of particles, the setup must be endowed with a collision operator to describe the action of particles on one another. For the right-hand side of Eq. , the Boltzmann operator, which is valid for binary, elastic collisions, is given by

C_{B} (f_{s}, f_{s^{'}}) = \iint (f_{s}^{'} f_{s^{'}}^{'} - f_{s} f_{s^{'}}) g σ d Ω d^{3} v_{s^{'}},

(21)

where

f^{'}

is the post-collision distribution function, subscripts s and

s^{'}

denote two different species, g =

| v_{s^{'}} - v_{s} |

is the magnitude of the relative velocity,

σ

is the collision cross section, and

d Ω

is the differential solid angle of the deflected particle, including polar and azimuthal angles. In general, this collision integral is difficult to evaluate because it is nonlinear in the distribution functions and analytic solutions do not exist. Numerical integration is also computationally costly because the integrand is a function of seven variables

{x_{i}, v_{i}, t}

⁠, and the integration is over five dimensions.¹⁹ Under certain assumptions, the Boltzmann operator can be approximately integrated to obtain the fluid source terms (C and c) using the theory by Chapman and Cowling.²⁰ Recent work in plasma simulations has evaluated these integrals of the Boltzmann operator in the calculation of collision frequencies of the various interactions extant in plasmas.^21,22

We employ an approximation to the Boltzmann operator, the Bhatnagar–Gross–Krook (BGK) operator:²³

C_{s s^{'}} (f_{s}, f_{s'}) = \sum_{r, s s^{'}} ν_{r, s s^{'}} (f_{r, s s^{'}}^{post} - f_{s}),

(22)

where the sum is over collision types r between species s and

s^{'}

⁠,

ν_{r, s s^{'}}

is the characteristic collision frequency, and

f_{r, s s^{'}}^{post}

is the post-collision VDF that

f_{s}

locally relaxes toward through collisions. The descriptions of

ν_{r, s s^{'}}

and

f_{r, s s^{'}}^{post}

are detailed in Sec. IV C. For fluid models, the collisional sources

C^{(n)}

are determined by taking the moments of Eq. .

For particle methods, collisions on particle p can be performed stochastically using a Monte Carlo collision (MCC) model, by updating the particle velocity at a rate $ν_{r, s s^{'}}$ for each collision type. The magnitude of the post-collision velocity $| v_{p} |^{post}$ is determined by conservation of energy, and its direction is assumed to be isotropic; a more accurate description of the post-collision particle VDF is reserved for future work. Equation (22) is consistent with the fluid model in terms of conservation of mass, momentum and energy. However, there may be significant discrepancies if the particles are significantly non-Maxwellian, which is not captured in the estimate of $f^{post}$ ⁠.

D. Heat flux closure model

The final piece required for the 10-moment model is the prescription of the heat flux,

q_{ijk}

⁠, in Eq. . The Chapman–Enskog closure method as described in textbooks²⁴ assumes that the VDF can be approximated as

f = f^{0} + f^{1}

⁠, where

f^{0}

is an equilibrium distribution, which is assumed to be Maxwellian based on the H-theorem, and

f^{1}

is a perturbation from the equilibrium distribution. Considering the BGK collision operator,

C = ν (f^{post} - f)

⁠, where

ν

is the collision frequency, through algebraic manipulations, the perturbation VDF can be written as^14,25

\begin{matrix} f^{1} = - \frac{1}{ν} [(\frac{ρ w_{i} w_{j}}{p} - \frac{p_{i j}}{p} \frac{ρ w^{2}}{3 p} + \frac{p_{i j}}{p} - δ_{i j}) \frac{\partial u_{i}}{\partial x_{j}} \\ + w_{j} (\frac{ρ w^{2}}{2 p} - \frac{5}{2}) \frac{\partial}{\partial x_{j}} ln (\frac{p}{ρ})] f^{0} . \end{matrix}

(23)

By taking the third moment of

f^{0} + f^{1}

⁠, one finds that

q_{ijk} = - \frac{p}{ν} [δ_{i j} \frac{\partial}{\partial x_{k}} (\frac{p}{ρ}) + δ_{k i} \frac{\partial}{\partial x_{j}} (\frac{p}{ρ}) + δ_{j k} \frac{\partial}{\partial x_{i}} (\frac{p}{ρ})] .

(24)

Note that the 5-moment heat flux, often known as Fick's law, is recovered by taking the trace

q_{iik}

⁠. Without loss of generality, the contracted heat flux can be obtained as

q_{x} = \frac{1}{2} (q_{xxx} + q_{yyx} + q_{zzx}) = - κ \frac{\partial T}{\partial x},

(25)

where the thermal conductivity can be obtained as

κ = 5 p k_{B} / (2 m ν)

when considering a BGK collision operator, where

k_{B}

is the Boltzmann constant.

However, there have been empirical models developed for the 5-moment

κ

which would be useful to generalize to the 10-moment system. Note that in Eq. , for instance,

q_{xxx} : q_{xyy} : q_{xzz} = 3 : 1 : 1

⁠. This

3 : 1 : 1

decomposition of the vector

q_{k}

to the on-diagonal elements of the tensor

q_{ijk}

is valid for a wide range of perturbation functions

f^{1}

⁠, as detailed in Appendix A. Thus, where in the 5-moment model

q_{k} = - κ_{k} \frac{\partial T}{\partial x_{k}}

(with no applied summation), the heat flux tensor for the 10-moment model can be written as

q_{ijk} = - \frac{2}{5} κ_{k} (δ_{i j} \frac{\partial T}{\partial x_{k}} + δ_{k i} \frac{\partial T}{\partial x_{j}} + δ_{j k} \frac{\partial T}{\partial x_{i}}),

(26)

where, for this equation, we do not imply summation over k, and we have allowed for a potentially anisotropic heat conductivity

κ_{k}

⁠.

E. Field equations

In this study, it is assumed that any induced magnetic fields due to currents in the fluid are small compared to applied magnetic fields from an external source. Thus, we employ the electrostatic assumption and only consider the induced electric field. Thus, Maxwell's equations can be reduced to Poisson's equation for electrostatic potential

\frac{\partial^{2}}{\partial x_{i} \partial x_{i}} ϕ = - \sum_{s} q_{s} n_{s},

(27)

where

n_{s} = ρ_{s} / m_{s}

is the species number density, and the electric field is then given by

E_{i} = - \frac{\partial ϕ}{\partial x_{i}} .

(28)

III. NUMERICAL METHODS

The continuum differential equations for the evolution of the plasma must now be put into a form that is solvable numerically. We employ a finite volume method, where the domain is discretized into a finite number of cells and we solve for the spatially averaged plasma quantities within each cell.

A. Timestepping

To numerically solve the time evolution equations, e.g., Eqs. (2)–(4), (19), and (20), they must be spatially and temporally discretized and numerically integrated forward in time. The fluid equations, i.e., Eqs. (2)–(4) and (9)–(11) for 10-moment and Eqs. (15)–(17) for 5-moment, are numerically explicitly integrated forward in time. By determining the conservative fluxes at the interfaces between volume elements, we can ensure that mass, momentum, and energy are conserved. Thus, the left-hand side of the fluid equations is treated conservatively.

Although working with the conservative form of the equations is numerically convenient for the fluxes, there is no preferred set of variables for all of the source terms. Consider the electromagnetic source terms in the energy equations: the conservative form in Eqs. (2)–(4) couples the electric field to the energy through bulk acceleration of the fluid momentum; however, when written in a primitive form as Eqs. (9)–(11), it is clear that the magnetic field only acts to rotate the pressure tensor in velocity space and the electric field does not directly affect the pressure tensor. In this way, a primitive formulation may mitigate numerical errors arising from the coupling between internal and bulk kinetic energy.⁴ However, certain inelastic collisions may cause a net drag as well as energy loss. In such a case, the conservative equations provide a convenient way to handle momentum and energy losses independently, without needing to predetermine its decomposition into internal and bulk kinetic energy losses.

In general, the primitive and conservative equations can be written as

\frac{\partial U}{\partial t} + \frac{\partial F_{i}}{\partial x_{i}} = S,

(29)

\frac{\partial Π}{\partial t} + A_{i} \frac{\partial Π_{i}}{\partial x_{i}} = Σ,

(30)

respectively, where

F_{i}

is the vector of conservative fluxes in the i direction,

S

is the vector of conservative sources,

A_{i}

is the Jacobian matrix for the fluxes transformed under a change-of-basis to the primitive variables, and

Σ

is the vector of primitive sources. Each of the flux terms is described in Ref. 14.

Second-order Strang splitting is employed,²⁶ wherein the temporal update due to the fluxes and sources are staggered by a half time step. Within a single time step, the variables are first pushed by a half time step using only the sources, then a full time step using only the fluxes, and then finally another half time step with only the sources. The numerical procedure of updating the primitive and conservative variables from time step n to

n + 1

is shown below:

\begin{matrix} [Π^{n}, U^{n}] \to_{Sources}^{Δ t / 2} [Π^{'}, U^{'}] \to_{Fluxes}^{Δ t} [Π^{″}, U^{″}] \to_{Sources}^{Δ t / 2} [Π^{n + 1}, U^{n + 1}], \end{matrix}

(31)

where

Δ t

is the time step and the variable vectors

Π, U

with the same superscript denote a consistent set of primitive and conservative variables, respectively. The first

Δ t / 2

source update includes updating the fluid variables of certain processes in a primitive or conservative manner; the update

[Π^{n}, U^{n}] \to_{Sources}^{Δ t / 2} [Π^{'}, U^{'}]

is performed as

{\begin{matrix} {(\frac{Π_{ℓ}^{⋆} - Π_{ℓ}^{n}}{Δ t / 2})}_{Prim . Sources} = Σ_{p} (Π_{ℓ}^{n}), \\ {(\frac{U_{ℓ}^{'} - U_{ℓ}^{⋆}}{Δ t / 2})}_{Cons . Sources} = S_{c} (U_{ℓ}^{n}), \end{matrix}

(32)

where

Σ_{p}

are the sources which are to be updated primitively, and

S_{c}

are the sources which are to be updated conservatively; the conservative

U_{ℓ}^{⋆}

is calculated from the primitive

Π_{ℓ}^{⋆}

for cell index

ℓ

⁠. Finally,

Π_{ℓ}^{'}

is calculated from

U_{ℓ}^{'}

⁠. For our simulations, electromagnetic interactions and self-collisions occur in the species frame of reference are handled primitively, while inter-species collisions are handled conservatively since they must respect a total energy balance.

Next, the second step in Eq. is the flux update. The discretized form of the one-dimensional flux update for variables at cell index

ℓ

due to both primitive and conservative sources can be written as

{(\frac{U_{ℓ}^{″} - U_{ℓ}^{'}}{Δ t})}_{Fluxes} + \frac{F_{x, ℓ + 1 / 2} - F_{x, ℓ - 1 / 2}}{Δ x} = 0,

(33)

where

Δ x

is the cell size and the flux function is given by

\begin{matrix} F_{x, ℓ + 1 / 2} = F_{x}^{+} [U^{L} (x_{ℓ + 1 / 2})] + F_{x}^{-} [U^{R} (x_{ℓ + 1 / 2})] \\ + F_{x, c} (U^{L}, U^{R}, \nabla U) . \end{matrix}

(34)

Here,

F_{x}^{\pm}

are the rightward and leftward moving Steger–Warming fluxes,²⁷

U^{L / R}

are the states reconstructed using monotonic upwind scheme for conservation laws (MUSCL) with the van Leer harmonic limiter²⁸ at the left

(L)

and right

(R)

sides of the interface at position

x_{ℓ + 1 / 2}

⁠, and

F_{x, c}

is the viscous flux due to closure, including gradients of

q_{ijk}

for 10-moment and gradients of

q_{k}

for 5-moment. The gradient values,

\nabla U

⁠, at cell interfaces are calculated using a second-order central differencing.

For increased temporal accuracy and to reduce spurious oscillations with larger timesteps, each time step is taken using the third-order strong stability-preserving Runge-Kutta (SSP-RK3) method.⁴ All of the species are updated at each substep of RK3.

The time step size

Δ t

is taken to satisfy the Courant–Friedrichs–Lewy (CFL) condition. In the present 10-moment model, it was observed that non-oscillatory (stable) results are obtained for the CFL condition of

C = \frac{| λ_{i} |_{max} Δ t}{Δ x_{i}} \leq 0.1,

(35)

where

| λ_{i} |_{max}

is the fastest wave speed in the i direction,

| λ_{i} |_{max} = | u_{i} | + \sqrt{3 p_{i i} / ρ}

⁠, where the index i is not summed over in this equation. In addition, the sources impose stability constraints on the time step

Δ t ν_{r, s s^{'}} < 1, Δ t ω_{p, e} < 1, Δ t ω_{c, e} < 1,

(36)

where

ω_{p, e} = \sqrt{n e^{2} / m_{e} ϵ_{0}}

is the electron plasma frequency and

ω_{c, e} = | q_{e} B / m_{e} |

is the electron cyclotron frequency. The first constraint ensures that we are not over-relaxing the VDF using the BGK operator as defined in Eq. , the second ensures that we resolve plasma oscillations, and the third ensures that we resolve cyclotron oscillations. The fixed time step size is chosen at the start of the simulation as the maximum that will satisfy each of the above constraints. All simulations are run to steady state, but the different models and data points may use different time step sizes.

B. Kinetic flux boundary condition for the 10-moment model

Boundary conditions play an important role in plasma simulations. While it is somewhat straightforward to determine how individual particles behave at the wall, the proper handling of the boundary condition for a fluid model is not trivial.

For each of the fluid models, we employ the so-called kinetic flux boundary conditions. The 5-moment description is detailed in Ref. 4. For benchmarking the 10-moment model with 5-moment and kinetic models, here we consider the kinetic boundary conditions for the 10-moment model.

The kinetic flux boundary condition is predicated on estimating the particle VDF at the wall from the fluid variables and determining the corresponding flux to the wall based on the integrating the VDF.^4,25 In this way, the kinetic flux formulation assigns a (e.g., mass, momentum, and energy) flux consistent with kinetic methods. Also, it provides a continuous and differentiable flux as a function of the macroscopic quantities (e.g.,

ρ, u_{i}, p_{i j}

⁠). For this study, we set values in two ghost cells outside of the computational domain, from which the fluid variables at the boundaries are determined using MUSCL reconstruction. If we consider a boundary interface in the

+ x

direction, i.e., to the right of the plasma, the rightward moving flux of a flow variable Q to the interface can be written as

F_{x, Q}^{+} = \int_{0}^{\infty} d v_{x} \int_{- \infty}^{\infty} d v_{y} \int_{- \infty}^{\infty} d v_{z} v_{x} Q {\tilde{f}}_{w},

(37)

where

{\tilde{f}}_{w}

is the VDF of the plasma at the interface. Likewise, the leftward moving flux,

F_{x, Q}^{-}

has bounds of integration for

v_{x}

from

- \infty

to 0 so that the net flux at the boundary interface is the sum of the two fluxes:

F^{net} = F^{+} + F^{-}

⁠. In the 10-moment framework,

{\tilde{f}}_{w}

may be determined from the reconstructed primitive variables at the boundary

{\tilde{f}}_{w} = \frac{ρ}{m} {[\frac{ρ^{3}}{{(2 π)}^{3} | ℘ |}]}^{1 / 2} exp (- \frac{ρ}{2} w_{i} ℘_{i j}^{- 1} w_{j}),

(38)

where

| ℘ |

and

℘_{i j}^{- 1}

are the determinant and the inverse, respectively of the pressure tensor

℘_{i j}

⁠. Given the 5-moment properties, the maximum-entropy VDF that matches the given moments is a Maxwellian VDF. When 10-moment properties are specified, the maximum-entropy VDF that matches the 10-moment variables in Eqs. and is the Gaussian VDF in Eq. . The explicit forms of the flux functions are shown in Appendix B.

C. Poisson's equation

For the fluid models, the number density is explicitly captured; however, for particle models, it is necessary to prescribe a way to gather particle information at grid points for calculation of a density. In this study, we use linear interpolation whence the density at a grid point

ℓ

is calculated as

n_{ℓ} = \sum_{p} w_{p} S_{1} (\frac{x_{p} - x_{ℓ}}{Δ x}),

(39)

where

S_{1} (x) = max (0, 1 - x)

is the unit linear shape function.

Given the value of the number density of each species at the cell-centers,

n_{ℓ}

⁠, we must solve the 1-dimensional Poisson's equation with Dirichlet boundary conditions to calculate the potential structure between an anode and cathode. The discretized form of the equation can be written as

\begin{matrix} \frac{ϕ_{0} + ϕ_{1}}{2} = ϕ_{A}, \\ \frac{ϕ_{ℓ - 1} - 2 ϕ_{ℓ} + ϕ_{ℓ + 1}}{{(Δ x)}^{2}} = - \sum_{s} q_{s} n_{s, ℓ}, \\ \frac{ϕ_{N} + ϕ_{N + 1}}{2} = ϕ_{B}, \end{matrix}

(40)

where L is the domain length,

ϕ_{ℓ}

is the potential at the

ℓ^{th}

cell center (⁠

ℓ = 1, \dots, N

⁠) in the domain, and

ϕ_{0}

and

ϕ_{N + 1}

are the left and right ghost cells introduced to satisfy Dirichlet boundary conditions:

ϕ (x = 0) = ϕ_{A}

and

ϕ (x = L) = ϕ_{B}

⁠. Equation is solved using the Thomas tridiagonal matrix solver.

IV. THE HALL-EFFECT THRUSTER MODEL

The computational domain in this study consists of a one-dimensional (1D) representation of a Hall-effect thruster as shown in Fig. 1, following Refs. 4 and 29. The left boundary is at $x = 0$ cm, and stands in for an annular anode of inner radius 3.45 cm and outer radius 5 cm, set to a potential of $+ 300$ V relative to the cathode. The channel length $L_{ch}$ is 2.5 cm.

FIG. 1.

View large Download slide

Simplified one-dimensional HET configuration used for numerical simulations, following Ref. 4.

Within the domain, we retain the full 10-moment equations for the electrons. Ions are treated using an isothermal 5-moment model, with (isotropic) ion temperature

k_{B} T_{i} = 0.025

eV, and neutrals are modeled using a linear advection equation for simplicity. The neutral density is updated as

\frac{\partial n_{n}}{\partial t} + u_{n} \frac{\partial n_{n}}{\partial x} = - {\dot{n}}_{iz},

(41)

where

{\dot{n}}_{i z}

is the ionization rate calculated consistently with the other species. A first-order upwind method is used for the flux term. The neutral velocity and temperature are assumed to be constant:

u_{n} = 270

m/s and

T_{n} = 0.025

eV. In reality, the neutral velocity may change due to selective ionization, collisions, and expansions,³⁰ but this simple model is chosen for benchmarking purposes.

A. Boundary conditions

At the anode, it is assumed that a constant flux of neutral xenon gas (5 mg/s mass flow rate) at 270 m/s flows into the system.^30,31 Xenon is assumed to be the only atomic species present and dominated by the ground-state neutral and singly charged ion. In addition, the ions that collide with the anode are assumed to recombine.

For electrons and ions, the ghost cells at the anode, which affect the reconstructed variables at the boundary and thus the calculated fluxes, have their density and bulk velocity set to zero, with a Neumann boundary condition for the pressure tensor. At the cathode, the incoming flux is assumed to consist of isotropic stationary thermal electrons at 3 eV. The density of the incoming electrons is calculated so that the plasma is quasineutral in the last cell at all times. Given the fluxes into the last cell from within the domain,

Γ_{N - \frac{1}{2}}

⁠, where N is the total number of cells, and assuming that there is no ion influx through the cathode, the incoming electron flux can be given as

\begin{matrix} Γ_{e, N + \frac{1}{2}}^{-} = Γ_{e, N - \frac{1}{2}} - Γ_{e, N + \frac{1}{2}}^{+} + \frac{m_{e}}{m_{i}} (Γ_{i, N + \frac{1}{2}} - Γ_{i, N - \frac{1}{2}}) + \frac{m_{e} Δ x}{Δ t} (n_{e, N} - n_{i, N}), \end{matrix}

(42)

where

n_{s, N}

is the number density of species s in the last cell. As

{u_{i}, T_{i j}}_{cathode}

are assumed to be constant, we can use the kinetic flux formula, i.e., Eq. , to determine the

n_{cathode}

required to achieve the desired electron mass flux as shown in Eq. . Using

{n, u_{i}, T_{i j}}_{cathode}

⁠, the momentum fluxes and fluxes for the anisotropic pressures can be evaluated using the kinetic flux formula, i.e., Eqs. .

B. Magnetic field

For the present cross field plasma benchmarking, a static, applied magnetic field is assumed in the radial (z) direction

B_{z} (x) = {\begin{array}{l} B_{max} exp [- \frac{1}{2} {(\frac{x - L_{ch}}{0.011})}^{2}], for x < L_{ch} \\ B_{max} exp [- \frac{1}{2} {(\frac{x - L_{ch}}{0.018})}^{2}], for x \geq L_{ch}, \end{array}

(43)

where the maximum magnetic field, achieved at the channel exit, is

B_{max} = 150

G.

C. Collisions

1. Collision types

We assume the following collisional processes for the electron dynamics for our benchmarking test: electron-neutral (en), electron-ion (ei), electron–electron (ee), electron-wall (ew). In addition, the momentum transfer that results from the so-called anomalous cross field transport is also included. It is further assumed that the effect of collisions on the heavy species (ie, ne, in, ni) is small and does not affect the dynamics.

The electron-neutral collisions can be broken down into two components: elastic and inelastic. The momentum transfer collision frequency due to electron-neutral collisions can be written as

ν_{en} = ν_{elas} + ν_{inel},

(44)

where

ν_{elas}

and

ν_{inel}

are the elastic and inelastic contributions, respectively, to the total electron-neutral collision frequency

ν_{en}

⁠. The elastic component is calculated using the model provided in Ref. 32

ν_{elas} = k_{0} n_{n},

(45)

where

k_{0} = 2.5 \times 10^{- 13} m^{3} / s

⁠. The inelastic collision frequency is given as

ν_{inel} = \sum_{c} k_{c} n_{n} = k_{eff} n_{n},

(46)

where the sum is over the various inelastic collision types c is condensed to an effective rate

k_{eff}

⁠.

The reaction rates for inelastic collisions, including ionization and excitation, are generated from BOLSIG+⁹ using data for electron–xenon cross section from the Biagi database³³ the local mean energy approximation of the two-term Boltzmann equation. Note that the Hayashi database from Ref. 34 has also been tested but showed negligible differences to the results that used the Biagi database. The inelastic reaction rate coefficients k are tabulated as a function of electron effective isotropic temperature,

k_{B} T_{eff} = (k_{B} T_{i i} + m u_{i} u_{i}) / 3

⁠, assuming that the neutral velocities are much smaller than the electron velocities and that the electron VDFs are Maxwellian. In Sec. VI, we describe how this approach may be generalized for non-Maxwellian electrons. While the collisions contribute as drag in the momentum equation, they also contribute as energy loss mechanisms. An energy loss rate coefficient

Q_{eff}

is calculated as

Q_{eff} = \sum_{c} Δ ϵ_{c} k_{c} n_{n},

(47)

where

Δ ϵ_{c}

is the excitation or ionization energy of the inelastic reaction c. Given the large electron-ion mass ratio (⁠

m_{Xe} / m_{e} ≃ 2.39 \times 10^{5}

⁠), the energy transferred from the electrons to the ions is small, so we do not consider any elastic energy transfer between the electrons and the heavy species.

The electron-ion momentum-transfer collision frequency is calculated using the model provided in Ref. 35, which assumes subsonic Maxwellian electrons and ions

ν_{ei} = α n_{e} λ {(k_{B} T_{e})}^{- 3 / 2},

(48)

where the coefficient

α = 2.9 \times 10^{- 12} [m^{3} {eV}^{3 / 2} s^{- 1}]

and the value of the Coulomb logarithm is taken to be

λ ≃ 23 - 0.5 ln [n_{e} {(k_{B} T_{e})}^{- 3} [{cm}^{3} {eV}^{- 3}]]

⁠.

In the 10-moment model, anisotropic pressure relaxes to an isotropic pressure through self-collisions for electrons, i.e., electron–electron collisions. While it is assumed to be effectively infinite for 5-moment fluid models, proper treatment of the intra-species collisions is important for the 10-moment model. While the physical electron–electron Coulomb collision frequency can be evaluated from plasma properties, we would like to numerically study the transition from the 10-moment model to the 5-moment model by artificially varying the self-collision frequency. The electron–electron collision frequency can be considered to consist of two components, classical and artificial

ν_{ee} = ν_{ee, cl} + ν_{ee, ar} .

(49)

The classical electron–electron frequency,

ν_{ee, cl} = ν_{ei}

⁠, while the artificial collision frequency,

ν_{ee, ar}

is a model parameter we vary. In the PIC-MCC, Coulomb collisions are neglected, i.e.,

ν_{ee} \to 0

⁠, whereas in the 5-moment model,

ν_{ee} \to \infty

⁠. For the 10-moment results presented, unless otherwise stated,

ν_{ee} = ν_{ei}

⁠, that is,

ν_{ee, ar} \to 0

⁠.

The electron-wall collision frequency is assumed to be

ν_{ew} = 10^{7} s^{- 1}

throughout the domain.^32,36 Furthermore, upon collisions with the wall, electrons are assumed to lose energy proportionally (not equally) from each component of

P_{i j}

⁠, according to the relationship

Δ ϵ_{w, i j} = P_{e, i j} exp (- \frac{U_{c}}{P_{e, i i} / 2 n_{e}}),

(50)

where the factor of two comes from the fact that

P_{i i}

is twice the total electron energy, and the parameter

U_{c} = 20

eV.

Finally, the anomalous momentum-transfer frequency is taken from the two-region model presented in Ref. 37,

ν_{ano} = α_{ano} ω_{c, e},

(51)

where

α_{ano}

is the anomalous transport coefficient. The coefficient is assumed to be

α_{in} = 1 / 160

within the channel and

α_{out} = 1 / 16

outside of the channel. To eliminate the sharp discontinuity at the channel exit,

α_{ano}

is linearly smoothed between the two values within the region between

L_{ch} - Δ L < x < L_{ch} + Δ L

where

Δ L

is taken to be 0.5 cm, following Refs. 4 and 29.

2. Collision model

For consistency between 10-moment, 5-moment, and PIC-MCC models, electrons are assumed to be isotropically scattered in the reference frame of the heavy species. For non-ionizing collisions, including electron-neutral (en), electron-ion (ei), and electron-wall (ew) collisions, the post-collision VDF

f_{c, s s^{'}}^{post}

in Eq. is assumed to be a Maxwellian with zero drift velocity

f_{c, en}^{post} = f_{c, i}^{post} = f_{c, ew}^{post} = \frac{n_{e}}{π^{3 / 2} v_{th, c}^{3}} exp (- \frac{v^{2}}{v_{th, c}^{2}}),

(52)

where the post-collision thermal velocity for collision type c is

m_{e} v_{th, c}^{2} = \frac{4}{3} [P_{e, i i} / (2 n_{e}) - Δ ϵ_{c}]

to account for the energy loss due to collisions. Note that for electron-neutral elastic collisions, it is assumed that

Δ ϵ_{c} = 0

as the energy exchange due to inelastic collisions dominate over the elastic collisions for electron temperatures above 5 eV, which is observed in this test case.

For electron-impact ionization reactions (iz), we assume that the primary and secondary electron VDFs are similar. Other works³⁸ consider the secondary electron VDF in more detail, but for the purpose of this study, the secondary electrons are assumed to have a negligible effect on the electron bulk velocity and temperature. In this case, the post-collision VDF can be given as

f_{iz}^{post} = \frac{2 n_{e}}{π^{3 / 2} v_{th, iz}^{3}} exp (- \frac{v^{2}}{v_{th, iz}^{2}}),

(53)

where

m_{e} v_{th, iz}^{2} = \frac{4}{3} [P_{e, i i} / (2 n_{e}) - Δ ϵ_{e, iz}]

⁠, i.e., each of the two post-collision electrons have the same energy. The factor of 2 in the numerator of Eq. arises from summing both the primary and secondary electron VDF.

For electron–electron Coulomb collisions, the post-collision VDF is assumed to be an isotropic Maxwellian VDF in the electron frame

f_{ee}^{post} = \frac{n_{e}}{π^{3 / 2} v_{th, e}^{3}} exp [- \frac{{(v - u_{e})}^{2}}{v_{th, e}^{2}}],

(54)

where

m_{e} v_{th, e}^{2} = 2 p_{e, i i} / (3 n_{e}) = 2 k_{B} T_{e}

⁠.

Finally, we substitute the post-collision VDFs for different collision processes, i.e., Eqs. , into the kinetic transport equation assuming a BGK collision operator, i.e., Eq. . Taking the moments of the kinetic equation results in the following collisional terms in Eqs. :

C_{e}^{(0)} = ν_{iz} ρ_{e},

(55)

C_{e, i}^{(1)} = - ν_{m} Γ_{e, i} = - ν_{m} ρ_{e} u_{e, i},

(56)

\begin{matrix} C_{e, i j}^{(2)} = ν_{m} (P_{e} δ_{i j} - P_{e, i j}) + ν_{ee} (p_{e} δ_{i j} - p_{e, i j}) \\ - (ν_{ew} Δ ϵ_{w, i j} + \frac{2}{3} δ_{i j} Q_{eff} ρ_{e}), \end{matrix}

(57)

where the total momentum transfer collision frequency

ν_{m} = ν_{en} + ν_{ei} + ν_{ew} + ν_{ano}

and the total collision frequency

ν_{tot} = ν_{m} + ν_{ee}

⁠. Equations describe the source term in conservation of mass due to ionization, the collisional drag in conservation of momentum, the time rate change of anisotropic pressure due to collisional heating and energy loss. Note that the factor of two-thirds in the second term of Eq. comes from the fact that

P_{e, i i} = 2 ϵ_{e}

and

δ_{i i} = 3

⁠. See Appendix C for a description of the primitive collisional terms.

3. Heat flux model

The thermal conductivity in the cross-field direction is given in Ref. 39 as

κ_{e, ⊥} = \frac{1}{1 + Ω_{H}^{2}} \frac{2.4}{1 + ν_{ei} / (\sqrt{2} ν_{m})} \frac{p_{e}}{ν_{m} ρ_{e}} \frac{k_{B}}{m_{e}},

(58)

where the Hall parameter

Ω_{H} = ω_{B} / ν_{m}

⁠. This formula is then used in Eq. to determine the heat flux components

q_{ijk}

⁠.

V. RESULTS

The 10-moment results are compared with the 5-moment results in Ref. 4 and the PIC-MCC results in Ref. 29. The PIC-MCC model evolves the particle evolution equations [Eqs. (19) and (20)] for the ion and electron species, and uses the same simplified neutral species update [Eq. (41)] as the fluid models. The PIC-MCC simulations take about 4–5 days using 30 processors on the Sherlock cluster at Stanford University to reach a steady state at about 1 ms, and the simulations are run for an additional 1 ms to average macroscopic quantities, which takes an additional 10 days.²⁹ The 5-moment results take about 5 h using 4 processors, while the 10-moment results take about 46 h on a single processor. The increased computational runtime for the 10-moment model as compared to the 5-moment model is due to both an increased number of fluid equations to solve and also requiring approximately $40 %$ smaller time step owing to the faster wave speed present in the 10-moment model. As discussed in Eq. (35), the ratio of the speed of sound between the 10-moment and 5-moment equations is $\sqrt{3 / (5 / 3)} ≃ 35 %$ near equilibrium.

A. Plasma profiles

Figure 2 shows a comparison of the steady-state results between the 10-moment, 5-moment, and PIC-MCC models. The kinetic model is solved using the algorithms described in Ref. 29. As described in Ref. 4, the values of $α_{ano}$ in Eq. (51) are chosen so that the system attains a steady-state profile and that there are no large-scale oscillations (cf. breathing mode⁴⁰) to provide a quantitative comparison. In general, there is good agreement between the 10-moment and 5-moment simulations, demonstrating that the first-order approximation of the electrons having a Maxwellian VDF can capture much of the dynamics in this regime. However, there are a number of key differences between the 10-moment and 5-moment results.

FIG. 2.

View large Download slide

Comparison between the results obtained from the 10-moment (10M), 5-moment (5M), and PIC-MCC models: (a) ion (dark) and electron (light) number densities, (b) neutral number density, (c) ion axial velocity, (d) axial electric field, (e) electron axial velocity, and (f) electron azimuthal velocity.

Figure 2(a) shows that the 10-moment and 5-moment estimates of electron and ion densities are higher than those of the PIC-MCC simulation. Furthermore, the 10-moment estimate within the channel is greater than that of the 5-moment. The anode sheath is resolved for all results, exhibiting an ion-attracting, electron-repelling sheath. Figure 2(b) shows that the neutral atom density obtained from the 5- and 10-moment models is lower than the PIC-MCC results in the downstream (⁠ $x \geq 1$ cm), which is attributable to the greater degree of ionization in the fluid models. Figure 2(c) shows a discrepancy of the ion bulk velocity between the PIC-MCC model and the two fluid models. As described in Ref. 29, this is most likely due to the PIC-MCC model being able to capture non-Maxwellian ion VDFs: both the fast ions accelerated from upstream and the slow ions generated from ionization events. Figure 2(d) shows the axial electric field. The PIC-MCC and 10-moment results agree better than the 5-moment results, especially near the channel exit where the $E \times B$ drift is largest, e.g., $x \in [2$ cm, 2.5 cm]. This demonstrates that the 10-moment model captures the electron kinetics near the acceleration region. We further investigate the anisotropic electron temperature and cross field electron transport in Sec. V B. Figure 2(e) shows the electron axial bulk velocity. Here again, good agreement between 5-moment and 10-moment is obtained. The anode discharge current, $I_{D} = \sum_{s} q_{s} Γ_{s, x}$ ⁠, is in good agreement between each of the three models, within 1%. The steady-state $I_{D}$ is 8.484, 8.457, and 8.415 A for the 10-moment, 5-moment, and PIC-MCC models, respectively. Because the anode current is similar between the two models, the bulk velocity follows a reciprocal trend to the density inside the channel, e.g., $x \leq 2$ cm. Figure 2(f) shows the electron azimuthal bulk velocity. It is clear that the azimuthal velocity consists predominantly of the $E \times B$ drift, as the profile of $u_{e, y}$ closely resembles the electric field profile. The maximum magnitude of $u_{e, y}$ agrees well between 10-moment and PIC-MCC due to the agreement in $E_{x}$ as shown in Fig. 2(d). A notable difference between 5- and 10-moment models is the azimuthal electron transport near the anode. The difference in ionization and transport near the anode can be understood by studying the electron temperature profiles.

Figure 3 presents a comparison of the anisotropic temperature profiles between 10-moment and PIC-MCC simulations. While the 5-moment model only captures the isotropic temperature, the 10-moment and PIC-MCC models capture all of the elements of the temperature tensor. Figure 3(a) shows the on-diagonal (normal-stress) components of the 10-moment temperature tensor compared to the isotropic 5-moment temperature. The isotropic temperatures obtained from the 10-moment and 5-moment models are in good agreement. It is also shown that the deviation from isotropy occurs most strongly near the anode and the channel exit, where there is a sharp change in the anomalous collision frequency.

FIG. 3.

View large Download slide

Comparison of the values of the on-diagonal and off-diagonal electron temperature tensor elements between [(a) and (c)] 10-moment and [(b) and (d)] PIC-MCC simulations. To help visualize the comparison, the 5-moment results are shown as dotted lines in both FIGS. The channel exit at $x = 2.5$ cm is denoted by the dotted teal vertical line.

Figure 3(b) shows a similar comparison between the PIC-MCC on-diagonal temperatures and the 5-moment isotropic temperature. It can be seen from Figs. 3(a) and 3(b) that the 10-moment and PIC-MCC results agree qualitatively, though PIC-MCC predicts a greater difference between $T_{e, z z}$ and the other on-diagonal temperatures near the channel exit. The temperature within the PIC-MCC simulation is about 2 eV lower than that of the fluid models. One factor in this discrepancy is that the magnitude of the axial electron velocity in the plume [ $x \geq 2.5$ cm in Fig. 2(e)] is greater in the fluid models; this leads to more Joule heating, whence the hotter electrons are advected toward the anode.

Except near the anode, the radial temperature $T_{e, z z}$ is the smallest. This is to be expected because the majority of heating arises from collisional drag (resistivity). Without any radial velocity, $u_{e, z}$ ⁠, the corresponding temperature will be lower compared to the other two. Furthermore, the radial magnetic field $B_{z}$ does not affect $T_{e, z z}$ ⁠, so there is no mechanism for it to equilibrate to the other temperatures except purely through self-collisions in the present 1D simulation. Just within the channel at $x ≃ 2.3$ cm, the collision frequency drops due to the change in $ν_{ano}$ ⁠, and the plasma sustains significantly more anisotropy than in the plume. As shown in Fig. 2(e), the radial speed $| u_{e, x} |$ also decreases at $x ≃ 2.3$ cm, leading to compression and axial heating. However, just afterwards (from the perspective of the electrons), at $x ≃ 2$ cm, the azimuthal speed, $| u_{e, y} |$ ⁠, increases rapidly as the electrons experience the strong $E \times B$ drift shown in Fig. 2(f). This increased velocity leads to azimuthal heating through collisions. As the azimuthal velocity is much larger than the axial velocity, the azimuthal temperature $T_{e, y y}$ remains the highest at $x \leq 2$ cm.

Another difference is the magnitude and profile of the electric field near the channel exit (i.e., $x \in [2 cm, 2.5 cm]$ ⁠) seen in Fig. 2(d). This discrepancy can be explained by recalling that in the 10-moment model $T_{e, x x} > T_{e}$ ⁠, and since axial transport is mediated by $T_{e, x x}$ ⁠, the electrons have increased thermal transport across the magnetic field; this allows them to more readily neutralize ion density perturbations, thus decreasing the magnitude of the electric field. This can also be thought of as a decrease in the effective resistivity and consequently, the electric field. As a result, the 10-moment electric field profile is in better agreement with PIC-MCC than the 5-moment model. This reduction in the electric field also manifests itself as a smaller peak azimuthal velocity [cf. Fig. 2(f)], also in better agreement with PIC-MCC.

At the anode, a sheath potential forms, which truncates the electron VDF in the axial direction, i.e., electrons with sufficiently negative velocity will be absorbed by the anode, therefore depleting the electrons in the $- x$ direction. This manifests itself as a decrease in the axial temperature, $T_{e, x x}$ ⁠, near the anode as seen in Figs. 3(a) and 3(b). This anode sheath region is larger in the PIC simulation than in the fluid, as can be seen in Fig. 3, most likely due to the kinetic simulation being able to sustain non-thermal populations while the fluid simulation necessarily relaxes back to equilibrium more quickly. The suppressed axial temperature leads to a smaller particle flux to the wall, cf. Eq. (B4).

Figures 3(c) and 3(d) show the only non-zero off diagonal (shear) electron temperature component, $T_{e, x y}$ ⁠, for the 10-moment model and the PIC-MCC simulations, respectively. Note that $T_{e, x y} = 0$ in the 5-moment model. We see quantitative agreement between the 10-moment and kinetic models in capturing $T_{e, x y}$ ⁠. There are shear temperatures near the wall and the channel exit, and notably, it is non-zero even in regions of high collision frequency. The origin of $T_{e, x y}$ will be elaborated in detail in Sec. V B.

B. Cross-field electron transport

It is observed from both the 10-moment and PIC-MCC results that the electron distribution functions within a Hall-effect thruster plasma are anisotropic, as shown in Fig. 3. In this section, we discuss the mechanisms of the electron anisotropy and their effects on cross field electron transport.

1. Mechanism of anisotropy

By subtracting the primitive equations for

p_{y y}

from that for

p_{x x}

in Eqs. , the evolution of pressure anisotropy,

Δ p \equiv p_{x x} - p_{y y}

⁠, can be obtained as

\begin{matrix} \frac{\partial Δ p}{\partial t} + u_{x} \frac{\partial Δ p}{\partial x} + (2 p_{x x} + Δ p) \frac{\partial u_{x}}{\partial x} - 2 p_{x y} \frac{\partial u_{y}}{\partial x} + \frac{\partial Δ q}{\partial x} \\ = \frac{4 e B_{z}}{m} p_{x y} - ν_{ee} Δ p + ν_{m} (ρ u_{x}^{2} - ρ u_{y}^{2} - Δ p) \\ + ν_{ew} (Δ ϵ_{w, x x} - Δ ϵ_{w, y y}) + ν_{iz} ρ (u_{x}^{2} - u_{y}^{2}), \end{matrix}

(59)

where

Δ q = q_{xxx} - q_{yyx}

⁠. At steady state, the first term vanishes and the off diagonal component

p_{x y}

can be isolated, and terms can be grouped and identified

\begin{matrix} (4 ω_{B} + 2 \frac{\partial u_{y}}{\partial x}) p_{x y} = \underset{Transport}{\underset{︸}{u_{x} \frac{\partial Δ p}{\partial x} + (2 p_{x x} + Δ p) \frac{\partial u_{x}}{\partial x}}} \\ + \underset{Heat Flux}{\underset{︸}{\frac{\partial Δ q}{\partial x}}} + \underset{Isotrop .}{\underset{︸}{ν_{tot} Δ p}} + \underset{Elastic Energy Exchange}{\underset{︸}{ν_{m} ρ (u_{y}^{2} - u_{x}^{2})}} + \underset{Inelastic Energy Exchange}{\underset{︸}{ν_{w} (Δ ϵ_{w, y y} - Δ ϵ_{w, x x}) + ν_{iz} ρ (u_{y}^{2} - u_{x}^{2})}}, \end{matrix}

(60)

where the sources and sinks of anisotropy are identified as (i) anisotropic inviscid transport, (ii) conductive heat flux, (iii) isotropizing collisions, (iv) elastic energy exchange, and (v) inelastic energy exchange. Dividing through by the first term on the left-hand side yields an expression for the steady-state off diagonal pressure

p_{x y}

⁠. The relative magnitude of each of the terms can then be analyzed to determine the effect of each on generating anisotropy.

Figure 4 shows that $p_{e, x y}$ is dominated by elastic energy exchange in the plume (⁠ $x \geq 3$ cm); i.e., the high anomalous transport frequency causes drag, which stretches the distribution function unequally in the x and y directions due to $u_{e, x}$ and $u_{e, y}$ ⁠, as shown in Figs. 2(e) and 2(f), respectively. In addition, in this region, the transport term is small because $Δ p ≃ 0$ due to the high collisionality.

FIG. 4.

View large Download slide

Decomposition of anisotropic electron pressure $p_{x y}$ into components as described in Eq. (60): (i) transport (gray solid line), (ii) heat flux (red solid line), (iii) isotropizing collisions (magenta dashed line), (vi) elastic energy exchange (green dotted line), (v) inelastic energy exchange (blue dotted line), and (vi) $p_{x y}$ (black solid line).

Within the channel (⁠ $x \leq 2.5$ cm), momentum transfer collision frequency, $ν_{m}$ ⁠, is significantly smaller than in the plume, and the inviscid transport has a greater effect on the anisotropy. In one spatial dimension, the transport-based anisotropy is primarily due to the compression of the fluid, $\partial_{x} u_{e, x}$ ⁠. It can be seen from Eq. (60) that a non-zero $p_{e, x y}$ can be generated by $p_{e, x x} \partial_{x} u_{e, x}$ even when $Δ p_{e} = 0$ ⁠. Consider the profile of $u_{e, x}$ in Fig. 2(e). From $x \geq 1.75$ cm, $\partial_{x} u_{e, x} < 0$ ⁠, i.e., the electrons are being compressed; this leads to an increase in $p_{e, x x}$ ⁠, while leaving $p_{e, y y}$ unchanged and is a source of $p_{e, x y}$ ⁠. Conversely, toward the anode, $x < 1.75$ cm, where $\partial_{x} u_{e, x} > 0$ ⁠, i.e., the electrons are accelerating and expanding, $p_{e, x x}$ decreases and there is a sink of $p_{e, x y}$ ⁠.

Inviscid transport dominates until the sheath boundary layer beside the anode, i.e., $x \leq 0.5$ cm. In this region, the heat flux is large (⁠ $Δ q = 2 q_{e, yyx} = - 0.8 κ_{⊥} \partial_{x} T_{e}$ ⁠) and leads to a reversal in the sign of $p_{e, x y}$ at the anode, cf. Fig. 3. The large gradients of shear pressure $p_{e, x y}$ and heat flux $q_{e, x}$ at the anode boundary indicate that proper models for capturing the strongly non-Maxwellian VDFs at phase boundaries are important for understanding plasma-wall interactions like electron absorption and erosion.

2. Drift velocity decomposition

By taking the cross product of the primitive momentum equation [Eqs. ] with the magnetic field, the velocity perpendicular to the magnetic field,

u_{⊥} = u - (u \cdot \hat{b}) \hat{b}

can be written as

\begin{matrix} u_{⊥} = \underset{E \times B drift}{\underset{︸}{\frac{E \times B}{B^{2}}}} + \underset{Hall effect}{\underset{︸}{(\frac{ν_{m} + ν_{iz}}{ω_{B}} u) \times \hat{b}}} \\ + \underset{Shear Drift}{\underset{︸}{(\frac{u \cdot \nabla}{ω_{B}} u) \times \hat{b}}} + \underset{Diamag . Drift}{\underset{︸}{\frac{(\nabla \cdot \bar{\bar{p}}) \times B}{e n B^{2}}}}, \end{matrix}

(61)

where B is the magnetic field strength and

\hat{b} = B / B

is the unit vector in the direction of the magnetic field. The cross field electron bulk velocity can be due to (i) the classical

E \times B

drift, (ii) the Hall effect drift due to collisional drag, (iii) a shear drift from electron inertial effects, and (iv) a generalized diamagnetic drift, including both the effects from the isotropic pressure (which corresponds to the conventional diamagnetic drift) and anisotropic pressure.²⁹ Each of these terms in Eq. , as well as

u_{⊥}

is presented in Fig. 5. Because there is no azimuthal electric field

E_{y}

in our simulation, there is no component of the

E \times B

drift in the axial direction.

FIG. 5.

Comparison of electron (a) axial and (b) azimuthal velocity profiles with drifts due to E×B (red line), Hall effect from collisional drag (green line), generalized diamagnetic drift (magenta line), and shear drift (dotted purple line).

View large Download slide

Comparison of electron (a) axial and (b) azimuthal velocity profiles with drifts due to $E \times B$ (red line), Hall effect from collisional drag (green line), generalized diamagnetic drift (magenta line), and shear drift (dotted purple line).

Figure 5(a) shows that within the domain, the electron axial velocity is mostly dominated by the Hall effect, given the large azimuthal velocity. It can be seen that the other the shear and diamagnetic drift effects become significant in the near-anode region (⁠ $x \leq 0.5$ cm), near the channel exit (⁠ $x ≃ 2.5$ cm), and at the cathode (⁠ $x \geq 4.5$ cm). Near the anode, the azimuthal velocity decreases and instead the electron inertial effects dominate, particularly the generalized diamagnetic drift, namely, the term proportional to $\partial_{x} p_{e, x y}$ ⁠. This finding is consistent with that in Ref. 29, reaffirming the capability of the 10-moment model in analyzing plasma-wall interactions where shear effects are non-negligible. Note that when comparing with the 5-moment model, the estimate of the axial velocity near the anode is similar, even though the 5-moment model does not account for the anisotropic pressure. This is achieved by the 5-moment model overpredicting the azimuthal velocity near the anode, as shown in Fig. 2(f). That is to say, to achieve the same axial transport properties and steady state flux of electrons to the anode, the 5-moment model substitutes $ρ u_{x} \partial_{x} u_{y}$ for $\partial_{x} p_{x y},$ leading to a different profile of $u_{e, y}$ near the anode.

Near the channel exit, e.g., $x \in [2$ cm, 2.6 cm], where both the azimuthal velocity and the anomalous collision frequency have strong gradients, inertial (shear and generalized diamagnetic) effects also play a significant role. Finally, inertial effects are also large near the cathode, where electrons are injected, e.g., $x \geq 4.8$ cm. This is most likely due to the fact the electrons are assumed to be non-drifting and thermally isotropic as they enter the domain, which is slightly inconsistent with the natural plasma response at $x = 5$ cm. Thus, this effect is most likely an artifact of the electron injection method. Previous works (Refs. 4 and 29) have found that the details of the injection method do not significantly affect the plasma profiles within the domain or global parameters like the discharge current. A more compatible injection scheme is reserved for future work.

Figure 5(b) shows the decomposition of the azimuthal velocity, $u_{e, y}$ ⁠. The profile is mainly dominated by the $E \times B$ drift as one might expect but inertial effects, in particular, the diamagnetic drift, $\partial_{x} p_{e, x x}$ ⁠, again become significant near the channel exit at $x = 2.5$ cm and near the anode, i.e., $x < 1$ cm, due to the steep gradients of $u_{e, x}$ and $p_{e, x x}$ ⁠. While they are not zero, the contributions from the Hall effect and shear terms to the azimuthal velocity in Eq. (61) are negligible.

C. Importance of self-consistent heat flux closure

The derivation of the heat conductivity, i.e., Eq. (58), from Ref. 39 makes the assumptions that (a) electrons are nearly Maxwellian and (b) the contributions of electron–electron collisions to the thermal conductivity are negligible. Clearly, this is valid in the limit of $ν_{ee, ar} \to \infty$ in Eq. (49). In this limit, the electron VDF immediately relaxes to an equilibrium Maxwellian and any associated heat flux, $\propto 1 / ν_{ee}$ ⁠, vanishes. For this reason, the use of this equation in the 5-moment model is valid. Conversely, for the previous 10-moment results shown, $ν_{ee, ar} \to 0$ and no additional heat flux due to ee collisions is considered. Although the first assumption of near-Maxwellian is not strictly true, the electron collisions are dominated by interspecies collisions and thus the second assumption is still satisfied. In the intermediate regime of $ν_{ee, ar}$ ⁠, however, neither assumption is valid. We can estimate the minimum $ν_{m}$ for which anisotropy becomes significant; i.e., when the dynamic length scale is on the order of the collisional mean free path. We can estimate the dynamic scale by the gradient-length scale, $max (| Q | / | \nabla Q |)$ ⁠. Taking $T_{e}$ to be the quantity of interest, we find that $L_{dyn} ≃ 0.01$ cm. So then, we would expect our formulation to be inconsistent with Eq. (58) around $ν_{ee} ≃ u_{th} / L_{dyn} \approx 3 \times 10^{10}$ s⁻¹.

Figure 6 shows the simulated anode current as a function of the artificial electron–electron collision frequency, $v_{ee, ar}$ in Eq. (49). As $ν_{ee, ar}$ decreases from infinity, the distribution function begins to be more anisotropic, but because in this regime $ν_{ee} ≫ ν_{m}$ ⁠, the heat flux and therefore the transport is grossly overestimated, leading to a rise in the discharge current estimate. Indeed, we see significant deviation from the equilibrium estimate near the transition value of $ν_{ee, ar} \approx 10^{10}$ s⁻¹ calculated above. Conversely, electron–electron collisions are negligible when $ν_{ee} < ν_{m} ≃ 5 \times 10^{7}$ s⁻¹, on average. Indeed, we see that by this point, the second assumption is satisfied. Therefore, these results demonstrate the importance of having a self-consistent heat flux closure in the 10-moment model and making sure that the assumptions of prior derivations are still satisfied when testing a new model. In addition, it shows that the derivation of a new closure model, consistent with anisotropic non-Maxwellian VDFs may be necessary when studying plasmas which span wide ranges of collisionality.

FIG. 6.

View large Download slide

Comparison of the calculated discharge current $I_{D}$ as a function of the artificial electron–electron collision frequency $ν_{ee, ar}$ ⁠.

VI. CONCLUSIONS

In this paper, we developed a 10-moment model applicable to low-temperature partially magnetized plasmas and applied the model to the discharge chamber of a simplified Hall-effect thruster geometry. The results show good agreement with 5-moment fluid and fully kinetic (PIC-MCC) results, both in spatial plasma profiles and global discharge current. The model was developed to study the ability of higher-order fluid moment models to accurately capture non-Maxwellian electron VDFs.

The anisotropic pressure in low-temperature magnetized plasmas is captured qualitatively well, and the nature of the fluid model provides insight into how electron anisotropy is generated within the cross field discharge plasma and the relative importance of each of the components on the overall profile. It was found that capturing a full pressure tensor was vital to understanding the dynamics present at plasma-wall interactions, as well as phase transition regions like that between the thruster chamber and plume. For this reason, the 10-moment model is in better agreement with PIC-MCC than the 5-moment model, not only in the estimates of the pressure tensor, but also in the behavior of transverse velocities near the wall and the maximum electric field. The importance of heat flux closure to the 10-moment model was also demonstrated, as well as the necessity of self-consistency when applying the closure.

Future work may study more general means of heat flux closure compatible with far-from-Maxwellian electron VDFs, and varying levels of inter- and intra-species collisionality. It may also look into generalized methods for calculating rate coefficients that incorporate the information of the full temperature tensor. While in general, the rate coefficient k may be a function of all six components of $T_{i j}$ and the drift speed u, finding low-rank models to accurately calculate k may be useful in future studies including the 10-moment model.

Regions of non-zero $p_{x y}$ indicate non-equilibrium shear that may be indicative of plasma turbulence. Future studies may look into relaxing the assumptions of $ν_{ano}$ ⁠, which were made to achieve a steady-state profile, and determining if the anomalous electron cross field transport may be explained in part by the non-equilibrium effects captured in the 10-moment model. Studies of three-dimensional plasma turbulence would be greatly accelerated by the use of a fluid model which can capture the relevant non-equilibrium physics.

ACKNOWLEDGMENTS

The authors thank A. R. Mansour for many helpful discussions on modeling the dynamics within Hall effect thrusters. This work was supported by a NASA Space Technology Graduate Research Opportunity, Grant No. 80NSSC21K1302, the U.S. Department of Energy National Nuclear Security Administration Stewardship Science Graduate Fellowship under Cooperative Agreement No. DE-NA0003960, NASA through the Joint Advanced Propulsion Institute, a NASA Space Technology Research Institute, under Grant No. 80NSSC21K1118, and the Office of Naval Research under Award No. N00014-21-1-2698. The authors would like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that contributed to these research results.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Derek Kuldinow: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Yusuke Yamashita: Software (supporting); Writing – review & editing (equal). Kentaro Hara: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Methodology (supporting); Project administration (lead); Visualization (supporting); Writing – review & editing (equal).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX A: JUSTIFICATION OF THE FORM OF THE 10-MOMENT HEAT FLUX MODEL

The perturbation function

f^{1}

is commonly written in the form

f^{1} = - τ [g (w) (w \cdot F) + E (w)] f^{M},

(A1)

for some even function

g (w)

⁠, another even function

E (w)

⁠, and a vector independent of w, F, where

f^{M}

is a Maxwellian distribution. One example is the typical Chapman–Enskog expansion

\begin{matrix} f^{1} = - τ [(\frac{ρ w_{i} w_{j}}{p} - \frac{p_{i j}}{p} \frac{ρ w^{2}}{3 p} + \frac{p_{i j}}{p} - δ_{i j}) \frac{\partial u_{i}}{\partial x_{j}} \\ + w_{j} (\frac{ρ w^{2}}{2 p} - \frac{5}{2}) \frac{\partial}{\partial x_{j}} ln (\frac{p}{ρ})] f^{M}, \end{matrix}

(A2)

where

\begin{matrix} g (w) = \frac{ρ w^{2}}{2 p} - \frac{5}{2}, F = \nabla [ln (\frac{p}{ρ})], \\ E (w) = (\frac{ρ w_{i} w_{j}}{p} - \frac{p_{i j}}{p} \frac{ρ w^{2}}{3 p} + \frac{p_{i j}}{p} - δ_{i j}) \frac{\partial u_{i}}{\partial x_{j}} . \end{matrix}

(A3)

Another example is the first term in the gyrophase-independent part of Braginskii's expansion of the VDF⁴¹

\begin{matrix} {⟨ f^{1} ⟩}_{ϑ} = - τ w_{| |} (\frac{5}{2} - \frac{w^{2}}{v_{e, th}^{2}}) \times (- 0.506 \nabla_{| |} ln T_{e} - 0.284 \frac{m u_{e, | |}}{τ T_{e}}) f^{M}, \end{matrix}

(A4)

where

{⟨ \cdot ⟩}_{ϑ}

denotes the gyroaverage over cyclotron orbits,

τ

is the electron-ion collision timescale, subscript

∥

denotes the direction along the magnetic field,

\hat{b}

⁠. In this case,

g (w) = \frac{5}{2} - \frac{w^{2}}{v_{th}^{2}}

⁠,

E = 0

⁠, and

F = \hat{b} (- 0.506 \nabla_{| |} ln T_{e} - 0.284 \frac{m u_{e, | |}}{τ T_{e}})

⁠. Likewise for the other components of the Braginskii closure.

The heat flux is a

3 \times 3 \times 3

tensor. Without loss of generality, let us consider the heat flux in the z direction. One of the heat flux elements in the z direction can be given by

q_{zzz} = \int w_{z}^{3} f^{1} d w = - τ F \cdot \int w_{z}^{3} w g (w) f^{M} d w .

(A5)

It can be seen that the even part,

E (w)

⁠, in Eq. does not contribute to the heat flux because the overall odd-order integral over

f^{M}

vanishes. Furthermore, the integrand vanishes for odd functions except for the component of w in the z direction. Thus,

q_{zzz} = - τ F_{z} \int w_{z}^{4} g (w) f^{M} d w .

(A6)

Then, converting to spherical coordinates yields

\begin{matrix} q_{zzz} = - τ F_{z} \int {(w cos ϕ)}^{4} g (w) f^{M} (w^{2} sin ϕ) dwd ϕ d θ \\ = - τ F_{z} \int sin ϕ {cos}^{4} ϕ d ϕ d θ \int w^{6} g (w) f^{M} d w . \end{matrix}

(A7)

Likewise, next we look at the flux of transverse energies in the z direction

q_{xxz} = \int w_{x}^{2} w_{z} f^{1} d w = - τ F \cdot \int w_{x}^{2} w_{z} w g (w) f^{M} d w .

(A8)

Again, the term we need to pair with to make even is

w_{z}

⁠, so again we extract the component in the z direction

q_{xxz} = - τ F_{z} \int w_{x}^{2} w_{z}^{2} g (w) f^{M} d w .

(A9)

This can be converted to spherical coordinates

q_{xxz} = - τ F_{z} \int {cos}^{2} θ {sin}^{3} ϕ {cos}^{2} ϕ d ϕ d θ \int w^{6} g (w) f^{M} d w .

(A10)

The equation for

q_{xxz}

is exactly the same as the equation for

q_{zzz}

except for the prefactor, which is a computable scalar, independent of the form of

f^{1}

\int sin ϕ {cos}^{4} ϕ d ϕ d θ = \frac{4 π}{5},

(A11)

\int {cos}^{2} θ {sin}^{3} ϕ {cos}^{2} ϕ d ϕ d θ = \frac{4 π}{15} .

(A12)

So, by comparing Eqs. and it can be seen that

q_{zzz} = 3 q_{xxz} = 3 q_{yyz}

⁠; and, if the contracted heat flux is written as

q_{z} = q_{zzz} + q_{yyz} + q_{xxz}

⁠, they must be similar forms but scaled in a ratio of

3 : 1 : 1

⁠.

APPENDIX B: 10-MOMENT KINETIC FLUXES

The definition of the kinetic fluxes is presented in Eq. (37), wherein the flux is calculated based on an integral of a reconstructed VDF at the boundary. For our 10-moment model, the fluid variables at the boundary are reconstructed using MUSCL and the associated VDF is as in Eq. (38).

From these equations, we need to determine the fluxes of each of the conservative quantities, ${ρ, ρ u_{i}, p_{i j} + ρ u_{i} u_{j}}$ ⁠. Performing these integrals yields the following equations for the rightward (+) and leftward (−) moving fluxes which are then used in Eqs. (33) and (42) to update the conservative variables.

To express the form of the fluxes in a succinct manner, let us introduce the normalized velocities as

\hat{u} = \frac{u_{x}}{\sqrt{2 p_{x x} / ρ}}, \hat{v} = \frac{u_{y}}{\sqrt{2 p_{y y} / ρ}}, \hat{w} = \frac{u_{z}}{\sqrt{2 p_{z z} / ρ}},

(B1)

the thermal speed based on

T_{x x}

as

u_{th} = \sqrt{\frac{8 p_{x x}}{π ρ}},

(B2)

and the normalized anisotropic pressures as

τ_{x y} = \frac{p_{x y}}{\sqrt{p_{x x} p_{y y}}}, τ_{x z} = \frac{p_{x z}}{\sqrt{p_{x x} p_{z z}}}, τ_{y z} = \frac{p_{y z}}{\sqrt{p_{y y} p_{z z}}} .

(B3)

Below we show the kinetic fluxes for each moment.

1. 0th velocity moment

Γ_{ρ}^{\pm} = \pm \frac{ρ u_{th}}{4} [exp (- {\hat{u}}^{2}) + \hat{u} \sqrt{π} [erf (\hat{u}) \pm 1]] .

(B4)

2. 1st velocity moments

Γ_{ρ v_{x}}^{\pm} = \pm p_{x x} [\frac{\hat{u}}{\sqrt{π}} exp (- {\hat{u}}^{2}) + (\frac{1}{2} + {\hat{u}}^{2}) [erf (\hat{u}) \pm 1]],

(B5)

\begin{matrix} Γ_{ρ v_{y}}^{\pm} = \pm \sqrt{p_{x x} p_{y y}} [\frac{\hat{v}}{\sqrt{π}} exp (- {\hat{u}}^{2}) \\ + (\frac{τ_{x y}}{2} + \hat{u} \hat{v}) [erf (\hat{u}) \pm 1]], \end{matrix}

(B6)

\begin{matrix} Γ_{ρ v_{z}}^{\pm} = \pm \sqrt{p_{x x} p_{z z}} [\frac{\hat{w}}{\sqrt{π}} exp (- {\hat{u}}^{2}) \\ + (\frac{τ_{x z}}{2} + \hat{u} \hat{w}) [erf (\hat{u}) \pm 1]] . \end{matrix}

(B7)

3. 2nd velocity moments

\begin{matrix} Γ_{ρ v_{x}^{2}}^{\pm} = \pm \frac{p_{x x} u_{th}}{2} [(1 + {\hat{u}}^{2}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} (\frac{3}{2} \hat{u} + {\hat{u}}^{3}) [erf (\hat{u}) \pm 1]], \end{matrix}

(B8)

\begin{matrix} Γ_{ρ v_{x} v_{y}}^{\pm} = \pm \frac{\sqrt{p_{x x} p_{y y}} u_{th}}{2} [(τ_{x y} + \hat{u} \hat{v}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} (\frac{\hat{v}}{2} + τ_{x y} \hat{u} (1 + \hat{u} \hat{v})) [erf (\hat{u}) \pm 1]], \end{matrix}

(B9)

\begin{matrix} Γ_{ρ v_{x} v_{z}}^{\pm} = \pm \frac{\sqrt{p_{x x} p_{z z}} u_{th}}{2} [(τ_{x z} + \hat{u} \hat{w}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} (\frac{\hat{w}}{2} + τ_{x z} \hat{u} (1 + \hat{u} \hat{w})) [erf (\hat{u}) \pm 1]], \end{matrix}

(B10)

\begin{matrix} Γ_{ρ v_{y} v_{y}}^{\pm} = \pm \frac{p_{y y} u_{th}}{2} [({\hat{v}}^{2} + \frac{1 + τ_{x y}^{2}}{2}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} (\frac{\hat{u}}{2} + \hat{v} (τ_{x y} + \hat{u} \hat{v})) [erf (\hat{u}) \pm 1]], \end{matrix}

(B11)

\begin{matrix} Γ_{ρ v_{y} v_{z}}^{\pm} = \pm \frac{\sqrt{p_{y y} p_{z z}} u_{th}}{2} [(\frac{τ_{x y} τ_{x z} + τ_{y z}^{2}}{2} + \hat{v} \hat{w}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} (\frac{τ_{y z} \hat{u} + τ_{x z} \hat{v} + τ_{x y} \hat{w}}{2} + \hat{u} \hat{v} \hat{w}) [erf (\hat{u}) \pm 1]], \end{matrix}

(B12)

\begin{matrix} Γ_{ρ v_{z} v_{z}}^{\pm} = \pm \frac{p_{z z} u_{th}}{2} [({\hat{w}}^{2} + \frac{1 + τ_{x z}}{2}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} (\frac{\hat{u}}{2} + \hat{w} (τ_{x z} + \hat{u} \hat{w})) [erf (\hat{u}) \pm 1]] . \end{matrix}

(B13)

Note that the total energy flux,

Γ_{\frac{1}{2} ρ v^{2}} = \frac{1}{2} Γ_{ρ v_{i} v_{i}}

can be derived as

\begin{matrix} Γ_{ρ v^{2}}^{\pm} = \pm \frac{u_{th}}{2} {(p_{x x} + \frac{p_{y y} + p_{z z}}{2} + \frac{τ_{x y}^{2} + τ_{x z}^{2}}{2} \\ + p_{x x} {\hat{u}}^{2} + p_{y y} {\hat{v}}^{2} + p_{z z} {\hat{w}}^{2}) exp (- {\hat{u}}^{2}) \\ + \sqrt{π} [\hat{u} (\frac{3}{2} p_{x x} + \frac{p_{y y}}{2} + \frac{p_{z z}}{2}) + τ_{x y} p_{y y} \hat{v} + τ_{x z} p_{z z} \hat{w} + \hat{u} (p_{x x} {\hat{u}}^{2} + p_{y y} {\hat{v}}^{2} + p_{z z} {\hat{w}}^{2})] [erf (\hat{u}) \pm 1]} . \end{matrix}

(B14)

When we consider a Maxwellian VDF where

p_{i j} = p δ_{i j}

⁠, the total energy flux reduces to

\begin{matrix} Γ_{\frac{1}{2} ρ v^{2}}^{\pm} = \pm \frac{u_{th}}{4} p [({\hat{u}}^{2} + {\hat{v}}^{2} + {\hat{w}}^{2} + 2) exp (- {\hat{u}}^{2}) + \sqrt{π} \hat{u} (\frac{5}{2} + {\hat{u}}^{2} + {\hat{v}}^{2} + {\hat{w}}^{2}) [erf (\hat{u}) \pm 1]], \end{matrix}

(B15)

which is consistent with the 5-moment equation.⁴

APPENDIX C: PRIMITIVE COLLISION QUANTITIES

Here, we show the detailed forms of the collisional terms. From Eqs. and Eqs. , the collisional terms in the primitive form can be written as

c_{e}^{(0)} = ν_{iz} ρ_{e},

(C1)

c_{e, i}^{(1)} = \frac{1}{ρ_{e}} (- ν_{m} ρ_{e} u_{e, i} - u_{e, i} ν_{iz} ρ_{e}) = - (ν_{m} + ν_{iz}) u_{e, i},

(C2)

\begin{matrix} c_{e, i j}^{(2)} = ν_{m} (P_{e} δ_{i j} - P_{e, i j}) + ν_{ee} (p_{e} δ_{i j} - p_{e, i j}) \\ - (ν_{ew} Δ ϵ_{w, i j} + \frac{2}{3} δ_{i j} Q_{eff} ρ_{e}) \\ + ν_{m} ρ_{e} u_{e, j} u_{e, i} + u_{e, j} ν_{m} ρ_{e} u_{e, i} + ν_{iz} u_{e, i} u_{e, j} ρ_{e} \\ = ν_{m} \frac{ρ u_{e}^{2} δ_{i j}}{3} + (ν_{m} + ν_{iz}) ρ_{e} u_{e, i} u_{e, j} + ν_{tot} (p_{e} δ_{i j} - p_{e, i j}) \\ - (ν_{ew} Δ ϵ_{w, i j} + \frac{2}{3} δ_{i j} Q_{eff} ρ_{e}) . \end{matrix}

(C3)

REFERENCES

1.

J.

Birn

,

J.

Drake

,

M.

Shay

,

B.

Rogers

,

R.

Denton

,

M.

Hesse

,

M.

Kuznetsova

,

Z.

Ma

,

A.

Bhattacharjee

,

A.

Otto

et al, “

Geospace environmental modeling (GEM) magnetic reconnection challenge

,”

J. Geophys. Res.: Space Phys.

106

,

3715

–

3719

, https://doi.org/10.1029/1999JA900449 (

2001

).

Google Scholar

Crossref

2.

J.

Büchner

,

C.

Dum

, and

M.

Scholer

,

Space Plasma Simulation

(

Springer Science & Business Media

,

2003

), Vol.

615

.

Google Scholar

Crossref

3.

L.

Wang

,

K.

Germaschewski

,

A.

Hakim

,

C.

Dong

,

J.

Raeder

, and

A.

Bhattacharjee

, “

Electron physics in 3-D two-fluid 10-moment modeling of Ganymede's magnetosphere

,”

J. Geophys. Res.: Space Phys.

123

,

2815

–

2830

, https://doi.org/10.1002/2017JA024761 (

2018

).

Google Scholar

Crossref

4.

R.

Sahu

,

A. R.

Mansour

, and

K.

Hara

, “

Full fluid moment model for low temperature magnetized plasmas

,”

Phys. Plasmas

27

,

113505

(

2020

).

https://doi.org/10.1063/5.0021474

Google Scholar

Crossref

5.

L.

Tonks

and

I.

Langmuir

, “

Oscillations in ionized gases

,”

Phys. Rev.

33

,

195

(

1929

).

https://doi.org/10.1103/PhysRev.33.195

Google Scholar

Crossref

6.

K.

Hara

, “

An overview of discharge plasma modeling for Hall effect thrusters

,”

Plasma Sources Sci. Technol.

28

,

044001

(

2019

).

https://doi.org/10.1088/1361-6595/ab0f70

Google Scholar

Crossref

7.

M.

Lieberman

and

A.

Lichtenberg

,

Principles of Plasma Discharges and Materials Processing

(

John Wiley & Sons

,

2005

).

Google Scholar

Crossref

8.

V. M.

Zhdanov

,

Transport Processes in Multicomponent Plasma

(

CRC Press

,

2002

).

Google Scholar

Crossref

9.

G.

Hagelaar

and

L. C.

Pitchford

, “

Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models

,”

Plasma Sources Sci. Technol.

14

,

722

(

2005

).

https://doi.org/10.1088/0963-0252/14/4/011

Google Scholar

Crossref

10.

M.

Bowden

,

T.

Okamoto

,

F.

Kimura

,

H.

Muta

,

K.

Uchino

,

K.

Muraoka

,

T.

Sakoda

,

M.

Maeda

,

Y.

Manabe

,

M.

Kitagawa

et al, “

Thomson scattering measurements of electron temperature and density in an electron cyclotron resonance plasma

,”

J. Appl. Phys.

73

,

2732

–

2738

(

1993

).

https://doi.org/10.1063/1.353046

Google Scholar

Crossref

11.

B.

Vincent

,

S.

Tsikata

, and

S.

Mazouffre

, “

Incoherent Thomson scattering measurements of electron properties in a conventional and magnetically-shielded hall thruster

,”

Plasma Sources Sci. Technol.

29

,

035015

(

2020

).

https://doi.org/10.1088/1361-6595/ab6c42

Google Scholar

Crossref

12.

J.

Simmonds

and

Y.

Raitses

, “

Ion acceleration in a wall-less Hall thruster

,”

J. Appl. Phys.

130

,

093302

(

2021

).

https://doi.org/10.1063/5.0062607

Google Scholar

Crossref

13.

A. H.

Hakim

, “

Extended MHD modelling with the ten-moment equations

,”

J. Fusion Energy

27

,

36

–

43

(

2008

).

https://doi.org/10.1007/s10894-007-9116-z

Google Scholar

Crossref

14.

D. A.

Kuldinow

,

Y.

Yamashita

,

A. R.

Mansour

, and

K.

Hara

, “

Ten-moment fluid model with heat flux closure for gasdynamic flows

,”

J. Comput. Phys.

508

,

113030

(

2024

).

https://doi.org/10.1016/j.jcp.2024.113030

Google Scholar

Crossref

15.

R.

Herdan

and

B.

Liley

, “

Dynamical equations and transport relationships for a thermal plasma

,”

Rev. Mod. Phys.

32

,

731

(

1960

).

https://doi.org/10.1103/RevModPhys.32.731

Google Scholar

Crossref

16.

S. T.

Miller

and

U.

Shumlak

, “

A multi-species 13-moment model for moderately collisional plasmas

,”

Phys. Plasmas

23

,

082303

(

2016

).

https://doi.org/10.1063/1.4960041

Google Scholar

Crossref

17.

T.

Ruggeri

,

Extended Thermodynamics

, Springer Tracts in Natural Philosophy (

Springer

,

New York

,

1993

).

Google Scholar

18.

S.

Boccelli

,

F.

Giroux

,

T. E.

Magin

,

C.

Groth

, and

J. G.

McDonald

, “

A 14-moment maximum-entropy description of electrons in crossed electric and magnetic fields

,”

Phys. Plasmas

27

,

123506

(

2020

).

https://doi.org/10.1063/5.0025651

Google Scholar

Crossref

19.

L.

Pareschi

and

G.

Russo

, “

Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator

,”

SIAM J. Numer. Anal.

37

,

1217

–

1245

(

2000

).

https://doi.org/10.1137/S0036142998343300

Google Scholar

Crossref

20.

S.

Chapman

and

T. G.

Cowling

,

The Mathematical Theory of Non-Uniform Gases: An Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases

(

Cambridge University Press

,

1990

).

Google Scholar

21.

A. A.

Laguna

,

B.

Esteves

,

J.

Raimbault

,

A.

Bourdon

, and

P.

Chabert

, “

Discussion on the transport processes in electrons with non-Maxwellian energy distribution function in partially-ionized plasmas

,”

Plasma Phys. Controlled Fusion

65

,

054002

(

2023

).

https://doi.org/10.1088/1361-6587/acc422

Google Scholar

Crossref

22.

A.

Alvarez Laguna

,

B.

Esteves

,

A.

Bourdon

, and

P.

Chabert

, “

A regularized high-order moment model to capture non-Maxwellian electron energy distribution function effects in partially ionized plasmas

,”

Phys. Plasmas

29

,

083507

(

2022

).

https://doi.org/10.1063/5.0095019

Google Scholar

Crossref

23.

P. L.

Bhatnagar

,

E. P.

Gross

, and

M.

Krook

, “

A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems

,”

Phys. Rev.

94

,

511

(

1954

).

https://doi.org/10.1103/PhysRev.94.511

Google Scholar

Crossref

24.

W.

Vincenti

and

C.

Kruger

,

Introduction to Physical Gas Dynamics

(

Wiley

,

New York

,

1965

).

Google Scholar

25.

I. D.

Boyd

and

T. E.

Schwartzentruber

,

Nonequilibrium Gas Dynamics and Molecular Simulation

(

Cambridge University Press

,

2017

), Vol.

42

.

Google Scholar

Crossref

26.

G.

Strang

, “

On the construction and comparison of difference schemes

,”

SIAM J. Numer. Anal.

5

,

506

–

517

(

1968

).

https://doi.org/10.1137/0705041

Google Scholar

Crossref

27.

J. L.

Steger

and

R.

Warming

, “

Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods

,”

J. Comput. Phys.

40

,

263

–

293

(

1981

).

https://doi.org/10.1016/0021-9991(81)90210-2

Google Scholar

Crossref

28.

B.

Van Leer

, “

Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method

,”

J. Comput. Phys.

32

,

101

–

136

(

1979

).

https://doi.org/10.1016/0021-9991(79)90145-1

Google Scholar

Crossref

29.

Y.

Yamashita

,

R.

Lau

, and

K.

Hara

, “

Inertial and anisotropic pressure effects on cross-field electron transport in low-temperature magnetized plasmas

,”

J. Phys. D: Appl. Phys.

56

,

384003

(

2023

).

https://doi.org/10.1088/1361-6463/acdb83

Google Scholar

Crossref

30.

K.

Hara

,

I. D.

Boyd

, and

V. I.

Kolobov

, “

One-dimensional hybrid-direct kinetic simulation of the discharge plasma in a Hall thruster

,”

Phys. Plasmas

19

,

113508

(

2012

).

https://doi.org/10.1063/1.4768430

Google Scholar

Crossref

31.

W.

Huang

,

A. D.

Gallimore

, and

R. R.

Hofer

, “

Neutral flow evolution in a six-kilowatt Hall thruster

,”

J. Propul. Power

27

,

553

–

563

(

2011

).

https://doi.org/10.2514/1.B34048

Google Scholar

Crossref

32.

J.

Boeuf

and

L.

Garrigues

, “

Low frequency oscillations in a stationary plasma thruster

,”

J. Appl. Phys.

84

,

3541

–

3554

(

1998

).

https://doi.org/10.1063/1.368529

Google Scholar

Crossref

33.

S.

Biagi

, “

Fortran program

,” MAGBOLTZ v.8.97,

2011

.

Google Scholar

34.

M.

Hayashi

, “

Bibliography of electron and photon cross sections with atoms and molecules published in the 20th century. Xenon

,”

Technical Report No. NIFS-DATA-79

(

National Institute for Fusion Science

,

Nagoya

,

2003

).

Google Scholar

35.

J. D.

Huba

,

NRL Plasma Formulary

(

Naval Research Laboratory

,

1998

).

Google Scholar

36.

K.

Hara

, “

Non-oscillatory quasineutral fluid model of cross-field discharge plasmas

,”

Phys. Plasmas

25

,

123508

(

2018

).

https://doi.org/10.1063/1.5055750

Google Scholar

Crossref

37.

J. W.

Koo

and

I. D.

Boyd

, “

Modeling of anomalous electron mobility in Hall thrusters

,”

Phys. Plasmas

13

,

033501

(

2006

).

https://doi.org/10.1063/1.2172191

Google Scholar

Crossref

38.

D.

Sydorenko

, “

Particle-in-cell simulations of electron dynamics in low pressure discharges with magnetic fields

,” Ph.D. thesis (

University of Saskatchewan

,

2006

).

Google Scholar

39.

M.

Mitchner

and

C. H.

Kruger

, Jr.,

Partially Ionized Gases

(

John Wiley and Sons, Inc

.,

New York

,

1973

).

Google Scholar

40.

K.

Hara

,

M. J.

Sekerak

,

I. D.

Boyd

, and

A. D.

Gallimore

, “

Mode transition of a Hall thruster discharge plasma

,”

J. Appl. Phys.

115

,

203304

(

2014

).

https://doi.org/10.1063/1.4879896

Google Scholar

Crossref

41.

S.

Braginskii

, “

Transport processes in a plasma

,” in

Reviews of Plasma Physics

(

Consultants Bureau

,

New York

,

1965

), Vol.

1

, p.

205

.

Google Scholar

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International (CC BY-NC-ND) license (https://creativecommons.org/licenses/by-nc-nd/4.0/).

This work is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License .

Ten-moment fluid model for low-temperature magnetized plasmas

I. INTRODUCTION

II. GOVERNING EQUATIONS

A. The fluid approximation

1. The 10-moment equations

2. The 5-moment equations

B. The particle-in-cell method

C. Collisions

D. Heat flux closure model

E. Field equations

III. NUMERICAL METHODS

A. Timestepping

B. Kinetic flux boundary condition for the 10-moment model

C. Poisson's equation

IV. THE HALL-EFFECT THRUSTER MODEL

A. Boundary conditions

B. Magnetic field

C. Collisions

1. Collision types

2. Collision model

3. Heat flux model

V. RESULTS

A. Plasma profiles

B. Cross-field electron transport

1. Mechanism of anisotropy

2. Drift velocity decomposition

C. Importance of self-consistent heat flux closure

VI. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX A: JUSTIFICATION OF THE FORM OF THE 10-MOMENT HEAT FLUX MODEL

APPENDIX B: 10-MOMENT KINETIC FLUXES

1. 0th velocity moment

2. 1st velocity moments

3. 2nd velocity moments

APPENDIX C: PRIMITIVE COLLISION QUANTITIES

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only