Combining the hybrid mimetic mixed method with the Scharfetter-Gummel scheme for magnetised transport in plasmas

1. Introduction

Magnetostatic fields play an important role in plasma physics. One useful property of a magnetic field is that it confines the plasma, limiting charged particle loss to the walls; hence, allowing magnetron discharges to operate at relatively low voltages. Some examples of these can be seen in Hall-effect thrusters (HETs) and electron-cyclotron-resonance (ECR) discharges [13]. These physical processes are expensive to set-up, and so we turn to mathematical models and perform numerical simulations in order to get an overview and general understanding of how particles behave in these situations.

Early numerical studies for plasmas involve particle-based models, which involve Monte Carlo simulations [3, 8]. However, the inclusion of a magnetic field complicates the physics and hence the numerical modeling of magnetised discharges. In this paper, we consider fluid models. These consist of drift-diffusion equations for the particles (electrons and ions), coupled to the Poisson equation for the electric field, see e.g.: [6, 13, 14, 17, 23]. This is an approximation for the real physics which is computationally more efficient and avoids some of the numerical complications and high computational costs associated to using Monte Carlo type methods on purely particle-based models [13].

The model system of equations is highly nonlinear due to the coupling between the Poisson equation and the drift-diffusion equations. In order to resolve this, we decouple the equations by using the transversal method of lines, where a semi-discrete system is obtained after discretising with respect to the time variable. Following this, we look at the spatial discretisation. The Poisson equation is elliptic, and can thus be solved by classical techniques, such as the central difference method. In the absence of a magnetic field, the drift-diffusion equations can be solved efficiently and precisely by the Scharfetter-Gummel (SG) scheme [25]. On the other hand, the presence of a magnetic field introduces anisotropy in the drift-diffusion equations; hence requiring a numerical scheme which captures the anisotropic diffusion precisely. Recent studies propose the use of high order finite difference gradient reconstruction methods, together with a magnetic field aligned mesh [24, 29] in order to resolve this.

The aim of this paper is to propose a numerical scheme which is computationally efficient, and is able to handle the anisotropy induced by the magnetic field. In order to do so, we propose, as in [2], to separate the discretisation of the diffusive and advective fluxes. We then use the hybrid-mimetic-mixed (HMM) method [10]

for discretising the diffusive flux, and a modification of the SG scheme for the advective flux. The main interest of considering the HMM method is its ability to handle the anisotropic diffusion tensor. Moreover, the method can be fully hybridised, leading to a purely local stencil for solving the system of equations, which allows a very efficient implementation of the scheme. As a remark, the HMM can also be viewed as a gradient reconstruction method or a gradient scheme

[9], and can thus be used in combination with the methods proposed in [24, 29].

The outline of this paper is as follows. We start by presenting the mathematical model for magnetised transport in Section 2. This model consists of a set of drift-diffusion equations, coupled to the Poisson equation. We then describe in Section 3 the choice of time discretisation, and how to decouple the system. Following this, we employ a finite volume type spatial discretisation in Section 4, using the HMM for the diffusive fluxes, and a modified SG for the advective fluxes. Details of the implementation are then given in Section 5. Numerical tests are provided in Section 6 to show the effectiveness of the proposed method. Finally, in Section 7 we provide a summary of the results and some recommendations for future research.

2. Mathematical model

Suppose that an electric field $E$ is applied between two parallel electrodes and that an electron flux is forced through a uniform medium, where an oblique magnetic field $B$ is applied. The electric field $E = - \nabla V$ in a spatial domain $Ω$ over a time interval $(0, T)$ is governed by the Poisson equation

	$\nabla \cdot (ε E) = ρ o n Ω \times (0, T),$	(1a)
where $ε$ is the permittivity of vacuum, $V$ is the potential, and $ρ$ is the space charge density given by
	$ρ = q (n_{i} - n_{e}) .$
Here, $q$ is the elementary charge, $n_{i}$ and $n_{e}$ are the ion and electron densities, respectively. Under the assumption that ions are not magnetised, the model which describes the evolution of the ion and electron densities is then given by the drift-diffusion equations:

	$\frac{\partial n_{i}}{\partial t} + \nabla \cdot Γ_{i} - k n_{e, l o c} n_{g} = 0 o n Ω \times (0, T),$		(1b)
	$\frac{\partial n_{e}}{\partial t} + \nabla \cdot M_{e} Γ_{e} - k n_{e, l o c} n_{g} = 0 o n Ω \times (0, T),$		(1c)

where $Γ_{i}$ and $Γ_{e}$ are the drift-diffusion flux densities for the ions and electrons, respectively, defined by

	$Γ_{i}$	$:= μ_{i} E n_{i} - D_{i} \nabla n_{i},$
	$Γ_{e}$	$:= μ_{e} E n_{e} - D_{e} \nabla n_{e} .$

Equations (1b) and (1c) are coupled to the Poisson equation (1a) due to the electric field $E$ and the space charge density $ρ$ . The parameters and constants used in the model are enumerated in Tables 1 and 2. For quantities which are common to both ions and electrons, we use the subscript $p = i, e$ , referring to a particle (either an ion or an electron).

$k$	ionisation rate
$μ_{p}$	mobility coefficient
$D_{p}$	diffusion coefficient
$M_{e}$	magnetic tensor
$n_{e, l o c}$	local electron density

Table 1. Parameters in the model.

$n_{g}$	density of the background gas
$¯ ϵ = 2$ eV	mean electron energy
$ε$	permittivity of vacuum
$k_{B}$	Boltzmann constant
$q$	elementary charge
$T_{g}$	temperature of the background gas

Table 2. Constants in the model.

We now discuss the parameters used in the model. Firstly, the ionisation rate $k$ and the mobility coefficients $μ_{p}$ are computed via local field approximations, and are hence functions of the reduced electric field $∥ E ∥ / n_{g}$ . In this paper, we consider helium gas, whose ion mobility $μ_{i}$ (Figure 2, left) is obtained from [11]. The ionisation rate $k$ (Figure 1) and electron mobility $μ_{e}$ (Figure 2, right) were computed in [15] using cross-sectional data for helium gas in [7, 20] for low energies and [16] for high energies. Here and throughout the paper, $∥ \cdot ∥$ refers to the $L^{2}$ -norm.

Figure 2. Mobility coefficients (Left: $μ_{i}$ , Right: $μ_{e}$ ).

Secondly, we discuss the diffusion coefficients. These are governed by the Stokes-Einstein relation

D_{p} = \frac{k_{B} T_{p} μ_{p}}{q_{p}},

where, $q_{i} = q$ and $q_{e} = - q$ . For ions, we assume that $T_{i} = T_{g}$ and so we have

	$D_{i} = \frac{k_{B} T_{g} μ_{i}}{q} .$	(2a)
We assume that electrons follow a Maxwellian distribution so that $k_{B} T_{e} = 2 / 3 ¯ ϵ$ and thus
	$D_{e} = - ¯ ϵ \frac{2 μ_{e}}{3 q} .$	(2b)

Thirdly, the the local electron density $n_{e, l o c}$ is computed by the relation [15]

n_{e, l o c} = \frac{∥ M_{e} Γ_{e} ∥}{μ_{e} ∥ M_{e} E ∥} .

(3)

Finally, the magnetic tensor $M_{e}$ is given by


	$M_{e}$
		$= \frac{1}{∥ β ∥^{2} + 1} (I + β β^{T} - C (β)),$	(4a)
	$C (β)$	$= ⎡ ⎢ ⎣ \begin{matrix} 0 & - β_{3} & β_{2} β_{3} & 0 & - β_{1} - β_{2} & β_{1} & 0 \end{matrix} ⎤ ⎥ ⎦,$
with Hall parameters $β_{i} = μ_{e} B_{i}$ , $i = 1, 2, 3$ and magnetic field $B$ . We now note that $C (β)$ is a cross product matrix such that for any vector $v \in R^{3}$ , $C (β) v = β \times v$ . Taking note that

	$M_{e} β$	$= \frac{1}{∥ β ∥^{2} + 1} (β + β (β^{T} β) - β \times β)$
		$= β,$
we then see that $β$ is an eigenvector of $M_{e}$ corresponding to the eigenvalue 1. Essentially, this means that $M_{e}$ preserves the component of the vector that is parallel to the magnetic field.

For dimension $d = 2$ , the magnetic field can be specified at an angle $θ$ with respect to the $x$ -axis by setting

	$B$	$= [B_{1}, B_{2}, 0]^{T}$		(4b)
		$= ∥ B ∥ [cos θ, sin θ, 0]^{T},$		(4b)

where $∥ B ∥$ is the strength of the magnetic field. Considering now a vector

β^{⊥} = μ_{e} ∥ B ∥ [- sin θ, cos θ, 0]^{T}

orthogonal to the magnetic field, we see that

M_{e} β^{⊥} = \frac{1}{∥ β ∥^{2} + 1} (β^{⊥} - ∥ β ∥^{2} e_{z}) .

For two-dimensional models, we neglect the $z$ -component, and hence we see that the components orthogonal to the magnetic field are scaled by a factor $1 / (∥ β ∥^{2} + 1)$ . These properties agree with the classical definition [18] of a magnetic tensor.

2.1. Boundary conditions

We now describe the boundary conditions for the model. In the electrodes (corresponding to either anodes or cathodes), Dirichlet boundary conditions are imposed for the Poisson equation (1a),


	$V$	$= V_{a} a n o d e,$	(5a)
	$V$	$= V_{c} c a t h o d e .$	(5a)
Then, flux boundary conditions are imposed for the drift-diffusion equations (1b) and (1c), i.e.,

	$Γ_{i} \cdot n$	$= (μ_{i} E \cdot n)^{+} n_{i} + \frac{1}{4} v_{t h, i} n_{i},$	(5b)
	$(M_{e} Γ_{e}) \cdot n$	$= (μ_{e} M_{e} E \cdot n)^{+} n_{e} + \frac{1}{4} v_{t h, e} n_{e} - γ_{i} Γ_{i} \cdot n,$	(5c)
where $v_{t h, p}$ is the thermal velocity, $f^{+} = max (f, 0)$ denotes the positive part of the function, and $n$ is the outward unit normal vector at the boundary. The first two terms in the right hand side of equations (5b) and (5c) describe the loss of electrons and ions, respectively, through the wall via the drift and thermal fluxes. On the other hand, the third term in (5c) is the secondary electron emission, which describes the production of electrons as a result of ion fluxes bombarding the wall [4], where $γ_{i}$ is the average number of electrons emitted per incident ion.

Finally, homogeneous Neumann boundary conditions are imposed over the other boundaries of the domain. That is, for (1a), (1b), and (1c) we impose:

$\nabla V \cdot n$	$= 0,$	(5d)
$\nabla n_{i} \cdot n$	$= 0,$
$\nabla n_{e} \cdot n$	$= 0.$

For the drift-diffusion equations, this means that there is no diffusion of ions and electrons across these boundaries. For the Poisson equation, this means that the electric field lines are parallel to the boundary.

Given the electron and ion densities $n_{e} (x, 0), n_{i} (x, 0)$ at time $t = 0$ , we are then interested in calculating the densities $n_{e} (x, T), n_{i} (x, T)$ at time $t = T$ . In this work, we consider a numerical scheme based on the transversal method of lines. That is, we start by presenting a semi-discrete version of the model system of equations with respect to the time variable. Following this, we then describe the discretisation in space.

3. Time discretisation

We form a partition of the time interval $(0, T)$ by taking $0 = t^{(0)} < t^{(1)} < \dots t^{(M)} = T$ . In order to avoid issues with the stability of the numerical solution, we use the implicit Euler time integration. Using a uniform time step of $Δ t = t^{(m + 1)} - t^{(m)}$ , for $m = 0, \dots M - 1$ , and given the densities $n_{e}^{(m)} \approx n_{e} (x, t^{(m)}), n_{i}^{(m)} \approx n_{i} (x, t^{(m)})$ , we seek $E^{(m + 1)}, n_{e}^{(m + 1)}, n_{i}^{(m + 1)}$ such that

(6a)

\frac{n_{i}^{(m + 1)} - n_{i}^{(m)}}{Δ t} + \nabla \cdot Γ_{i}^{(m + 1)} - k^{(m + 1)} n_{e, l o c}^{(m + 1)} n_{g} = 0,

(6b)

\frac{n_{e}^{(m + 1)} - n_{e}^{(m)}}{Δ t} + \nabla \cdot M_{e}^{(m + 1)} Γ_{e}^{(m + 1)} - k^{(m + 1)} n_{e, l o c}^{(m + 1)} n_{g} = 0.

(6c)

Here,

Γ_{p}^{(m + 1)} = μ_{p}^{(m + 1)} E^{(m + 1)} n_{p}^{(m + 1)} - D_{p}^{(m + 1)} \nabla n_{p}^{(m + 1)}

is the drift-diffusion flux density of particle $p$ (ions or electrons) at time $t^{(m + 1)}$ . As a remark, we note that even though the scheme was presented with a uniform time stepping, it can be easily adapted into a scheme with an adaptive time stepping $Δ t^{(m)} = t^{(m + 1)} - t^{(m)}$ .

We note that due to the electric field, the Poisson equation (6a) is coupled to the continuity equations (6b) and (6c). We now describe a few ways for decoupling the system of equations (6a)-(6c).

3.1. First-order approximation of the charge density

The simplest way to decouple the Poisson equation from the continuity equations is to use a first-order approximation of the charge density. This leads to seeking, for $m = 0, \dots M - 1$ , $E^{(m + 1)}$ such that

\nabla \cdot (ε E^{(m + 1)}) = q (n_{i}^{(m)} - n_{e}^{(m)}) .

(7)

The advantage of using (7) is that the charge density will then only depend on the electron and ion densities at the previous time step. On the other hand, this is a first-order approximation, which means that the leading order of error in time is $O (Δ t)$ . Hence, a very small time step, at most equal to the dielectric relaxation time [27, Section II.C.], is needed in order to ensure that the electric field does not change sign during one time step, and that the approximation of the electric field is good enough to be used for computing the electron and ion densities in the drift-diffusion equations.

3.2. Second-order approximation of the charge density

We now present a second-order approximation by performing a Taylor expansion on the charge density. Centering the Taylor expansion at time $t^{(m + 1)}$ , we may write

ρ (x, t^{(m + 1)}) = ρ (x, t^{(m)}) + Δ t \frac{\partial ρ (x, t^{(m + 1)})}{\partial t} + O (Δ t^{2}) .

(8)

We then note that the time derivative $\frac{\partial ρ}{\partial t}$ can be computed from the drift-diffusion equations, and can be written as

\frac{\partial ρ}{\partial t} = - q (\nabla \cdot Γ_{i} - \nabla \cdot M_{e} Γ_{e}) .

(9)

Using now (8) for the charge density in the Poisson equation and evaluating $\frac{\partial ρ}{\partial t}$ in (9) at time $t^{(m + 1)}$ , we then seek, for $m = 0, \dots M - 1$ , $E^{(m + 1)}$ such that

	$\nabla \cdot (ε E^{(m + 1)}) =$	$q (n_{i}^{(m)} - n_{e}^{(m)}) - q Δ t \nabla \cdot Γ_{i}^{(m + 1)}$		(10)
		$+ q Δ t \nabla \cdot (M_{e}^{(m + 1)} Γ_{e}^{(m + 1)}) .$		(10)

In this case, the quantities $E, M_{e}, μ_{i}, μ_{e}, D_{i}, D_{e}, n_{i},$ and $n_{e}$ at time $t^{(m + 1)}$ are unknown and hence the Poisson equation is still coupled to the drift-diffusion equations. One way to deal with this is to consider a semi-implicit treatment of these quantities [27, 28]. That is, we keep $E$ implicit and take first-order approximations to $M_{e}, μ_{i}, μ_{e}, D_{i}, D_{e}, n_{i}, n_{e}$ , which leads to

	$\nabla \cdot (ε E^{(m + 1)}) =$	$q (n_{i}^{(m)} - n_{e}^{(m)}) - q Δ t \nabla \cdot (μ_{i}^{(m)} E^{(m + 1)} n_{i}^{(m)} - D_{i}^{(m)} \nabla n_{i}^{(m)})$
		$+ q Δ t \nabla \cdot M_{e}^{(m)} (μ_{e}^{(m)} E^{(m + 1)} n_{e}^{(m)} - D_{e}^{(m)} \nabla n_{e}^{(m)}),$

or equivalently,

	$\nabla \cdot ((ε I + q Δ t$
		$= q (n_{i}^{(m)} - n_{e}^{(m)}) + q Δ t \nabla \cdot (D_{i}^{(m)} \nabla n_{i}^{(m)})$
		$- q Δ t \nabla \cdot (M_{e}^{(m)} D_{e}^{(m)} \nabla n_{e}^{(m)}) .$

Denoting the correction flux at time level $m + 1$ by

Γ_{c}^{(m + 1)} := (μ_{i}^{(m)} n_{i}^{(m)} I - μ_{e}^{(m)} n_{e}^{(m)} M_{e}^{(m)}) E^{(m + 1)},

(11)

the above equation can be rewritten as

	$\nabla \cdot (ε E^{(m + 1)}) + q Δ t \nabla \cdot Γ_{c}^{(m + 1)}$	$= q (n_{i}^{(m)} - n_{e}^{(m)}) + q Δ t \nabla \cdot (D_{i}^{(m)} \nabla n_{i}^{(m)})$		(12)
		$- q Δ t \nabla \cdot (M_{e}^{(m)} D_{e}^{(m)} \nabla n_{e}^{(m)}) .$		(12)

4. Spatial discretisation

For the spatial discretisation, we start by forming $Ω_{h}$ , a partition of the domain $Ω$ into $N_{K}$ control volumes. Each of the control volumes $K$ has edges $σ \in E_{K}$ such that ${max}_{σ \in E_{K}, K \in Ω_{h}} | σ | = h$ , and we denote by $N_{σ}$ the total number of edges contained in this partition. We then denote by $X_{D}$ the collection of unknowns (one on each cell, and one on each edge). For this work, we focus on finite volume type schemes.

4.1. Finite volume discretisation

To write the equations (6a),(6b), and (6c) in their finite volume forms, we take their integrals over all control volumes $K \in Ω_{h}$ . This results to, in each control volume $K$ ,

\int_{K} \nabla \cdot (ε E^{(m + 1)}) d A = q \int_{K} (n_{i}^{(m + 1)} - n_{e}^{(m + 1)}) d A,

\frac{1}{Δ t} \int_{K} (n_{i}^{(m + 1)} - n_{i}^{(m)}) d A + \int_{K} \nabla \cdot Γ_{i}^{(m + 1)} d A - n_{g} \int_{K} k^{(m + 1)} n_{e, l o c}^{(m + 1)} d A = 0,

\frac{1}{Δ t} \int_{K} (n_{e}^{(m + 1)} - n_{e}^{(m)}) d A + \int_{K} \nabla \cdot M_{e}^{(m + 1)} Γ_{e}^{(m + 1)} d A - n_{g} \int_{K} k^{(m + 1)} n_{e, l o c}^{(m + 1)} d A = 0,

respectively. By applying the divergence theorem and using $E^{(m + 1)} = - \nabla V^{(m + 1)}$ , these can be rewritten as

$- \sum σ \in E_{K} \int_{σ} ε \nabla V^{(m + 1)} \cdot n_{K, σ} d s = q \int_{K} (n_{i}^{(m + 1)} - n_{e}^{(m + 1)}) d A,$		(13a)

$\frac{1}{Δ t} \int_{K} (n_{i}^{(m + 1)} - n_{i}^{(m)}) d A$	$+ \sum σ \in E_{K} \int_{σ} Γ_{i}^{(m + 1)} \cdot n_{K, σ} d s$	(13b)
		(13b)

$\frac{1}{Δ t} \int_{K} (n_{e}^{(m + 1)} - n_{e}^{(m)}) d A$	$+ \sum σ \in E_{K} \int_{σ} M_{e}^{(m + 1)} Γ_{e}^{(m + 1)} \cdot n_{K, σ} d s$	(13c)
		(13c)

Here, $n_{K, σ}$ is the outward unit normal vector of $K$ at $σ$ . Due to the definition (3) of $n_{e, l o c}$ , the term $n_{e, l o c}^{(m + 1)}$ is nonlinear in $n_{e}$ . Here, we choose a first-order approximation and take the local electron density at time $m$ so that

n_{e, l o c}^{(m)} := \frac{∥ M_{e}^{(m)} Γ_{e}^{(m)}) ∥}{μ_{e}^{(m)} ∥ M_{e}^{(m)} E^{(m)} ∥} .

(14)

Key to the definition of a finite volume scheme is the choice of how to discretise the fluxes. In order to be able to treat the anisotropy, we propose to split the fluxes into a diffusive and an advective component. Diffusive fluxes are of the form $- \int_{σ} Λ \nabla c \cdot n_{K, σ} d s$ for some diffusion tensor $Λ$ and scalar $c$ ; on the other hand, the advective fluxes take the form $\int_{σ} c V \cdot n_{K, σ} d s$ , where $V$ is a velocity field. In equations (13a)-(13c), three diffusive fluxes and two advective fluxes are present. We then denote the discrete diffusive and advective fluxes for each edge $σ$ in a cell $K \in Ω_{h}$ by

	$F_{K, σ}^{D_{V}} \approx - \int_{σ} ε \nabla V \cdot n_{K, σ} d s,$		(15a)
	$F_{K, σ}^{D_{i}} \approx - \int_{σ} D_{i} \nabla n_{i} \cdot n_{K, σ} d s,$		(15b)
	$F_{K, σ}^{D_{e}} \approx - \int_{σ} M_{e} D_{e} \nabla n_{e} \cdot n_{K, σ} d s,$		(15c)
	$F_{K, σ}^{A_{i}} \approx \int_{σ} n_{i} μ_{i} E \cdot n_{K, σ} d s,$		(15d)
	$F_{K, σ}^{A_{e}} \approx \int_{σ} n_{e} μ_{e} M_{e} E \cdot n_{K, σ} d s .$		(15e)

Replacing the fluxes in (13a)-(13c) with their discrete counterparts (15a)-(15e) at the appropriate time level and under the assumption that $n_{i}^{(m)}$ and $n_{e}^{(m)}$ are piecewise constant with values $n_{i, K}^{(m)}$ and $n_{e, K}^{(m)}$ respectively in a cell $K$ , we then obtain our numerical scheme. That is, for $m = 0, \dots, M$ , given the densities $n_{e}^{(m)} \approx n_{e} (x, t^{(m)}), n_{i}^{(m)} \approx n_{i} (x, t^{(m)})$ , we seek $E^{(m + 1)}, n_{e}^{(m + 1)}, n_{i}^{(m + 1)}$ such that for each $K \in Ω_{h}$ ,

	$\sum σ \in E_{K} F_{K, σ}^{D_{V}, (m + 1)} = q \| K \| (n_{i, K}^{(m)} - n_{e, K}^{(m)}),$	(16a)
	$\| K \| \frac{n_{i, K}^{(m + 1)}}{Δ t} + \sum σ \in E_{K} F_{K, σ}^{D_{i}, (m + 1)} + \sum σ \in E_{K} F_{K, σ}^{A_{i}, (m + 1)} = \| K \| \frac{n_{i, K}^{(m)}}{Δ t} + n_{g} \| K \| k_{K}^{(m + 1)} n_{e, l o c, K}^{(m)},$	(16b)
	$\| K \| \frac{n_{e, K}^{(m + 1)}}{Δ t} + \sum σ \in E_{K} F_{K, σ}^{D_{e}, (m + 1)} + \sum σ \in E_{K} F_{K, σ}^{A_{e}, (m + 1)} = \| K \| \frac{n_{e, K}^{(m)}}{Δ t} + n_{g} \| K \| k_{K}^{(m + 1)} n_{e, l o c, K}^{(m)},$	(16c)
where $k_{K}^{(m + 1)}$ is the average ionisation rate in cell $K$ at time $t^{(m + 1)}$ , and $n_{e, l o c, K}^{(m)}$ is the value of $n_{e, l o c}^{(m)}$ computed at cell $K$ .

Starting with $m \geq 1$ , we then use the second-order approximation for the charge density in solving the Poisson equation. That is, instead of using (16a), we use a semi-implicit discretisation based on (12), given by

	$\sum σ \in E_{K} F_{K, σ}^{D_{V}, (m + 1)}$	$+ q Δ t \sum σ \in E_{K} F_{K, σ}^{C_{V}, (m + 1)}$		(16d)
		$= q \| K \| (n_{i, K}^{(m)} - n_{e, K}^{(m)}) + q Δ t (\sum σ \in E_{K} (F_{K, σ}^{D_{e}, (m)} - F_{K, σ}^{D_{i}, (m)})),$		(16d)

where

F_{K, σ}^{C_{V}, (m + 1)} \approx - \int_{σ} Γ_{c}^{(m + 1)} \cdot n_{K, σ} d s

is a discretisation of the semi-implicit correction flux $Γ_{c}$ defined in (11).

We note here that we cannot immediately use (16d) at the first time step due to the yet unknown mobility coefficients, correction and diffusive fluxes at time $t = 0$ . Hence, we choose to use (16a) for solving the Poisson equation at the first time step when $m = 0$ .

4.2. Boundary conditions and conservation of fluxes

We now see that there are $N_{K}$ equations (corresponding to the number of cells) involved in each of (16a)-(16d). However, we have $N_{K} + N_{σ}$ unknowns for each of these equations: $N_{K}$ corresponding to the number of cells, and $N_{σ}$ corresponding to the number of edges. This is standard for hybrid finite volume methods, and the equations for the $N_{σ}$ unknowns along the edges are obtained by first imposing conservation of fluxes for interior edges. That is, if $σ$ is an edge shared by two cells $K$ and $L$ , we should have that the total (diffusive and advective) flux $F_{K, σ}$ going from $K$ to $L$ via $σ$ is equal to the negative total flux $F_{L, σ}$ going from $L$ to $K$ via $σ$ , i.e.

(F_{K, σ}^{D} + F_{K, σ}^{A}) + (F_{L, σ}^{D} + F_{L, σ}^{A}) = 0.

(17)

Finally, we impose boundary conditions on the remaining edges, which are located along the boundary of the domain. The Dirichlet and Neumann boundary conditions are straightforward to impose, so we only describe here how to impose the flux boundary conditions (5b) and (5c). This is done by setting, for the corresponding boundary edges $σ \in E_{K}$ ,

$F_{K, σ}^{D_{i}, (m + 1)} + F_{K, σ}^{A_{i}, (m + 1)} = \| σ \| ((μ_{i}^{(m + 1)} E^{(m + 1)} \cdot n_{K, σ})^{+} n_{i, σ}^{(m + 1)} + \frac{1}{4} v_{t h, i} n_{i, σ}^{(m + 1)}) .$		(18a)

$F_{K, σ}^{D_{e}, (m + 1)} + F_{K, σ}^{A_{e}, (m + 1)} =$		(18b)
	$- γ_{i} (F_{K, σ}^{D_{i}, (m + 1)} + F_{K, σ}^{A_{i}, (m + 1)}) .$	(18b)

Here, $n_{e, σ}$ and $n_{i, σ}$ are the values (to be determined) of $n_{e}$ and $n_{i}$ at $σ$ , respectively. To complete the definition of the numerical scheme, we are then left to describe the choice of how to compute the discrete diffusive and advective fluxes.

4.3. Diffusive fluxes

In this section, we discuss the discretisation of the diffusive fluxes $F_{K, σ}^{D} \approx - \int_{σ} Λ \nabla c \cdot n_{K, σ} d s$ . Here, we choose to use the hybrid mimetic mixed (HMM) method, due to its ability to handle anisotropic diffusion tensors on generic meshes, see e.g. [10, 12]. For this paper, we only present the HMM method for Cartesian meshes, since the tests and models are performed on square cells. Denoting by $X_{D}$ the collection of all discrete values, one on each cell, and one on each face (see Figure 3, left), we then define the diffusive fluxes by the relation:

		$\forall K \in Ω_{h}, \forall v \in X_{D},$		(19)
		$\sum σ \in E_{K} F_{K, σ}^{D} (v_{K} - v_{σ}) = \int_{K} Λ_{K} \nabla_{D} c (x) \cdot \nabla_{D} v (x) d A,$		(19)

where $\nabla_{D}$ is a discrete gradient that is defined as follows: For $K \in Ω_{h}$ , we set

	$\forall w \in X_{D}, \forall x \in K,$
	$\nabla_{D} w (x)$	$= {¯ ¯¯¯ ¯ \nabla}_{K} w + S_{K} w (x),$		(20)

	${¯ ¯¯¯ ¯ \nabla}_{K} w = \frac{1}{\| K \|} \sum σ \in E_{K} \| σ \| (w_{σ} - w_{K}) n_{K, σ} .$	(21a)
We see in (20) that the discrete gradient consists of two components. The first component, ${¯ ¯¯¯ ¯ \nabla}_{K} w$ in (21a), is a linearly exact reconstruction of the gradient. The second component, $S_{K} w (x)$ , is a stabilisation term, which is piecewise defined such that for each convex hull $(T_{K, σ})_{σ \in E_{K}}$ of $σ$ and $x_{K}$ (see Figure 3, right), $S_{K} w (x) = S_{K, σ} w$ , with
	$S_{K, σ} w = \frac{\sqrt{2}}{d_{K, σ}} [w_{σ} - w_{K} - {¯ ¯¯¯ ¯ \nabla}_{K} w \cdot ({¯ ¯ ¯ x}_{σ} - x_{K})] n_{K, σ} .$	(21b)

Figure 3. Notation for a Cartesian cell in dimension $d = 2$ .

Since we work on square cells, the expressions (21a) and (21b) can be simplified. Firstly, we have for every edge $σ \in E_{K}$ , $| σ | = h$ . By adopting the compass notation and denoting by $σ_{N}$ , $σ_{S}$ , $σ_{E}$ , $σ_{W}$ the north, south, east and west edges of cell $K$ , respectively, we find that

{¯ ¯¯¯ ¯ \nabla}_{K} w = \frac{1}{h} [\begin{matrix} w_{σ_{E}} - w_{σ_{W}} w_{σ_{N}} - w_{σ_{S}} \end{matrix}] .

Moreover, we obtain the following simplified expression for the stabilisation term.

S_{K, σ} = ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} \frac{\sqrt{2}}{h} (w_{σ_{E}} - 2 w_{K} + w_{σ_{W}}) n_{K, σ} f o r σ = σ_{E} o r σ = σ_{W}, \frac{\sqrt{2}}{h} (w_{σ_{N}} - 2 w_{K} + w_{σ_{S}}) n_{K, σ} f o r σ = σ_{N} o r σ = σ_{S} . \end{matrix}

Finally, owing to (19), the diffusive flux $F_{K, σ}^{D}$ along the edge $σ$ of cell $K$ can be uniquely determined by choosing $v \in X_{D}$ such that $v_{K} = 0, v_{σ} = - 1$ and $v_{σ^{'}} = 0$ for $σ^{'} \neq σ$ .

Example 4.1 (Computing the diffusive flux along an edge).

As an example, to compute the diffusive flux along the eastern edge $σ_{E}$ of a cell $K$ , we choose $v \in X_{D}$ such that $v_{K} = 0$ , $v_{σ_{E}} = - 1$ , $v_{σ_{N}} = v_{σ_{S}} = v_{σ_{W}} = 0$ . Substituting into (19) and denoting by $Λ_{K}$ the value of the diffusion tensor at cell $K$ , we then have

F_{K, σ_{E}}^{D} = \int_{K} Λ_{K} \nabla_{D} c \cdot \nabla_{D} v d A = \sum σ \in E_{K} \int_{T_{K, σ}} Λ_{K} \nabla_{D} c \cdot \nabla_{D} v d A,

with

Λ_{K} = [\begin{matrix} λ_{K, 11} & λ_{K, 12} λ_{K, 21} & λ_{K, 22} \end{matrix}] .

We then compute

\nabla_{D} v = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} - \frac{1}{h} e_{x} & f o r σ = σ_{N} o r σ = σ_{S}, (- \frac{1}{h} - \frac{\sqrt{2}}{h}) e_{x} & f o r σ = σ_{E}, (- \frac{1}{h} + \frac{\sqrt{2}}{h}) e_{x} & f o r σ = σ_{W} . \end{matrix}

Using definition (20) of the discrete gradient and denoting the values of $c$ at the cell $K$ and its edges by $c_{K}, c_{σ_{N}}, c_{σ_{S}}, c_{σ_{E}}, c_{σ_{W}}$ , respectively, we then compute the expression corresponding to $\nabla_{D} c$ . Using the appropriate values of $\nabla_{D} c$ and $\nabla_{D} v$ for each of the regions $T_{K, σ}$ , we obtain

F_{K, σ_{E}}^{D} = - 1 (2 λ_{K, 11} (c_{σ_{E}} - c_{K}) + λ_{K, 12} (c_{σ_{N}} - c_{σ_{S}})) .

Here, we see that $F_{K, σ_{E}}^{D}$ is an approximation of

- \int_{σ_{E}} Λ \nabla c \cdot n_{K, σ_{E}} d s = - \int_{σ_{E}} (λ_{K, 11} c_{x} + λ_{K, 12} c_{y}) d s

in the sense that $c_{x}$ , $c_{y}$ are approximated by $\frac{2}{h} (c_{σ_{E}} - c_{K})$ and $\frac{1}{h} (c_{σ_{N}} - c_{σ_{S}})$ , respectively. ∎

4.4. Advective fluxes

We now discuss the discretisation of the advective fluxes $F_{K, σ}^{A} \approx \int_{σ} c V \cdot n_{K, σ} d s$ . The definition of the discrete advective fluxes is motivated by the Scharfetter-Gummel or exponential scheme [25, 26]. Here, we briefly recall the Scharfetter-Gummel flux, and for simplicity of exposition, we consider an isotropic diffusion tensor $Λ = ε I$ and velocity field $V = [v_{1}, 0]^{T}$ . Defining the Péclet number

P = \frac{h v_{1}}{ε},

(22)

the Scharfetter-Gummel flux across the eastern edge $σ_{E}$ of cell $K$ shared with cell $E$ is given by

F_{K, σ} = - \frac{ε | σ_{E} |}{h} (B (P) c_{E} - B (- P) c_{K}),

(23)

where $B$ is the Bernoulli function with

B (z) = \frac{z}{e^{z} - 1} .

We now note that the definition of the Scharfetter-Gummel flux (23) is based on the Péclet number $P$ . Hence, for a good definition of advective fluxes $F_{K, σ}^{A}$ to be used in (16b) and (16c), we need a proper definition of a grid-based Péclet number for anisotropic advection-diffusion problems. In particular, for each $K \in Ω_{h}, σ \in E_{K}$ , we define as in [5, Section 3.2.2] a local grid-oriented Péclet number

P_{K, σ} = \frac{h V_{K, σ}}{2 λ_{K, σ}},

(24)

with

	$V_{K, σ}$	$= \frac{1}{\| σ \|} \int_{σ} V \cdot n_{K, σ} d s,$
	$λ_{K, σ}$	$= n_{K, σ}^{T} Λ_{K} n_{K, σ} .$

Here, the quantities $V_{K, σ}$ and $λ_{K, σ}$ measure the strength of advection and diffusion across the outward normal $n_{K, σ}$ , respectively; hence allowing the Péclet number $P_{K, σ}$ to capture the relative strength of advection over diffusion along $n_{K, σ}$ properly. As a remark, we note that for isotropic diffusion, the expression (24) of the Péclet number boils down to the classical definition (22) of the Péclet number.

Since the diffusive fluxes have already been defined in Section 4.3, we consider only the advective component of the Scharfetter-Gummel flux. This is done by following the ideas in [1, 2], where instead of the Bernoulli function, we extract the advective component of the Scharfetter-Gummel fluxes by defining

A_{s g} (z) = B (- z) - 1.

The advective fluxes are then given by

F_{K, σ}^{A}

(25)

5. Summary and implementation of numerical scheme

In this section, we present a summary of the numerical scheme for the magnetised transport problem described by the model system of equations (1a), (1b), and (1c). Starting with $m = 0$ and initial ion and electron density profiles of $n_{i}^{(0)}$ and $n_{e}^{(0)}$ , respectively, the idea is to sequentially solve, the potential $V^{(m + 1)}$ , followed by the ion density $n_{i}^{(m + 1)}$ and then the electron density $n_{e}^{(m + 1)}$ . This process is repeated until we reach a time $T = t^{(m + 1)}$ such that the solution is already at steady state. Numerically, the solution is said to be at steady state whenever

∥ V

Combining the hybrid mimetic mixed method with the Scharfetter-Gummel scheme for magnetised transport in plasmas

Authors

A Generalised Complete Flux scheme for anisotropic advection-diffusion equations

Flux correction for nonconservative convection-diffusion equation

A Hybrid High-Order scheme for the stationary, incompressible magnetohydrodynamics equations

Mixed methods for large-strain poroelasticity/chemotaxis models simulating the formation of myocardial oedema

Lattice Boltzmann method for simulation of diffusion magnetic resonance imaging physics in heterogeneous tissue models

A hybrid WENO method with modified ghost fluid method for compressible two-medium flow problems

Visualizing the world's largest turbulence simulation

1. Introduction

2. Mathematical model

2.1. Boundary conditions

3. Time discretisation

3.1. First-order approximation of the charge density

3.2. Second-order approximation of the charge density

4. Spatial discretisation

4.1. Finite volume discretisation

4.2. Boundary conditions and conservation of fluxes

4.3. Diffusive fluxes

Example 4.1 (Computing the diffusive flux along an edge).

4.4. Advective fluxes

5. Summary and implementation of numerical scheme

Combining the hybrid mimetic mixed method with the Scharfetter-Gummel scheme for magnetised transport in plasmas

Authors

Related Research

A Generalised Complete Flux scheme for anisotropic advection-diffusion equations

Flux correction for nonconservative convection-diffusion equation

A Hybrid High-Order scheme for the stationary, incompressible magnetohydrodynamics equations

Mixed methods for large-strain poroelasticity/chemotaxis models simulating the formation of myocardial oedema

Lattice Boltzmann method for simulation of diffusion magnetic resonance imaging physics in heterogeneous tissue models

A hybrid WENO method with modified ghost fluid method for compressible two-medium flow problems

Visualizing the world's largest turbulence simulation

This week in AI

1. Introduction

2. Mathematical model

2.1. Boundary conditions

3. Time discretisation

3.1. First-order approximation of the charge density

3.2. Second-order approximation of the charge density

4. Spatial discretisation

4.1. Finite volume discretisation

4.2. Boundary conditions and conservation of fluxes

4.3. Diffusive fluxes

Example 4.1 (Computing the diffusive flux along an edge).

4.4. Advective fluxes

5. Summary and implementation of numerical scheme