Numerical Modeling of Partial Discharge Development Process

Abstract

Partial discharge (PD), a type of low-temperature plasma, indicates a discharge event that does not bridge the electrodes of an electrical insulation system under high voltage stress. It is common in power equipment, such as transformers, cables, gas-insulated switchgears, and so on. The occurrence of PD could deteriorate the insulation performance of the equipment, but, meanwhile, it is often used to diagnose the insulation status. Therefore, it is very necessary to clarify the PD mechanism, and through modeling the PD process, a better understanding of the phenomenon could be attained. Although PD is essentially a gas discharge phenomenon, it possesses some distinctive features, for example, very narrow discharge channel, short time duration, and stochastic behavior, which determine the simulation method of PD different from that for the other types of plasmas. This chapter seeks to propose a simulation method that could reflect the physical processes of PD development after introducing some background knowledge about PD and analyzing the shortcomings of existent models.

1. Introduction

Partial discharge (PD) is usually observed in power equipment, such as transformers, cables, gas insulated switchgears, and so on, which indicates a gas breakdown in essence induced by a local electric field distortion. It should be noted that it does not bridge the electrodes, differing from the gas breakdown across conductors. The remaining component of the whole insulation which does not suffer from PD could be oil, solid, or gas. On one hand, during the PD process, the heat energy, the charges with high velocity and chemical-active substances are released to erode and change the composite of the remaining component, leading to the deterioration of insulation performance and even the insulation failure. For example, as for high voltage power cable, PD resulting from the insulation defects could induce degradation of the solid dielectric due to chemical effect and physical attack by charge bombardment, and electric trees will be present after long-term service [1]. When the solid dielectric is across by the trees, an insulation fault takes place. On the other hand, PD parameters, such as discharge magnitude, discharge time, and so on, are determined by the characters of the gas and the remaining insulation. In terms of this, the PD measurement is often employed to diagnose the insulation status of power equipment. Whether understanding the negative effect of PD on insulation or equipment condition maintenance in the usage of PD measurement, it is based on the clear PD mechanism.

In essence, PD is a gas breakdown phenomenon. Similar to the other types of low-temperature plasmas, the temperature of electrons during a PD is much higher than that of ions, which is equivalent to the neutral gas molecules. However, PD also shows some distinctive features. For example, because PD always results from the local defect with a high electric field, the discharge channel is very narrow (the radius may be 100 μm) and the duration time is very short (several to tens of nanoseconds). During a PD sequence, once previous PD is terminated, and the subsequent one may take place after several milliseconds or even several days [2]. This phenomenon indicates that PD has a stochastic behavior, due to not only the effect of gas itself but also the interaction between gas breakdown and the remaining insulation. Therefore, as for the PD, the mere investigation of gas breakdown is meaningless. On the contrary, the interaction between PD and the remaining insulation should be considered. More importantly, a large number of PD data should be obtained to seek for its statistical characters because of its stochastic behavior.

According to the type of the remaining insulation and electrode configuration, PD could be divided into three categories [3]: internal discharge, surface discharge, and corona, as in Figure 1. Internal discharge indicates a gas breakdown taking place in a cavity embedded in solid or liquid dielectric. Generally, the former is more common. It consists of the streamer development and the interaction between streamer and cavity walls. A surface charge usually occurs along the solid dielectric surface due to a large tangential component of electric field, during which the interaction between streamer development and dielectric dominates. Corona often takes places in the local region around a conductor, which mainly involves the streamer development. Therefore, internal discharge could best represent PD, because it includes the two processes. In fact, the majority of PD simulations are concentrated on the internal discharge (also called cavity discharge) [4, 5, 6, 7]. And in this chapter, we also focus on it.

Figure 1.

Three categories of PD: (a) internal discharge, (b) surface discharge, and (c) corona.

There are many factors that could affect PD characters, such as the applied voltage (voltage waveform, amplitude, and frequency), electrode configuration, cavity (transportation parameters of gas, location, and size), remaining dielectric (permittivity, conductivity, and surficial parameters, e.g., morphology, surface trap distribution), and so on. To sum it up, two intrinsic factors behind them determine the evolution of PD behavior, that is, electric field and seed electrons. Generally speaking, two conditions must be simultaneously satisfied in order that a gas breakdown can take place: there must be at least one free electron in the gas, and the electric field must be of sufficient strength and duration time to ensure that this electron generates a sequence of avalanches [8]. Based on the conditions, it is inferred that the supply of free electrons and electric field affect not only the occurrence of PD but also its characters. Actually, the electric field is related to the applied voltage, electrode configuration, residual charges within the cavity, the cavity size, and the permittivity of remaining dielectric, while the supply of free electron depends on the gas status and surficial conditions of dielectric, corresponding to the volume generation and surface emission, respectively [9].

Looking back at the evolution of PD simulation methods, the a-b-c model was initially proposed [10, 11, 12], in which the discharge process was considered as charging-discharging of capacitors. Subsequently, some researchers held that the discharge could be represented by the increase of gas conductivity, and the current continuity equation was used to calculate discharge parameters [13, 14, 15]. On the contrary, others thought that a discharge was actually the deployment of charges in the cavity, and Poisson’s equation was enough [16, 17, 18]. Obviously, these models could represent the transient phenomenon of a discharge, but not reflect its physical processes. In recent years, a plasma model was employed to simulate single PD [19, 20, 21], in which the impact ionization, drift, diffusion, recombination, and other processes were quantitatively described by fluid equations. This model successfully obtained microscopic physical processes of a PD, but did not take the stochastic characters into account.

In this chapter, we firstly reviewed PD simulation models in brief, which consisted of the a-b-c model, Pedersen’s model, conductance model, Niemeyer’s model and plasma model, and analyzed their merits and drawbacks. Then, an advanced model was constructed to obtain physical processes, including the streamer propagation and surface charge dynamics, and macroscopic parameters, for example, discharge magnitude and moment of continuous PDs, so that a comprehensive analysis was available.

2. Review of PD simulation models

Since a-b-c model was proposed, numerical modeling of PD has been developed for decades of years. During this period, many kinds of simulation models have been constructed, which could be roughly divided into two categories: based on the point of view of circuit and based on the point of view of field. The former indicates a-b-c model and the latter consists of Pedersen’s model, conductance model, and Niemeyer’s model.

2.1. a-b-c model

The a-b-c model or the three-capacitor model is the original one to interpret the PD mechanism [3], and then it is usually employed to simulate the stochastic characters of PD [10, 11]. In the model, the dielectrics between electrodes, including the gas and solid insulation, are considered as capacitors, as in Figure 2. In detail, C₁ indicates cavity capacitance, C₂ is the capacitance of dielectric in series with the cavity, and C₃ is the capacitance of solid dielectric in parallel with the cavity. Besides, R₁, R₂, and R₃ indicate the resistance of corresponding part, respectively.

Figure 2.

a-b-c model.

The occurrence and termination of PD depend on the potential difference across the cavity, U₁. When U₁ exceeds the inception voltage, a discharge will take place and will stop when it is less than the extinction voltage. If a discharge occurs, C₁ is short-circuited, leading to a fast transient current to flow in the circuit due to a voltage difference between the voltage source and across C₂. Based on the analysis of capacitor charging-discharging processes, the apparent charge magnitude, which reflects PD intensity, could be calculated.

It could be found that this model is very simple, but it can represent the transient process related to a discharge event and is often used to explain some experimental results. However, it could not describe the discharge process physically, and the concept, capacitor, is not strictly valid, because the interface between the cavity and the solid dielectric is not equipotential when a discharge takes place [22].

2.2. Pedersen’s model

There are two important parameters of PD, that is, physical charges and apparent charges. The former indicate the charges generated during a discharge process, while the latter are measured charges through external circuit. In order to establish the link between physical charges and apparent charges, Pedersen proposed a model to describe the transient process [23]. Without considering the charge exchange between solid dielectric and the adjacent electrode, the amount of apparent charges equals the induced charges at an electrode surface due to charge generation, recombination, and movement during a discharge process. Therefore, if the physical charge distribution is known, the apparent charges could be calculated [24]

$$Q_{\mathrm{app}}=-\iiint \mathit{\lambda \rho }\mathrm{d}V-\iint \mathit{\lambda \sigma }\mathrm{d}s$$E1

where ρ and σ indicate volume and surface charge density within the cavity, respectively. λ, a dimensionless function, depends on the charge location, which satisfies Laplace equation

$$\nabla ·\left({\epsilon }_0{\epsilon }_{\mathrm{r}}\nabla \lambda \right)=0$$E2

where ε₀ is the vacuum permittivity, and ε_r the relative permittivity.

Pedersen’s model is helpful to understand the measured results by using the pulse current method. However, the apparent charges depend on physical charge distribution which results from the discharge process and keeps unknown in this model.

2.3. Conductance model

When PD takes place, a plasma region with a high charge concentration in the cavity is formed, so the gas conductivity largely increases in comparison with the initial state. Based on this fact, the discharge process is simplified by the variation of gas conductivity [13], which can be described by the following equations:

$$\nabla ·\mathit{D}=\rho$$E3

$$\nabla ·\mathit{J}+\frac{\partial \rho }{\partial t}=0$$E4

where D is the electric displacement field, J the free current density. At the initial state, the gas conductivity is set to be zero. When a discharge takes place, it is set to be γ_gd and hence the electric field distribution within the cavity changes. In terms of the electric field evolution, some PD parameters are obtained, for example, apparent charges and physical charges.

Forssen compared the simulation results with the experimental data, and they were in general agreement but with a slight difference. Furthermore, Illias developed the simulation model by taking the surface emission and temperature variation during the discharge into account [14]. However, in any case, the increment of gas conductivity could not represent the PD process.

2.4. Niemeyer’s model

Niemeyer considered PD within the cavity as a streamer-type discharge, because only this type could be detected and has engineering significance [9]. After analyzing the physical processes of PD, he proposed several equations to describe PD, as follows:

$$\underset{0}{\overset{x_{\mathrm{cr}}}{\int }}\overline{\alpha }\left[E\left(x\right)\right]\mathit{dx}\ge K_{\mathrm{cr}}$$E5

$$\Delta U_{\mathrm{res}}\approx E_{\mathrm{ch}}{l}_{\mathrm{str}}$$E6

$$q=\pm g{\mathit{\pi \epsilon }}_0\mathit{l\Delta }{U}_{\mathrm{PD}}$$E7

Eq. (5) is actually the well-known critical avalanche criterion, in which $\overline{\alpha }$, the function of electric field, indicates the effective ionization coefficient, K_cr the logarithm of a critical number of electrons that has to accumulate in the avalanche head to make the avalanche self-propagating by its own space charge field, and x_cr the distance within $\overline{\alpha }$ which exceeds zero. In terms of it, the inception field of PD occurrence could be obtained. Eq. (6) simply describes the streamer propagation, where E_ch is the electric field in the discharge channel, U_res the residual voltage instantaneously after discharge, and l_str the distance to which streamer could propagate. Eq. (7) establishes the relationship between physical charges and potential difference before and after a PD, in which g is a dimensionless proportionality factor and l the cavity scale.

Based on the model, Niemeyer simulated PD behaviors within a spherical cavity by considering the stochastic supply of free electrons, which agreed with experimental data qualitatively and quantitatively although there was a slight disagreement in the phase and magnitude distributions of PD. However, there is a significant shortcoming that the electric field distribution was assumed to be uniform within the cavity. Considering this point, Illias developed the simulation model in which the deployed charges were not uniform and Poisson’s equation was employed to calculate PD parameters [16, 17].

2.5. Plasma model

In terms of physical processes, a cavity PD is similar to the filamentary dielectric barrier discharge (DBD) [25]. As for the latter, fluid equations are widely used to simulate gas discharge process [26, 27], which describe the impact ionization, charge drift, diffusion, recombination, and some secondary effects. In recent years, several researchers employed them to simulate the PD occurring in a cavity [18, 19, 20]. For example, Novak and Bartnikas established a two-dimensional breakdown model based on the continuity equations for electrons and ions to examine the influence of surface charges upon the partial discharge behavior [19]. In terms of it, the evolution of electric field and charge concentration distribution within the cavity during the discharge process was obtained, as well as the discharge current pulse.

However, the behaviors of single PD could not represent that of continuous PDs due to the memory effect. On one hand, residual charges generated by previous discharge land on the cavity surface and affect the electric field distribution within the cavity, leading to the change of subsequent PD characters. On the other hand, the accumulated surface charges may provide free electrons for the next PD occurrence. The interaction between adjacent PDs could not be represented by singe PD. Therefore, it is necessary to establish a simulation model which could present the discharge development process and take the memory effect into account to obtain the stochastic characters of PD sequences.

3. Numerical modeling of PD sequences using fluid equations

As for the PD simulation, on one hand, the model should reflect physical processes as much as possible, and on the other hand, a large number of data should be obtained to get the statistical parameters of repetitive PDs due to the stochastic characters. There is a contraction that taking too much physical processes into account must result in the model complexity and a large calculation consumption which is not beneficial to statistical analysis. Therefore, some important processes should be considered in the simulation model, while others are abandoned.

By reviewing the PD simulation models, it is found that two processes are crucial to cavity PD characters, that is, streamer development and surface process. Obviously, the apparent charges that could be detected by pulse current method are determined by streamer development in the cavity. Surface process mainly consists of charge accumulation on the interface and surface emission of charge. After the streamer lands on the dielectric surface, charges accumulate and will affect the subsequent PD behavior. Besides, surface emission could provide free electrons for the next PD. It should be noted that the distribution of surface charges generated by previous discharge does not keep unchanged until subsequent one takes place. Due to the surface or bulk conductivity of dielectric, the accumulated charges may decay. To sum it up, the streamer development and surface charge accumulation reflect a single PD process, while surface charge accumulation, decay, and emission represent the interaction of adjacent discharges during a PD sequence, which should be considered in the simulation model.

3.1. Simulation model construction

Because sandwich-type samples are widely used in the experimental researches on PD, a cylindrical cavity with a diameter of 2 mm and a height of 0.25 mm is employed in our simulation model, as in Figure 3. The cavity, full of atmospheric pressure air, is embedded within the solid dielectric, of which the relative permittivity equals 2.3. The thickness of dielectric barriers is set to be identical to the cavity height. Although during the discharge process, the temperature of cavity may slightly increase due to the joule heating from discharges, the temperature variation is neglected in our model, which means that the pressure in the cavity keeps unchanged.

Figure 3.

Configuration of simulation model.

The streamer development is quantitatively described by fluid equations, as follows:

$$\frac{\partial N_{\mathrm{e}}}{\partial t}=N_{\mathrm{e}}\alpha \left|{\mathit{W}}_{\mathrm{e}}\right|-N_{\mathrm{e}}\eta \left|{\mathit{W}}_{\mathrm{e}}\right|-N_{\mathrm{e}}{N}_{\mathrm{p}}\beta -\nabla ·\left(N_{\mathrm{e}}{\mathit{W}}_{\mathrm{e}}-D\nabla N_{\mathrm{e}}\right)$$E8

$$\frac{\partial N_{\mathrm{p}}}{\partial t}=N_{\mathrm{e}}\alpha \left|{\mathit{W}}_{\mathrm{e}}\right|-N_{\mathrm{e}}{N}_{\mathrm{p}}\beta -N_{\mathrm{n}}{N}_{\mathrm{p}}\beta -\nabla ·\left(N_{\mathrm{p}}{\mathit{W}}_{\mathrm{p}}\right)$$E9

$$\frac{\partial N_{\mathrm{n}}}{\partial t}=N_{\mathrm{e}}\eta \left|{\mathit{W}}_{\mathrm{e}}\right|-N_{\mathrm{n}}{N}_{\mathrm{p}}\beta -\nabla ·\left(N_{\mathrm{n}}{\mathit{W}}_{\mathrm{n}}\right)$$E10

where N indicates the bulk charge concentration within the cavity, e, p, and n the symbols for electron, positive ion, and negative ion, respectively, t discharge time, α, η, β, and D denote the ionization, attachment, recombination, and electron diffusion coefficients, respectively, and W the drift velocity. Eqs. (8)–(10) reflect the transportation processes of electrons, positive and negative ions, which includes impact ionization, drift, diffusion, attachment, and recombination. However, the secondary processes, for example, photoionization, are neglected due to two reasons: (1) photoionization is crucial to the streamer development in long gaps but not so important for short gaps [28] and (2) the calculation of the secondary effect is extremely complicated, especially for the photoionization [29], which would bring about great difficulties of the PD sequence simulation. The detailed expressions of the above transport parameters come from Morrow’s paper [30], and we list them in Appendix A.

After the streamer arrives at the interface between the cavity and the dielectric, the charges will accumulate on the dielectric surface. We use the following equation to describe the transition from volume charges to surface charges:

$$\mathit{\sigma \Delta S}=\left(N_{\mathrm{p}}-N_{\mathrm{e}}-N_{\mathrm{n}}\right)\mathit{e\Delta V}$$E11

where ΔS and ΔV represent the area and volume of unit grid after meshing, respectively. Surface charge distribution is assumed to keep unchanged during the discharge process.

During the streamer development, the influence of space charges on the electric field should not be neglected, so Poisson’s equation is employed to obtain the electric field within the cavity:

$${\nabla }^2\phi =-\frac{e}{{\epsilon }_{\mathrm{r}}{\epsilon }_0}\left(N_{\mathrm{p}}-N_{\mathrm{e}}-N_{\mathrm{n}}\right)$$E12

At the upper and lower surfaces of cavity, the boundary conditions for Poisson’s equation are

$${\epsilon }_0{E}_{\mathrm{z}}\left(d_{\mathrm{g}}^{-}\right)-{\epsilon }_{\mathrm{r}}{\epsilon }_0{E}_{\mathrm{z}}\left(d_{\mathrm{g}}^{+}\right)={\sigma }_{\mathrm{u}}$$E13

$${\epsilon }_{\mathrm{r}}{\epsilon }_0{E}_{\mathrm{z}}\left(0^{-}\right)-{\epsilon }_0{E}_{\mathrm{z}}\left(0^{+}\right)={\sigma }_{\mathrm{d}}$$E14

where d_g is the cavity height, $E_z\left(d_{\mathrm{g}}^{-}\right)$ and $E_z\left(d_{\mathrm{g}}^{+}\right)$ indicate the z-component of electric field at both sides of the upper surface, while $E_z\left(0^{-}\right)$ and $E_z\left(0^{+}\right)$ represent the z-component of electric field at both sides of the lower surface, σ_u and σ_d denote the surface charge density at the upper and lower surfaces.

An initial electron-positive ion pair with a concentration of 10¹³ cm⁻³ is placed near the upper or lower surface to induce the streamer and avoid Townsend phase of gas discharge [28]. It should be noted that this assumption differs from the consideration of free electrons, which will be described in the later text. During the streamer development, charge concentration varies quickly, and an area with a steep concentration gradient appears at the head of the streamer. Meanwhile, the value of charge concentration should maintain positive, which cannot be guaranteed by the traditional finite difference method. So, the flux-corrected transport (FCT) algorithm is used to solve the convection term of charge continuity equations to overcome the two problems [31, 32, 33], which is listed in Appendix B.

In general, the time step for FCT is chosen based on the electron drift velocity, however, which may not apply to the circumstance in our simulation model. It is because apart from the streamer development, its extinguishment process also needs to be obtained which is responsible for the accumulation of electrons and ions. However, the drift velocity of electrons is about 100 higher than that of ions, and the choice of time step must lead to the large increase of calculation consumption at the later stage of discharge when ion drift dominates. Instead, if it is chosen based on the ion drift velocity, the accuracy of the calculation cannot be guaranteed at the initial stage of discharge. Therefore, as a compromise, the time step is set according to whether there are any electrons within the cavity volume. In detail, during the initial stage of streamer development, it is determined by the electron drift velocity. After electrons completely accumulate at the interface, it depends on the drift velocity of a positive ion or a negative one (both are the same). The expression for the time step is

$$\Delta t_{\mathrm{e},\mathrm{p}}=0.1\frac{\mathit{\Delta z}}{{\left|{\mathit{W}}_{\mathrm{e},\mathrm{p}}\right|}^{\max }}$$E15

where $∆t$ is the time step, $∆z$ the grid length along z-direction, and ${\left|\mathit{W}\right|}^{\max }$ the maximum value of charge drift velocity within the cavity.

According to Pedersen’s model, the apparent charges are determined by charge transportation within the cavity, which could be detected by pulse current method. However, due to the effect of dielectric barriers, the pulse obtained at the external circuit may not reflect the streamer propagation. So, we use Sato’s equation to calculate the current due to free charge movement [34], as follows:

$$I=\frac{e}{U_{\mathrm{a}}}{\iiint }_V\left(N_{\mathrm{p}}{\mathit{W}}_{\mathrm{p}}-N_{\mathrm{e}}{\mathit{W}}_{\mathrm{e}}-N_{\mathrm{n}}{\mathit{W}}_{\mathrm{e}}\right)·{\mathit{E}}_{\mathrm{a}}\mathit{dV}$$E16

where U_a indicates the applied voltage, E_a the applied field, and V the discharging volume.

On one hand, the field within the cavity should exceed a critical value so that a discharge may take place. Based on the ignition condition of streamer, the critical field is expressed as follows [35]:

$$E_b=\frac{24410\left(\frac{P}{760}{d}_{\mathrm{g}}\right)+6730\sqrt{\frac{P}{760}{d}_{\mathrm{g}}}}{d_{\mathrm{g}}}$$E17

where P is in Torr. After a discharge takes place, electrons and ions accumulate on the dielectric surface. Due to the recombination of charges from gas, surface, and bulk conduction of dielectric, the accumulated surface charges will decay until the next discharge occurs. It is found from our previous experiments that the decaying discipline of surface charges could be expressed as [36].

$$\frac{{\sigma }_{\mathrm{p}}}{{\sigma }_{\mathrm{p}0}}=e^{\frac{-t}{{\eta }_{\mathrm{p}}}}$$E18

$$\frac{{\sigma }_{\mathrm{e}}}{{\sigma }_{\mathrm{e}0}}=e^{\frac{-t}{{\eta }_{\mathrm{e}}}}$$E19

where σ_p0 and σ_e0 indicate initial positive charge and electron density at dielectric surfaces, respectively. η_p and η_e equal 312.5 and 568.8 ms, both of which represent the surface charge decay time for positive ions and electrons. The negative ion is neglected because its concentration is much lower in comparison with electron and positive ions.

On the other hand, although free electrons from the volume ionization and surface emission are formulated, their supply shows a strong scholastic behavior. Hence, there is usually a time delay between the instant of application of an electric field in excess of the critical field and the onset of breakdown, which is called a discharge time lag (strictly speaking, it is a statistical time lag, but the formative time lag is very short for cavity discharge and could be neglected). In order to simplify the physical process of free electron production, the discharge time lag is introduced to our model. Some experimental and simulation results show that the discharge time lag is not completely random, but is subject to exponential distribution [37, 38], which is expressed as

$$\tau =\left\{\begin{array}{l}\infty \\ -\text{ln}\left(1-P_{\mathrm{d}}\right)/\zeta \end{array}\right.\begin{array}{l}{E}_{\mathrm{z}}<E_{\mathrm{b}}\\ E_{\mathrm{z}}\ge E_{\mathrm{b}}\end{array}$$E20

where P_d indicates the discharge probability, which belongs to [0, 1) and is random, ζ the rate parameter of exponential distribution.

In terms of Eq. (17), the critical field for gas breakdown is calculated, and it equals 67,000 V/cm. In this case, the potential difference across the electrodes is 3130 V. Because the PD mechanism at AC voltage has been studied by many authors [2, 4, 5, 6, 7], and a comprehensive understanding about it has been obtained, the PD mechanism under DC voltage needs to be clarified. In this chapter, the DC voltage with an amplitude of 3200 V is applied to the anode, and the cathode is grounded all the time. Of course, this model is also applied to the circumstance of AC voltage application.

3.2. Simulation results

3.2.1. A PD development process

The process of PD development in the cavity consists of two stages: the streamer propagation and surface charge accumulation. Figures 4 and 5 show the temporal and spatial distribution of electrons and positive ions during this process, respectively. After discharge conditions are satisfied, the streamer is initiated near the lower surface of dielectric. With the help of applied field, electrons propagate toward the anode. At 0.72 ns, the head of streamer arrives at the upper surface of dielectric. Based on this, the streamer development velocity could be calculated, which equals 3.5 × 10⁷ cm/s and the order is in accordance with other researcher’s simulation result [39]. Then, electrons begin to accumulate on the upper surface of dielectric, and the density of surface charges reaches a saturation value after 1.4 ns. During this period, positive ions almost maintain stationary because the drift velocity is approximately 1/100 of the electron. However, positive ions seem to move according to Figure 5, and the distribution appearance looks like a ladle, which are attributed to the impact ionization of electrons. At 11.9 ns, a large number of positive ions land on the lower surface of dielectric, and the accumulation is terminated at 147.8 ns. Therefore, the accumulation time of electrons is much shorter than that of positive ions.

Figure 4.

Evolution of electron concentration distribution during the first PD (a) within the cavity volume (unit: cm⁻³) and (b) on the upper surface of the cavity (unit: cm⁻²).

Figure 5.

Evolution of positive ion concentration distribution during the first PD (a) within the cavity volume (unit: cm⁻³), (b) on the lower surface of cavity (unit: cm⁻²).

Based on the simulation results, it is found that the distribution of surface charges appears as a spot, and the maximum charge density locates at the middle of a spot. Compared with the experimental results [36], the distribution shape and surface density level (0.1 nC/mm²) are identical, which show that the simulation results are reasonable. However, there are some slight differences due to the simplification of model.

Charge transportation within the cavity will induce a current pulse, as in Figure 6, which could reflect the streamer development. The peak value of pulse appears at 0.72 ns; at this moment, the streamer head arrives at the upper surface of the cavity. The pulse width lasts for 1.4 ns; during this period, the accumulation of electrons is terminated. On the contrary, positive ions still move in the cavity volume. It is inferred that positive ions have a minor contribution to the current pulse because of their low drift velocity. A low-inductance resistor connected to the cathode is usually employed to detect a current, but this current slightly differs from that in Figure 6 [40].

Figure 6.

Current pulse waveform of the first PD obtained by simulation.

3.2.2. A PD sequence

A PD sequence consisting of 100 continuous discharges is obtained by the simulation (Figures 4–6 show the first discharge development process). Figure 7 shows the discharge time and the peak value of current of each discharge. In terms of this information, some statistical parameters of PDs, for example, discharge frequency and average discharge magnitude, could be calculated, and discharge patterns could be depicted.

Figure 7.

A PD sequence with 100 continuous discharges.

Besides, by analyzing the PD sequence, the interaction between adjacent discharges is obtained. Figure 8 shows the temporal evolution of surface charges and electric field within the cavity of first eight discharges. The first PD does not take place immediately after the voltage application due to the existence of a discharge time lag. After the discharge is terminated, the electric field within the cavity is dramatically reduced (as in Figure 8c), which is attributed to the effect of surface charge accumulation. Then, the surface charges begin to decay, and the electric field within the cavity gradually recovers. After it exceeds the critical value, and the condition for discharge time lag is satisfied, the next PD takes place.

Figure 8.

The first eight discharges during the PD sequence: (a) discharge time and magnitude, (b) surface charge decaying process, and (c) evolution of electric field within the cavity.

During the process of surface charge decaying (as in Figure 8b), the initial concentration of electrons and positive ions is approximately identical, but residual charges are completely distinct at the moment when a next discharge occurs. Due to the decay rate of positive charges faster than that of electrons, the concentration of residual negative surface charges is much higher. Therefore, compared with positive ions, residual electrons resulting from previous discharge have a larger influence on the subsequent one during a PD sequence.

4. Conclusions

PD, a type of low-temperature plasma, has some distinctive features, which determines its simulation method different from that of other types. In detail, as for the most representative PD type, cavity PD, it is necessary to take the streamer propagation, surface charge accumulation and decay, free electron supply into account so that the PD mechanism could be clarified. Besides, due to the stochastic character of PD, a large number of PD data must be obtained with the help of simulation.

Traditional simulation models about PD could be mainly divided into two categories: based on the point of view of circuit and based on the point of view of field. The former indicates a-b-c model, in which the discharge process is replaced by capacitor charging and discharging. The latter consists of Pedersen’s model, conductance model, and Niemeyer’s model, in which the discharge process is modeled by the variation of gas volume conductivity or significant simplification of discharge process. Anyway, these models could not reflect the PD development process physically.

Based on the simulation method for a single PD, we develop it by using fluid equations combined with Poisson’s equation. In terms of the model, microscopic physical processes, that is, streamer development and surface charge accumulation, could be obtained, as well as macroscopic parameters, that is, discharge current and discharge time, and the interaction between adjacent discharges. It is found that electrons and positive ions, respectively, land on the two surfaces of the cavity, and the accumulation time of positive ions is much longer than that of electrons. During a PD sequence, the decay of surface charges resulting from previous discharge could be considered to be the key factor, contributing to the occurrence of the subsequent one.

Acknowledgments

We would like to thank Prof. George Chen from the University of Southampton for the inspiring discussion and appreciate the financial support from the National Natural Science Foundation of China (51607128).

The transportation parameters for air are expressed by the following equations:

N = 2.69 × 10¹⁹ cm⁻³ indicates the number of gas molecules per unit volume, and E is the local field in V/cm

$$W_{\mathrm{e}}=\left\{\begin{array}{l}-\left(E/\left|E\right|\right)\left(7.4\times {10}^{21}\left|E\right|/N+7.1\times {10}^6\right)\left|E\right|/N>2.0\times {10}^{-15}\\ -\left(E/\left|E\right|\right)\left(1.03\times {10}^{22}\left|E\right|/N+1.3\times {10}^6\right){10}^{-16}\le \left|E\right|/N\le 2.0\times {10}^{-15}\\ -\left(E/\left|E\right|\right)\left(7.2973\times {10}^{21}\left|E\right|/N+1.63\times {10}^6\right)2.6\times {10}^{-17}\le \left|E\right|/N<{10}^{-16}\\ -\left(E/\left|E\right|\right)\left(6.87\times {10}^{22}\left|E\right|/N+3.38\times {10}^4\right)\left|E\right|/N<2.6\times {10}^{-17}\end{array}\right.$$EA1

$$W_{\mathrm{n}}=\left\{\begin{array}{l}-2.7E\left|E\right|/N>5.0\times {10}^{-16}\\ -1.86E\left|E\right|/N\le 5.0\times {10}^{-16}\end{array}\right.$$EA2

$$W_{\mathrm{p}}=2.34E$$EA3

$$\frac{\alpha }{N}=\left\{\begin{array}{l}2.0\times {10}^{-16}{e}^{\frac{-7.248\times {10}^{-15}}{\left|E\right|/N}}\left|E\right|/N>1.5\times {10}^{-15}\\ 6.619\times {10}^{-17}{e}^{\frac{-5.593\times {10}^{-15}}{\left|E\right|/N}}\left|E\right|/N\le 1.5\times {10}^{-15}\end{array}\right.$$EA4

$$\frac{{\eta }_2}{N}=\left\{\begin{array}{l}8.889\times {10}^{-5}\left|E\right|/N+2.567\times {10}^{-19}\left|E\right|/N>1.05\times {10}^{-15}\\ 6.089\times {10}^{-4}\left|E\right|/N-2.893\times {10}^{-19}\left|E\right|/N\le 1.05\times {10}^{-15}\end{array}\right.$$EA5

$${\eta }_3/N^2=4.7778\times {10}^{-59}{\left(\left|E\right|/N\right)}^{-1.2749}$$EA6

$$\eta ={\eta }_2+{\eta }_3$$EA7

where η₂ and η₃ are the two-body and three-body attachment coefficients, respectively.

$$\beta =2.0\times {10}^{-7}$$EA8

$$D=0.3341\times {10}^9{\left(\left|E\right|/N\right)}^{0.54069}\left|W_e/E\right|$$EA9

Based on the axisymmetric character of sample configuration in our model, the cylindrical coordinate system is employed, so the convection term could be rewritten as

$$\frac{\partial w}{\partial t}=-\frac{\partial f}{\partial r}-\frac{\partial g}{\partial z}$$EA10

where w = rN, f = rNW_r, g = rNW_z, W = W_re_r + W_ze_z, e_r, and e_z are the unit vectors along r and z directions, respectively. To solve this equation, six steps are needed:

(1) to obtain the low-order flux ${F^L}_{i+\frac{1}{2},j}$, ${G^L}_{i,j+\frac{1}{2}}$

$${F^L}_{i+\frac{1}{2},j}=\left(2i+1\right)\mathit{\text{πΔzΔrΔt}}{\left(W_r\right)}_{i+\frac{1}{2},j}w\left\{\begin{array}{l}w=\mathit{i\Delta }{\mathit{rN}}_{i,j}{\left(W_r\right)}_{i+\frac{1}{2},j}\ge 0\\ w=\left(i+1\right)\Delta {\mathit{rN}}_{i+1,j}{\left(W_r\right)}_{i+\frac{1}{2},j}<0\end{array}\right.$$EA11

$${G^L}_{i,j+\frac{1}{2}}=\Delta {\mathit{tS}}_i{\left(W_z\right)}_{i,j+\frac{1}{2}}w\left\{\begin{array}{l}w=\mathit{i\Delta }{\mathit{rN}}_{i,j}{\left(W_z\right)}_{i,j+\frac{1}{2}}\ge 0\\ w=\mathit{i\Delta }{\mathit{rN}}_{i,j+1}{\left(W_z\right)}_{i,j+\frac{1}{2}}<0\end{array}\right.$$EA12

$$S_i=2\mathit{i\pi \Delta }{r}^2$$EA13

$${\left(W_r\right)}_{i+\frac{1}{2},j}=\frac{{\left(W_r\right)}_{i,j}+{\left(W_r\right)}_{i+1,j}}{2}$$EA14

$${\left(W_z\right)}_{i,j+\frac{1}{2}}=\frac{{\left(W_z\right)}_{i,j}+{\left(W_z\right)}_{i,j+1}}{2}$$EA15

where i and j are the sequence number of node along r and z directions, respectively.

(2) to obtain high-order flux ${F^H}_{i+\frac{1}{2},j}$, ${G^H}_{i,j+\frac{1}{2}}$

$$\begin{array}{l}{F^H}_{i+\frac{1}{2},j}=\left(2i+1\right)\mathit{\text{πΔzΔrΔt}}[\frac{533}{840}\left(f_{i+1,j}+f_{i,j}\right)-\frac{139}{840}\left(f_{i+2,j}+f_{i-1,j}\right)\\ +\frac{29}{840}\left(f_{i+3,j}+f_{i-2,j}\right)-\frac{1}{280}\left(f_{i+4,j}+f_{i-3,j}\right)]\end{array}$$EA16

$$\begin{array}{l}{G^H}_{i,j+\frac{1}{2}}=S_i\mathit{\Delta t}[\frac{533}{840}\left(g_{i,j+1}+g_{i,j}\right)-\frac{139}{840}\left(g_{i,j+2}+g_{i,j-1}\right)+\frac{29}{840}\left(g_{i,j+3}+g_{i,j-2}\right)\\ -\frac{1}{280}\left(g_{i,j+4}+g_{i,j-3}\right)]\end{array}$$EA17

$$f_{i,j}=\mathit{i\Delta }{\mathit{rN}}_{i,j}{\left(W_r\right)}_{i,j}$$EA18

$$g_{i,j}=\mathit{i\Delta }{\mathit{rN}}_{i,j}{\left(W_z\right)}_{i,j}$$EA19

(3) to define antidiffusion flux

$$A_{i+\frac{1}{2},j}\equiv {F^H}_{i+\frac{1}{2},j}-{F^L}_{i+\frac{1}{2},j}$$EA20

$$A_{i,j+\frac{1}{2}}\equiv {G^H}_{i,j+\frac{1}{2}}-{G^L}_{i,j+\frac{1}{2}}$$EA21

(4) to obtain the temporary solution

$$w_{i,j}^{\mathit{td}}=w_{i,j}^n-\Delta V_{i,j}^{-1}\left[F_{i+\frac{1}{2},j}^L-F_{i-\frac{1}{2},j}^L+G_{i,j+\frac{1}{2}}^L-G_{i,j-\frac{1}{2}}^L\right]$$EA22

$$w_{i,j}=\mathit{i\Delta }{\mathit{rN}}_{i,j}$$EA23

$$\Delta V_{i,j}=\Delta {\mathit{zS}}_i$$EA24

(5) to restrict the antidiffusion flux

If $A_{i+\frac{1}{2},j}\left(w_{i+1,j}^{\mathit{td}}-w_{i,j}^{\mathit{td}}\right)<0\&\left\{A_{i+\frac{1}{2},j}\left(w_{i+2,j}^{\mathit{td}}-w_{i+1,j}^{\mathit{td}}\right)<0\Vert A_{i+\frac{1}{2},j}\left(w_{i,j}^{\mathit{td}}-w_{i-1,j}^{\mathit{td}}\right)<0\right\}$

$$A_{i+\frac{1}{2},j}=0$$EA25

If $A_{i,j+\frac{1}{2}}\left(w_{i,j+1}^{\mathit{td}}-w_{i,j}^{\mathit{td}}\right)<0\&\left\{A_{i,j+\frac{1}{2}}\left(w_{i,j+2}^{\mathit{td}}-w_{i,j+1}^{\mathit{td}}\right)<0\Vert A_{i,j+\frac{1}{2}}\left(w_{i,j}^{\mathit{td}}-w_{i,j-1}^{\mathit{td}}\right)<0\right\}$

$$A_{i,j+\frac{1}{2}}=0$$EA26

$$A_{i+\frac{1}{2},j}^C=C_{i+\frac{1}{2},j}{A}_{i+\frac{1}{2},j}0\le C_{i+\frac{1}{2},j}\le 1$$EA27

$$A_{i,j+\frac{1}{2}}^C=C_{i,j+\frac{1}{2}}{A}_{i,j+\frac{1}{2}}0\le C_{i,j+\frac{1}{2}}\le 1$$EA28

$$C_{i+\frac{1}{2},j}=\{\begin{array}{l}\text{min}\left(R_{i+1,j}^{+},R_{i,j}^{-}\right)A_{i+\frac{1}{2},j}\ge 0\\ \text{min}\left(R_{i,j}^{+},R_{i+1,j}^{-}\right)A_{i+\frac{1}{2},j}<0\end{array}$$EA29

$$C_{i,j+\frac{1}{2}}=\{\begin{array}{l}\text{min}\left(R_{i,j+1}^{+},R_{i,j}^{-}\right)A_{i,j+\frac{1}{2}}\ge 0\\ \text{min}\left(R_{i,j}^{+},R_{i,j+1}^{-}\right)A_{i,j+\frac{1}{2}}<0\end{array}$$EA30

$$w_{i,j}^a=\text{max}\left(w_{i,j}^n,w_{i,j}^{\mathit{td}}\right)$$EA31

$$w_{i,j}^{\max }=\text{max}\left(w_{i-1,j}^a,w_{i,j}^a,w_{i+1,j}^a,w_{i,j-1}^a,w_{i,j+1}^a\right)$$EA32

$$w_{i,j}^b=\text{min}\left(w_{i,j}^n,w_{i,j}^{\mathit{td}}\right)$$EA33

$$w_{i,j}^{\min }=\text{min}\left(w_{i-1,j}^b,w_{i,j}^b,w_{i+1,j}^b,w_{i,j-1}^b,w_{i,j+1}^b\right)$$EA34

$$P_{i,j}^{+}=\text{max}\left(0,A_{i-\frac{1}{2},j}\right)-\text{min}\left(0,A_{i+\frac{1}{2},j}\right)+\text{max}\left(0,A_{i,j-\frac{1}{2}}\right)-\text{min}\left(0,A_{i,j+\frac{1}{2}}\right)$$EA35

$$Q_{i,j}^{+}=\left(w_{i,j}^{\max }-w_{i,j}^{\mathit{td}}\right)\Delta V_{i,j}$$EA36

$$R_{i,j}^{+}=\{\begin{array}{l}\text{min}\left(1,Q_{i,j}^{+}/P_{i,j}^{+}\right)P_{i,j}^{+}>0\\ 0P_{i,j}^{+}=0\end{array}$$EA37

$$P_{i,j}^{-}=\text{max}\left(0,A_{i+\frac{1}{2},j}\right)-\text{min}\left(0,A_{i-\frac{1}{2},j}\right)+\text{max}\left(0,A_{i,j+\frac{1}{2}}\right)-\text{min}\left(0,A_{i,j-\frac{1}{2}}\right)$$EA38

$$Q_{i,j}^{-}=\left(w_{i,j}^{\mathit{td}}-w_{i,j}^{\min }\right)\Delta V_{i,j}$$EA39

$$R_{i,j}^{-}=\{\begin{array}{l}\text{min}\left(1,Q_{i,j}^{-}/P_{i,j}^{-}\right)P_{i,j}^{-}>0\\ 0P_{i,j}^{-}=0\end{array}$$EA40

(6) to solve the charge concentration

$$w_{i,j}^{n+1}=w_{i,j}^{\mathit{td}}-\Delta V_{i,j}^{-1}\left[A_{i+\frac{1}{2},j}^C-A_{i-\frac{1}{2},j}^C+A_{i,j+\frac{1}{2}}^C-A_{i,j-\frac{1}{2}}^C\right]$$EA41

where n indicates $n∆t$ and n+1 indicates $\left(n+1\right)∆t$.