Numerical Modelling of Electrical Discharges

Numerical Modelling of Electrical Discharges F. J. Dur´an-Olivencia1,∗, F. Pontiga1,∗∗ and A. Castellanos2,∗∗∗ 1Departamento de F´ısica Aplicada II, University of Seville, Spain 2Departamento de Electr´onica y Electromagnetismo, University of Seville, Spain E-mail: ∗[email protected], ∗∗ pontig[email protected], ∗∗∗castell[email protected] Abstract. The problem of the propagation of an electrical discharge between a spherical electrode and a plane has been solved by means of finite element methods (FEM) using a fluid approximation and assuming weak ionization and local equilibrium with the electric field. The numerical simulation of this type of problems presents the usual difficulties of convection-diffusion-reaction problems, in addition to those associated with the nonlinearities of the charged species velocities, the formation of steep gradients of the electric field and particle densities, and the coexistence of very different temporal scales. The effect of using different temporal discretizations for the numerical integration of the corresponding system of partial differential equations will be here investigated. In particular, the so-called θ-methods will be used, which allows to implement implicit, semi-explicit and fully explicit schemes in a simple way. 1. Introduction It is a well known fact that implicit Galerkin FEM may introduce spurious oscillations in convection dominated problems, which may lead to severe errors in the results. These oscillations appear especially when the solutions exhibit strong gradients, discontinuities and shock fronts. The source of this instability is directly related to the symmetry of the weighting functions used in the Galerkin method. In an attempt to eliminate such problems, formulations such as streamline upwind Petrov-Galerkin methods add a consistent artificial diffusion by using modified weighting functions, in which the element in the upwind node is more heavily weighted than that in the downwind node [1]. In contrast, the approach adopted by the θ-method is to evaluate the temporal derivative using a combination of implicit and explicit procedures [2]. Therefore, the computed numerical solution is not smoothed by the addition of any artificial diffusion. In principle, different temporal discretizations could be applied within the θ-method to the various terms (convection, diffusion and reaction) that appear in the governing equations. However, finding the optimal time-weighting for each term is a challenging task, due to the relationship of these terms with the electric field, which undergoes strong variations during the development of discharge. Thus, a global time discretization is more appropriated for the simulation of electrical discharges. Here, the optimal combination between explicit and implicit procedures in the θ-method will be investigated. 2. Problem formulation 2.1. Electrode configuration The electrode configuration consists of a spherical cathode, with radius R= 5 mm, situated at 20 mm from a grounded plane. The gas filling the gap between electrodes is pure oxygen, at a 2nd International Conference on Mathematical Modeling in Physical Sciences 2013 IOP Publishing Journal of Physics: Conference Series 490 (2014) 012209 doi:10.1088/1742-6596/490/1/012209 Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Published under licence by IOP Publishing Ltd 1 pressure of 50 Torr. The sphere is subjected to negative high voltage (−2900 V) and a certain initial density of electrons (400 seed electrons) is introduced in the vicinity of the cathode to trigger the discharge. The spatial distribution of these seed electrons is Gaussian, with is maximum at 55 µm from the cathode surface. 2.2. Governing equations One of the most significant numerical studies about corona discharge in oxygen was carried out by Morrow [3], using a flux-corrected transport (FCT) algorithm in finite differences. In his work, he studied the spatiotemporal development of the first discharge pulse (known as Trichel pulse) and he succeeded in computing the drastic distortion of the electric field under the effect of the spatial charge. Here, Morrow’s model of corona discharge will be solved by means of finite element methods implemented through the software Comsol Multiphysics. Therefore, the reader is referred to the work of Morrow for a detailed description of the physical model. Assuming a weak ionization and a local equilibrium of the transport coefficients with the electric field, the governing equations can be expressed as a set of continuity equations, one for each species, coupled with Gauss equation, ∂Ne ∂t +∂ ∂x −De ∂Ne ∂x +NeWe=α|We|Ne−η|We|Ne−βNeNp,(1) ∂Np ∂t +∂ ∂xNpWp=α|We|Ne−βNp(Ne+Nn),(2) ∂Nn ∂t +∂ ∂xNnWn=η|We|Ne−βNpNn,(3) ∇ · ε0E=ρ(x),(4) where Nis the number density; Wis the drift velocity; Eis the electric field; α,η,βand D are the ionization, attachment, recombination and diffusion coefficients; and subscripts e,pand nrefer to electron, positive ions and negative ions, respectively. Following Morrow’s model [3], the continuity equations (1)–(3) are solved along the axis of symmetry, x. However, the electric field must be evaluated in two dimensions. Boundary conditions for the species densities and the electric field can be written as follows, Cathode (x= 0): Ne(t) = γiNp|Wp/We|+Np e(t), Nn(t) = 0, φ =−2900 V.(5) Anode (x=d): Ne(t) = 0, Np(t) = 0, φ = 0,(6) where γiis the ion secondary ionization coefficient, Np e(t) is the number of secondary electrons released at cathode by photons created in discharge volume, and φis the electrical potential (E=−∇φ). 2.3. Temporal differentiation The temporal derivative of a variable Tis approximated in the θ-method as Tn+1 −Tn ∆t=θ˙ Tn+1 + (1 −θ)˙ Tn,(7) where 0 ≤θ≤1. The essence of this method is that by varying θit is possible to conduct different time integration schemes. The limits θ= 0 and θ= 1 correspond to explicit Euler (backward difference) and implicit Euler (forward difference) integration schemes, respectively. 2nd International Conference on Mathematical Modeling in Physical Sciences 2013 IOP Publishing Journal of Physics: Conference Series 490 (2014) 012209 doi:10.1088/1742-6596/490/1/012209 2 Figure 1. Electric field distribution at different times obtained by FEM with Crank- Nicolson time weighting (θ=1/2). Figure 2. Electron density distribution at different times obtained by FEM with Crank- Nicolson time weighting (θ=1/2). 3. Results and discussion Figure 1, 2 and 3 show, respectively, the spatial distribution of electric field, electron density and negative ion density computed at different times by means of FEM using a Crank-Nicolson time weighting (θ= 1/2). Clearly, the space charge generated by electron avalanches modifies the initial Laplacian electric field (fig. 1, t= 0 ns), and generates a sharp gradient in the vicinity of the cathode. The electric field is reinforced on the cathode surface, and a quasi-neutral region (or plasma region) starts to develop in front of this electrode (see fig. 1). The dramatic evolution of the electric field originates further numerical difficulties. Firstly, Figure 3. Negative density distribution at different times obtained by FEM with Crank- Nicolson time weighting (θ=1/2). 2nd International Conference on Mathematical Modeling in Physical Sciences 2013 IOP Publishing Journal of Physics: Conference Series 490 (2014) 012209 doi:10.1088/1742-6596/490/1/012209 3 -6 -4 -2 0 2 0.3 0.5 0.7 electric field, E (104 V/m) x (mm) 51ns, θ = 0.5 51ns, θ = 1 60ns, θ = 0.5 60ns, θ = 1 Figure 4. Spatial distribution of the electric field at times t= 51 ns and t= 60 ns obtained by using a fully implicit scheme (θ= 0) and a semi-explicit Crank-Nicolson scheme (θ=1/2). the intense drop of the electric field intensity decelerates the electrons that migrate towards the anode. Therefore, electrons start accumulating in the quasi-neutral region and a sharp peak of electron density is formed (fig. 2). Then, negative ions are generated by attachment of electrons to oxygen molecules, which produces in turn a maximum of its number density (fig. 3). All these processes confront the numerical schemes with very stiff problems, and the solutions they give depends crucially on the numerical implementation of the temporal differentiation. In figure 4, the results obtained by using implicit differentiation (θ= 0) and a hybrid implicitexplicit differentiation (θ= 1/2) are compared in these critical situations. Clearly, the abrupt change of the electric field in front of the cathode is dealt much more satisfactorily by the semiexplicit scheme than by the fully implicit scheme, which exhibits spurious oscillations. However, the stability of the semi-explicit scheme imposes a constraint on the computational step, which must be smaller than a certain critical value. 4. Conclusion Under the circumstances concurring in electrical discharges and in other similar problems, where the variables undergo abrupt changes and fast variations in time, semi-explicit methods can help to reduce or avoid the spurious oscillations that implicit methods may exhibit. In particular, in the problem of a Trichel pulse propagating between a spherical cathode and a grounded plane, the use of a Crank-Nicolson time weighting (θ= 1/2) was able to accurately reproduce the sharp variations of electric field near the cathode without the need to introduce any artificial diffusion. Furthermore, the improvement of the stability shown by the semi-explicit method was not at the expenses of a longer computational time, compared to the fully implicit method. Acknowledgments This work has been supported by the European Regional Development Fund, under contracts FIS2011-25161 (Ministerio de Ciencia e Innovaci´on) and P09-FQM-4983 (Junta de Andaluc´ıa). References [1] Brooks A N and Hughes T J R 1982 Comput. Meth. Appl. Mech. Eng. 32 199. [2] Donea J and Huerta A 2003 Finite Element Methods for Flow Problems (Chichester, West Sussex: John Wiley & Sons, Ltd.) [3] Morrow R 1985 Phys. Rev. A 32 1799. 2nd International Conference on Mathematical Modeling in Physical Sciences 2013 IOP Publishing Journal of Physics: Conference Series 490 (2014) 012209 doi:10.1088/1742-6596/490/1/012209 4