A Finite-Difference Time-Domain Formulation for Modeling Air Ionization Breakdown Streamers

Abstract

A numerical model based on the finite-difference time-domain (FDTD) method is developed for calculating discharge currents and fields produced due air ionization between two electrodes. Electric conductivity σ of the ionized channel between two electrodes is calculated using non-linear equations, i.e., σ is obtained using functions of time. The functions for calculating σ are adapted and proposed in this work for each stage of the plasma formation and disrupting processes. The proposed numerical model is validated by making comparisons among the obtained FDTD results and experimental data published in literature. In order to perform the validation procedures, an experimental setup from literature is numerically reproduced using the proposed FDTD methodology.

Index Terms
Air Ionization; Discharge Streamer; FDTD method; Plasma Channels

I. INTRODUCTION

When two electrodes are set up such as in the needle-plane configuration [¹]–[³], immersed in a fluid and they are subjected to a significantly high electric voltage, a discharge is started in the vicinity of the thinner electrode. Discharges take place because the atoms of the fluid do not support a certain critical electric field threshold level E_c, which may be reached and exceed by electric field produced in the vicinity of the electrodes, so that the electrons closest to the surface of the thinnest electrode start to detach due to Coulomb forces [⁴]. These electrons collide with adjacent atoms in the fluid, causing more electron releases, and so on. This process, which evolves with successive ionization of atoms in the fluid, is known as the Townsend avalanche [⁵]–[⁷]. Thus, the atoms become electrically charged, i.e., they are transformed into ions. Townsend avalanche exponentially increase amounts of electrons with kinetic energy, which in turn constitute the electric currents between the electrodes.

When the voltage applied between the electrodes is sufficient to generate electric field strength greater than E_c, the increasing amount of free electrons may reach and exceed the critical number of 10⁸ around the anode electrode tip. At this stage, there is a sharp deformation of electric field spatial distribution with respect to the field produced by the external source (in the space between the conductors), giving rise to the formation of the so called streamer, which is a plasma channel of poor conductivity also referred to as streamer channel. At the beginning of the streamer channel formation, there is an exponential increasing of the channel current. The streamer then increases in length, with its tip being pushed towards the cathode electrode. Once the streamer channel reaches the cathode, the conducting channel electrically connects the high voltage electrodes, constituting a partial short-circuit, thus giving rise to a velocious exponential intensification of the electric current due to fast-increasing number of releases of electrons from atoms to the channel [¹], [²]. During this stage, expressive particle velocity increases and substantial and fast electric conductivity elevation are observed, along with plasma channel heating up and limited volume expansion. As a consequence, there is a redistribution of electric field levels over the gap space between anode and cathode. Because of the massive increase of the number of electrons in the channel, electron-ion re-associations are favored. Those phenomena cause progressive deionization of the plasma channel, in which electric conductivity and electric current decay exponentially, caused mainly by the re-associations of electrons to ionized oxygen molecules [¹]. A current function is thus formed over time, with specific profile depending on the just described non-linear phenomena.

The interest of the scientific community on this matter has been growing recently in literature [³], [⁸], [⁹] because predictive methods can be developed for antecipating the occurence of electric faults in high voltage systems [¹⁰]–[¹⁴]. In [³], it is presented a two-dimensional finite element method for calculating discharge electrostatic fields using Poisson’s equation, the parameter number of electrons and energy. The method described in [³] requires fine discretization around the tip of the needle electrode and around the streamer head, demanding considerable computational resources especially if three-dimensional simulations need to be executed. Since time is not explicitly considered, although advance of the streamer head is computed through a series of field-independent simulations, magnetic field is not calculated nor the wave equation is solved. Thus, electromagnetic induction is not taken into account in [³]. In [⁹], air ionization was modeled by using a recursive least squares (RLS) method for corona discharges. The plasma channel is represented by a resistor and a capacitor connected in series, describing the gap space. Resistance and capacitance are functions of time obtained by performing RLS fitting from experimental data. Since the approach proposed in [⁹] is not a full-wave method, electromagnetic induction and mutual coupling cannot be studied. Finally, a two-dimensional finite-difference time-domain (FDTD) algorithm, using equations of diffusion motion of charge carrier particles was developed in [¹⁵] for modelling ionization discharges between a spherical (two-dimensionally circular) anode and a planar (two-dimensionally linear) cathode in corona electrostatic field analysis.

In this paper, a three-dimensional finite-difference time-domain (FDTD) algorithm [¹⁶] was developed for modeling air breakdown. Differently from [¹⁵], the developed FDTD formulation is based on [¹⁷], which shows that resistance of streamer filament in nitrogen-rich gases decreases exponentially over time, influenced by the Townsend coefficient and the drift velocity. The associated equation is adapted in this paper for properly representing the progress of breakdown streamer electric conductivity over time during ionization in FDTD simulations. Deionization is described in this paper by an adapted conductivity decrease profile also described in [¹⁷], which is given as an exponentially decreasing function of time (regulated by attachment frequency). Additionally, we identify and model the following ionization stages: phases 1 and 2 of prebreakdown streamer, breakdown streamer and deionization. Gap size influence has been treated by detecting the moment in which the streamer tip reaches the anode plate (initiating the breakdown streamer phase). Notice that electromagnetic coupling among the circuit parts (including the ionized channel) for any geometry is treated by solving Maxwell’s equations via the FDTD method. This way, the coupling influence of simple or geometrically complex objects placed around the ionized channel into temporal evolution of conductivity of the ionized filament can be taken into account by employing regular FDTD modelling procedures. We, thus, present simple formulas and procedures for numerically representing the ionized region by means of a time-evolving electric conductivity under concurrent influence of voltage source and electromagnetic coupling of neighboring objects. Distinctively from [³], the proposed method does not require fine discretization around the tip of the needle electrode nor around the streamer head. Thus, uniform FDTD grid is used and the ionized channel is represented by a subcellular filament. The transient FDTD results in this paper, obtained for air breakdown streamer currents, agree very well with experimental data from [²].

II. BRIEF REVIEW OF THE FDTD METHOD

A. Discretization of Maxwell’s Equations

The well-known time-domain Maxwell’s differential equations describing the propagation of electromagnetic waves in nondispersive-isotropic-lossy materials are given by

(1) $$\nabla \times \overrightarrow{E}=-\mu \frac{\partial \overrightarrow{H}}{\partial t}$$

and

(2) $$\nabla \times \overrightarrow{H}=\overrightarrow{J}+ϵ\frac{\partial \overrightarrow{E}}{\partial t},$$

where $\overrightarrow{H}$, $\overrightarrow{E}$ and $\overrightarrow{J}=\sigma \overrightarrow{E}$ are the magnetic field (A/m), electric field (V/m) and electric current density (A/m²), respectively. The parameters ϵ, µ and σ are the electric permittivity (F/m), magnetic permeability (H/m) and electric conductivity (S/m), and t is time (s). When (1) and (2) are written in three-dimensional Cartesian space, six scalar equations are obtained for the time derivatives of the scalar components of $\overrightarrow{E}$ and $\overrightarrow{H}$. For instance, for ∂E_z/∂t and ∂H_x/∂t, we have, respectively,

(3) $$\frac{\partial E_z}{\partial t}=\frac{1}{ϵ}(\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}-{\sigma }_z{E}_z)$$

and

(4) $$\frac{\partial H_x}{\partial t}=\frac{1}{\mu }(\frac{\partial E_y}{\partial z}-\frac{\partial E_z}{\partial y}).$$

In 1966, Kane Yee [¹⁸] developed a numerical method based on central finite differences [¹⁶] to solve the system of scalar differential equations from (1) and (2): the FDTD (finite-difference time-domain) method. The space is numerically represented by a set of Yee cells (Fig. 1), where the vector components of $\overrightarrow{E}$ and $\overrightarrow{H}$ are positioned in order to satisfy curl operators in (1) and (2). The edges of the Yee cell measure Δx, Δy and Δz and the computational coordinates of the discrete space are given by the integer indexes i, j and k. The physical coordinates of the cell reference point (i, j, k) in Fig. 1 are given by x = i · Δx, y = j η Δy and z = k Δz. For sake of exemplification, space coordinates of E_z, in Yee Cell of Fig. 1, are x = i η Δx, y = j η Δy and $z=(k+\frac{1}{2})\cdot \mathrm{\Delta }z$. Electric conductivity is represented by ${\sigma }_x(i+\frac{1}{2},j,k)$, ${\sigma }_{y(i,j+\frac{1}{2},k)}$ and ${\sigma }_z(i,j,k+\frac{1}{2})$, which are respectively defined at the physical coordinates of E_x, E_y and E_z. In addition, ${\sigma }_x(i+\frac{1}{2},j,k)$, ${\sigma }_{{}_y(i,j+\frac{1}{2},k)}$ and ${\sigma }_{{}_z(i,j,k+\frac{1}{2})}$ are used to calculate E_x, E_y and E_z, respectively.

Fig. 1
The three-dimensional Cartesian Yee cell.

By applying central finite differences to space and time derivatives in (3) and (4), one obtains

(5) $$\begin{array}{r}{E}_{z(i,j,k+\frac{1}{2})}^{n+1}=\left(\frac{1-{\sigma }_{z(i,j,k+\frac{1}{2})}\frac{{\mathrm{\Delta }}_t}{2ϵ}}{1+{\sigma }_{z(i,j,k+\frac{1}{2})}\frac{{\mathrm{\Delta }}_t}{2ϵ}}\right)\cdot E_{z(i,j,k+\frac{1}{2})}^n+\\ +\frac{{\mathrm{\Delta }}_t}{ϵ(1+{\sigma }_{z(i,j,k+\frac{1}{2})}\frac{{\mathrm{\Delta }}_t}{2ϵ})}\cdot [\frac{H_{y(i+\frac{1}{2},j,k+\frac{1}{2})}^{n+\frac{1}{2}}-H_{y(i-\frac{1}{2},j,k+\frac{1}{2})}^{n+\frac{1}{2}}}{{\mathrm{\Delta }}_x}-\frac{H_{x(i,j+\frac{1}{2},k+\frac{1}{2})}^{n+\frac{1}{2}}-H_{x(i,j-\frac{1}{2},k+\frac{1}{2})}^{n+\frac{1}{2}}}{{\mathrm{\Delta }}_y}]\end{array}$$

and

(6) $$\begin{array}{r}{H}_{x(i,j+\frac{1}{2},k+\frac{1}{2})}^{n+\frac{1}{2}}=H_{x(i,j+\frac{1}{2},k+\frac{1}{2})}^{n-\frac{1}{2}}+\\ +\frac{{\mathrm{\Delta }}_t}{\mu }[\frac{E_{y(i,j+\frac{1}{2},k+1)}^n-E_{y(i,j+\frac{1}{2},k)}^n}{{\mathrm{\Delta }}_z}-\frac{E_{z(i,j+1,k+\frac{1}{2})}^n-E_{z(i,j,k+\frac{1}{2})}^n}{{\mathrm{\Delta }}_y}],\end{array}$$

respectively. Because central approximations are applied to calculate time derivatives, time coordinate of all scalar components of $\overrightarrow{E}$ is t = Δt and time coordinate of all scalar components of $\overrightarrow{H}$ is $t=(n+\frac{1}{2})\mathrm{\Delta }t$, where n is the time index used to control the time loop of FDTD computer implementations.

B. Stability and Accuracy

The FDTD method has criteria to guarantee stability and accuracy of numerical solutions. Accuracy is understood as the convergence of the numerical solution to the analytic or experimental data, and stability is associated to decreasing error over time [¹⁶]. These features are guaranteed by limiting spatial and temporal steps. In order to guarantee stability, the time step Δt must satisfy the Courant condition [¹⁶], which is given by

(7) $$\mathrm{\Delta }t\le \frac{1}{v_{max}\sqrt{\frac{1}{\mathrm{\Delta }{x}^2}+\frac{1}{\mathrm{\Delta }{y}^2}+\frac{1}{\mathrm{\Delta }{z}^2}}},$$

where v_max is the maximum wave propagation velocity in the problem space. In order to restrict errors associated with dispersion and accuracy, all spatial increments, represented by Δ_x,y,z, must satisfy

(8) $${\mathrm{\Delta }}_{x,y,z}\le \frac{{\lambda }_{min}}{10},$$

where λ_min is the shortest propagating wavelength in the space domain.

The analysis region needs to be truncated to avoid errors associated with non-physical reflections of waves at the domain boundaries. For this goal, it is necessary to use absorption layers at the boundaries of the computational domain. The formulation used in this work is the Uniaxial Perfectly Matched Layer (UPML) [¹⁶] [¹⁹], designed with ten absorbing layers in this work.

III. THE PROPOSED AIR IONIZATION MODEL

A. Ionization Fundamental Equations

In [¹⁷], using experimental curves obtained for discharges in nitrogen-rich gases, it was deduced that streamer resistance between electrodes can be expressed by

(9) $$R_s=R_{so}+R_{ss}exp(-t/{\tau }_s),$$

in which t is measured from the time at which the ionization is initiated, R_so and R_ss are obtained from experiments and τ_s = 1/(αv_d), where α is the first Townsend coefficient and v_d is the drift velocity [²⁰], [²¹]. During deionization, it is shown in [²²] that

(10) $$n_e=n_{0}exp(-{\nu }_{a}t),$$

where n_e is the electron density and ν_a is the attachment frequency. As it is well-known [²²], electric conductivity σ is given in a general way by

(11) $$\sigma =e{\mu }_e{n}_e,$$

in which e is the elementary charge and µ_e is the electron mobility.

The developed FDTD formulation is based on (9), which is adapted in this work for properly representing the progress of breakdown streamer electric conductivity over time during ionization in FDTD simulations. Deionization is described by a conductivity decrease profile as function of time, adapted from (10) and (11), and by a proposed time-evolving critical electric field.

B. The Proposed Description of the Physical Phenomena

The formation of plasma in the air region between ionizing-high-voltage electrodes is governed by three fundamental stages: no ionization, ionization and deionization. Ionization is divided into stages called prebreakdown streamer (PBS), which in turn is subdivided into phases 1 and 2, and breakdown streamer (BS). Fig. 2 shows the air ionization discharge current experimentally obtained in [²]. Additionally, we have indicated all the stages of plasma formation using vertical strips along with respective ionization names or acronyms of stages.

Fig. 2
Current measured in [²] and the stages of ionization process.

No ionization is the stage in which, in a given region of space, the electric field has not yet reached the critical electric field E_c, i.e., the air has not yet been ionized. In this stage, the electric conductivity remains constant and independent on the electric field intensity since the excitation source level is not yet able to release electrons from the air atoms in the referred region. Therefore, at this stage, conductivity is simply given by

(12) $${\sigma }_{\xi }(t)={\sigma }_0,$$

where ξ can be x, y or z, and σ₀ is the initial conductivity of the air, which is typically between 10⁻¹³ and 10⁻⁹ S/m [²³].

When electric field exceeds E_c, there is an electrical rupture of the air and, as a consequence, the amount of free electrons gradually increases, initiating the ionization process of the channel. Thus, the phase 1 of the prebreakdown streamer is triggered. During this phase, considering excitation sources of tens of kilovolts and milimetric gaps separating the electrodes, the conductivity increases typically six orders of magnitude, beginning from σ₀, i.e., it can reach the order of 10⁻³ S/m.

For this first ionization phase, we adapt (9) from [¹⁷] for describing the exponential temporal evolution of electric conductivity, which, in this work, is given by

(13) $${\sigma }_{\xi }(t)={\sigma }_{\iota ,1}{e}^{(t-t_{\iota ,1})/{\tau }_{\iota ,1}}$$

where σ_ι,1 and t_ι,1 are, respectively, the conductivity and the time at which phase 1 of the prebreakdown streamer starts. The parameter τ_ι,1 is the time constant that governs the conductivity increase over time during the PBS phase 1, of which current can be seen in Fig. 2 during the approximate time range from 2.5 ns to 41.87 ns, for the experimental current obtained in [²].

As described in [¹], once a certain level of free charges is reached, and since the streamer tip has not yet reached the cathode, the rate of free charge release is reduced due to interactions among free electrons and ions. Therefore, the plasma current (and channel conductivity) increases less sharply over time, as it is shown by Fig. 2 in the time range from 41.87 to 73.83 ns. This is designated as phase 2 of the prebreakdown streamer. In this phase, conductivity is obtained from

(14) $${\sigma }_{\xi }(t)={\sigma }_{\iota ,2}{e}^{(t-t_{\iota ,2})/{\tau }_{\iota ,2}},$$

where σ_ι,2 and t_ι,2 are, respectively, the initial conductivity of phase 2 (last value of σ_ξ(t) in phase 1) and the time at which phase 2 of PBS starts. The parameter τ_ι,2 is the time constant governing σ during phase 2 of PBS.

The third stage of ionization is the breakdown streamer. This stage is initiated when the discharge channel reaches the cathode and it is characterized by a rapid increase of electric conductivity since the discharge channel electrically connects the high-voltage electrodes. In this stage, the discharges reach their current peak. The streamer conductivity evolves over time according to

(15) $${\sigma }_{\xi }(t)={\sigma }_{\iota ,3}{e}^{(t-t_{\iota ,3})/{\tau }_{\iota ,3}},$$

where σ_ι,3 is the last conductivity value reached in phase 2 of prebreakdown streamer, t_ι,3 the starting time of the streamer breakdown and τ_ι,3 is the time constant of conductivity during the breakdown streamer ionization phase.

The progressive increase of channel diameter due to temperature rise and the subsequent reduction of electronic mobility due to the large number of free electrons, as described in [²] and [²²], are represented here by a progressive increase of critical electric field during the breakdown streamer process, which follows the linear equation we propose in this work, given by

(16) $$E_c(t)=E_{c1}+\left(\frac{t-t_{\iota ,3}}{m_t}\right)(E_{c2}-E_{c1}),$$

where E_c1 is the initial critical electric field of air, E_c2 is the maximum critical electric field reached during the breakdown streamer and m_t is the parameter determining the critical electric field increasing rate over time.

The increase of channel diameter and the reduction of electronic mobility, represented by the increase of critical electric field governed by (16), favor the beginning of the last stage: deionization. The electrons gradually return to the atomic orbits and consequently the conductivity progressively is reduced towards its initial value σ₀ (conductivity of non-ionized air) as (10) and (11) suggest. The conductivity, during this stage, is calculated in this work by

(17) $${\sigma }_{\xi }(t)={[\frac{1}{{\sigma }_d}+(\frac{1}{{\sigma }_0}-\frac{1}{{\sigma }_d})(1-e^{-\left(\frac{t-t_d}{{\tau }_d}\right)})]}^{-1}.$$

In (17), σ_d is the last instantaneous conductivity obtained during the previous ionization phase and t_d is the time at which deionization begins (when |E_ξ| < E_c(t)). Finally, the parameter τ_d is a function of time, introduced in this paper to mathematically describe non-linear temporal evolution of charge re-association rates [¹], strongly influenced by µ_e over time, which is given by

(18) $${\tau }_d=({\tau }_{d,1}-{\tau }_{d,2})\{1-\frac{1}{1+e^{-c_1(t-t_d-c_2)}}\}+{\tau }_{d,2},$$

where τ_d,1 and τ_d,2 are, respectively, the maximum and minimum of τ_d. The parameter c₁ is defined to control the sigmoid slope and c₂ defines an adjusting time shift for the function.

C. The proposed algorithm

The physical processes described in the previous section are computationally defined by the proposed algorithms, which are illustrated by Figs. 3 and 4. Fig. 3 shows an overview flowchart of the FDTD algorithm for calculating electromagnetic fields in ionized air and surroundings. Modeling of multiplestage-air-ionization streamer, proposed in this work, is contained in the block “Calculation of σ”, which contains the main differences with respect to the traditional FDTD algorithm: establishing the ionizable air path in which the discharge channel is to be formed and updating the conductivity within the discharge channel(s) using (13)-(18) during the time loop. The parameters of ionization processes are specified prior to the time loop.

Fig. 3
FDTD Algorithmic flowchart: overview of the proposed numerical model for calculating electromagnetic fields in ionized air and surroundings. Electric and magnetic fields are calculated all over the numerical space.

Fig. 4
The proposed algorithm to calculate electric conductivity over time in ionized air.

In Fig. 4, the developed procedure for updating the streamer electric conductivity is presented, considering a single ionizable linear path parallel to z-axis. In this case, the calculation of σ_z is performed in a spatial loop over the z-direction along the channel line, in which streamer cells are indexed by k from the cathode to the anode. For z-aligned streamer channels, Yee indexes i and j are constant and set to i_e and j_e, respectively (see Fig. 1). The procedure is performed before the calculation of electric field all over the numerical space, as Fig. 3 indicates. Time-related indexes n and n_max in Fig. 3 are, respectively, the time index as previously defined and the number of time steps necessary for properly analyzing the physical problem. In Fig. 4, k_f is the index of the streamer cell closest to the anode. The steps related to the ionization stages, i.e, prebreakdown streamer (phases 1 and 2), breakdown streamer and deionization, are enveloped by closed dashed lines, containing steps of each stage separately, which are identified by their respective stage names.

As discussed in subsection III-B, the critical electric field E_c is considered to be a function of time t. Before starting the time loop, the state E_c = E_c1 is established, which is the initial value of critical electric field. As long as $\left|E_{z(i_e,j_e,k)}^n\right|$ at a point in the streamer is smaller than E_c, the ionization process does not occur and, therefore, ${\sigma }_z^n$ remains equal to σ₀. This procedure corresponds to the dashed box labeled as no ionization in Fig. 4.

When $\left|E_{z(i_e,j_e,k)}^n\right|$ exceeds E_c, ionization of air begins with phase 1 of the prebreakdown streamer stage. During this stage, identified in Fig. 4 as PBS - phase 1, the conductivity ${\sigma }_z^n$ increases governed by (13) while it does not exceed the threshold σ₁.

Once the conductivity exceeds the threshold σ₁ in points of streamer channel, reduction of increase rate of ${\sigma }_z^n$ over time takes place governed by (14). This slower increasing profile of conductivity corresponds to phase 2 of the prebreakdown streamer stage. The steps for calculating ${\sigma }_z^n$ during this stage is labeled as PBS - phase 2 in Fig. 4.

As soon as the conductivity at all points in the ionizable path becomes greater than the threshold σ₁, the plasma channel is considered to be formed and then the conductivity starts to increase very quickly, since the plasma channel is subjected to the high voltage from the source, as described in subsection III-B. This stage is modeled by the breakdown streamer phase, in which ${\sigma }_z^n$ follows (15). Once the breakdown streamer phase begins, the increasing of the critical electric field E_c takes place governed by (16) for the streamer cell closest to the anode. For the other streamer cells, the increasing of E_c begins at (i_e, j_e, k) once ${\sigma }_{{}_z(i_e,j_e,k)}\ge {\sigma }_{z(i_e,j_e,k+1)}$.

As ${\sigma }_z^n$ increases, |E_z| is reduced at the points in the plasma channel as a natural consequence of Maxwell’s equations. The deionization phase occurs when |E_z| is reduced to a level where $\left|E_z^n\right|\le E_c$. At this stage, conductivity decreases governed by (17). Additionally, the deionization parameter τ_d evolves driven by (18), which is a decreasing sigmoid used in order to properly represent interactions and re-associations of charges in the channel [¹], which substantially affect the current decay profile during deionization. Finally, notice that the developed methodology for streamer conductivity calculation can be straightforwardly adapted to model plasma channels formed parallelly to x and y axes.

IV. RESULTS

Experiments carried out by Eichwald et al. in [²] have been numerically reproduced in order to validate the multiple-stage air ionization FDTD model we developed. The experimental setup used in [²] and [⁴] for producing plasma channels is based on a high-voltage air gap between a metallic needle and a metallic plate acting as anode and cathode, respectively. The plasma channel is formed in a 7 mm air gap. As broadly detailed in [⁴], the setup also has a set of metallic conductors, an 1 kΩ resistor and a high voltage source v(t), such as illustrated by Fig. 5. The circuit and electrodes are immersed in air. Fig. 6 shows a three-dimensional representation of FDTD boundary conditions produced using the in-house-made FDTD simulator SAGS (Synthesis and Analysis of Grounding Systems) [²⁴], which we adapted for performing the streamer simulations. The electrodes and wires are modelled as perfect electrical conductors (PECs), on which tangential electric field components are forced to zero during all time. The PEC metal plate is positioned in such a way that the vertical electrode above the plate (the needle anode) points to the center of the flat cathode structure, as shown in Fig. 6.

Fig. 5
The experimental setup described in [⁴] for plasma generation.

Fig. 6
SAGS three-dimensional representation of FDTD boundary conditions for the needle-plate setup.

For all simulated cases, the computational domain contains 144 × 140 × 79 cells, with Δx = Δy = Δz = 0.5 mm and Δt = 0.953 ps, respecting the FDTD Courant and accuracy criteria [¹⁶] given by (7) and (8), respectively. By using a single 64-bit eight-core 3.7 GHz processor running the Linux operating system, our parallel-FDTD C code simulator requires approximately two hours for running n_max = 160,000 FDTD time steps by employing 16 threads. The simulations presented here needed 4112 MB of RAM. We used the thin-wire formulation described in [²⁵] for modeling the thin cylindrical metallic wire conductors, whose radii are set to 0.02 mm. The resistor with resistance R is included in FDTD mesh by simply setting σ_z = Δ_z/(RΔ_xΔ_y) in (5) solely at the Yee cell in which the resistor is contained [¹⁶]. Voltage source v(t) is modelled by setting $E_z^n(i_s,j_s,k_s)=v\left(n{\mathrm{\Delta }}_t\right)/{\mathrm{\Delta }}_z$, i.e., it acts as a hard source [¹⁶] at the source cell (i_s, j_s, k_s).

The 7.2 kV and 8.2 kV time-dependent high voltage source waveforms from the experimental work [²] are reproduced here using Padé approximants [²⁶] and the Least Squares Method [²⁷], which are then used in our FDTD simulations. Fig. 7 shows a comparison between the 7.2 kV excitation voltage from [²] and the obtained Padé approximant. The 8.2 kV voltage function used to excite the second FDTD validation case is also presented in Fig. 7. All the parameters used in this work for calculating σ in streamer channels during ionization and deionization processes, with the proposed algorithm illustrated in Fig. 4, are given in Table I for the 7.2 kV and 8.2 kV excitation sources.

Fig. 7
Voltage source signals: the experimental 7.2 kV signal from [²] and Padé approximants for 7.2 kV and 8.2 kV signals. The Padé approximants are used in FDTD validation cases as excitation functions in this work.

Table I
STREAMER PARAMETERS USED IN FDTD SIMULATIONS OF THIS WORK

		Values
Parameter	Description	Excitation: 7.2 kV	Excitation: 8.2 kV
σ 0	Initial electric conductivity of air	1.0 nS/m	1.0 nS/m
σ 1	Electric conductivity threshold	0.5 mS/m	0.06 mS/m
Ec 1	Minimum critical electric field	0.1 MV/m	0.1 MV/m
E c2	Maximum critical electric field	3.0 MV/m	3.0 MV/m
mt	Parameter defining the temporal change rate of critical electric field	3.0 ns	2.0 ns
τ ι,1	Time constant of PBS - phase 1	3.0 ns	4.0 ns
τ ι,2	Time constant of PBS - phase 2	8.0 ns	4.0 ns
τ ι,3	Time constant of BS	2.0 ns	1.0 ns
τ d,1	Maximum time constant of deionization phase	1.0 s	4.0 s
τ d,2	Minimum time constant of deionization phase	0.2 s	0.6 s
c 1	Constant defining the temporal change rate of τd	60×106	100×106
c 2	Constant defining the temporal shift of τd	30 ns	20 ns

Fig. 8 illustrates spatial distributions of electric conductivity σ in the ionized channel over time, for the 7.2 kV excitation voltage. As one can see, the streamer tip moves from the anode to the cathode at a speed of approximately 2 × 10⁵ m/s, which can be estimated by observing the conductivity distributions calculated at 60 and 70 ns. This velocity obtained in our FDTD simulation is compatible to the speed range given in [¹] (between 2 × 10⁵ m/s and 2 × 10⁶ m/s), demonstrating the physical consistency of the proposed model. As the streamer tip reaches the anode plate, at approximately 75 ns, the breakdown streamer phase begins and it is visible in Fig. 8 that conductivity increases rapidly over the channel extent, which displays more evenly spatially distributed high levels of conductivity as soon as at 80 ns.

Fig. 8
Electric conductivity (S/m) in the plasma channel, for the 7.2 kV excitation, obtained using the proposed FDTD ionization algorithm.

The deionization phase begins and conductivity levels decrease as observed in Fig. 8 from 100 ns to 200 ns. Illustratively, Fig. 9 shows the z-component of electric field on x-z plane for several moments of the ionization process. Observe that radiated wave electromagnetically couples all the parts of the circuit, affecting its transient responses.

Fig. 9
Spatial distributions of |E_z | (dB) on x-z plane at: (a) 0 ns, (b) 1 ns, (c) 3 ns, (d) 20 ns, (e) 70 ns and (f) 110 ns.

Considering the problem excited with the 7.2 kV signal, Fig. 10 shows electric conductivity at three points in the streamer channel, as functions of time. The observation points are 0.25 mm below the needle tip, at the center of the channel and 0.25 mm above the plate electrode. The peak of σ is approximately 0.046 S/m and it is observed in the cell adjacent to the tip of the anode needle. In all cells along the plasma channel, σ increases progressively over time governed by (12)-(18) and by local electric field calculated using Maxwell’s equations, as described by Fig. 4. The obtained peaks of conductivity σ increase from the cell closest to the metal plate to the Yee cell closest to the needle. Furthermore, the time at which peaks occur increases from the needle to the plate due to the expected delay caused by the electromagnetic wave propagation finite speed, from anode to cathode, as governed by Maxwell’s equations.

Fig. 10
Electric conductivity at three points of the plasma channel for the 7.2 kV excitation source.

During the beginning of the ionization process (up to approximately 40 ns), the conductivity value is negligible (Fig. 10). Therefore, the current in the initial moments is small (practically zero). Then, σ reaches higher values due to (13), implying the exponential increasing of the current. In a second stage, the conductivity rises moderately according to (14). At t ≈ 75 ns, the breakdown streamer phase begins and, as a result, conductivity increases considerably in a short period of time (Fig. 10), governed by (15). In this process, the current reaches its peak of 14 mA. After breakdown streamer phase (see Fig. 8), streamer channel heats up, it goes through dilation and electric field diminishes over the plasma channel due to a massive increase of electron density n_e. Thus, after some time, $\left|E_z^n\right|$ becomes less than E_c and therefore the deionization phase begins (re-associations of electrons and ions), in which σ reduction is governed by (17) and (18). Observe the inset of Fig. 10, in which ionization phases are identified in the conductivity curve, over time, obtained right below the needle tip.

In Fig. 11, we see that the current calculated using the proposed FDTD method presents improved conformity to experimental results with respect to current calculated in [²]. Electric current is calculated 0.25 mm below the tip of the needle electrode. The quality of numerical results in Fig. 11 is calculated using the Mean Square Error (MSE) [²⁸], which is applied over time. Numerical results obtained in this work with FDTD and those obtained numerically in [²] are compared to experimental results from [²]. Mean Square Error is defined as

(19) $$MSE=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}{[I_{exp}(n)-I_{\text{num}}(n)]}^2},$$

where n is the time index of signal sample, N is the total number of samples, I_exp(n) is the n-th sample of the experimental current obtained in [²] and I_num (n) is the n-th sample of the current obtained numerically. The calculated value of MSE for the current obtained with the proposed FDTD model is 0.177 while MSE computed for the numerical results in [²] is 6.643, which quantitatively confirms what is seen in Fig. 11: the proposed numerical model has substantially higher conformity to experimental current curve.

Fig. 11
Plasma channel currents obtained by applying the 7.2 kV excitation: experiment and calculations of [²] and the current calculated using the proposed FDTD method.

In Fig. 12, the current signals for the 8.2 kV excitation are shown. Similarly to the previous case, the obtained FDTD result shows excellent agreement with the signal measured in [²]. Fig. 13 shows conductivity obtained at the same observation points of the previous case. As expected, conductivity reaches higher levels due to the higher applied voltage. As one can see by comparing Figs. 10 and 13, σ increases non-linearly with voltage, as expected. As voltage source goes from 7.2 kV and 8.2 kV, conductivity peak increases from approximately 0.046 S/m to 0.105 S/m, respectively.

Fig. 12
Comparison of plasma channel currents obtained by applying the 8.2 kV excitation: experimental curve from [²] and current calculated in this paper using the proposed FDTD method.

Fig. 13
Electric conductivity at three points of the plasma channel for the 8.2 kV excitation source.

Similar behavior is observed for time constants, which change non-linearly with applied source voltage. As observed in Table I, when voltage source is increased from 7.2 kV to 8.2 kV, the time constant of PBS phase 2 decreases non-linearly from 8.0 ns to 4.0 ns and the time constant of breakdown streamer is reduced from 2.0 ns to 1.0 ns, as one would expect, since as ionizing voltage is increased, more electrons are extracted from atoms per second [²²]. In other words, avalanche tends to evolve faster. This also means that exciting the gap using the 8.2 kV signal tends to release more electrons than using the 7.2 kV source, as the higher plasma conductivity indicates (see Figs. 10 and 13). Finally, on the other hand, the deionization constants τ_d,1 and τ_d,2 increase with source voltage level because it takes more time to re-associate the larger number of free charges available when the 8.2 kV signal is used.

V. FINAL REMARKS

In this work, an electromagnetic air ionization formulation is developed, based on the finite- difference time-domain method (used for calculation of electromagnetic fields) and on the streamer particle physics, defining instantaneous space-localized electric conductivity. In the proposed model, electric conductivity σ evolves over time and it is calculated for each proposed ionization stage: phases 1 and 2 of prebreakdown streamer, breakdown streamer and deionization. The calculation of σ, which is used in time domain Maxwell-Ampère FDTD equation for updating electric field over time, is performed using functions of time adapted from literature, for each ionization stage. Specifically, new functions describing conductivity reduction over time during air streamer deionization and time-evolving critical electric field are proposed. Therefore, we have presented general formulas and procedures that can be used to numerically represent the ionized filament. Important gap size effect has been treated by detecting the moment in which the streamer tip reaches the anode plate and, thus, the breakdown streamer phase begins. These are the main contributions of this paper. Electromagnetic coupling among the circuit parts (including the ionized channel) for any geometry is treated by solving Maxwell’s equations using the FDTD method itself. Thus, changing geometric features of the problem should not be a formation limitation.

The proposed streamer modeling was validated by performing comparisons among our FDTD results and experimental data from literature for an anode needle over a metallic sheet setup. Between the vertical needle and the horizontal cathode sheet, voltages of 7.2 kV and 8.2 kV were stablished for each analyzed case, where it is possible to observe excellent agreements between numerical and experimental data for each applied excitation voltage.

In order to have a quantitative assessment of agreement between the literature experimental current and the FDTD-simulated current signal, the well-known Mean Square Error (MSE) indicator was used for comparing the signals over time. A MSE of 0.177 was obtained between the experimental current and the FDTD current signal, where both of which were obtained when the excitation source is 7.2 kV. For the same excitation case, when both experimental and numerical currents obtained in the reference paper are compared, the calculated MSE is 6.643, i.e., it is 37.53 times the MSE obtained using the FDTD model proposed in this work, demonstrating a substantial precision improvement over the simulated current obtained previously in literature.

For future works, the streamer parameters can be calculated as functions of the instantaneous local field information for updating σ over time for each cell of the ionized channel. Furthermore, it is suggested possible extension and application of the developed method to model (partial) discharges in high voltage equipment, such as those in insulating parts of hydroelectric generators and transformers.

REFERENCES

• ^[1]
  R. S. Sigmond, “The residual streamer channel: Return strokes and secondary streamers,” Journal of Applied Physics, vol. 56, no. 5, pp. 1355–1370, 1984.
• ^[2]
  O. Eichwald, O. Ducasse, D. Dubois, A. Abahazem, N. Merbahi, M. Benhenni, and M. Yousfi, “Experimental analysis and modelling of positive streamer in air: towards an estimation of O and N radical production,” Journal of Physics D: Applied Physics, vol. 41, no. 23, p. 234002, nov 2008.
• ^[3]
  M. Talaat, A. El-Zein, and A. Samir, “Numerical and simulation model of the streamer inception at atmospheric pressure under the effect of a non-uniform electric field,” Vacuum, vol. 160, pp. 197–204, 2019.
• ^[4]
  E. Moreau, P. Audier, and N. Benard, “Ionic wind produced by positive and negative corona discharges in air,” Journal of Electrostatics, vol. 93, pp. 85 – 96, 2018.
• ^[5]
  E. Kuffel, W.S. Zaengl, J. Kuffel, High Voltage Engineering Second edition, British Library Cataloguing in Publication Data, New Delphi, 2000.
• ^[6]
  R Arora, W Mosch, High Voltage Insulation Engineering 1 ed., New Age International Publishers, New Delphi, 2004.
• ^[7]
  Y. P. Raizer, Gas Discharge Physics Springer-Verlag, Berlim, 1991.
• ^[8]
  A. N. Esfahani, S. Shahabi, G. Stone, and B. Kordi, “Investigation of corona partial discharge characteristics under variable frequency and air pressure,” in 2018 IEEE Electrical Insulation Conference (EIC), pp. 31–34, 2018.
• ^[9]
  D. Raouti, S. Flazi, and D. Benyoucef, “Modeling and Identification of Electrical Parameters of Positive DC Point-to-Plane Corona Discharge in Dry Air Using RLS Method,” IEEE Transactions on Plasma Science, vol. 44, no. 7, pp. 1144–1149, 2016.
• ^[10]
  R. M. S. Oliveira, J. F. M. Modesto, V. Dmitriev, F. S. Brasil, and P. R. M. Vilhena, “Spectral method for localization of multiple partial discharges in dielectric insulation of hydro-generator coils: Simulation and experimental results,” Journal of Microwaves, Optoelectronics and Electromagnetic Applications, vol. 15, pp. 170 – 190, 09 2016.
• ^[11]
  Institute of Electrical and Electronics Engineers - IEEE, “IEEE Guide to the Measurement of Partial Discharges in Rotating Machinery,” IEEE Std 1434-2000, pp. 1–64, 2000.
• ^[12]
  G. C. Stone, I. Culbert, E. A. Boulter, and H. Dhirani, Electrical Insulation for Rotating Machines, 2nd ed. Wiley, 2004.
• ^[13]
  C. Hudon and M. Belec, “Partial discharges signal interpretation for generator diagnostics,” pp. 297–319, 11 2005.
• ^[14]
  R. C. F. Araujo, R. M. S. de Oliveira, F. S. Brasil, and F. J. B. Barros, “Novel Features and PRPD Image Denoising Method for Improved Single-Source Partial Discharges Classification in On-Line Hydro-Generators,” Energies, vol. 14, no. 11, p. 3267, Jun 2021.
• ^[15]
  B. B. Alagoz, H. Z. Alisoy, S. Alagoz, and F. Hansu, “A space charge motion simulation with FDTD method and application in negative corona electrostatic field analysis,” Applied Mathematics and Computation, vol. 218, no. 17, pp. 9007–9017, 2012.
• ^[16]
  A. Taflove and S. C. Hagness, Computational Electrodynamics, The Finite-Difference Time-Domain Method, 3rd ed. Artech House, 2005.
• ^[17]
  M. J. Kushner, A. L. Pindroh, C. H. Fisher, T. A. Znotins, and J. J. Ewing, “Multidimensional modeling of transverse avalanche laser discharges: Applications to the HgBr laser,” Journal of Applied Physics, vol. 57, no. 7, pp. 2406–2423, 1984.
• ^[18]
  K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” Antennas and Propagation, IEEE Transactions on, vol. 14, pp. 302–307, 1966.
• ^[19]
  S. D. Gedney, “An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Transactions on Antennas and Propagation, vol. 44, no. 12, pp. 1630–1639, 1996.
• ^[20]
  P. Persephonis, V. Giannetas, A. Ioannou, J. Parthenios, and C. Georgiades, “The time dependent resistance and inductance of the electric discharges in pulsed gas lasers,” IEEE Journal of Quantum Electronics, vol. 31, no. 10, pp. 1779–1784, 1995.
• ^[21]
  S. Panahibakhsh and A. Hariri, “Calculation of electrical conductivity of fast discharges in nitrogen gas using the performance of a transversely excited N2 laser,” Physics of Plasmas, vol. 24, no. 9, p. 093112, 2017.
• ^[22]
  E. M. Bazelyan and Y. P. Raizer, Spark Discharge, 1st ed. CRC Press, 1998.
• ^[23]
  E. Seran, M. Godefroy, E. Pili, N. Michielsen, and S. Bondiguel, “What we can learn from measurements of air electric conductivity in 222Rn-rich atmosphere,” Earth and Space Science, vol. 4, no. 2, pp. 91–106, 2017.
• ^[24]
  R. M. S. de Oliveira and C. Sobrinho, “Computational environment for simulating lightning strokes in a power substation by finite-difference time-domain method,” IEEE Transactions on Electromagnetic Compatibility, vol. 51, no. 4, pp. 995–1000, 2009.
• ^[25]
  Y. Taniguchi, Y. Baba, N. Nagaoka, and A. Ametani, “An Improved Thin Wire Representation for FDTD Computations,” IEEE Transactions on Antennas and Propagation, vol. 56, no. 10, pp. 3248–3252, 2008.
• ^[26]
  C. Brezinski, “Extrapolation algorithms and Padé approximations: a historical survey,” Applied Numerical Mathematics, vol. 20, no. 3, pp. 299–318, 1996.
• ^[27]
  H. Oraizi, “Application of the method of least squares to electromagnetic engineering problems,” IEEE Antennas and Propagation Magazine, vol. 48, no. 1, pp. 50–74, 2006.
• ^[28]
  R. J. Hyndman and A. B. Koehler, “Another look at measures of forecast accuracy,” International Journal of Forecasting, vol. 22, no. 4, pp. 679 – 688, 2006.

Publication Dates

• Publication in this collection
  23 Sept 2022
• Date of issue
  Sept 2022

History

• Received
  17 Dec 2021
• Reviewed
  11 Jan 2022
• Accepted
  28 June 2022