Sheath governing equations in computational weakly-ionized plasmadynamics

Abstract

To date, fluid models of plasma sheaths have consisted of the coupling of the electric field potential equation obtained through Gauss’s law to the charged species conservation equations obtained through the drift–diffusion approximation. When discretized using finite-difference stencils, such a set of equations has been observed to be particularly stiff and to often require more than hundreds of thousands of iterations to reach convergence. A new approach at solving sheaths using a fluid model is here presented that reduces significantly the number of iterations to reach convergence while not sacrificing on the accuracy of the converged solution. The method proposed herein consists of rewriting the sheath governing equations such that the electric field is obtained from Ohm’s law rather than from Gauss’s law. To ensure that Gauss’s law is satisfied, some source terms are added to the ion conservation equation. Several time-accurate and steady-state cases of dielectric sheaths, anode sheaths, and cathode sheaths (including glow and dark discharges) are considered. The proposed method is seen to yield the same converged solution as the conventional approach while exhibiting a reduction in computational effort varying between one-hundred-fold and ten-thousand-fold whenever the plasma includes both quasi-neutral regions and non-neutral sheaths.

Introduction

In plasmadynamics, a “sheath” refers to a plasma region that is located adjacent to surface boundaries and that is significantly non-neutral. That is, a sheath always exhibits a substantial difference between its positive and negative charge densities. Several different types of sheaths are recognized, depending on the type of boundary material and also depending on the direction of the electric field near the surface. When the boundary is a dielectric, the electric field generally points towards the surface. In such a case, the plasma near the dielectric has an excess of positive charge, is a few Debye lengths thick, and is denoted as a “dielectric sheath”. When the boundary is an electrode and the electric field points away from the surface, the plasma close to the electrode has an excess of negative charge and is denoted as an “anode sheath”. When the boundary is an electrode and the electric field points towards the surface, the plasma close to the electrode has an excess of positive charge and is denoted as a “cathode sheath”.

An accurate solution of the cathode sheath is generally crucial to predict correctly the current and electric field distribution through the rest of the plasma. This follows from the plasma near the cathode being characterized by an electron number density that is so low that the current in that region is mostly ionic rather than electronic. Because the ion mobility is substantially less than the electron mobility, it is not uncommon for the conductivity in the vicinity of the cathode to be less than one thousandth of the conductivity within the rest of the plasma. Such a low conductivity near the cathode results in a large voltage drop for a desired current, and this in turn leads to a large amount of power dissipated within the sheath in form of heat. Because of this, accurately modeling sheaths is generally deemed essential to assess adequately the performance of plasma-based devices such as MHD generators and accelerators [1], [2], or plasma actuators [3], [4], [5].

Because the sheath thickness is typically very small compared to the size of the plasma, its accurate resolution may necessitate a substantial refinement of the grid near the boundaries. To overcome this problem, the sheath and the plasma bulk can be solved through independent methods and “patched” to each other through some adequate boundary conditions [6], [7]. While the patching approach is advantaged by not requiring the sheath and the rest of the plasma to be solved in coupled form (hence providing substantial savings in computational effort) it is questionable how well it can model more complex situations where a strong coupling may occur, such as those involving motion of the neutrals and magnetic field effects.

To overcome the limitations of the patching method, a coupled solution of the plasma bulk and sheath is necessary. This can be accomplished through a kinetic simulation [8], a particle-in-cell simulation [9], or through a fluid model [10], [11], [12]. Also known as the “drift–diffusion” model, the fluid model is currently the method of choice in simulating sheaths that would occur within plasma aerodynamic applications (see for instance Refs. [13], [14], [15]). Indeed, when solving the relatively high density flows typical of plasma aerodynamics problems, the more detailed physical models associated with kinetic simulations or particle-in-cell simulations entail computing efforts that are beyond the capabilities of present computers.

To date, fluid models of plasma sheaths have consisted of the coupling of the electric field potential obtained through Gauss’s law to the electron and ion conservation equations obtained through the drift–diffusion approximation. While such a strategy has had considerable success, it suffers from very slow convergence when solving plasmas that include quasi-neutral regions. When a quasi-neutral region forms within the plasma, the stiffness of the governing equations is such that typically one hundred thousand iterations or more are needed to reach steady-state. Such slow convergence has been observed not only for explicit relaxation schemes but also for fully-coupled implicit solvers. For this reason, quasi-neutral plasmas that do not include sheaths are typically solved by obtaining the electric field from Ohm’s law rather than from Gauss’s law. The stiffness of the governing equations is then relieved substantially and the computing effort can be reduced one-hundred-fold or even more. Such a strategy cannot be used for a plasma that includes non-neutral sheaths, however, because an electric field that is obtained from Ohm’s law may not necessarily satisfy Gauss’s law, and because the correct solution of Gauss’s law is crucial to model the non-neutral regions of the plasma correctly.

It is here argued that the last statement is not necessarily correct. That is, it is not mandatory to obtain the electric field from Gauss’s law to ensure that the latter is satisfied. As will be demonstrated in this paper, it is possible to use Ohm’s law to obtain the electric field within the non-neutral sheath regions while satisfying Gauss’s law. In so-doing, remarkable gains in computational efficiency are realized when solving plasmas that include both sheaths and quasi-neutral regions.

This paper does not constitute the first attempt at solving sheaths using Ohm’s law. For instance, in Refs. [16], [17], [18], a thermionic cathode sheath is solved by determining the electric field from a form of Ohm’s law. The method outlined in the latter references, however, does not ensure that Gauss’s law is satisfied and is hence only applicable to thermionic sheaths, and not to cold cathode sheaths, dielectric sheaths, or anode sheaths. In contrast, we here propose a set of governing equations that is such that the electric field is determined from Ohm’s law while ensuring that Gauss’s law is satisfied. The novel set of governing equations proposed herein results in the same converged solution as the conventional set (in which the electric field is determined directly through Gauss’s law) while yielding as much as a ten-thousand-fold reduction in computing effort. This is attributed to the proposed governing equations being considerably less stiff and hence permitting the use of significantly larger time steps when being integrated.

It is noted that a “stiff” system of equations here does not necessarily denote a system in which there is a large discrepancy between the eigenvalues or the physical time scales. Rather, we here adopt Lambert’s definition of stiffness in the more general sense [19, p. 220] with a slight modification to make it applicable to non-linear systems of equations:

Definition If a numerical method, applied to a system with any initial conditions, is forced to use in a certain interval of integration a steplength which is excessively small in relation to the smoothness of the exact solution in that interval, then the system is said to be stiff in that interval.

When defined in this manner, stiffness is not an intrinsic property of the physical model: stiffness can originate from the disparity between the physical time scales, but can also originate from other attributes of the system of equations (see Chapter 6 in Ref. [19]). For instance, when a physical model involves several conservation laws dependent on each other, there is more than one way that the governing equations can be expressed, and the equations may exhibit more or less stiffness depending on how they are written, despite solving the same physical model and sharing the same physical time scales. Rewriting the governing equations differently while keeping the physical model the same was the strategy used to reduce the stiffness of the system associated with a quasi-neutral multicomponent plasma model in Ref. [20], with a currentless plasma model in Ref. [21], and with a two-fluid Euler–Poisson plasma model in Ref. [22]. The methods outlined in the latter references, however, cannot be readily extended to the drift–diffusion model of sheaths. In this paper, we show that the stiffness of the system of equations associated with the drift–diffusion model of sheaths can also be alleviated substantially through a recast in different form, despite the recast set of equations solving the same physical model as the original set.