Space and Time Domain Finite Volume Method for Numerical Simulation of Negative Corona Discharge in Air

1. Introduction

Corona discharge occurs when the voltage applied to the corona electrode exceeds a critical value [1]. With the rapid development of extra high-voltage transmission lines, a corona discharge in air has become one of the critical problems, which leads to deterioration of insulation systems, power loses, radio noises, etc. [2, 3]. In the meantime, corona discharge has been employed as an effective tool in many industrial applications such as electrostatic precipitators, disinfection, ozone generators, surface treatment, ignition, and electrostatic painting [4-6]. The features of a corona discharge vary considerably among different gases. A corona discharge in electronegative gas can be classified into three phases [7]: (1) auto-stabilization, (2) Trichel pulse, and (3) continuous current discharge. Trichel pulses [8] are normally regular pulses with a short duration (tens of ns) and a pulse period of several ms. Trichel pulses only occur in electronegative gases, such as air [9], whereas in electropositive gas, such as nitrogen, a corona discharge exists in continuous current without a pulse. Research on the characteristics of negative corona discharge in electronegative gas has been promoted.

Trichel [8] found the existence of the Trichel pulse in oxygen. Pulse frequency is related with corona current and cathode diameter, whereas it is unaffected by gap length. Trichel also explained the shielding effect produced by the positive ions near the cathode. Several studies on the Trichel pulse have been conducted. Loeb et al. [9] affirmed the existence of the Trichel pulse only in electronegative gases. Lama and Gallo [10] designed several experiments to determine the dependence of pulse frequency and corona current on applied voltage, cathode diameter, and gap length.

For the convenience of numerical simulation, some models have been established to study the Trichel pulse. In 1985, Morrow [11] proposed a 1-D model that combines finite difference method and flux-corrected transport technology to simulate Trichel pulses. This model presents theoretical explanations for the corona discharge phenomena and the different stages of a Trichel pulse. However, this method can only simulate the first Trichel pulse and ignore secondary electron emissions. Stattari [12] introduced a new 2-D numerical model that considers three ionic species. An optimized technique is utilized to raise calculation efficiency and decrease calculation time. The electrical field and densities of three ionic species during different Trichel pulse phases are presented. The properties of Trichel pulse under different applied voltages are analyzed. The present paper provides references and an optimized technique that can be used to decrease calculation time. Yin [13] proposed a line-to-plane geometry to simulate Trichel pulses in a negative corona. A hybrid numerical algorithm is applied. Finite element method (FEM) is used for Poisson’s equations, and finite volume method (FVM) is used to solve a charge transport equation. The distribution of three ionic species is calculated only in the region near the line. However, an artificial boundary is set and the calculation area is only limited in the region near the transmission line, which results in errors. The space charges outside the artificial boundary are neglected and have effects on the distribution of ionic species in the calculation region.

Previous studies [14-18] proposed numerical methods using FEM and FVM to calculate the ion flow under a high-voltage current transmission line. Yin [19] combines FVM with a time step to simulate the formation process of the ion-flow field of the ±800 kV HVDC test line. Liao and Wu [20, 21] proposed an improved hybrid numerical model that considers six kinds of plasma chemical reactions, photo ionization, and secondary electron emission. This model presents the variations of electrical field and space charge density and focuses on the discussion of heavy ion transport properties. However, the chemical reactions involved in the present study are rather complicated and very time consuming.

In this study, a space and time domain FVM is proposed to simulate the negative corona discharge in the air. FEM is applied for Poisson’s equation and FVM is used in the charge transport equation. In diverse pulse phases and discharge areas, different time steps are applied to obtain the distribution of electrical field, electrons, and ion densities. The corona current is discussed and compared with the testing results to justify the proposed model.

 

2. Mathematical Model

2.1 Governing equations

Fig. 1 shows the configuration of the bar-to-plate system. A hemisphere bar with radius of curvature R and length L is mounted perpendicularly to an infinitely large plate at a distance d. A negative DC potential V is produced by a high-voltage power supply and applied to a bar electrode. The corona currents through the plate electrode are measured by the current acquisition board, and then showed on the oscilloscope.

Fig. 1.Calculation model of the bar-to-plate geometry

The governing equations in the gap between the bar and the plate are as follows [13]:

where φ is the electric potential, E is the electrical field, ε0 is the air permittivity, and e is the electron charge, which is 1.6×10−19 C; Ne, Np, Nn are the densities of electron, positive ions, and negative ions, respectively. ve, vp, and vn are the drift velocities of electrons, positive ions, and negative ions, which can be obtained by [22]:

where μp and μn are the mobility of positive and negative ions, taking as 1.4×10−4 and 1.8×10-4m2·V−1·S−1. α, β, and η are the ionization, recombination, and attachment of coefficients, which can be calculated as bellow [23]:

Eq. (1) is a Poisson’s equation on electric potential. Eqs. (2), (3), and (4) are charge continuity equations. They show the transport properties of electrons, positive ions, and negative ions.

2.2 Boundary condition

The boundary conditions in solving Poisson’s equation on voltage are:

φ=V on the corona electrode φ=0 on the ground plate

where V is the applied voltage on corona electrode.

The boundary conditions of the charge continuity equation are:

Np=0 on the ground plate Nn=0 on the corona electrode

For electrons, the secondary emission considered to update the electron density on corona electrode is

where Nec and Npc are the densities of electrons and positive ions on the corona electrode, and γ is the secondary electron emission coefficient, which is 0.01.

2.3 Initial condition

A Gaussian distribution of seed electrons with a maximum value of 1016 m−3and width of 25 µm is set at the tip according to the following Eq. [24]:

where Nmax = 1016 m−3, (x0,y0) is the position of bar tip, s0 = 25 μm. Ne(t=0) is the electron density at t=0. The initial condition only accelerates the pulse formation and does not change the discharge characteristics [25, 26].

 

3. Algorithm Overview

Poisson’s equation describes the variation of electric potential, whereas the charge continuity equation shows the distribution of three ionic species. In this paper, FEM is used to solve the Poisson’s equation, and FVM is used to solve the charge continuity equation. To improve calculation efficiency, an optimized technique based on space and time domain is utilized.

3.1 Solution of poisson’s equation

The typical triangle meshes (elements) are used as shown in Fig 1. The solution of Equation (1) is equal to minimize the following function [14]:

where the subscript i represents the ith element and Ω is the region corresponding to the ith element. Through the derivative of Fi in all elements, a linear equation can be obtained. The electric potential can be solved by combining with the boundary condition.

3.2 Solution of charge continuity equation

Charge continuity equations describe the transport properties of three ionic species under the effect of an electrical field. They also explain the interactions between particles, including ionization, attachment, and recombination. The auxiliary meshes (control volume) are constructed by connecting the center of each triangle and midpoint of each edge. These three equations are coupled by the interaction effect between the three ionic species. When solving Eq. (2), the densities of negative ion and positive ion in the last time step are used to decouple the equation. The same solution is used to the other two equations.

First, the integral of Eq. (2) is obtained in the control volume around the ith node, and then the equation is discretized in time domain, turning differential equation to a linear equation.

When integrating the second terms on the left of the equal sign in Eq. (2) around the control volume, divergence theorem can be used. Integrating in the control volume can be changed into integrating around the edge of the control volume. A number of interpolation algorithms can be applied to calculate charge densities on the edges. When dealing with the convective problem or the charge transport problem in this paper, central difference shows low accuracy and poor convergence. Thus, upwind difference is used to calculate the charge densities on the edges. The equation for the calculation is as follows [24]:

where ρi is the charge density of element i; and ρx is the charge density of element which shares edge with element i. When integrating in the time domain, three accessor methods, namely, the explicit format, the Crank–Nicolson format, and the fully implicit format, can be mainly applied to determine the density of ionic species as shown in Eq. (16).

Considering convergence and stability, the fully implicit format is chosen. The density of ionic species is taken as the value of the latter time step when integrating in the time domain. Finally, the following equation can be obtained:

where the subscript i represents the ith node, im represents the edge of the control volume around the ith node, M is the number of edges, Lm is the length of edges, S is the area of control volume around the ith node, t represents current time step, and ∆t is the time step. The equations of the positive ion and the negative ion can be solved similarly.

3.3 Space domain and time domain for time step

The time step used in the model is usually quite short which causes the entire simulation process to be time consuming. The time step mainly depends on the charge velocity v and satisfies the following formula:

where ∆d demonstrates the mesh size. For electrons, a time step of 1 ns is suitable. However, for positive and negative ions, 100 ns is suitable. The period of a single pulse is in the order of 1 ms. If 1 ns is applied as time step, the calculation steps should be in the order of 106 to 107, which is not advisable. As a result, the optimization technique should be used to reduce simulation time.

At different pulse phases, the properties of three ionic species are discrepant. Positive ions and electrons are active during the phase of half pulse rising slope, whereas negative ions exist in almost the whole pulse.

At different gap areas, various active ionic species are present. Positive ions are active near the corona electrode, whereas negative ions are almost filled in the whole gap. Electrons are mainly located in the area between positive and negative ions, and their density remains at a low level during the two pulses. In the later period of a single pulse, negative ions are far away from positive ions. As a result, negative ions have a little influence on the electrical field near the corona electrode, and positive ions have little influence on the electrical field away from the corona electrode.

As shown in Fig. 2, different calculation models can be utilized according to space and time domain. Mode A for calculation in the area near the corona electrode (area A1 in Fig. 1), and model B for calculation in the area far away from the corona electrode (area A2 in Fig. 1). In model A, three ionic species are involved and time step is 1 ns, whereas only negative ion is involved in model B, in which the time step is 100 ns. When the density of electrons in the gap decreases to zero, electrons can be neglected. Model C is used now for calculation, in which only negative and positive ions are considered. The calculation area includes the whole gap (area A1 and area A2), and the time step is 100 ns. The time domain is divided according to the electron density. Models A and B are applicable when the total charge of electrons in the gap is high. When the total charge of electrons in the gap has reduced to zero, model C is applicable. The space domain is divided according to the degree of ionization. In model A, ionization is greater than attachment (α>η). In model B, attachment is greater than ionization (α<η).

Fig. 2.Flow chart of the algorithm

To test this optimized algorithm, corona discharge in line-cylinder geometry is simulated in two method (using optimized algorithm and not using this algorithm). The results shows that this optimized technique is effective.

 

4. Validation

4.1 Trichel pulse

A corresponding experiment model is established in the laboratory. In this study, tip radius is 0.3 mm, and gap distance is set as 3.3 mm.

The experimental and simulated results of repeatable pulses and single pulse under −1 kV are shown in Fig. 3 and Fig. 4; the simulated results are the Trichel pulses in steady state. A considerable difference in the first several pulses is found, which is removed and not shown in Fig. 3. The amplitude of the first pulses is much larger than the rest of the pulses which agrees with previous literature [12]. The first few pulses are produced in a charge-free space. The negative ion cloud in the bar-to-plate gap reduced the electrical field between the negative ion cloud and the corona electrode. As a result, ionization becomes weaker, whereas negative ions do not exist during the gap in the first pulse. Thus, ionization is much more intense and additional electrons are generated by avalanche ionization, which leads to a larger corona pulse in the first pulse. To compare the simulation and experiments results, the first few pulses in the simulated results are removed and the pulses in steady state are shown. In addition, the corona currents under different voltages are calculated and measured as shown in Fig. 5. Corona current increases with increase of applied voltage, which agrees with previous literatures.

Fig. 3.Repeatable pulse

Fig. 4.Curve of a single Trichel pulse

Fig. 5.Corona current under different applied voltage

4.2 Axial distribution of three ionic species and electric field

To study the pulse in detail, four classical time points are marked to show the feature of a single Trichel pulse at four different phases, as shown in Fig. 4, where t1 indicates the beginning of pulse, t2 is for half-pulse rising slope of the pulse, t3 is for half-pulse decreasing slope of the pulse, and t4 is for end of the pulse. In Fig. 6-Fig. 9, the horizontal coordinate demonstrates the distance away from corona electrode along the axis, 0 mm represents the tip of corona electrode, whereas 3.3 mm represents the ground electrode.

Fig. 6.Density of positive ions along the axis of symmetry of the four points.

Fig. 7.Density of negative ions along the axis of symmetry of the four points.

Fig. 8.Electron density along the axis of symmetry of the four points

Fig. 9.Electrical field along the axis of symmetry of the four points.

The feature of each phase is closely related to the behaviour of three kinds of ionic species and the distribution of electrical field strength. To study the Trichel pulse further, determining the distribution of ionic species and electrical field strength is essential. Fig. 6 - Fig. 9 show the distributions of the electron density, positive ion density, negative ion density, and electrical field strength along the axis of symmetry of the four time points under −1 kV.

At t1, negative and positive ions remain at a relatively low level. The positive and negative ions are relatively few at the tip of corona electrode. With a short distance away from corona electrode, the number of positive and negative ions increases rapidly. Negative ions have the same polarity as the corona electrode; thus, the electric repulsion moves the negative ions away from the surface of corona electrode. Closer distance to the corona electrode entails greater repulsion; a more rapid movement of ions means lower ion density. Finally, negative ions gather dozens of micrometer away from the surface of the corona electrode surface to reach a relatively high order, whereas for positive ions, their polarity is opposite to the corona electrode. The intensity of the electrical field near the surface of the corona electrode then increases, which accelerates the movements of positive ions towards the corona electrode, and finally force them to impact the surface of the corona electrode. This secondary electron emission leads to the neutralization of positive ions in the surface of the corona electrode. As a result, positive and negative ions remain at a low level in the surface of corona electrode.

The number of negative ions generated in the last pulse decreases to a low level, which leads to the increase of the electrical field around the corona electrode, as shown in Fig. 9. The enhanced electrical field near the corona electrode can promote the generation of electrons and electron avalanches. With the rapid increase of electrons, a new pulse is developed and the drift velocity of electrons under the effect of electrical field increases. The corona currents seem to step to its amplitude, as shown in Fig. 4.

At time point t2, the corona currents reach their amplitude. As seen in the electric field distortion in Fig. 9, two extreme values at time point t2 exist along the axis of the symmetry. This observation is closely related to negative ions in the gap. Electrons accelerate under electrical field and impact air molecules to produce new electrons between t1 and t2. As a result, electron density reaches its maximum in point t2. With increasing distance from the corona electrode, the electrical field drops and ionization is weaker, whereas attachment intensifies. Negative ions are generated under the attachment of electrons. The intensity of the attachment is proportional to the number of negative ions. The density curve of negative ions is similar to that of electron density. Numerous negative ions gather to form a negative ion cloud in the area near the corona electrode. The negative ion cloud moves away from the corona electrode under applied electrical field. Having the same polarity as the corona electrode, the negative ion cloud reduces the electrical field between them. As a result, electric field distortion appears.

By comparing the curves at the four time points in Fig. 7 and Fig. 8, the amplitude electron density reaches a maximum at time point t2 corresponding to the amplitude of corona currents at time point t2, which then declines gradually. While the intensity of ionization weakens, the intensity of attachment strengthens. This trend leads to a sustainable growth of the amplitude of negative ions produced by the electron attachment. Consequently, the existence of negative ions reduces the electrical field strength. As a result, the amplitude of negative ions continues to increase after time point t2. Positive ions are produced according to the impact of electrons. During the half-pulse rising slope of the pulse, ionization is more intense than attachment, especially in the area near corona electrode. Thus, positive ions are almost two orders of magnitude higher than negative ions.

At time point t3, the electric field distortion intensifies with the growth of negative ions, as shown in Fig. 9. The electrical field strength and the densities of electrons, negative ions, and positive ions in the whole gap are relatively well-distributed, especially in the region away from the corona electrode (Fig. 6-Fig. 9). The amplitudes of electrons and positive ions decrease with respect to the intensity of ionization and attachment. In addition, the density curves of the negative and positive ions are similar to that of the electron density.

The electrons move much faster than positive and negative ions. During the movement through the gap, the electrical field is at a relatively high level, i.e., ~30 kV/cm, which is the breakdown field strength of air. As a result, ionization and attachment continue to proceed and produce a mass of new negative and positive ions in the motion path of electron. Ionization and attachment are closely related to the electrical field, which can be described by the ionization and attachment coefficients in Eqs. (9) and (11). The electrical field that is well-distributed in the gap affects the intensity of the ionization and attachment. Accordingly, positive and negative ions are distributed uniformly.

The time phase t4 is the end of the Trichel pulse, all of the three ionic species decrease to a low level. Negative ions diffuse under electrical field. With the increase in distance of the negative ion cloud from the corona electrode, the electric field distortion becomes weaker and the electrical field has a tendency of recovering to time point t1. When it satisfies a critical condition, a new pulse begins. In addition, two-dimensional distribution of three kinds of particles are presented in Fig. 10. As you can see, spatiotemporal distributions of electron, negative ion and positive ion agree with the change rule as mentioned above.

Fig. 10.Two-dimensional distribution

 

5. Conclusion

An improved hybrid model that combines FEM and FVM is proposed to simulate the Trichel pulse with a bar-to-plate geometry under negative DC and high voltage. FEM is used for calculation of the Possion’s equation, and the FVM is for calculation of the charge transport equation. An optimized method based on space and time domain is utilized. The optimized method on time step according to the time and space domain can be used to improve the numerical model.

Density distributions of electrons, positive ions, and negative ions, as well as electrical filed strength along the axis of the symmetry in the Trichel pulse are determined. Four classical time points are marked to focus on the feature of the Trichel pulse, especially the characteristic of the three ionic species. Electrons are the most active species, which are directly related to ionization and attachment. The ionization of electrons produces an electron avalanche, which initiates the pulse. The attachment of electron produces a negative ion cloud, which reduces the electrical field and terminates development of the pulse. The feature of the three ionic species is impressed by the electrical field. However, this feature affects the distribution of the electrical field. These effects are coupled by the interaction of ionization and attachment under the electrical field. A corresponding model is established in the laboratory. A good agreement between simulated and experimental results verifies the validity of numerical model.