*2.2. Field-to-Line Coupling*

The field-to-line coupling computation is obtained considering the well-known Agrawal model [26], which is here presented in its extended version taking into account the presence of a finite-conducting ground and a multi-conductor line.

$$\begin{cases} \frac{\partial V\_i^s}{\partial \mathbf{x}}(\mathbf{x}, t) + \sum\_{j=1}^M L\_{ij} \frac{\partial I\_j}{\partial t}(\mathbf{x}, t) + V\_i^\%(\mathbf{x}, t) = E\_{\text{inc}, \mathbf{x}, j}(\mathbf{x}, t) \\\frac{\partial I\_i}{\partial \mathbf{x}}(\mathbf{x}, t) + \sum\_{j=1}^M \mathbb{C}\_{ij} \frac{\partial V\_j^s}{\partial \mathbf{f}}(\mathbf{x}, t) = 0 \end{cases} \tag{1}$$

with

$$V\_i^{\mathcal{S}}(\mathbf{x}, t) = \int\_0^t \xi\_{\mathcal{S}}^i(t - s) \frac{\partial I\_i}{\partial s}(\mathbf{x}, s) ds \tag{2}$$

where *Vsi* (*<sup>x</sup>*, *t*), *Ii*(*<sup>x</sup>*, *t*) and *Einc*,*x*,*<sup>i</sup>*(*<sup>x</sup>*, *t*) are the scattered voltage, the current and the tangential component of the exciting electric field (computed in the previous subsection) on the *i*th conductor at distance *x* from the beginning of the line. As expressed in Equation (1), the knowledge of the inductance and capacitance matrices (*L* and *C*) is required. Please note that *ξgi*is the time-domain expression for the ground impedance [27].

The total voltage occurring on the *i*-th conductor at the point *x* can be then expressed as the sum of the scattered voltage and the incident voltage, whose value depends on the vertical electric field (computed in the previous subsection).

The proposed methodology is adapted to an EMT-type software (in this framework Simulink-Simscape is used), through the finite-difference time-domain (FDTD) technique. In this case, a second-order scheme is adopted with *dt* = 10 ns and *dx* = 9 m, which satisfies the well-known Courant stability condition. Further details can be found in [22].

### **3. Tower Modeling**

The modeling of the tower is usually neglected in lightning-induced voltages studies. However, in this framework, the tower is included in the model according to [28,29] and is modeled as a lossless transmission line, whose characteristic impedance is:

$$Z\_{\varepsilon} = 60 \left[ \ln \left( \sqrt{2} \frac{2h}{r} \right) - 1 + \frac{r}{4h} + \left( \frac{r}{4h} \right)^2 \right] \tag{3}$$

*h* being the tower height and *r* the tower radius.

### **4. Electrical Grounding Modeling**

In this paper, the grounding transient behavior is modeled by HEM [18,19]. This model is an electromagnetic computational method developed for the numerical solution of lightning problems and, according to CIGRÉ [30], it is classified as a hybrid electromagneticcircuit approach. The main motivations for using HEM are as follows: (i) it is accurate and flexible, i.e., it can be used in different types of grounding configuration; (ii) its results have been extensively validated experimentally, such as measurements in TS [18,19], horizontal electrodes [18,31], vertical rods [31,32], and typical substation grounding grids [33] and (iii) it is faster than traditional full-wave methods (without losing accuracy). It is worth mentioning that the usage of HEM has increased significantly recently [30,34,35].

Basically, HEM consists of subdividing the actual system (in this case, electrical grounding) into N small conductive cylindrical segments and, for each segment, the electromagnetic theory is applied. After that, by using the circuit theory, it is possible to obtain a matrix system that computes the wideband response of the electrical grounding. For the sake of clarity, we present a brief overview of HEM below. More details about HEM are described in [18,19].

It is worth commenting on the fact that HEM corresponds to an electromagnetic model developed in the frequency domain. Thus, it is necessary first to determine the frequency spectrum (depending on the phenomenon of interest). After that, the electrical grounding is divided into N segments, where the length of each segmen<sup>t</sup> is equal to 10 times its radius (thin wire approximation). A discussion about segmentation length is presented in [36].

Each segmen<sup>t</sup> is considered a source of two currents, one longitudinal that flows along the electrode (*IL*) and another transversal that flows from the electrode to the surrounding soil (*IT*). It is worth noting that *IL* generates a non-conservative electric field and *IT* a conservative one. With the aid of the magnetic vector and electric scalar potentials, both voltage drops ( Δ *V*) and electric potentials ( *V*) in each pair of segments (transmitter and receiver) are determined. Additionally, double integral equations are established for Δ *V* and *V*. These integrals depend on the frequency, geometry, soil parameters and *IT* and *IL* distributions. However, the distributions of *IT* and *IL* are not known and are integrands of the integrals. From this point on, the Method of Moments (MoM) is applied to solve these integral equations [37]. The effect of the air-soil interface is included using the method of images, similar to [38,39].

The *IT* and *IL* distributions considered in this paper are of the piecewise-constant function type [18,19]. MoM makes it possible to transform integral equations into algebraic ones, the solution of which allows determining all the quantities of interest (in the frequency domain), such as *IT*, *IL*,Δ *V* and *V* distributions; transverse (capacitive and conductive couplings) and longitudinal (resistive and inductive couplings) impedances (self and mutual); electromagnetic field; harmonic grounding impedance ( *<sup>Z</sup>*(*ω*)); low-frequency grounding resistance ( *RLF*), etc.

Also, it has been documented in the literature, over almost one hundred years [40], that the soil is a dispersive medium, i.e., the response is not instantaneous. Several researchers have presented a numerical solution to consider this dispersivity in the frequency domain, such as [31,40–49]. According to [50], for the values of conductivity considered in this paper, the frequency-dependent soil electric parameters can be neglected in the computation of the EM fields that illuminate the line. On the other hand, according to [17], the impact of

considering the frequency-dependence of the soil in the grounding modeling generates sensible differences. Thus, in this paper, the frequency-dependence of the soil parameters is considered in the grounding modeling. Use is made by the formulations proposed in [31], since a CIGRE Brochure has suggested them [51]. Equations (4) and (5) illustrate the formulation. 

$$
\sigma(f) = \sigma\_0 + \sigma\_0 h(\sigma\_0) \left( f \times 10^{-6} \right)^{\gamma} \tag{4}
$$

$$
\epsilon\_r(f) = \epsilon\_{r\infty} + \frac{\tan(\pi\gamma/2) \times 10^{-3}}{2\pi\epsilon\_0(10^6)^{\gamma}} \sigma\_0 h(\sigma\_0) f^{(\gamma - 1)} \tag{5}
$$

where *σ*(*f*) is the frequency-dependent soil conductivity (in mS/m), *σ*0 is the low-frequency soil conductivity (in mS/m),  *r*(*f*) is the frequency-dependent soil permittivity and 0 is the vacuum permittivity. To obtain mean results (more details about it in [31]), one should use *h*(*<sup>σ</sup>*0) = 1.26*σ*<sup>−</sup>0.73 0 , *γ* = 0.54 and  *r*∞ = 12.

As specified before, the calculation of induced voltages is performed directly in the time domain; however, *<sup>Z</sup>*(*ω*) is a frequency domain quantity. Thus, the well-known Vector Fitting (VF) approach is used for fitting the calculated frequency domain grounding response with rational function approximations [20]. The passivity is enforced by perturbation [21].

Finally, based on the obtained rational function, it is possible to synthesize an electric network that can be promptly included in the time-domain simulation. It is important to note that this electric circuit generates the same frequency response as the harmonic grounding impedance provided by HEM. Thus, it includes reactive and electromagnetic wave propagation effects.

### **5. Test Cases**

This section presents the test cases related to the comparison between two different grounding modeling, i.e., the low-frequency grounding resistance (*RLF*) and the harmonic grounding impedance (*Z*(*ω*)).

Let us consider a 1.2 km matched three-phase DS (Figure 1). The three-phase conductors' heights are 10, 11 and 12 m, respectively, while the shield wire height is 14 m. The horizontal distance between each conductor and the shield wire is 2.4 m. The conductors' diameter is 1.83 cm, while the shield wire diameter is 0.72 cm.

The span between each tower is 50 m. To consider the influence of the adjacent towers, a total of five towers are modeled in detail, because the towers in longer distance have little impact on the overvoltage. According to [52], for lightning-related phenomena adjacent towers place a moderate influence. Each tower is 14 m high and with a base diameter of 0.5 m. According to Equation (3), a value of *Zc* = 244.17 Ω is considered. The propagation velocity along the tower is considered to be 0.8*c* [53,54]. Moreover, the insulators are modeled taking into account their parasitic capacitance. The value of parasitic capacitance of each insulator is 7.68 nF, and they were calculated considering the information in [55–57].

Each tower is grounded with a grounding system as shown in Figure 2. This is a typical configuration for grounding distribution networks in the State of Minas Gerais, Brazil. It consists of three vertical rods 2.5 m long interconnected by a horizontal galvanized steel cable 6 m long. The vertical rods are copper-plated steel, with a diameter of 15 mm.

The equivalent circuit of the system composed of a three-phase distribution line, tower and grounding system is shown in Figure 3.

**Figure 1.** Line configuration.

**Figure 2.** Grounding grid of the distribution tower.

Two cases for the soil parameters will be considered. The soil conductivity and permittivity will be frequency-dependent according to [31], where *σ*0 = 10 mS/m and 1 mS/m, respectively. These cases correspond to two different grounding harmonic responses according to the grounding modeling proposed in Section 4. Figures 4 and 5 show *<sup>Z</sup>*(*ω*) and *RLF* of the two considered cases. Based on the behaviors described in these figures, it is possible to verify that: (i) grounding can only be represented by *RLF* in the low-frequency range, where *<sup>Z</sup>*(*ω*) tends to *RLF*; (ii) the limit frequency of the low-frequency range increases with a reduction in conductivity; (iii) in the intermediate-frequency range there is a predominance of capacitive behavior of the grounding, verified by the decrease of *<sup>Z</sup>*(*ω*) in relation to the *RLF*; (iv) the limit frequency of the intermediate-frequency range also increases with the reduction in conductivity and (v) only in the high-frequency range inductive effect is predominant, mainly for higher conductivity values. Thus, the response of the system under study (DS and grounding) will be a direct function of the frequency spectrum of the electromagnetic signal that requests it. As a consequence, it is expected that the overvoltages in the insulator string are sensitive to grounding modeling.

**Figure 3.** Equivalent circuit of the power system, tower and grouding.

**Figure 4.** Grounding harmonic impedance *<sup>Z</sup>*(*ω*) with *σ* = 10 mS/m. Adjusted Z refers to the impedance obtained with the equivalent circuit.

When we consider a grounding model described by *RLF*, the implementation in the EMT-type software is trivial, while when we consider the harmonic grounding impedance, it is possible to obtain the synthesis of the electric circuit to be implemented in the EMT-type software by using the approach presented in Section 4.

The general layout of the circuit obtained from the Vector fitting approach is described in Figure 6, while the values of the passive circuit are proposed in Tables 1 and 2 for *σ* = 10 mS/m and in Tables 3 and 4 for *σ* = 1 mS/m. It is worth mentioning that these equivalent circuits are mathematical models that have a frequency response very close to *<sup>Z</sup>*(*ω*), but their electrical parameters do not have physical consistency, it is also important mentioning that although some elements may have negative values, the circuit is passive in overall. Hence the existence of negative values for resistance, inductance and capacitance in Tables 1–4 do not mean that the circuit is not passive. The i-index appearing in Tables 1–4 refers to the electrical branch.

**Figure 5.** Grounding harmonic impedance *<sup>Z</sup>*(*ω*) with *σ* = 1 mS/m. Adjusted Z refers to the impedance obtained with the equivalent circuit.

**Figure 6.** General layout of the grounding circuit for the harmonic grounding impedance.


**Table 1.** Elements of Real Poles of the grounding circuit. *σ* = 10 mS/m.

**Table 2.** Elements of Complex Poles of the grounding circuit. *σ* = 10 mS/m.



**Table 3.** Elements of Real Poles of the grounding circuit. *σ* = 1 mS/m.

**Table 4.** Elements of Complex Poles of the grounding circuit. *σ* = 1 mS/m.


To compare the grounding modeling, 12 different tests have been implemented (Table 5), each one differing for the soil conductivity, stroke location and stroke type (first or subsequent). The stroke location is always placed in front of the middle of the line, Figure 7 illustrates the distance between the lightning-channel and the tower under study. It is important to highlight that if the closest point of the DS to the stroke location is in the mid-span, it would, naturally, change the maximum overvoltage but not the perceptual differences between the approaches. The lightning return stroke channel is characterized by a height of 8 km and a speed equal to one-half the speed of light in a vacuum. The channel-base current is modeled as a sum of two Heidler's functions as in Equation (6), with parameters reported in Table 6. The representation of the channel-base current is proposed in Figures 8 and 9.

$$I\_0(t) = \frac{I\_{01}}{\eta\_1} \frac{\left(\frac{t}{\tau\_{11}}\right)^{\eta\_1}}{1 + \left(\frac{t}{\tau\_{11}}\right)^{\eta\_1}} e^{-\frac{t}{\tau\_{12}}} + \frac{I\_{02}}{\eta\_2} \frac{\left(\frac{t}{\tau\_{21}}\right)^{\eta\_2}}{1 + \left(\frac{t}{\tau\_{21}}\right)^{\eta\_2}} e^{-\frac{t}{\tau\_{22}}} \tag{6}$$

being

$$\eta\_i = \exp\left(-\frac{\tau\_{i1}}{\tau\_{i2}} \left(n\_i \frac{\tau\_{i2}}{\tau\_{i1}}\right)^{\frac{1}{n\_i}}\right) \tag{7}$$

**Table 5.** Test details.


**Figure 7.** Sketch that illustrates the distance between the tower under study and the lightning stroke Location.


**Table 6.** Heidler's current parameters.

**Figure 8.** Channel-base current: Heidler's current-first stroke.

**Figure 9.** Channel-base current: Heidler's current-subsequent stroke.
