Mitigation of Insulator Lightning-Induced Voltages by Installing Parallel Low-Voltage Surge Arresters

Abstract

In this paper, we propose a mitigation method for reducing lightning-induced insulator voltages based on the installation of low-voltage surge arresters aligned parallelly to the insulator. The three-dimensional finite-difference time-domain (FDTD) method is applied to numerically model a real surge arrester residual voltage evaluation system. The application of a transient current pulse, typical of lightning discharges, is considered in our numerical model. We considered cases with one or two surge arresters installed per phase, in three different geometric and parametric configurations for installing distribution surge arresters. In addition to the Kirchhoff current division, which reduces both the absorbed energy and the thermal stress, the results associated with the installation of two surge arresters parallelly aligned to the insulator show that the interaction of magnetic fields generated by the surge arresters’ currents can produce an additional strong reduction in lightning-induced voltage over the insulator, as presented in this paper. Conditions for maximum voltage reduction are also identified. A brief cost-effectiveness analysis is also provided.

1. Introduction

Lightning discharges produce strong transient electromagnetic fields associated with fast high-current flows (tens of kiloamperes during a few microseconds, in general). In Brazil alone, the occurrence of 50 to 100 million lightning events is estimated per year [1,2], which is the main cause of service interruptions linked to electric power distribution overhead system issues [3,4,5]. In some Brazilian regions, for instance, 30% to 40% of unavailability events of such distribution systems are caused by lightning occurrences [6,7] and approximately 47% of power transformers are damaged by atmospheric discharges [8], leading to high financial losses to power companies and various damage to society due to electric damage to equipment and interruption in electric power supply.

The proposal of efficient remedial measures applied to power distribution networks is not a simple issue due to the associated complex physics. Therefore, the use of adequate numerical methods and computational tools for modelling the physics of such networks is necessary. The Electromagnetic Transient Program (EMTP) [9], which is based on circuit and transmission line theories for calculating overvoltages and overcurrents in electric power distribution networks, has been frequently employed [10,11,12,13]. However, over the past few years, the use of full-wave numerical approaches for calculating and analyzing transient electromagnetic fields over distribution networks has produced physically complete solutions and, consequently, comprehensive conclusions and new possibilities to solve engineering problems. Among the methods for numerically solving Maxwell’s equations, the finite-difference time-domain (FDTD) method [14,15] has been widely adopted for analyzing transient fields over distribution networks, since it enables modelling complex three-dimensional structures, such as grounding structures [16,17,18,19], power substations [20] and distribution networks [21,22], among others. It naturally considers wave phenomena, such as electromagnetic coupling and induction, reflections and diffractions, allowing one, for example, to perform analyses involving lightning-induced voltages in geometrically and parametrically complex distribution networks [23,24,25,26].

The literature reports the use of two zinc oxide (ZnO) surge arresters (SAs) per phase and set up parallelly to each other as one of the corrective measures for mitigating the undesirable effects of electromagnetic transients in electric power distribution networks. Sugimoto et al. [27] conducted tests in an experimental electric power distribution network by installing two SAs per phase. He showed the effectiveness of that configuration regarding energy absorption by the arresters since they are of the same voltage class and similar voltage-current characteristics. He et al. [28] performed laboratory experiments, statistic and numerical analysis in order to examine the nonuniformity among electrical parameters of ZnO varistors (commercial low-voltage) as a useful tool to select and coordinate the ZnO varistors in parallel operation and improve the protection capability of a surge-protection device (SPD).

Haddat et al. [29] conducted laboratory tests for investigating the disparity in currents of the same voltage class SAs and coincident characteristics, which were set up parallelly, regarding three aspects: (1) the use of ballast resistors, minimizing unbalances and improving the energy absorption rates, (2) study of influence of the arrangement on the current parity and (3) verification of the temperature influence on the residual voltage of the SAs. Tuczec et al. [30] also evaluated the connection of SAs in parallel towards the implications on the surge current division in laboratory tests by analyzing the safety of replacing defective SAs of an arrester bank developed for protecting series compensation capacitors. The experimental results achieved by Sabiha et al. in [31] and their numerical analysis with EMTP indicate that the installation of a spark gap in parallel with the surge arrester reduced the thermal stress on the device and, hence, the risks of failures in SAs. Leal and Oliveira et al. [1] computed electric field distribution in a small riverboat due to lightning strikes using the FDTD method and showed that the electric field produced inside the boat is large enough to cause air breakdown. For this reason, they proposed a cost-effective method (using more than one rod, arranged vertically, parallel to each other) to improve lightning safety for people that use small boats in the Amazon region.

In this work, a theoretical review and FDTD-based numerical analyses are conducted in order to extend the comprehension of SA operation when set up parallel to the insulator, with three geometric and parametric configurations for installing SAs, operating under fast electromagnetic transient lightning fields. For the first time, the three-dimensional analysis region was inspired in a real-world residual voltage test environment for surge arresters. Results showed that the main effect of the magnetic field interactions generated by parts of lightning current flowing through SAs is the significant reduction in the induced electric field over the insulator, being highlighted as a scientific contribution in this paper. Conditions when the insulator is subjected to maximum voltage reduction are also identified and described and a brief cost-effectiveness analysis is provided.

The paper is organized as follows: Section 2 provides a theoretical review regarding quasi-static magnetic fields, on the full-wave FDTD method, and on the typical lightning current mathematical representation; in Section 3, the main research proposal of this paper (installation of SAs parallel to insulators) is presented and technically discussed; details regarding the conceived FDTD model are provided in Section 4; the obtained numerical results are shown and discussed in Section 5; and, finally, conclusions are drawn in Section 6.

2. Theoretical Review

2.1. Interaction of Quasi-Magnetostatic Fields

To develop an effective method for reducing induced voltages over the insulator, Ampère’s law must be revisited. Maxwell’s equation for quasi-magnetostatic fields [32] can be approximated to

$$\nabla \times \overline{H}\approx \overline{J}.$$

(1)

Considering current density $\overline{J}$ in the conductor of radius $a$ as in Figure 1 implies the circulation of the quasi-magnetostatic field $\overline{H}$, so we may integrate (1) on a surface S limited by the amperian contour ℓ, providing

$${{\displaystyle \iint }}_S^{}\nabla \times \overline{H}·d\overline{S}\approx {{\displaystyle \iint }}_S^{}\overline{J}·d\overline{S}.$$

(2)

We may assume that

$$\overline{J}=\frac{I}{\pi a^2}{\overline{a}}_z\mathrm{and}d\overline{S}=\rho d\varphi d\rho {\overline{a}}_z.$$

(3)

Since ρ, $\varphi$ and z are the coordinates of the cylindrical system and I is the current flowing through the conductor, by applying Stokes’ theorem to the left-hand side of (2), we may write

$$I={{\displaystyle \iint }}_S^{}\overline{J}·d\overline{S}\approx {{\displaystyle \oint }}_{\ell }^{}\overline{H}·d\overline{l}.$$

(4)

The right-hand side of (4) can be calculated by

$${{\displaystyle \oint }}_{\ell }^{}\overline{H}·d\overline{l}=H_{\varphi }2\pi \rho .$$

(5)

Since $\overline{H}=H_{\varphi }{\overline{a}}_{\varphi }$, where $H_{\varphi }$ is the scalar $\varphi$ component of the magnetic field $\overline{H}$, using (4) and (5), we see that, around the cylindrical conductor at distance $\rho >a$ from its center, one has

$$\overline{H}\approx \frac{I}{2\pi \rho }{\overline{a}}_{\varphi }.$$

(6)

For illustrative purposes, we numerically solved (6) for two cases: (1) the first problem is illustrated by Figure 2a, which contains a section of a cylindrical conductor with radius measuring $0.02\mathrm{m}$ (transverse section $S_1$) in a $0.8\mathrm{m}\times 0.8\mathrm{m}$ vacuum analysis region. We established a z-component current density $J_{Z1}=1000\mathrm{A}/\mathrm{m}²$; (2) Figure 2b illustrates the second problem, which consists of two parallel conductors immersed in vacuum, of which transverse sections $S_1$ and $S_2$ are represented (transverse sections $S_1$ and $S_2$ have identical areas and $J_{Z2}=J_{Z1}=1000\mathrm{A}/\mathrm{m}²$).

Figure 3 shows solutions for both cases. Magnetic field $\overline{H}$ in Figure 3a was obtained as expected: it obeys the right-hand rule and the field magnitude decreases with ρ. In Figure 3b, because the vector current densities $J_{Z1}{\widehat{a}}_z$ and $J_{Z2}{\widehat{a}}_z$ are identical, $\overline{H}$ intensity is reduced in the region between the parallel conductors, as predicted by the superposition principle [32] and by directions of magnetic fields produced by the currents, as predicted by (1) and (6).

Let us consider the case in which the insulator is parallelly placed between the z-aligned conductors, in the region of interest in Figure 3b. Voltage between its terminals is produced by the associated electric field $\overline{E}$, which is governed by

$$\frac{\partial \overline{E}}{\partial t}=\frac{1}{\epsilon }\nabla \times \overline{H}.$$

(7)

where $\epsilon$ is the medium electric permittivity and t is time. Therefore, if circulation of quasi-magnetostatic $\overline{H}$ is reduced between conductors, temporal derivative of the electric field $\overline{E}$ is also reduced, limiting its evolution over time. Consequently, since voltage $V_{AB}$ induced between insulator terminals A and B is related to the electric field line integral

$$V_{AB}={{\displaystyle \int }}_A^B\overline{E}·d\overline{l}={{\displaystyle \int }}_A^B{E}_{z}dz,$$

(8)

the induced voltage is also reduced across the insulator. Notice that (8) is calculated along the z-axis, i.e., $d\overline{l}=dz {\widehat{a}}_z$.

Although (1) can be used as an approximation, full-wave electromagnetic methods are the most rigorous approaches for modelling and analyzing transient fields associated with overhead lines [21,22] and grounding systems [16,17,18,19] since Maxwell’s equations are directly solved numerically. Among the most adopted methodologies, we have the finite element method (FEM) [33], the method of moments (MoM) [34] and the finite-difference time-domain method (FDTD) [14,15]. The methodology adopted in this work is the FDTD approach.

2.2. The FDTD Method and Computational Modelling of Residual Voltage Test Environment

2.2.1. Review of the FDTD Method

FDTD is a technique developed by Kane Yee in 1966 [14] for numerically solving Maxwell’s equations in the time domain. It is based on the following two main aspects: (1) space discretization using cells with discrete spatial distribution of electric and magnetic field components, satisfying Faraday’s and Ampère’s laws in differential and integral forms, and (2) approximation of partial derivatives with respect to time and space in Maxwell’s equations by using second-order central finite differences. This way, approximate explicit equations are obtained for the scalar components of $\overline{E}$ and $\overline{H}$ [15].

A computational mesh based on Yee cell (Figure 4) is defined for representing space in the analysis region. Considering a general three-dimensional region, the dimensions of cell edges are ${\mathsf{\Delta }}_x$, ${\mathsf{\Delta }}_y$ and ${\mathsf{\Delta }}_z,$ as shown in Figure 4. The cell spatial positions in the computational domain are defined in terms of array integer indexes $i$, $j$ and $k$ defined for each cell. The physical coordinates of the corner of the cell ($i$,$j$,$k$) are given by $x=i{\mathsf{\Delta }}_x$, $y=j{\mathsf{\Delta }}_y$ and $z=k{\mathsf{\Delta }}_z$.

Temporal coordinates are provided by index $n$ and time increment ${\mathsf{\Delta }}_t$, such that the physical instants at which electric and magnetic fields are calculated are, respectively, given by $t=n{\mathsf{\Delta }}_t$ and $t=(n+0.5){\mathsf{\Delta }}_t$ for a given value of $n$ [15]. Thus, a scalar function $f$, depending on $x$, $y$, $z$ and $t$, is denoted, in its approximate discrete form, by

$$f\left(x,y,z,t\right)\approx f_d^n\left(i,j,k\right).$$

(9)

Maxwell’s equations, for Faraday’s and Ampère’s laws, are, respectively, given by

$$\nabla \times \overline{E}=-{\overline{M}}_i-\frac{\partial \overline{B}}{\partial t},$$

(10)

and

$$\nabla \times \overline{H}={\overline{J}}_i+{\overline{J}}_c+\frac{\partial \overline{D}}{\partial t}.$$

(11)

Since no magnetic current density ${\overline{M}}_i$ is present in our analysis and the media parameters do not depend on frequency (i.e., $\overline{B}=\mu \overline{H}$, $\overline{D}=\epsilon \overline{E}$ and ${\overline{J}}_c=\sigma \overline{E}),$ Equation (10) for the three-dimensional Cartesian space provides the scalar equations

$$\frac{\partial H_x}{\partial t}=\frac{1}{\mu }\left(\frac{\partial E_y}{\partial z}-\frac{\partial E_z}{\partial y}\right),$$

(12a)

$$\frac{\partial H_y}{\partial t}=\frac{1}{\mu }\left(\frac{\partial E_z}{\partial x}-\frac{\partial E_x}{\partial z}\right),$$

(12b)

and

$$\frac{\partial H_z}{\partial t}=\frac{1}{\mu }\left(\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x}\right).$$

(12c)

When electric current density is not imposed, ${\overline{J}}_i=0$ and (11) produces

$$\frac{\partial E_x}{\partial t}=\frac{1}{\epsilon }\left(\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}-\sigma E_x\right),$$

(13a)

$$\frac{\partial E_y}{\partial t}=\frac{1}{\epsilon }\left(\frac{\partial H_x}{\partial z}-\frac{\partial H_z}{\partial x}-\sigma E_y\right),$$

(13b)

and

$$\frac{\partial E_z}{\partial t}=\frac{1}{\epsilon }\left(\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}-\sigma E_z\right).$$

(13c)

The concept of central finite differences can be applied, for instance, to (12a), (12b) and (13c), respectively, producing the FDTD field equations

$$H_{x\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)}^{n+\frac{1}{2}}=H_{x\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)}^{n-\frac{1}{2}}+\frac{{\mathsf{\Delta }}_t}{\mu }\left[\frac{E_{y\left(i,j+\frac{1}{2},k+1\right)}^n-E_{y\left(i,j+\frac{1}{2},k\right)}^n}{{\mathsf{\Delta }}_z}-\frac{E_{z\left(i,j+1,k+\frac{1}{2}\right)}^n-E_{z\left(i,j,k+\frac{1}{2}\right)}^n}{{\mathsf{\Delta }}_y}\right],$$

(14)

$$H_{y\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)}^{n+\frac{1}{2}}=H_{y\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)}^{n-\frac{1}{2}}+\frac{{\mathsf{\Delta }}_t}{\mu }\left[\frac{E_{z\left(i+1,j,k+\frac{1}{2}\right)}^n-E_{z\left(i,j,k+\frac{1}{2}\right)}^n}{{\mathsf{\Delta }}_x}-\frac{E_{x\left(i+\frac{1}{2},j,k+1\right)}^n-E_{x\left(i+\frac{1}{2},j,k\right)}^n}{{\mathsf{\Delta }}_z}\right],$$

(15)

and

$$E_{z\left(i,j,k+\frac{1}{2}\right)}^{n+1}=E_{z\left(i,j,k+\frac{1}{2}\right)}^n\left(\frac{\frac{1-\sigma {\Delta }_t}{2\epsilon }}{\frac{1+\sigma {\Delta }_t}{2\epsilon }}\right)+\frac{{\Delta }_t}{\epsilon \left(1+\frac{\sigma {\Delta }_t}{2\epsilon }\right)}\left[\frac{H_{y\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)}^{n+\frac{1}{2}}-H_{y\left(i-\frac{1}{2},j,k+\frac{1}{2}\right)}^{n+\frac{1}{2}}}{{\mathsf{\Delta }}_x}-\frac{H_{x\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)}^{n+\frac{1}{2}}-H_{x\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)}^{n+\frac{1}{2}}}{{\Delta }_y}\right].$$

(16)

2.2.2. Accuracy and Stability

Since numerical methods are defined from approximations, numerical errors are present, which are cumulative and must be minimized so that the iterative processes converge to an exact solution. In FDTD, errors are dependent on spatial increments ${\Delta }_x$, ${\Delta }_y$ and ${\Delta }_z$ and also on time step ${\Delta }_t$. A numerical phenomenon produced by errors from FDTD discretization schema is numerical dispersion (phase velocities of wave depending on frequency and discretization parameters) [20]. The space discretization criterion

$${\Delta }_{x,y,z}\le \frac{{\lambda }_{min}}{10},$$

(17)

is adopted for minimizing numerical dispersion [15] (minimum wavelength must be represented by at least ten cells). On the other hand, for guaranteeing stability for the FDTD temporal marching, ${\Delta }_t$ is limited by the Courant’s condition [15], given by

$${\Delta }_t\le \frac{1}{v_{max}\sqrt{{\left(\frac{1}{\Delta x}\right)}^2+{\left(\frac{1}{\Delta y}\right)}^2+{\left(\frac{1}{\Delta z}\right)}^2}} .$$

(18)

If Yee cells are cubic with edges measuring Δ, (18) is reduced to

$${\Delta }_t\le \frac{\Delta }{v_{max}\sqrt{3}} .$$

(19)

In this work, cubic Yee cells are used, with edges measuring Δ = 0.1 m.

2.2.3. The Excitation Source: Lightning Current

The lightning current $I_s\left(t\right)$ exciting our numerical model is described by (20) [16], given in amperes by

$$I_s\left(t\right)=\{\begin{array}{c}\left(\frac{I_{\max }}{A_0}\right)\left(e^{a_{1}n\mathsf{\Delta }t}-e^{a_{2}n\mathsf{\Delta }t}\right){\sin }^2\left({\omega }_{0}n\mathsf{\Delta }t\right), \mathrm{if} n\mathsf{\Delta }t\le 1.5T_f\\ \left(\frac{I_{\max }}{A_0}\right)\left(e^{a_{1}n\mathsf{\Delta }t}-e^{a_{2}n\mathsf{\Delta }t}\right), \mathrm{if} n\mathsf{\Delta }t>1.5T_f,\end{array}$$

(20)

where $T_f=1.66\mathsf{\mu }\mathrm{s}$ is the front time, $T_t=13\mathsf{\mu }\mathrm{s}$ is the tail time, $A_0=e^{a_1{T}_0}-e^{a_2{T}_0}$, $a_1=0.08147180/T_f$, $a_2=2.358427881/T_t$, $T_0=\ln \left[\left(a_1/a_2\right)/\left(a_1-a_2\right)\right]$, ${\omega }_0=\mathsf{\pi }/3T_f$ and $I_{\max }=10 \mathrm{kA}$ is the current peak in (20). The function $I_s\left(t\right)$, given by (20) and its parametric values produce the waveform illustrated by Figure 5a (truncated at 10 $\mathsf{\mu }\mathrm{s}$), which is used for calculating the Fourier transform ${\overline{I}}_s\left(f\right)$, shown in Figure 5b in its absolute normalized form $\left|{\overline{I}}_s\left(f\right)\right|/{\left|{\overline{I}}_s\left(f\right)\right|}_{\max }$. The frequency domain lightning current has appreciable energy up to approximately f = 1 MHz, which implies that ${\lambda }_{min}\approx 106 \mathrm{m}$ and ${\mathsf{\Delta }}_{x,y,z}\le 10.6 \mathrm{m}$, in such a way that the cubic Yee cell edge dimension $\mathsf{\Delta }=0.1$ m used in this work satisfies the condition (17). Finally, in this work, time step ${\mathsf{\Delta }}_t$ is calculated by ${\Delta }_t=0.99\Delta /(v_{max}\sqrt{3})$, where $v_{max}$ is the light speed in a vacuum.

3. Research Proposal: Geometric and Parametric Configurations for Installing Distribution Surge Arresters and Insulators

By observing Figure 3b and Ampère’s law, as given by (7) or (13c), we see that it is possible to diminish insulator-induced voltage $V_{AB}$, given by (8), by installing the insulation device in the region of interest in Figure 3b. Observe from (7), (13c) and (9) that reducing circulation of $\overline{H}$ on planes parallel to the xy-plane, in the region of interest in Figure 3b, implies reducing $E_z$, consequently decreasing $V_{AB}$ if the insulator is installed in the region of interest. This effect can be achieved by allowing magnetic fields, which are generated by same-direction parallel electric currents, to superpose in the region of interest. During a lightning event, the currents should flow through two (or more) lightning arresters placed equidistantly and parallelly to the insulator (see Sub-item 2.1). Although surge arresters have been installed parallel to insulators before [27,28,29,30], full-wave electromagnetic analysis has not been performed. This way, based on full-wave FDTD modeling and analysis, we propose three geometrical configurations for positioning the insulator and electric power distribution class SAs at the low-voltage (LV) level, aiming at minimizing $V_{AB}$.

Figure 6a shows a typical Single-Phase Earth Return (SPER) system, often adopted for powering rural houses [35]. Those networks are frequently exposed to incidence of lightning discharges due to low shielding factor [4,10]. In SPER systems, each voltage line is supported by an insulator and it is protected by a single SA. If a line is hit by lightning current, electrons flow through the connected surge arrester to ground, producing magnetic fields similar to those illustrated by Figure 3a around the SA (see Figure 7a), which penetrate the insulator. Figure 6b shows one of the proposed configurations, in which each voltage line is protected by two SAs placed equidistantly and parallel to the corresponding line insulator. Once a lightning current hits a voltage line protected by two SAs, part of the surge current flows across each SA and, consequently, a magnetic field distribution similar to that in Figure 3b is produced, in such a way that the total magnetic field penetrating the insulator (which is mainly a result of superposed magnetic fields produced by currents flowing across SAs) is minimized (see also Figure 7b), thus, reducing the induced voltage across the insulator. The surge arresters’ impedances should ideally be equal and SA current phases should ideally be the same. Notice that, as previously addressed in [27,28,29], a strong decline in arrester current (such as the approximately 50% of proposed cases) decreases the dissipated energy and thermal stress in each device, mitigating occurrences of arrester failures. Additionally, reducing voltage $V_{AB}$ induced across insulator terminals has the effect of increasing the level of availability in the electric power system.

4. FDTD Computational Modelling

The computational modelling of a typical electric experimental setup for determining residual voltages in surge arresters is the starting point for the investigation of our proposal. The geometries represented in the FDTD-3D analysis region were inspired by the surge arrester residual voltage facility of the High- and Extra-High-Voltage Laboratory (LEAT) at the Federal University of Pará (UFPA). Figure 8a shows a photograph of the referred laboratory facility and Figure 8b illustrates the representation of the corresponding three-dimensional FDTD model produced using the software Synthesis and Analysis of Grounding Systems (SAGSs) [20]. In Figure 8b, the region where the SAs are modelled is detailed, in which $y_a$ is the distance between SAs (as indicated by red arrows). In this work, $y_a$ assumes the values 0.8 m, 0.6 m and 0.4 m.

Parameters of the Conceived FDTD Model

The experimental setup illustrated in Figure 8a is represented by conceiving a $28 \mathrm{m}\times 18.4\mathrm{m}\times 22.4\mathrm{m}$ analysis region using SAGS [20], which employs FDTD for performing realistic full-wave modeling. The Yee cell dimensions were set to ${\mathsf{\Delta }}_x={\mathsf{\Delta }}_y={\mathsf{\Delta }}_x=0.1 \mathrm{m}$, respecting (18) and dimensions and distances among the objects, in such a way to properly represent objects, electromagnetic waves and the associated wave phenomena (couplings and inductions, reflections, diffractions and refractions). Modeled elements and their respective identifications are shown in Figure 9a. Figure 9b,c depict the region of interest for typical and proposed arrangements, respectively.

Table 1. Electromagnetic parameters of materials represented in FDTD model.

Identification (Figure 9a)	Material	${\mathit{\epsilon }}_{\mathit{r}}$	$\mathit{\sigma }$$(\mathbf{S}/\mathbf{m})$	${\mathit{\mu }}_{\mathit{r}}$
C and F	Epoxy resin	4	0.0005	1
${\mathrm{L}}_1$$\mathrm{and}{\mathrm{L}}_3$	Concrete	8	${10}^{-6}$	1
${\mathrm{L}}_2$$\mathrm{and}{\mathrm{L}}_4$	Soil	4	0.022	1
K	Air	1	0	1

provides the electromagnetic properties of the modeled materials [32,36,37,38].

The analysis region was limited by using the Uniaxial Perfectly Matched Layer (UPML) technique [15]. The thin-wire technique developed in [17] was employed for performing sub-cellular representations of cylindrical conductors with transverse section areas much smaller than faces of Yee cells.

The free space among the horizontal SAs’ interconnection conductors and the fixation metallic base (between the SAs in Figure 9c) is where the insulator should be installed. Voltage is evaluated by a line-integrating vertical component in the electric field from the metallic base to the higher horizontal SA interconnection conductor, working as insulator terminals (see terminals A and B in Figure 9a,b, indicated using red color). Thus, the aim of the research is the evaluation of overvoltages between the insulator terminals caused by the direct incidence of atmospheric discharge on a phase line. The evaluation of the insulator performance during lightning hit is not on the scope of this paper.

The zinc oxide surge arresters were represented in our FDTD model as resistive elements. SA resistances were calculated from their $V\times I$ curve points (for $5 \mathrm{kA}$ and $10 \mathrm{kA}$), available in technical documentation of a commercialized device [39]. Such modelling is sufficient for performing investigations on induced overvoltages across insulators excited by magnetics fields produced by discharge currents flowing through surge arresters. When subjected to currents above $10\mathrm{kA}$, such SA behaves nearly as a resistor due to its fundamentally linear relation between voltage and current, as addressed in [40]. The representation of the grounding system in the model aimed at fidelity to the grounding project of LEAT-UFPA, encompassing the three grounding grids (${\mathrm{J}}_1$, ${\mathrm{J}}_2$ and ${\mathrm{J}}_3$) with distinct depths. Notice that the deep rods $\mathrm{I}$ interconnect the grounding grids and conduct reinforced concrete grid currents to soil.

5. Numerical Results and Discussion

Eight FDTD simulations were defined and executed, as detailed in

Table 2. Reference and proposed cases, which are numerically analyzed.

Identification	${\mathit{I}}_{\mathbf{max}} \left(\mathbf{kA}\right)$	$\mathbf{Resistance}\mathbf{of}\mathbf{SAs}\left(\mathbf{m}\mathbf{\Omega }\right)$	${\mathit{y}}_{\mathit{a}}$$\left(\mathbf{m}\right)$
Reference Case 1 (RC1)	5	280	-
Reference Case 2 (RC2)	10	170	-
Proposed Case 1 (PC1)	5	280	0.8
Proposed Case 2 (PC2)	10	170	0.8
Proposed Case 3 (PC3)	5	280	0.6
Proposed Case 4 (PC4)	10	170	0.6
Proposed Case 5 (PC5)	5	280	0.4
Proposed Case 6 (PC6)	10	170	0.4

. A typical line protection arrangement using a single SA connected parallelly to the insulator (see Figure 9b), subjected to lightning current peaks $I_{\max }$ set to $5\mathrm{kA}$ and $10\mathrm{kA}$, was modeled and referred to as reference cases RC1 and RC2. Furthermore, the proposal of using two SAs connected parallelly to the insulator (as Figure 9c illustrates) was also analyzed for lightning current peaks of $5\mathrm{kA}$ and $10\mathrm{kA}$, with $y_a$ set to $0.4\mathrm{m}$, $0.6\mathrm{m}$ and $0.8\mathrm{m}$. The cases associated with the proposed configuration are referred to as proposed cases PC1 to PC6 (see Table 2). For all cases, lightning current $I_s\left(t\right)$ given by (20) was used as an FDTD excitation source. Since the time derivative of lightning current $I_s\left(t\right)$ is more relevant during its few initial microseconds (see Figure 5), a total simulation time of $10\mathsf{\mu }\mathrm{s}$ was set for executing the FDTD simulations.

Figure 10a,b show the calculated transient surge arrester currents obtained for $I_{\max }=5\mathrm{kA}$ and $I_{\max }=10\mathrm{kA}$, respectively. Figure 10 shows that cases 1 and 5 (Figure 10a) and 2 and 6 (Figure 10b) are in accordance with Kirchhoff’s current law. As one would expect, the difference between current levels of reference and proposed cases is approximately 50% for all time (since resistances of SAs ${\mathrm{A}}_{\mathrm{a}}$ and ${\mathrm{A}}_{\mathrm{b}}$ are equal for the proposed cases). This also agrees with the results verified in [27,28,29]. Notice that as $y_a$ is reduced, arrester current is slightly lessened since electromagnetic coupling between SAs increases.

Figure 11 shows the calculated voltage curves obtained in the region of interest, across insulator terminals A and B, as indicated in Figure 9b,c. At the initial instants of the transient voltages (from 0 to 2.5 μs), the reference case displays peak values of approximately 2.5 kV (Figure 11a) and 4.5 kV (Figure 11b). In contrast, the implementation of the proposed strategy provides reductions in maximum residual voltages to 770 V and 1.3 kV, respectively, corresponding to decreases of 68.96% and 70.92%, as depicted in Figure 12. Those strong voltage attenuations contribute to minimizing electric stresses across the insulator. As previously pointed out, the obtained voltage reductions are due to the fact that the discharge currents conducted by the SAs installed parallelly to the insulator create opposite-direction magnetic fields that superpose in the region of interest, thus, resulting in a reduced voltage across insulator terminals. Additionally, this phenomenon is more noticeable at the initial instants of the transient phenomena.

Furthermore, by analyzing Figure 11, we see that the configurations with $y_a=0.4\mathrm{m}$ (PC5 and PC6) show, approximately, a 27.3% reduction in residual voltages with respect to those in which $y_a=0.8\mathrm{m}$ (PC1 and PC2). Such a result agrees with the theoretically expected increase in the level of electromagnetic coupling associated with reducing distances between SAs. The calculated voltage reductions remain above 50% in all proposed cases with respect to reference cases, suggesting that not only the Ohmic law is responsible for the decreases, but also the interactions between the transient magnetic fields reduce the induced voltage on the insulator, as illustrated in Figure 12.

For the sake of analysis clarity regarding the phenomena under consideration, the spatial distributions of $\left|\overline{H}\right|$ and $\left|\overline{E}\right|$ on planes xy and yz, respectively, obtained with the FDTD method at 3.33 µs, are shown in Figure 13. Figure 13a shows that, in fact, $\left|\overline{H}\right|$ is reduced in the region between surge arresters ${\mathrm{A}}_{\mathrm{a}}$ and ${\mathrm{A}}_{\mathrm{b}}$, minimizing $\left|\overline{E}\right|$ in the same region, as clearly revealed by Figure 13b. The above agrees, qualitatively, with the results seen in [1], since it is a different application of this proposal.

The promising benefits of the proposed SA arrangement, presented so far from a technical point of view, give rise to a brief cost-effectiveness analysis. Assuming the high exposure of rural networks to the incidence of lightning strikes [4] and the impacts/damage of these transients on the power system [3,4,5], mainly in tropical areas (as is the case in Brazil) [1], by themselves, already suggest the implementation of some remedial measures. In this sense, the total cost of adding two more LV SAs per pole (see Figure 6b) was estimated to be between USD 80 and 90 as of January 2023. However, considering the lightning density in the Brazilian Amazon region [1] and the measured parameters of lightning currents in Brazil [41], for example, 50% of lightning currents can reach 45 kA, this surpasses the maximum current withstand threshold of LV SAs [39]. In this context, balancing the currents in the SAs of the proposed arrangement mitigates: (1) the lightning-induced voltages across network insulators (or across bushings on the LV side of medium-voltage transformers, for instance); (2) the possibility of SA failures due to thermal avalanche and also, if one of the devices fails, the network would still be protected by a surge arrester, such as in the reference single-SA case; (3) the costs involved with customer claims for damages caused by power supply interruptions, which reached approximately USD 703,000.00 to a single power distribution company in Brazil in 2019 [42]; and (4) priceless losses, such as disservice of life support in hospitals, in addition to damage to the public image of power companies.

6. Conclusions

In this work, it is demonstrated that reductions in induced voltages across the terminals of overhead electric power line insulators can be achieved with the use of surge arresters (SAs) installed parallelly to the insulator. Once a lightning current hits the line, SAs provide the division of the surge electric current, producing interaction of opposite-direction transient magnetic fields in the region in which an insulator is to be installed. As observed, mainly during the initial instants of a lightning current, voltage induced across an insulator can be substantially reduced, as demonstrated by the results in this study. This analysis has not been addressed in previous research.

The numerical evaluation of the induced voltages at the insulator terminals between a reference configuration and a proposed configuration with two surge arresters (SAs) per phase installed parallelly to the insulator and equidistant from the insulating device showed reductions of up to 71%, the highest reduction being obtained when the SAs were placed 0.4 m apart. This reduction can be attributed to the influence of high electromagnetic coupling between the SAs at the considered distance. More specifically, the obtained spatial distributions of magnitudes of electric field $\overline{E}$ and magnetic field $\overline{H}$, at a given time instant immediately after the occurrence of the peak of the fast electromagnetic transient current, show that horizontal components in $\overline{H}$ are reduced at the insulator installation location, resulting in a corresponding decrease in vertical components of $\overline{E}$ in the same region, producing, as a result, a reduction in induced voltage across insulator terminals.

Lessening the voltages developed across insulator terminals has the effect of increasing the level of availability of electric power supply systems, mitigating damage caused by atmospheric transient currents to devices of power companies and customers, being a remedial measure with the capability to be highly cost-effective. Although installation of SAs with very different $V\times I$ characteristics is possible, premature aging and/or failure of one of the SAs in the array and installation mistakes are some risks in the proposed arrangement under real-world installation conditions (which can lead to unbalanced SA discharge currents), and such circumstances still tend to produce responses with advantages over the single-SA configuration, as long as horizontal components in the magnetic field are reduced in insulators.

Future studies involving the evaluation of the proposed arrangement at medium-, high- and extra-high-voltage levels (including in DC systems) and the numerical analysis in standardized bus test feeders (or based on real-world networks) can be carried out. Regardless of the advances achieved in this paper, studies involving laboratorial experiments and refined numerical modelling of three or more SAs are also proposed and are presently in progress.