Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 31, N o 3, September 2018, pp. 487-500 https://doi.org/10.2298/FUEE1803487A ELECTROMAGNETIC ANALYSIS OF SINGLE/MULTIPLE GROUNDING RODS * Vesna Arnautovski-Toseva, Leonid Grcev University SS Cyril and Methodius University, Faculty of Electrical Engineering and Information Technologies, Skopje, Macedonia Abstract. This paper presents electromagnetic modeling of multiple driven grounding rods in homogeneous/two-layer soil. The mathematical model is formulated by mixed potential integral equation (MPIE) on the basis of Sommerfield integrals. Several configurations of multiple driven rods located in homogeneous or two-layer soil are analyzed. The authors are focused on the calculation of the current density along the rods in wide frequency range from 100Hz to 1MHz. Key words: electromagnetic theory, grounding rod, high frequencies, homogeneous soil, two-layer soil. 1. INTRODUCTION Ground rods are simplest and the most often used means used for earth termination of different electrical systems, providing a conducting connection, whether intentional or accidental between an electrical circuit or equipment and the earth. Their behavior at DC (50 or 60 Hz) is well understood [3-4], but their high-frequency (HF) and transient performance is also of interest in different fields, such as, lightning protection, power and telecommunication systems, power system transients, electromagnetic compatibility, etc. The safety criteria based on “a minimum rise in the potential” are taken from the power systems analysis. In such cases the usual DC approximation leads to rather straightforward computations. The simulation studies show that grounding systems, even the simplest ones, behave quite differently at low and high frequencies [5]. The survey of the literature shows that high frequency analysis of grounding systems is realized by using lumped circuit equivalents [6-7], quasi-static method of images [8-9], rigorous electromagnetic model  Received February 19, 2018; received in revised form May 24, 2018 Corresponding author: Vesna Arnautovski-Toseva Faculty of Electrical Engineering and Information Technologies, Skopje, Macedonia (E-mail: atvesna@feit.ukim.edu.mk) * An earlier version of this paper was presented at the 13th International Conference on Applied Electromagnetics (ПЕС 2017), August 31 - September 01, 2017, in Niš, Serbia [1] 488 V. ARNAUTOVSKI-TOSEVA, L. GRCEV [10-11], or hybrid approaches [12-13]. However, there is a lack of papers that threat the problem of multiple driven rods at high frequencies, except at DC [3]. Our objective in this paper is to give sight into the problem of high frequency performance of multiple driven ground rods in homogeneous or two-layer soil with respect to the case of a single rod. The main interest is the current density along the rods of various configurations for which the DC behavior is known [4]. Preliminary results of authors’ research in this field are presented in [1]. The authors have recently presented similar analysis of a single rod at high frequencies [2]. 2. ELECTROMAGNETIC MODEL The rigorous treatment of the air/two-layer soil interfaces in electromagnetic models is based on the exact solution for the field of a Hertz dipole near a conducting half space/two-layer soil. This approach involves Green’s functions formulated by Sommefeld integrals that need numerical integration. This model is confirmed as theoretically most accurate since it is based on minimum approximations. The detailed description of the mathematical model is given in our previous work [10, 11, 15]. In this analysis the electromagnetic model is extended to take into account multiple rods geometry and the corresponding excitations. 2.1. Geometry of the problem In Fig. 1 we consider grounding system consisting of k identical parallel rods, each of length L and radius a penetrating homogeneous/two-layer soil. The upper layer (medium 1) is of finite depth d characterized by permittivity ε1, permeability μ0 and resistivity 1. When the soil is homogeneous it is assumed that d . In the case of two-layer soil, the bottom layer (medium 2) is characterized by permittivity ε2, permeability μ0 and resistivity 2. The air (medium 0) is characterized by permittivity ε0 and permeability μ0. Corresponding rod(s) lengths in the upper/bottom layer are L1 and L2 respectively. In the case of homogeneous soil L2=0. Fig. 1 Geometry of multiple rods in two-layer soil Electromagnetic Analysis of Single/Multiple Grounding Rods 489 2.2. Mathematical model Grounding rod may be considered in a circuit with an ideal harmonic current source of magnitude IS with one terminal connected to the ground electrodes and the other terminal to the remote earth theoretically at infinite distance. The influence of the connecting leads is ignored. Multiple rods excitation is assumed by using respectively k harmonic current sources of equal magnitude IS, leading to total excitation current of kIS. Fig. 2 Approximation of the current with triangle dipoles Following thin-wire approximation, the physical model of the system of ground rods is based on fictitious segmentation into n+1 straight tubular segments, Fig. 2. The segmentation is done in a way so that no segment penetrates through the boundary between the two soil layers. To solve current distribution along each rod the method of moments is applied by using thin wire approximation [10]. Following Galerkin formulation the current distribution in the system of rods is approximated by n overlapped triangular dipoles, each extended over two neighboring segments i = li-1 + li (i=2, 3, ... , n+1). One of the triangular dipoles passes through the interface between the two soil layers with its top point located just at the interface between the two soil layers (z = d). The total excitation current kIS is approximated by k additional triangular monopoles sk with length k = lk, each of magnitude IS, that are positioned at the top of each rod. The weighting functions are also triangular dipoles. The following matrix equation yields the current distribution along the grounding rods 1 2 1 2 1 2 1 1 1 11 12 1 1 2 2 221 22 2 2 1 2 k k k s S s S s S n s S S S S Sn n n nn n ns S ns S ns S z I z I z Iz z z I z I z I z Iz z z I z z z I z I z I z I                                  , (1) where matrix column [I] contains the coefficients In of unknown currents; [Z] is generalized impedance matrix related to self and mutual impedances between all triangular dipoles that represent all electromagnetic influences between all rods in the 490 V. ARNAUTOVSKI-TOSEVA, L. GRCEV configuration; [ZSIS] is excitation matrix where the corresponding multiple rods excitations and their influences are taken into account. Once the currents in the dipoles are computed, the current density along each rod in the configuration is easily estimated 1i i i i I I I l     . (2) The elements zij and the elements kjs z correspond to mutual impedances between the source dipole i and observation dipole j (i, j=1, 2, ... , n); and respectively between the source dipole i (i=1, 2, ... , n) and each of the k excitation monopoles sk as following 1 1 k k k ij ij zj i i j is is zs i i k u z E dz I I u z E dz I I         . (3) Here Ezj is tangential electric field at the surface of the observation dipole j over the length of a dipole segment lj due to current Ii in the dipole i. In this analysis, MPIE formulation for the electric field is used z z z E j A V   , (4) 1 i z Azz i z Vz i i dI A G I dz V G q dz q j dz       . (5) In (5) GAzz is z-component of the dyadic Green's function for the magnetic vector potential at observation point (x,y,z) due to a Hertzian vertical electric dipole (VED) of unit strength at source point (x',y',z'). Respectively, GVz is the corresponding scalar potential Green's function due to a single point charge associated with VED [14]. The exact formulation of spatial domain Green’s functions are formulated by Sommerfeld integrals that are solved by direct numerical integration  , , 0 0 , 0 1 1 ( ) ( ) 2 2 mn mn mn Azz Vz Azz Vz Azz Vz G G k J k k dk S G           . (6) where J0(k) is zero-order Bessel function of the first kind. The corresponding spectral domain Green’s functions for the magnetic vector potential are given below [15] Electromagnetic Analysis of Single/Multiple Grounding Rods 491 1 1 1 2 1 1 2 1 1 2 ( ) ( )11 0 1 ( ) ( ) ( )12 0 12 10 1 ( ) ( ) ( )21 0 21 10 2 22 0 2 e 2 e e e 2 e 2 e 2 z z z z z z z z z z jk z z jk z z jk z z Azz z jk z d jk d z jk d z Azz z jk z d jk d z jk d z Azz z jk Azz z G A e Be j k G T M R j k G MT e R e j k G j k                                         21 ( 2 )2 12 10 ( ) zz z z jk z z djk d M R R e e          . (7) 1 1 1 1 1 1 1 ( ) ( ) ( ) 12 10 ( 2 ) 10 12 1 2 12 10 ( ) ( ) 1 jkz d z jkz d z jkz d z jkz z jkz z jkz d z jk d A e R e R e M B e R e R e M M R R e                         . (8) In above relations, R21, R10, T21 and T10 are reflection and transmission coefficients of a TM wave incident on both interfaces between mediums 11 2 0 1 02 1 21 10 12 21 1 2 12 1 0 1 0 1 2 21 12 21 1 22 1 2 2 2 2 2 0 0 0 0 2 1 ( ) , 1, 2 z zz z z z z z z z z i zi i i i k kk k R R R R k k k k k T T T k k j k k k k k i                                  . (9) The corresponding spectral expressions for the scalar potential Green's functions are derived by using the following relation 2 2 1 mn mn Azz V zm G G z zk     . (10) 3. NUMERICAL RESULTS In this section frequency domain behavior of single/multiple grounding rods located in homogeneous/two-layer soil is analyzed. Three grounding test configurations are assumed as shown in Fig. 3; single rod (R1), 5-rods configuration (R5), and 9-rods configuration (R9). The geometry of each rod is identical, characterized by length 10m (extending from 0.05m to 10.05m) and radius 0.01m. The outer dimensions of the “mesh” of R5 and R9 configurations are: a) 20m20m, (R5a and R9a); and b) 10m10m (R5b and R9b) as given in [3]. The main objective of this analysis is to compare the behavior of the given grounding multiple rod structures at high frequencies with respect to the corresponding DC results. 492 V. ARNAUTOVSKI-TOSEVA, L. GRCEV Fig. 3 Test configurations: Single rod (R1) and multiple rods (R5) and (R9) In the case of two-layer soil the permittivity of both layers is 1 = 2 = 100. The resistivity of the upper layer with depth d=5m is fixed at 1=100Ωm, while the bottom layer resistivity is: 2=33.33Ωm (reflection factor K=0.5), and 2=300Ωm (reflection factor K=+0.5). The homogeneous soil corresponds to reflection factor K=0. The total excitation current in the corresponding analysis is k1kA, i.e. the excitation current applied to each rod is 1kA. 3.1. Homogeneous soil In this section, all test configurations are analyzed in homogeneous soil. The main objective is to investigate the influence of the number of multiple rods in R5 and R9 configurations, as well as how their mutual distance affects the current distribution. As a result, current density along the central rod and outer/corner rods is compared to the corresponding results obtained for single rod R1. The analysis is preformed in frequency range from 100Hz to 1MHz. In Fig. 4 it may be observed respectively the current density along the central and the outer/central rods of: a) R5a and R5b configurations, and b) R9a and R9b configurations at 100Hz with respect to single rod R1 behavior. As it may be observed, the current density along R1 is generally uniform except at the bottom end of the rod where much higher current density is observed. The results show that in multiple rods configurations (R5 and R9) the current density in the central rod is lower in the upper part, and higher in the bottom part of the rod as compared to the outer/corner rods. This effect is emphasized with smaller distance between the rods and with larger number of rods in the grounding configuration (R9b). The results are in good accordance with the reference results obtained at DC [3]. In Fig. 5 a) and b) it may be observed current density in the specific rods of R5 and R9 configurations respectively at 100kHz. Again, the injected current is almost uniformly discharged along each rod, i.e. the grounding system performance is quasi-static up to 100kHz. The differences in the current density in outer/corner rods are lower as compared to single rod. Next, in Fig. 6 it is shown the current density along the central and outer/corner rods of R5 and R9 configurations obtained at 1MHz with respect to single rod R1. The results show that each rod of R5 and R9 configurations acts as isolated rod, since the current density along each rod is identical and equal to the corresponding single rod behavior. As it may be observed, most of the injected current is discharged from the upper part of the rod. The results show that the distance between the rods has reduced influence on the current density at higher frequencies. Summarizing the results obtained for the analyzed configurations it may be expected that in the case when the distance between the rods is smaller higher differences in the current single rod R1 central corner side R9 central corner R5 20 (10) m 2 0 ( 1 0 ) m Electromagnetic Analysis of Single/Multiple Grounding Rods 493 densities would arise between the central rod and the outer rods. Also it may be assumed that such differences would decrease at higher frequencies it can be assumed that the differences will be reduced so that the results will converge with those given in Fig. 6. 0 2 4 6 8 10 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 Rod length (m) C u rr e n t d e n s it y ( k A /m ) a) Current density along R1, R5a and R5b at 100Hz R1 central rod of R5a corner rod of R5a central rod of R5b corner rod of R5b 0 2 4 6 8 10 0.08 0.1 0.12 0.14 0.16 0.18 0.2 b) Current density along R1, R9a and R9b at 100Hz Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b Fig. 4 Current densities along R1, R5 and R9 in homogeneous soil at 100Hz 494 V. ARNAUTOVSKI-TOSEVA, L. GRCEV 0 2 4 6 8 10 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 a) Current density along R1, R5a and R5b at 100kHz Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R5a corner rod of R5a central rod of R5b corner rod of R5b 0 2 4 6 8 10 0.08 0.09 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 b) Current density along R1, R9a and R9b at 100kHz Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b Fig. 5 Current densities along R1, R5 and R9 in homogeneous soil at 100kHz Electromagnetic Analysis of Single/Multiple Grounding Rods 495 0 2 4 6 8 10 0.05 0.1 0.15 0.2 0.25 0.3 Current density along R1, R5 and R9 at 1MHz for K=0 Rod length (m) C u rr e n t d e n s it y ( A /m ) R1 central rod of R5b corner rod of R5b central rod of R9b side rod of R9b corner rod of R9b Fig. 6 Current densities along R1, R5 and R9 in homogeneous soil at 1MHz 3.2. Two-layer soil In this section, the analysis is focused on R9 configuration since highest variations in the current density are observed in the central rod with respect to the outer/corner rods or single rod case, as observed in previous section. As may be seen in Fig.7 a) K=+0.5 and b) K=0.5 respectively, significant differences in the current density along ground rods are observed at 100Hz. The current density is much higher in the partition of the rod located in the lower resistivity layer. Also, current density is practically uniform along the rod partition located in one layer. This result is in accordance with the rod DC behavior, when the injected current is practically uniformly discharged in the surrounding soil [3]. However, the differences in current density between the upper and the bottom partition of the rod lead to large jump that occur at the interface between both soil layers (at depth d). When K=+0.5, the current density along the central rod is lower along the upper rod partition, but much higher at the bottom rod partition as compared to the outer/single rod. However for K=0.5, such differences in the current density are much less observed. In Fig. 8, the corresponding results obtained at 100kHz are shown. The results are generally similar with those shown in Fig. 7. However, higher deviations, expressed by small peaks, lead to larger discontinuity in the current density that occurs at the interface between both soil layers. 496 V. ARNAUTOVSKI-TOSEVA, L. GRCEV 0 2 4 6 8 10 0.04 0.06 0.08 0.1 0.12 0.14 0.16 a) Current density along R1, R9a and R9b at 100Hz for K=+0.5 Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b 0 2 4 6 8 10 0.05 0.1 0.15 0.2 0.25 b) Current density along R1, R9a and R9b at 100Hz for K=-0.5 Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b Fig. 7 Current densities along R1, R9a and R9b in two-layer soil at 100Hz Electromagnetic Analysis of Single/Multiple Grounding Rods 497 0 2 4 6 8 10 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 a) Current density along R1, R9a and R9b at 100kHz for K=+0.5 Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b 0 2 4 6 8 10 0 0.05 0.1 0.15 0.2 0.25 b) Current density along R1, R9a and R9b at 100kHz for K=-0.5 Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b Fig. 8 Current densities along R1, R9a and R9b in two-layer soil at 100kHz 498 V. ARNAUTOVSKI-TOSEVA, L. GRCEV As may be observed in Fig. 9 and in Fig. 10, the current density obtained at high frequencies, 1MHz and 10MHz respectively, differs significantly from the corresponding low frequency, quasi-static, behavior. 0 2 4 6 8 10 0 0.05 0.1 0.15 0.2 0.25 a) Current density along R1, R9a and R9b at 1MHz for K=+0.5 Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b 0 2 4 6 8 10 0.05 0.1 0.15 0.2 0.25 0.3 b) Current density along R1, R9a and R9b at 1MHz for K=-0.5 Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 central rod of R9a side rod of R9a corner rod of R9a central rod of R9b side rod of R9b corner rod of R9b Fig. 9 Current densities along R1, R9a and R9b in two-layer soil at 1MHz Electromagnetic Analysis of Single/Multiple Grounding Rods 499 0 2 4 6 8 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Current density along R1 and R9 at 10 MHz Rod length (m) C u rr e n t d e n s it y ( k A /m ) R1 K=0 R1 K=+0.5 R1 K=-0.5 central rod of R9b K=0 central rod of R9b K=+0.5 central rod of R9b K=-0.5 side rod of R9b K=0 side rod of R9b K=+0.5 side rod of R9b K=-0.5 corner rod of R9b K=0 corner rod of R9b K=+0.5 corner rod of R9b K=-0.5 Fig. 10 Current densities along R1, R9a and R9b in two-layer soil at 10MHz At high frequencies the influence of the number of multiple rods and their mutual distance is practically negligible. The differences in the current density are mainly due to the two-layer soil parameters. At 1MHz, the jump in the current density that occurs between the upper and the bottom rod partitions is especially high when the bottom layer is less resistive. In this case, due to high resistivity of the upper layer, most of the injected current is discharged from the bottom part of the rod surrounded by less resistive bottom layer. However, as the frequency increases this effect vanishes. As may be seen in Fig. 10, the results obtained at 10MHz for all rods in R9b configuration converge to the corresponding behavior of a single rod (R1). At high frequencies above 1MHz, practically, the behavior of the grounding multiple rods is not affected by the number of rods in the configuration or by the distance between the rods, since each rod in the configuration acts as isolated rod in homogeneous soil. Here also, the bottom layer and the corresponding reflection factor has almost no influence on the current density since all injected current is quickly discharged from the upper part of each rod into the upper soil layer. It may be expected that when the distance between the rods is smaller, higher differences in the current densities would arise at low frequencies together with higher discontinuities at the interface between the two soil layers. Again, it may be expected that such differences would decrease at higher frequencies leading to results similar as given in Fig. 10. 500 V. ARNAUTOVSKI-TOSEVA, L. GRCEV 4. CONCLUSION In this paper, the frequency domain behavior of single/multiple grounding rods in homogeneous/two-layer soil is analyzed. The results for the current density obtained for various multiple rods configurations show that grounding rods performance at high frequencies differs significantly from their low frequency performance. In case of homogeneous soil, the current in the central rod varies from the corresponding distribution in the outer/corner rods. At higher frequencies, this effect vanishes, i.e. at 1MHz the current density for all test cases converge to the distribution obtained for single rod. In case of two- layer soil, current density along the upper and the bottom rod partitions differs significantly and lead to large jump in the current distribution that occurs at the interface between the both soil layers. This effect is noticeable at higher frequencies also, while the influence of other parameters such as the number of rods and their distances are negligible. At very high frequencies as 10MHz, the effect of the two distinct soil layers also vanishes. The results obtained for rods all test configurations show that their performance is identical to the corresponding case of single rod in homogeneous soil. In the future work the authors will analyze in more details the influence of the multiple rods geometry on the current density and will extend their analysis to h grounding impedance. REFERENCES [1] V. Arnautovski-Toseva, L. Grcev and S. Cundeva, "High Frequency Behaviour of Ground Rods", In Proc. of the 13th International Conference on Applied Electromagnetics, Nis, Serbia, August 30 - Sep. 01, 2017, pp. 1–4. [2] V. Arnautovski-Toseva, L. Grcev and K. El Khamlichi Drissi, "High Frequency Performance of Ground Rod", In Proc. of the 17th IEEE International Conference on Smart Technologies, Ohrid, Macedonia, 6–8 July 2017, pp. 914–918. [3] F. Dawalibi, D. Mukhedekar, "Influence of Ground Rods on Grounding Grids", IEEE Trans. Power App. Syst, vol. 98, pp. 2089–2097, 1979. [4] J. M. Nahman, "Digital Calculation of Earthing Systems in Nonuniform Soil", Arch. Elektrotech., vol. 62, pp. 19–24, 1980. [5] R. G. Olsen, M. C. Willis, "A Comparison of Exact and Quasi-Static Methods for Evaluating Grounding Systems at High Frequencies", IEEE Trans. Power Del., vol. 11, pp. 1071–1081, 1996. [6] C. T. Mata, M. I. Fernandez, V. A. Rakov, M. A. Uman, "EMTP Modeling of a Triggered-Lightning Strike to the Phase Conductor of an Overhead Distribution Line", IEEE Trans. Power Del., vol. 15, pp. 1175–1181, 2000. [7] L. Grcev, M. Popov, "On High Frequency Circuit Equivalents of a Vertical Ground Rod", IEEE Trans. Power Del., pp. 15981603, 2005. [8] T. Takashima, T. Nakae, R. Ishibashi, "High Frequency Characteristics of Impedances to Ground and Field Distributions of Ground Electrodes", IEEE Trans. Power App. Syst, vol. 100, pp. 18931900, 1980. [9] S. Bourg, B. Sacepe et al., "Deep Earth Electrodes in Highly Resistive Ground: Frequency Behaviour”, In Proc. of the IEEE Int. Symp. on Electromagnetic Compatibility, Atlanta GA, USA, 14-18 Aug. 1995, pp. 584 –589. [10] L. Grcev, F. Dawalibi, "An Electromagnetic Model for Transients in Grounding Systems", IEEE Trans. Power Del., No. 4, pp. 17731781, 1990. [11] L. Grcev, V. Arnautovski-Toseva, "Grounding Systems Modeling for High Frequencies and Transients: Some Fundamental Considerations", In Proc. IEEE Bologna Power Tech. Bologna, Italy, June 23-26, 2003. [12] A. F. Otero, J. Cidras, J. L. del Alamo, "Frequency-dependent Grounding System Calculation by Means of a Conventional Nodal Analysis TechNique", IEEE Trans. Power Del., vol. 14, pp. 873–878, 1999. [13] Z-X Li, W. Chen, J-B Fan, J. Lu, "A Novel Mathematical Modeling of Grounding System Buried in Multilayer Earth", IEEE Trans. Power Del., vol. 21, pp.1267–1272, 2006. [14] G. Dural, M. I. Aksun, "Closed-form Green's Functions for General Sources and Stratified Media", IEEE Trans. Microwave Theory Techn., vol. 43, July 1995, pp. 1545-1552. [15] V. Arnautovski-Toseva, L. Grcev, "Image and Exact Models of a Vertical Wire Penetrating a Two- Layered Earth", IEEE Trans. on Electromagnetic Compatibility, vol. 9, 2011, pp. 1–9.