Transient electromagnetic field radiated by grounding systems caused by lightning strike in a dissipative half-space

Abstract

In this work, a method is proposed for the calculation of the electromagnetic field radiated by a grounding system in a semi-infinite dissipative medium subjected to a current excitation generated by the lightning strike. The assessment of the induced current distribution on the grounding system due to a transient excitation is obtained by solving Pocklington’s equation using the method of moments. A new method for the calculation of the transient electromagnetic field following the injection of a lightning current or short circuit is also developed. This method is based on the concept of the hertzian dipoles and the asymptotic expansion of a single dipole. A MATLAB computer program has been developed for the determination of the induced currents. The program is validated by comparing its results with those obtained with the commercial software NEC 4.

Keywords:

• lightning
• electromagnetic field
• Hertzian dipole
• NEC4
• transient
• Sommerfeld integrals
• grounding

Public Interest Statement

When lightning strikes an object on the ground, a large amount of electric current flows from the cloud to the object in less than one thousandth of a second, causing a multitude of effects ranging from mechanical damage, overheating, biological effects on humans and livestock, to indirect effects in the form of induced voltages in nearby metallic structures. This paper intends to contribute in better assessing such effects on metallic grounding systems buried underneath electrical installations, by proposing a method to compute the electromagnetic (EM) radiation generated by the flow of such large current. Evaluation of the EM radiation is of value in assessing the safety risks due to exposure of working persons and people to induced lightning effects, and also in evaluating the impact on sensitive electronic equipment. It can also help electric utility companies design better mitigating measures to protect the community.

1. Introduction

The increasing use of sensitive electronics and digital technologies in power system equipment has raised significant interest in power system electromagnetic compatibility studies.

When lightning hits an electrical substation, large currents generated by the stroke flow in the above ground structures and dissipate in the soil through the substation grounding system. The electromagnetic fields generated by such high currents may cause damage to equipment and may be dangerous to personnel working nearby.

For this purpose, many researchers have developed different models for analyzing the transient behavior of grounding systems under lighting strikes. These models can be classified as circuit approach (Araneo et al., 2015; Celli, Ghiani, & Pilo, 2010; Gatta, Geri, Lauria, & Maccioni, 2014; Geri, 1999), transmission-line approach (Gazzana et al., 2014, 2017; Kherif, Chiheb, Teguar, Mekhaldi, & Harid, 2017; Liu, Theethayi, & Thottappillil, 2005; Lorentzou, Hatziargyriou, & Papadias, 2003; Theethayi, Baba, Rachidi, & Thottappillil, 2008), and electromagnetic field approach (Karami, Sheshyekani, & Rachidi, 2017; Li, Li, Fan, & Yin, 2011; Poljak, Sesnic, Cavka, & Drissi, 2012; Qi, Cui, Zhao, & Li, 2007; Zhang, Cui, Zhao, Yin, & Li, 2003).

The electromagnetic field approach is the most rigorous method for modelling the transient behaviour of a grounding system. It solves the full set of Maxwell’s equations with minimum approximations. This approach can be implemented either by the Method of Moments (MOM) or by the Finite Element Method (FEM). In all above-mentioned models, the effect of lossy ground is taken into account using the modified image theory or by approximation in Sommerfeld integrals. This implies that if the effect lossy ground is ignored it leads to erroneous results in predicted electromagnetic fields.

The electromagnetic fields radiated by a grounding system above a semi-infinite dissipative medium subjected to a current excitation generated by the lightning strike are evaluated using the concept of the Hertzian dipoles. In this method, the effect of lossy ground is taken into account by the Sommerfeld integrals (Sommerfeld, 1949; Ye & Jin, 2010). This reflects the effect of finite conductivity of the interface, which requires a correction term, due to the air-soil interface, for the complete solution of the electric field. In addition, the Sommerfeld integrals depend on the horizontal distance r, measured parallel to the interface between the source and observation point and on the vertical distance between the observation point and the source image. There are many published numerical and analytical techniques dealing with the integrals arising in the Sommerfeld half-space problem, which have been reviewed in Chatterjee, Rao, and Kluskens (2016), Hochman and Leviatan (2010), Rahmat-Samii, Mittra, and Parhami (1981). Subsequently, a comprehensive formulation on the development of the lossy half-space problem, using Hertz potentials, was given by Banos (1966). The solution was expressed in terms of the solutions of the two basic Sommerfeld integrals V₁ and U₁, which are evaluated by numerical integration along contours in the complex plane and used subsequently to obtain the integral values needed for the integral equation solution using the method of moments.

Modeling this problem for the quantification of radiated electromagnetic fields resulted in the proposal of some computational approaches, especially in the frequency domain. The best known in this area is the one that uses the antenna theory.

In this paper, we analyze the performance of grounding system subjected to a lightning strike using a two-stage approach. In the first stage, a model of the grounding system is developed, and in the second, the electromagnetic field radiated by the grounding system under impulse current is computed. The solution is obtained using antenna theory and the method of moments, in the case of direct injection of a lightning current. In this procedure, the lightning current waveform is described using Heidler’s formula, and is injected at one point of the grid. The grid’s transient response is computed in the frequency domain by solving the electromagnetic field integral equation. The temporal and spatial distributions of the electromagnetic field near the grounding system are obtained by the inverse Fourier transform (Song & Chen, 1993). The proposed approach is applied to the case of a horizontal conductor buried in uniform soil and a grounding grid. To validate the approach, the results are compared with those obtained using the commercial software NEC-4 (Burke, 1992).

2. Analytical expressions for dipole in infinite conducting medium

The method consists of subdividing the wire structure element into small elements called dipoles (Figure 1), whose size is chosen in order to mask the propagation.

Figure 1. Segmented structure in dipoles.

The field at any point of the medium is obtained by superimposing the contributions of all dipoles.The size of the dipoles “dl” must meet the following two conditions:(1) $$\text{(a)}dl\le \mathit{\lambda }/20$$(1)

where λ is the smallest significant wavelength. This condition can mask the propagation along the dipole; the current amplitude and phase along the dipole are independent of z.(2) $$\text{(b)}dl\le R/10$$(2)

where R is the distance between the dipole under consideration and the observation point.

This condition can take into account small variations in current from a point very close to the wire structure.

In his work, Sivaprasad and King (1963) provides compact expressions for calculating the field radiated by a dipole located in a semi-dissipative medium. In harmonic regime and cylindrical coordinates, the electromagnetic field created by a horizontal dipole in a semi-infinite medium is expressed in terms of Sommerfeld integrals.

The electromagnetic fields of a horizontal dipole in a semi-infinite medium, due to a source are expressed in terms of the solutions of the two basic integrals V₁ and U₁. The derivations of the expressions of the field are given in Harrigton (1968) using cylindrical coordinates (r, θ, z).

The rigorous expressions for the r, θ and z-components of the electric and magnetic field at a near point due to a source current, in a local cylindrical coordinate system illustrated in Figure 2, are (Harrigton, 1968):(3) $$E_{rl}=\frac{-j{P}_z}{4\mathit{\pi }{\mathit{\epsilon }}_1}cos\mathit{\theta }\left\{\frac{{\mathit{\partial }}_{}^2}{\mathit{\partial }{r}^2}\left[\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}+V_1\right]+{\mathit{\beta }}_1^2\left[\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}+U_1\right]\right\}$$(3)

Figure 2. Coordinate system for one dipole.

(4) $$E_{\mathit{\theta }l}=\frac{-j{P}_z}{4\mathit{\pi }{\mathit{\epsilon }}_1}sin\mathit{\theta }\left\{\frac{1}{r}\frac{\mathit{\partial }}{\mathit{\partial }r}\left[\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}+V_1\right]+k_1^2\left[\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{j{k}_{1}R{}_2}}{R_2}+U_1\right]\right\}$$(4)

(5) $$E_{\mathit{zl}}=\frac{j{P}_z}{4\mathit{\pi }{\mathit{\epsilon }}_1}cos\mathit{\theta }\left\{\frac{\mathit{\partial }{}^2}{\mathit{\partial }z\mathit{\partial }r}\left[\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}+V_1\right]-n^2\frac{{\mathit{\partial }}^2{V}_1}{\mathit{\partial }r\mathit{\partial }z}\right\}$$(5) (6) $$H_{\mathit{rl}}=\frac{P_{z}sin\mathit{\theta }}{4\mathit{\pi }{k}_1^2}\frac{\mathit{\partial }}{\mathit{\partial }z}\left\{\left(1+n^2\right)\frac{1\mathit{\partial }{V}_1}{r\mathit{\partial }r}+k_1^2{U}_1-\frac{1\mathit{\partial }}{r\mathit{\partial }r}\left(\frac{2e^{jk{}_{1}R{}_2}}{R_2}\right)++k_1^2\left(\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}\right)\right\}$$(6) (7) $$H_{\mathit{\theta }l}=\frac{P_{z}cos\mathit{\theta }}{4\mathit{\pi }{k}_1^2}\frac{\mathit{\partial }}{\mathit{\partial }z}\left\{\left(1+n^2\right)\frac{{\mathit{\partial }}^2{V}_1}{\mathit{\partial }{r}^2}+k_1^2{U}_1-\left(\frac{2e^{jk{}_{1}R{}_2}}{R_2}\right)++k_1^2\left(\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}\right)\right\}$$(7) (8) $$H_{\mathit{zl}}=-\frac{P_{z}sin\mathit{\theta }}{4\mathit{\pi }}\frac{\mathit{\partial }}{\mathit{\partial }r}\left\{\frac{e^{jk{}_{1}R{}_1}}{R_1}-\frac{e^{jk{}_{1}R{}_2}}{R_2}+U_1\right\}$$(8)

where(9) $$V_1=\underset{-\infty }{\overset{+\infty }{\int }}\frac{k{}_1^2{H}_0^{\left(1\right)}\left(\mathit{\lambda }r\right)}{k{}_2^2\sqrt{{\mathit{\lambda }}^2-k{}_1^2}+k{}_1^2\sqrt{{\mathit{\lambda }}^2-\mathit{\beta }{}_2^2}}exp\left(-\sqrt{{\mathit{\lambda }}^2-k{}_1^2}\left(z+z_0\right)\right)\mathit{\lambda }d\mathit{\lambda }$$(9)

(10) $$U_1=\underset{-\infty }{\overset{+\infty }{\int }}\frac{2H_0^{\left(1\right)}\left(\mathit{\lambda }r\right)}{k{}_2^2\sqrt{{\mathit{\lambda }}^2-k{}_1^2}+k{}_1^2\sqrt{{\mathit{\lambda }}^2-\mathit{\beta }{}_2^2}}exp\left(-\sqrt{{\mathit{\lambda }}^2-k{}_1^2}\left(z+z_0\right)\right)\mathit{\lambda }d\mathit{\lambda }$$(10)

$H_0^{\left(1\right)}\left(\mathit{\lambda }r\right)$: Hankel function.

The electric and magnetic field at a point near the grounding system is the sum of contributions from the overlapped sinusoidal current sources in all segments. Asymptotic solutions obtained in Lytle and Lager (1974) and Sivaprasad (1962) may be used if the following conditions are satisfied (Lytle & Lager, 1974):

(1)	The source and the receiver are close to the interface.
(2)	The refractive index $\left|n^2\right|=\left|\frac{k_2^2}{k_1^2}\right|<<1$.
(3)	k1 is complex and has a large attenuation constant.
(4)	Re(k1r) >> 1 and k2r > 1.

3. Computation of the transient current distribution in a grounding system

The geometry of interest, shown in Figure 3, is a grounding grid system composed of interconnected conductors.

Figure 3. Geometry of a grounding system.

The grounding system is assumed to be a network of interconnected straight cylindrical metallic conductors with arbitrary orientation and finite conductivity. The conductors are subject to the thin-wire approximation, assuming that their radiuses are much smaller than the wavelength and their lengths are much greater than the radiuses.

The integral-differential equation of the Pocklington type for the unknown current I(r′) induced along the thin wire of an arbitrary shape buried in a lossy half-space is given by:(11) $$E^s=-\frac{1}{4\mathit{\pi }j\mathit{\omega }{\mathit{\epsilon }}_{\text{eff}}}\underset{l}{\int }{t}^{\prime }.I(r^{\prime }).\overline{r}.{\overline{r}}^{\prime }\left[\frac{{\mathit{\partial }}^2}{\mathit{\partial }r\mathit{\partial }{r}^{\prime }}-k^2\right]G(r,r^{\prime })d{r}^{\prime }$$(11)

where G(r, r′) is the corresponding Green’s function:(12) $$G(r,r^{\prime })=g_0(r,r^{\prime })-g_i(r,r^{\prime })+g_{\text{su}}(r,r^{\prime })$$(12)

E^s: the total scattered electric field along the surface of the conductor.

t′: the unit vector along wire axis, t′

I(r′): the current flowing along the conductor.

G(r, r′): the dyadic Green’s function for an unbounded medium with that corresponding to a half-space conducting medium.

g₀(r, r′) and g_i(r, r′) are Green’s dyadic functions for the electric field at r due to the current element and its image.

The third and last term g_su(r, r′) is a correction term necessary for the complete solution of the electric field.

4. Validation: Buried horizontal electrode

In order to validate the approach we propose for calculating the transient electromagnetic field emitted by a ground electrode, we consider a 20 m horizontal electrode of 0.07 m radius excited by a lightning current. For this application, a computer program was developed in MATLAB and the results are compared with those obtained with the Numerical Electromagnetic Code (NEC-4). Figure 4 shows the horizontal conductor in homogeneous soil with a resistivity 500 Ωm, buried at a depth h = 1 m. The soil’s relative permittivity and relative permeability are assumed equal to 10 and 1 respectively. The electrode is energized at one end by a lightning current having a waveform described by Heidler’s formula, as recommended by international standards on lightning protection (IEC 61312-1, 1995; ITU-T Recommendation K.40 (10/96), 1996):(13) $$i(t)=\frac{I_{max}}{\mathit{\eta }}\frac{{(t/{\mathit{\tau }}_1)}^{10}}{1+{(t/{\mathit{\tau }}_1)}^{10}}exp\left(-\frac{t}{{\mathit{\tau }}_2}\right)$$(13)

Figure 4. Buried horizontal electrode.

where I = 10.7 kA current peak value, η = 0.93 correcting factor, τ₁ = 19 μs, τ₂ = 0.8 μs – front and tail time constants.

The results illustrated in Figure 5 show the current at the injection point and at 10 m along the electrode length obtained by the proposed model and by NEC-4. For the two points of the electrode, we note a perfect agreement in amplitude and shape, particularly on the rising front (i.e. at the highest frequencies) with a slight difference on the impulse tail. This can be explained by the fact that our approach uses an impulse function as test function (point-matching method) while the software NEC-4 uses sinusoidal functions as test functions. Both results show that half-way through the electrode, more than 50% of the current have dissipated into the surrounding soil.

Figure 5. Current at different points along electrode.

Figures 6 and 7 show the total radiated electric field E_t and y- components of electric fields E_y computed at points P₁ (in the soil) and P₂ (in the air) respectively. As can be seen in Figure 6, there is satisfactory agreement between the proposed approach and the results obtained with NEC-4, with a slight over estimation of the peak values compared to NEC-4 results. Above ground (point P₂), a similar difference is observed with magnitudes but a somewhat faster rising time is predicted by the proposed approach. The difference may be due to the choice of the segment length introducing inaccuracies. The magnetic fields show better concordance with the NEC method on the rising impulse front, and significant differences on the tail, as shown in Figure 8.

Figure 6. Variation of the electric field at point P₁.

Figure 7. Variation of the electric field at point P₂.

Figure 8. Variation of the magnetic field at point P₁.

Figures 9 and 10 show the time-domain y-component of electric field and total electric field, respectively, at the earth’s surface, along the selected profile. For the typical lightning impulse current, the electric field on the earth surface can reach up to 1.5 kV/cm or higher at the point of energization. Figure 11 shows the temporal variation of the magnetic field on the earth’s surface, along the y-direction profile shown in Figure 4. Magnitudes of more than 140 A/m can be generated as a result of the lightning impulse current, corresponding to a peak magnitude field density of approximately 176 μT. The computed values of electric and magnetic field are in excess of the reference limits defined by ICNIRP (ICNIRP, 1998) for the frequencies that are within the lightning frequency spectrum.

Figure 9. Y-component of electric field along Profile 1.

Figure 10. Total electric field along profile 1.

Figure 11. Magnetic field along profile 1.

5. Case of a grounding grid

The substation grounding system consists of an above-ground structure bonded to a 40 × 40 m, 16-mesh square ground grid buried at a depth of 0.6 m. The soil is assumed to be homogeneous with a 500 Ωm resistivity, a relative permitivity of 10, and relative permiability of 1. The grid is energized by a lightning current waveform described with Heidler’s formula given in Equation (13(13) $$i(t)=\frac{I_{max}}{\mathit{\eta }}\frac{{(t/{\mathit{\tau }}_1)}^{10}}{1+{(t/{\mathit{\tau }}_1)}^{10}}exp\left(-\frac{t}{{\mathit{\tau }}_2}\right)$$(13) ). Two injection points were considered: first at one corner of the grid and then at the centre point of the grid. This surge current gradually dissipates in the soil while travelling along the conductors of the structure and ground grid. Simultaneously, strong electromagnetic fields are generated both in the air and in the soil. The grid configuration shown in Figure 12 is considered, and the electromagnetic field is computed at points P₃ and P₄.

Figure 12. 40 × 40 m susbstation grid in a homogenous soil.

Figure 13 shows the transient current at four points “A”, “B”, “C” and “D” in a 40 × 40 m grid with 10 m long ground rods at the corners, subjected to current impulse injected at point “A”. As expected, a large proportion of the current dissipates to ground as it travels along the grid conductors, and only a very small proportion reaches point D, furthest from the injection point. Close agreement is obtained between the present approach and the NEC-4 Code computations at the feeding point, with slight differences between at points “B”, and “D”.

Figure 13. Transient current at different points of the grid.

The radiated electric field at point P₃ (in the soil) and P₄ (in the air) are shown in Figures 14 and 15 respectively. The magnitudes of the radiated fields are comparable with those given by the NEC-4 Code, although the wave shapes show some differences. For assessing the level of exposure of the personnel and the public, the magnitude is the most relevant parameter that is considered for a given impulse current. Figure 16 shows the temporal evolution of the y-component of magnetic field intensity H_y and total field intensity H_t, computed at points P₃ (in the soil) using the two computation methods. Reasonable agreement is obtained for the total magnetic field, whereas a time delay is observed between the two methods in the computed y-component of the magnetic field.

Figure 14. Variation of the electric field at point P₃.

Figure 15. Variation of the electric field at point P₄.

Figure 16. Variation of the magnetic field at point P₃.

5.1. Current injection at corner point A

Figures 17–20 show the temporal evolution of the electric and magnetic field intensity respectively on the earth’s surface when the lightning surge current is injected at corner point A. The induced field is maximal at the injection point, falling rapidly with distance away from the grid’s corner. The fall in magnetic field with distance is much sharper, reducing to zero withinn the first 10 meters from the corner point.

Figure 17. Y-component of electric field along profile 1.

Figure 18. Total electric field along profile 1.

Figure 19. Y-component of magnetic field along profile 1.

Figure 20. Total magnetic field magnitude along profile 1.

5.2. Current injection at center point B

Figures 21^{Figure 22Figure 23}–24 show the results obtained when the impulse current is injected at the middle point on the grid’s periphery (point B).

Figure 21. Y-component of electric field along profile 2.

Figure 22. Total electric field along profile 2.

Figure 23. Y-component of magnetic field along profile 2.

Figure 24. Total magnetic field along profile 2.

As can be seen from these figures, much higher electromagnetic fields are radiated near the centre injection point compared with the corner injection point. This result shows that, for given site conditions, the emitted electromagnetic fields due to lightning currents flowing in grounding systems depend not only on current magnitude, but also on the location of the point of lightning strike. Note that the field radiated by the grounding grid is lower than that radiated by a single electrode (Figure 3), which is expected since the flow of fault between the two is different. For the grid, the branching of current through the mesh causes a reduction in the induced electromagnetic field.

6. Conclusion

The fundamental concepts and analytical development for the quantification of the electromagnetic field radiated by a grounding system in semi-infinite dissipative medium subjected to a current excitation generated by the lightning strike in a conducting half-space are described and discussed. A method is proposed based on the computation by Hertzian dipole concept in frequency domain of the scattered electric and magnetic field in the vicinity of the grid. The solution was given in terms of the solutions of two basic Sommerfeld integrals. These are evaluated by numerical integration along contours in the complex plane and used subsequently to obtain the Sommerfeld integral values needed for the moment method solution of an integral equation. The latter requires knowledge of the current distribution. To this end, the integral equation is solved for the electric field by the method of moments in the frequency domain.

The proposed method was applied to a horizontal electrode and a grounding grid to address the problem of computing the transient electromagnetic radiation of a grounding system. The results are compared with those computed using the NEC-4 code with acceptable agreement.

The results demonstrate that the proposed method can be used as a tool to estimate the electromagnetic fields and associated EM energy radiated in the vicinity of grounding systems, and to assess safety thresholds for persons and equipment. It can also be used to determine the sub-surface distribution of electromagnetic fields, useful for estimating the impulse performance of grounding systems and the occurrence of soil ionisation phenomena in grounding system design studies.

This approach can be refined to improve accuracy. It can be further validated by performing further comparisons with other established software packages. It can also be applied to compute the radiated electromagnetic field of a wide range of grounding system configurations.

Funding

The authors received no direct funding for this research.

Additional information

Notes on contributors

N. Harid

N. Harid has been conducting research on grounding systems for the past 15 years, during which he and his colleagues at Cardiff University have secured funding of several million pounds to develop new measurement methods and modelling techniques for assessing the safety risks associated with short-circuit currents and lightning currents flowing in electrical substations’ grounding systems. His research interests include condition monitoring of overhead lines, demand management, application of wireless communications in power systems, high-voltage insulation and power system transients. The collaboration with co-authors has resulted in the work reported in this paper, which, like several other publications by the author on the same topic, relates very well to the ongoing issue of safety and electrocution risks for people living near or working in high-voltage electrical installations. The results of part of the previous research are proposed in a relevant international standard.