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

by   Hanz Martin Cheng, et al.
TU Eindhoven

In this paper, we propose a numerical scheme for fluid models of magnetised plasmas. One important feature of the numerical scheme is that it should be able to handle the anisotropy induced by the magnetic field. In order to do so, we propose the use of the hybrid mimetic mixed (HMM) scheme for diffusion. This is combined with a hybridised variant of the Scharfetter-Gummel (SG) scheme for advection. The proposed hybrid scheme can be implemented very efficiently via static condensation. Numerical tests are then performed to show the applicability of the combined HMM-SG scheme, even for highly anisotropic magnetic fields.



page 17

page 19

page 20


A Generalised Complete Flux scheme for anisotropic advection-diffusion equations

In this paper, we consider separating the discretisation of the diffusiv...

Flux correction for nonconservative convection-diffusion equation

Our goal is to develop a flux limiter of the Flux-Corrected Transport me...

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

We propose and analyse a hybrid high-order (HHO) scheme for stationary i...

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

In this paper we propose a novel coupled poroelasticity-diffusion model ...

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

We report the first implementation of the Lattice Boltzmann method (LBM)...

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

In this paper, we develop a simplified hybrid weighted essentially non-o...

Visualizing the world's largest turbulence simulation

In this exploratory submission we present the visualization of the large...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 is applied between two parallel electrodes and that an electron flux is forced through a uniform medium, where an oblique magnetic field is applied. The electric field in a spatial domain over a time interval is governed by the Poisson equation

where is the permittivity of vacuum, is the potential, and is the space charge density given by
Here, is the elementary charge, and 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:

where and are the drift-diffusion flux densities for the ions and electrons, respectively, defined by

Equations (1b) and (1c) are coupled to the Poisson equation (1a) due to the electric field 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 , referring to a particle (either an ion or an electron).

ionisation rate
mobility coefficient
diffusion coefficient
magnetic tensor
local electron density
Table 1. Parameters in the model.
density of the background gas
eV mean electron energy
permittivity of vacuum
Boltzmann constant
elementary charge
temperature of the background gas
Table 2. Constants in the model.

We now discuss the parameters used in the model. Firstly, the ionisation rate and the mobility coefficients are computed via local field approximations, and are hence functions of the reduced electric field . In this paper, we consider helium gas, whose ion mobility (Figure 2, left) is obtained from [11]. The ionisation rate (Figure 1) and electron mobility (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, refers to the -norm.

Figure 1. Ionisation rate
Figure 2. Mobility coefficients (Left: , Right: ).

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

where, and . For ions, we assume that and so we have

We assume that electrons follow a Maxwellian distribution so that and thus

Thirdly, the the local electron density is computed by the relation [15]


Finally, the magnetic tensor is given by

with Hall parameters , and magnetic field . We now note that

is a cross product matrix such that for any vector

, . Taking note that
we then see that

is an eigenvector of

corresponding to the eigenvalue 1. Essentially, this means that

preserves the component of the vector that is parallel to the magnetic field.

For dimension , the magnetic field can be specified at an angle with respect to the -axis by setting


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

orthogonal to the magnetic field, we see that

For two-dimensional models, we neglect the -component, and hence we see that the components orthogonal to the magnetic field are scaled by a factor . 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),

Then, flux boundary conditions are imposed for the drift-diffusion equations (1b) and (1c), i.e.,
where is the thermal velocity, denotes the positive part of the function, and 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 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:


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 at time , we are then interested in calculating the densities at time . 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 by taking . 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 , for , and given the densities , we seek such that



is the drift-diffusion flux density of particle (ions or electrons) at time . 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 .

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 , such that


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 . 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 , we may write


We then note that the time derivative can be computed from the drift-diffusion equations, and can be written as


Using now (8) for the charge density in the Poisson equation and evaluating in (9) at time , we then seek, for , such that


In this case, the quantities and at time 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 implicit and take first-order approximations to , which leads to

or equivalently,

Denoting the correction flux at time level by


the above equation can be rewritten as


4. Spatial discretisation

For the spatial discretisation, we start by forming , a partition of the domain into control volumes. Each of the control volumes has edges such that , and we denote by the total number of edges contained in this partition. We then denote by 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 . This results to, in each control volume ,

respectively. By applying the divergence theorem and using , these can be rewritten as


Here, is the outward unit normal vector of at . Due to the definition (3) of , the term is nonlinear in . Here, we choose a first-order approximation and take the local electron density at time so that


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 for some diffusion tensor and scalar ; on the other hand, the advective fluxes take the form , where 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 by


Replacing the fluxes in (13a)-(13c) with their discrete counterparts (15a)-(15e) at the appropriate time level and under the assumption that and are piecewise constant with values and respectively in a cell , we then obtain our numerical scheme. That is, for , given the densities , we seek such that for each ,

where is the average ionisation rate in cell at time , and is the value of computed at cell .

Starting with , 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



is a discretisation of the semi-implicit correction flux 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 . Hence, we choose to use (16a) for solving the Poisson equation at the first time step when .

4.2. Boundary conditions and conservation of fluxes

We now see that there are equations (corresponding to the number of cells) involved in each of (16a)-(16d). However, we have unknowns for each of these equations: corresponding to the number of cells, and corresponding to the number of edges. This is standard for hybrid finite volume methods, and the equations for the 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 and , we should have that the total (diffusive and advective) flux going from to via is equal to the negative total flux going from to via , i.e.


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 ,


Here, and are the values (to be determined) of and 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 . 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 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:


where is a discrete gradient that is defined as follows: For , we set

We see in (20) that the discrete gradient consists of two components. The first component, in (21a), is a linearly exact reconstruction of the gradient. The second component, , is a stabilisation term, which is piecewise defined such that for each convex hull of and (see Figure 3, right), , with
Figure 3. Notation for a Cartesian cell in dimension .

Since we work on square cells, the expressions (21a) and (21b) can be simplified. Firstly, we have for every edge , . By adopting the compass notation and denoting by , , , the north, south, east and west edges of cell , respectively, we find that

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

Finally, owing to (19), the diffusive flux along the edge of cell can be uniquely determined by choosing such that and for .

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

As an example, to compute the diffusive flux along the eastern edge of a cell , we choose such that , , . Substituting into (19) and denoting by the value of the diffusion tensor at cell , we then have


We then compute

Using definition (20) of the discrete gradient and denoting the values of at the cell and its edges by , respectively, we then compute the expression corresponding to . Using the appropriate values of and for each of the regions , we obtain

Here, we see that is an approximation of

in the sense that , are approximated by and , respectively. ∎

4.4. Advective fluxes

We now discuss the discretisation of the advective fluxes . 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 and velocity field . Defining the Péclet number


the Scharfetter-Gummel flux across the eastern edge of cell shared with cell is given by


where is the Bernoulli function with

We now note that the definition of the Scharfetter-Gummel flux (23) is based on the Péclet number . Hence, for a good definition of advective fluxes 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 , we define as in [5, Section 3.2.2] a local grid-oriented Péclet number



Here, the quantities and measure the strength of advection and diffusion across the outward normal , respectively; hence allowing the Péclet number to capture the relative strength of advection over diffusion along 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

The advective fluxes are then given by


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 and initial ion and electron density profiles of and , respectively, the idea is to sequentially solve, the potential , followed by the ion density and then the electron density . This process is repeated until we reach a time such that the solution is already at steady state. Numerically, the solution is said to be at steady state whenever