Sheath governing equations in computational weakly-ionized plasmadynamics
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:
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.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.
Section snippets
Physical model
The physical model considered herein determines the electron, ion, and neutral properties through macroscopic-scale transport equations (i.e. the so-called drift–diffusion model). Such can accurately predict using the same set of differential equations both the non-neutral plasma sheaths and the quasi-neutral plasma bulk, including ambipolar diffusion and ambipolar drift phenomena. It is emphasized that the physical model presented in this section is commonly used to model sheaths taking place
Conventional sheath governing equations
The term “conventional governing equations” here denotes the set of partial differential equations that are generally used to solve the physical model outlined in the previous section. In the conventional set of governing equations, the electric field is determined from a form of the potential equation obtained from Gauss’s law. This can be found starting from Eq. (6) and assuming that a potential exists such that :The conventional governing equations hence consists
Proposed sheath governing equations
To reduce the computational effort needed to solve plasmas that include both non-neutral and quasi-neutral regions, a new method is presented here. The method proposed consists of obtaining the electric field potential from a form of Ohm’s law instead of Gauss’s law. Then, to guarantee that Gauss’s law is satisfied both in the neutral and quasi-neutral regions, some source terms are added to the ion conservation equation.
To yield a valid solution to the physical model presented in Section 2,
Discretization of the governing equations
Noting that the residual for the conventional set of sheath governing equations (see Eq. (8)) has the same form as the residual for the proposed sets of sheath governing equations (see Eq. (18), (26)), we can express both sets of governing equations in discrete form as follows:where is the discretized residual which we seek to minimize, and where and are some discretization operators which vary depending on the term being discretized. The discrete
Pseudotime relaxation of the discretized equations
To attain a converged solution, the discrete residual outlined in Eq. (28) must be reduced to a small quantity at every node. This is here accomplished through the use of pseudotime relaxation combined with a block-implicit method. A block-implicit pseudotime relaxation strategy is chosen because it is the preferred relaxation technique in various compressible fluid dynamics and plasma aerodynamics codes (see NASA’s OVERFLOW and CFL3D codes for instance or the recent papers in plasma
Boundary conditions
At the interface between the plasma and a solid surface, the boundary conditions for the electron and ion number densities depend on the direction of the electric field. For this purpose, it is convenient to consider a coordinate which is oriented perpendicular to the surface and which points away from the surface. Then, when the electric field vector points away from the surface, the boundary conditions take the form:On the other hand, when the electric field vector
Test cases
Several test cases are now considered to assess the performance of the proposed governing equations over the conventional governing equations in solving sheaths. In all cases, the medium is a three-component air plasma including electrons, neutrals, and one type of positive ions. The electron and ion mobilities are as specified in Table 1 while the chemical reactions taking place within the plasma are listed in Table 2. In the chemical model used herein, the electron creation mechanisms are
Conclusions
To date, the numerical simulation of plasma sheaths using macroscopic-scale transport equations has involved the coupling of the electric field potential equation obtained through Gauss’s law to the electron and ion transport equations obtained through the drift–diffusion model. When discretized using finite difference stencils, this set of equations (referred herein as the “conventional governing equations”) is particularly stiff and typically requires hundreds of thousands of iterations to
Acknowledgment
This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant #2010-0023957).
References (35)
- et al.
Three-dimensional simulation of the electric field and magnetohydrodynamic power generation during reentry
AIAA Journal
(2009) - et al.
Numerical study of an electron-beam-confined Faraday accelerator
Journal of Propulsion and Power
(2007) - et al.
Periodic electrodynamic field of dielectric barrier discharge
AIAA Journal
(2011) - et al.
Magnetoaerodynamic actuator for hypersonic flow control
AIAA Journal
(2005) - et al.
Modeling of dielectric barrier discharge plasma actuators driven by repetitive nanosecond pulses
Physics of Plasmas
(2007) - et al.
A suitable boundary condition for bounded plasma simulation without sheath resolution
Journal of Computational Physics
(1993) - et al.
Sheath criterion for a collisional sheath
Physics of Plasmas
(2003) - et al.
Numerical study of a direct current plasma sheath based on kinetic theory
Physics of Plasmas
(2002) - et al.
A PIC based procedure for the integration of multiple time scale problems in gas discharge physics
Journal of Computational Physics
(2009) - et al.
A continuum model of DC and RF discharges
IEEE Transactions on Plasma Science
(1986)
A one-dimensional model of DC glow discharges
Journal of Applied Physics
Modeling low pressure collisional plasma sheath with space-charge effect
Physics of Plasmas
Modeling of air plasma generation by repetitive high-voltage nanosecond pulses
IEEE Transactions on Plasma Science
Two-component plasma model for two-dimensional glow discharge in magnetic field
Journal of Computational Physics
Numerical simulation of direct current glow discharges for high-speed flow control
Journal of Propulsion and Power
A unified theory of free burning arcs, cathode sheaths and cathodes
Journal of Physics D: Applied Physics
A one-dimensional theory for the electrode sheaths of electric arcs
Journal of Physics D: Applied Physics
Cited by (21)
An implicit time integration approach for simulation of corona discharges
2024, Computer Physics CommunicationsDouble distribution function lattice Boltzmann modeling for energy transport in the DC argon arc plasma
2016, International Communications in Heat and Mass TransferCitation Excerpt :Therefore, numerical simulations of thermal behavior and fluid dynamics in plasma are necessary. Several numerical methods have previously been used to model plasma, such as the finite element method (FEM) [1–3], finite volume method (FVM) [4,5], finite difference method (FDM) [6,7] and Monte-Carlo method [8,9]. The Lattice Boltzmann method (LBM) is a robust numerical technique based on kinetic theory, and because of its kinetic nature and relatively easy calculation procedure, the LBM has a very wide application from single-phase flows to magneto-hydrodynamic flows [10–12].
AETHER: A simulation platform for inductively coupled plasma
2015, Journal of Computational PhysicsModeling weakly-ionized plasmas in magnetic field: A new computationally-efficient approach
2015, Journal of Computational PhysicsCitation Excerpt :Such stiffness has been observed to be independent of the type of integration strategy used, either explicit or fully-implicit. As first outlined in Ref. [15], the stiffness originates from the potential equation based on Gauss's law being particularly sensitive to small errors in the electron or ion densities whenever the plasma becomes quasi-neutral. One approach that has been shown successful in relieving the stiffness is by rewriting the governing equations such that the electric field is obtained from a potential based on Ohm's law rather than Gauss's law (see Refs. [15] and [17]).
Electron and ion transport equations in computational weakly-ionized plasmadynamics
2014, Journal of Computational PhysicsCitation Excerpt :This is problematic because the redefined boundary condition is approximate and it is uncertain whether it would remain valid in the general case. Further, the approach presented in Ref. [6] is limited to a three-component plasma (one type of positive ions, electrons, and neutrals) in one dimension, and it is not clear how it could be extended to multidimensional and multicomponent plasmas (plasmas with several types of positive ions and negative ions). The goal of this paper is to build upon the ideas presented in Ref. [6] and to derive a new set of electron and ion transport equations that is more computationally efficient than the conventional set and that is generally applicable to multidimensional and multicomponent weakly-ionized plasma flows.
Numerical Analysis for Argon Arc Plasma Jet Flow by Three-Dimensional Thermal Lattice Boltzmann Model
2023, Bulletin of the Lebedev Physics Institute
- 1
Current address: Lockheed Martin Aeronautics Company, 1011 Lockheed Way, Palmdale, CA 93599-0160, United States.