Thin Conductor Modelling Combined with a Hybrid Numerical Method to Evaluate the Transferred Potential from Isolated Grounding System

Grounding systems are essential parts of substations and power generation stations. The evaluation of transferred potentials from an active grounding system to other passive ones or to any near conductors is an important aspect to be considered, because transferred potentials may cause serious and fatal events. Moreover, it is an intrinsic issue of the Smart Grid where the ground systems of the power and ICT systems could be close to each other. Therefore, the estimation of the transferred potential is necessary at grounding system design stage for people safety and electric components safeguard. Numerical methods are the best choice to perform a truthful estimation, especially when large and complex grounding systems have to be designed. However, this task is complicated by the “unbounded” nature of the electromagnetic field and by the presence of components of extremely different size in the analysis domain. In this paper, an efficient hybrid finite element method is applied for the accurate and fast computation of transferred earth potentials from grounding systems. Moreover, the small dimensions of the components in the analysis domain are taken into account by the use of one-dimensional finite elements inserted in the tetrahedral mesh. It is worth mentioning the additional advantage of obtaining the electric potential on the earth surface without any post-processing operation.


Introduction
Grounding systems (GSs) are very important components of electric systems in substations and power plants, as well as in any civil, commercial, and industrial building, because they are essential for human safety [1] and equipment safeguard [2].
IEEE Standard 80 and IEC 479-1 establish the body current limits to be met when designing a GS [3].Such limits involve permissible touch-and-step voltages [4], which must be accurately evaluated.
The presence of passive GSs, over ground metallic elements (e.g., train tracks, fences), as well as underground ones (e.g., pipelines), in the vicinity of an active GS induces a modification of the electromagnetic field due to the active GS, as well as transferred potentials (TPs) that, in turn, may cause safety issues.The values of TPs must also be evaluated to verify the limit satisfaction outside the active GS in order to take suitable countermeasures.In this regard, the standard EN 50522 provides some qualitative indications [5].In this paper, TPs from isolated GSs are considered, while the TPs through conductors connected to the active GS [6,7] are not.
TPs from an active GS to a passive one or to metallic components in the proximity of the GS always occur, except for in a few cases [8].Very often, metal rails are present in large substations Energies 2019, 12, 1210 2 of 11 and power plants to facilitate the transportation of transformer equipment and fuel.These rails may transfer hazardous potentials around the switchyard area when a ground fault occurs [9].
Moreover, TPs may damage equipment in the neighborhood of windfarm turbines, which have increasing diffusion as green energy power generators [10].TPs may lead to undesired interfering effects in special power systems, such as hospital and data centers [11].
More in general, the study of the TPs is an intrinsic issue in smart grids due to the high integration of ICT and power systems [12].In urban areas, TPs in buried metallic parts, such as water and gas pipes, tramway tracks, and building foundations, are of particular concern because they involve worsening safety when there are floating conductors near an active GS [13].
During line maintenance, a de-energized line may be subject to uncontrollable transfer of the grid potential from other energized lines.Grounding the line at source end for personnel protection during maintenance involves greater hazard due to TP into this temporary GS from the active GS of other energized lines when they are faulted [14].
Commercial software tools usually used for GS analysis present limitations in evaluating the safety risk due to the transfer of potential out of a substation grounding system.Numerical methods based on the finite element method (FEM) are the best choice for the estimation of GS behavior, since they enable to account for more physical phenomena and to provide accurate results [15].In such a context, the effectiveness of hybrid numerical methods in GS analysis has already been proved [16,17].
The advantages that derive from the use of hybrid numerical methods for TP evaluation are presented in this paper.More specifically, the FEM-DBCI (Dirichlet boundary condition iteration) is effectively applied to TP evaluation.Furthermore, the real dimension and conductivity of the GS grid is considered by properly treating conductor slenderness, then the accuracy worsening due to some simplifications that are usually adopted in GS modelling is discussed.Finally, it is worth to note that, due to the possibility of decomposing the analysis domain in separated subdomains, the FEM-DBCI outperform other numerical methods, especially when the TP in large domains has to be evaluated.

Overview of Finite Element Analysis for the Evaluation of the Potential Transfer from Isolated Grounding Systems
The influence of buried metallic objects on the performance of power system grounds and the transferred earth potentials in power systems were initially performed by means of empirical methods.More accurate results have only been obtained by means of analytical methods implemented in computer programs since the 1970s [5,9,[18][19][20].
In this context, circuit method based on simplified formulas has been also adopted [21].On the other hand, these methods involve a high computational costs or unrealistic results when the section of conductors increases, as well as unpredictable margin of error [22,23].
The FEM provides better results for the computation of GS behavior, it enables to take into account for stratified soils, and it can consider different types and configurations of the GS, regardless of the size, length, type and shapes of the electrodes [12].An additional benefit of FEM is that the electric potential distribution on the earth surface is directly provided as part of the solution.
FEM has been used to compute the main quantities related to the GS by using the domain truncation combined with analytical expressions [24,25] or by means of spatial transformation techniques [26].An optimal size of the 3D finite elements near the GS conductors (modeled as 1D elements) is determined as a function of their radius [27].The distribution of the potential of a GS belonging to an L-shape building has been simulated by means of finite difference method (FDM) in [28].
FEM and FDM are useful for solving close-boundary value problems.When the problem exhibits open boundaries, as in the case of a GS analysis, these methods imply very poor accuracy to computational cost ratios since the discretization of a large domain is necessary in order to obtain acceptable results [29,30].These drawbacks are exacerbated when evaluating the TPs due to the very Energies 2019, 12, 1210 3 of 11 wide analysis domain to be considered.Moreover, the use of FDM is almost complex when the GS grid is irregular.
A numerical technique based on the transformation of the differential equations underlying the physical phenomena onto an equivalent boundary integral equation and the subsequent application of the boundary element method (BEM) has been proposed in References [31,32].The method of moments (MoM) has also been effectively adopted for obtaining a solution for the integral equations [33].
These approaches overcome the aforementioned limitations of empirical and analytical methods as well as the difficulties of FEM and FDM.On the other, the soil is considered to be homogeneous, or with limited inhomogeneity, and the size and resistivity of the conductors are neglected.The first limitation has been overcome by using a layered soil model in the application of BEM while maintaining the other simplifications [34,35].In the case of multiple layer models, the application of the "method of images" enables relevant improvement in convergence rate while managing the maximum admissible error [36], but the computation of the boundary integrals remains a very hard task due to the presence of singularity and of Green functions expressed by series.
Moreover, although the BEM well describes the GS geometry and is suitable to cope with the unbounded nature of the problem, it suffers from a cumbersome post-processing phase with high computational effort and limited accuracy.
The hybrid FEM-DBCI numerical method combines the benefits of the FEM and boundary integral methods.It has been used for evaluating the GS behavior [17] without considering TPs.Moreover, in Reference [17] it was presented as an accurate approach to properly consider both the conductivity and the size of each conductor of a GS.Such an approach has the very welcome feature of avoiding the thickening of the mesh of the regions around and inside the conductors.Such an advantage is obtained by modelling the conductors with 1D finite elements and merge the unknowns of these elements with those of the 3D ones adopted to discretize the soil.This approach is very useful in TP analysis, considering the presence of other conductors.For this reason, it has been combined in this work with the aforesaid hybrid numerical method.

The FEM-DBCI Method
In Figure 1, two GSs are shown buried in unbound soil.For the sake of simplicity, we assume that the soil is satisfactorily modelled as two-layer medium, in which the overlying layer has finite depth h and conductivity σ 1 whereas the underlying layer is an infinite half space of conductivity σ 0 (with σ 0 > σ 1 ).With the aim of studying the TP phenomenon, we consider that one GS (active) is energized at voltage V 0 , whereas the other GS (passive) is floating.
Energies 2018, 11, x FOR PEER REVIEW 5 of 11 until convergence takes place [38].FEM-DBCI has been implemented in ELFIN, a numerical code developed by the authors for electromagnetic research [43].

Numerical Results
This section firstly presents a simple example to verify the accuracy for numerical solution, then the computation of TP from an active GS to a near passive one, then the problem of TP from a GS to a metal rail has been treated.

A Single Rod Vertically Buried in a Homogeneous Soil
The first example is a very basic one and consists of a single rod vertically buried in a homogeneous soil.For this simple system, an analytical formula (9) for the grounding resistance is available: In order to apply the FEM-DBCI method for the computation of the electrical potential v and the current density J around the GSs, a truncation boundary Γ T is selected enclosing all the GSs and the non-homogeneous objects localized around them, if any.Since the two GSs are separated by a great distance, the boundary Γ T is constituted by two detached surfaces (see the red dashed lines in Figure 1).

The Laplace equation
holds in the bounded domain D, whose boundaries are Γ T and the soil surface Γ 0 .On Γ 0 , a homogeneous Neumann conditions ∂v/∂n=0 holds, whereas an unknown Dirichlet boundary condition is assumed to hold on Γ T , so that the FEM can be applied to the bounded domain D.
The soil inside D is discretized by means of tetrahedral finite elements, whereas the conductors of the GSs are modelled as 1D finite elements in addition to the 3D ones.A modified conductivity (in Sm) is attributed to these 1D elements, where σ c is the metal conductivity (in S/m), S c is its cross section area (in m 2 ), and η s is a numerical coefficient, which is set to 1, 1/2 or 1/4, according to the position of the conductor with respect to symmetry planes.For a 1D element, say E, the triple integral relative to its stiffness matrix is transformed into a line integral: where α m and α n are the shape functions of the m-th and n-th nodes, respectively, and L is the segment relative to the 1D finite element E. By introducing the following universal values [37]: where ξ is the local (non-dimensional) coordinate along the segment, the line integral in Equation ( 3) can be actually computed by means of such values as: In this way, along the conductors of the GSs, we can observe a variation of the potential.For the energized GS the Dirichlet condition v = V 0 is imposed on its point at the soil surface.
By applying the FEM to the bounded domain, the following algebraic system is obtained: where v is the array of the potential values at the nodes lying inside D and along Γ 0 , v T is the array of the potential values at the nodes lying on Γ T ; A, and A T , are sparse matrices of geometrical coefficients; b 0 is the known term array due to the assigned voltages.
The unknown Dirichlet condition on the truncation boundary is expressed by means of the integral equation [38]: ∂G(P, P ) ∂n − G(P, P ) ∂v(P ) ∂n dS (7) where Γ M is a closed surface (see the blue dotted line in Figure 1) enclosing all the conductors and non-homogeneities, but strictly enclosed by Γ T .The integration surface Γ M is conveniently selected homologous to Γ F and constituted by triangular faces of tetrahedral elements.Often, two layers of elements are inserted in between Γ M and Γ T .Other more suitable ways to choose the Γ M surface are described in Reference [39].In Equation ( 7) P and P are points lying on Γ T and Γ M , respectively, n' is the outward normal unit vector to Γ M and G(P, P ) is the two-layer soil Green's function [17].In this case, the use of the two-layer (or multi-layer) soil Green's function does not complicate the evaluation of the integral equation (Equation ( 7)) because the integration surface Γ M is separated from the truncation one Γ T, and no singularities arise.For this reason, in the presence of multi-layered soil, it is preferable to avoid the use of the BEM, even if it is combined with FEM [40], while the FEM-DBCI method is more suitable [41].
In numerical form, the integral equation (Equation ( 10)) reads: where F is a dense rectangular matrix having non null columns only for the nodes involved in the integral equation, that is, for the nodes of the tetrahedra external to Γ M and having a face lying on it.
The FEM-DBCI global algebraic system is constituted by Equations ( 6) and (8).This system can be solved by means of the generalized minimal residual (GMRES) solver, used as described in Reference [42].Alternatively, by taking advantage of the fact that the integral equation (Equation ( 8)) is explicit with respect to the array v T , the following iterative procedure may be used: (1) start from an arbitrary initial guess for v T ; (2) solve the FEM equation (Equation ( 6)) for v by means of the conjugate gradient (CG) solver; (3) use the DBCI equation (Equation ( 8)) to improve v T ; (4) iterate these steps until convergence takes place [38].FEM-DBCI has been implemented in ELFIN, a numerical code developed by the authors for electromagnetic research [43].

Numerical Results
This section firstly presents a simple example to verify the accuracy for numerical solution, then the computation of TP from an active GS to a near passive one, then the problem of TP from a GS to a metal rail has been treated.

A Single Rod Vertically Buried in a Homogeneous Soil
The first example is a very basic one and consists of a single rod vertically buried in a homogeneous soil.For this simple system, an analytical formula (9) for the grounding resistance is available: where L is the rod length, r its radius and σ the soil conductivity.Considering a rod with L = 10 m and r = 5 mm, buried in a soil with σ = 0.01 S/m, Equation ( 9) gives a value of grounding resistance R g = 12.71 Ω.
There are many factors that affect the solution accuracy when 3D computer aided design is adopted [44].In case of GS analysis, it is known that the accuracy of the finite element solution is strongly affected by the grid size in the vicinity of the rod [27].Hence, in order to verify the accuracy of the FEM-DBCI method with 1D elements, different numerical evaluations of the grounding resistance are performed, modifying the discretization of the rod and of its immediate neighborhood.The rod is discretized by means of 1D elements while the first-order tetrahedra have been adopted in the 3D domain.The parameter s = L/N s , where N s is the number of subdivisions of the rod, is varied from 0.1 m to 0.0125 m, doubling the number N s each simulation.The unbounded soil is truncated by means of a fictitious boundary Γ F placed at a distance d = 3 m from the rod.Symmetry is exploited by reducing the analysis domain to the octant x > 0, y > 0, and z < 0. Homogeneous Neumann conditions on the symmetry planes are assumed.The integration surface Γ M is selected at distance s from the rod.Table 1 shows the results, putting in evidence the correct behavior of the proposed numerical procedure: In fact, the error on the value of the computed grounding resistance decreases if a finer mesh is adopted in the proximity of the rod.Moreover, the use of 1D elements allows to obtain a very accurate value of the ground resistance with a low computational cost.Other simulations were carried Energies 2019, 12, 1210 6 of 11 out moving away the fictitious boundary but no significant difference was observed in the results, besides an expected growth of the computing time.

Transferred Potential from an Active to a Passive Grounding System
Figure 2 shows the geometry of the active GS where there are two equal squares, with dimension 1 m × 1 m, buried parallel with the soil surface at 0.6 m in depth.They present a rod of length 0.6 m that connects them with the soil surface and four other vertical rods, 1 m in length, at its corners.They are steel (σ c = 7.0•10 6 S/m) conductors and radius, r c , of their circular cross section is equal to 1 cm.The voltage of the active GS is set to V 0 = 100 V, hence the voltage across the passive GS may be considered directly as the percentage of TP to such a GS since the problem under study is linear.

Transferred Potential from an Active to a Passive Grounding System
Figure 2 shows the geometry of the active GS where there are two equal squares, with dimension 1 m x 1 m, buried parallel with the soil surface at 0.6 m in depth.They present a rod of length 0.6 m that connects them with the soil surface and four other vertical rods, 1 m in length, at its corners.They are steel (c = 7.010 6 S/m) conductors and the radius, rc, of their circular cross section is equal to 1 cm.The voltage of the active GS is set to V0 = 100 V, hence the voltage across the passive GS may be considered directly as the percentage of TP to such a GS since the problem under study is linear.
The intersection between the ground surface and the rod connecting the GS with such a surface is adopted as origin of the reference coordinate system used for the analysis.The z axis has been chosen orthogonal with the soil surface while the y axis has been chosen parallel with the common side between the two squares.
A cross with dimension 1 m  1 m buried at 0.6 m in depth and belonging to a plane parallel to the ground surface constitutes the passive GS.It is connected to the ground surface through a vertical rod of 0.6 m placed at the cross intersection.The arms of the cross are aligned along the x-and y-axes, far away 15 m from the center of the first system in the x direction.Figure 3 shows the passive GS Both xz and xy coordinate planes are symmetry planes for the active and passive GSs according to the chosen reference coordinate system and GSs position.The intersection between the ground surface and the rod connecting the GS with such a surface is adopted as origin of the reference coordinate system used for the analysis.The z axis has been chosen orthogonal with the soil surface while the y axis has been chosen parallel with the common side between the two squares.
A cross with dimension 1 m × 1 m buried at 0.6 m in depth and belonging to a plane parallel to the ground surface constitutes the passive GS.It is connected to the ground surface through a vertical rod of 0.6 m placed at the cross intersection.The arms of the cross are aligned along the x-and y-axes, far away 15 m from the center of the first system in the x direction.Figure 3  The soil is modelled as an unbounded medium with two strays: the upper one with depth h = 20 cm and conductivity 1 = 3.3310 −4 S/m while the bottom has conductivity 0 = 0.05 S/m.
The symmetry in the analysis domain enable the consideration of only a quarter, thus strongly reducing the computational effort necessary to perform the simulation.More specifically, the considered analysis sub-domain is the portion with y > 0, z < 0. The reduction is obtained by setting homogeneous Neumann conditions on the symmetry planes y = 0 and z = 0.The unbounded soil in such sub-domain is truncated by means of a fictitious boundary F constituted of two parallelepipedal boxes of sizes 2.8 m  1.4 m  2.0 m and 1.8 m  0.9 m  1.0 m, respectively for the active and passive GS.The integration surface M is selected homologous to F, at a distance of 0.2 m.
Second-order tetrahedra have been adopted to discretise the two detached octants, that is, the upper and lower part of the sub-domain.More specifically, 23,660 tetrahedra have been considered for the former while 56,160 for the latter.The analysis converges in four iterations being the stopping criteria a tolerance of 10 −4 per cent for the CG solver and 10 −2 per cent for the GMRES.Figure 4 shows the contours of the electrical potential in the xz symmetry plane and on F around the active GS.The floating potential at the centre of the cross of the passive GS is about VF= 4.32 V.In other words, the TP to the passive GS is the 4.32% of the voltage appearing at the ground surface of the active GS.

Transferred Potential from an Active to a Metal Rail
The active GS of Figure 3 is considered while the passive system is a massive steel rail with a length of 4 m aligned along the y-axis with cross section of 10 cm  20 cm.The bar is buried in the soil being the upper side aligned with the soil surface.The middle point of the bar is far away 15 m from the center of the GS in the x direction.The soil is modelled as an upper layer of depth h = 20 cm and 1 = 3.3310 −4 S/m and a below unbounded medium of conductivity 0 = 0.05 S/m.The two detached subdomains are discretized by means of tetrahedral finite elements of the second order, for a total 52301 and 81154 nodes.In Figure 5, the contours of the potential are drawn on the soil surface around the bar.The mean value of the electric potential of the bar is about 5.7 V. Once again, this means the TP at the bar is 5.7% of the active GS.Both xz and xy coordinate planes are symmetry planes for active and passive GSs according to the chosen reference coordinate system and GSs position.
The soil is modelled as an unbounded medium with two strays: the upper one with depth h = 20 cm and conductivity σ 1 = 3.33•10 −4 S/m while the bottom has conductivity σ 0 = 0.05 S/m.
The symmetry in the analysis domain enable the consideration of only a quarter, thus strongly reducing the computational effort necessary to perform the simulation.More specifically, the considered analysis sub-domain is the portion with y > 0, z < 0. The reduction is obtained by setting homogeneous Neumann conditions on the symmetry planes y = 0 and z = 0.The unbounded soil in such sub-domain is truncated by means of a fictitious boundary Γ F constituted of two parallelepipedal boxes of sizes 2.8 m × 1.4 m × 2.0 m and 1.8 m × 0.9 m × 1.0 m, respectively for the active and passive GS.The integration surface Γ M is selected homologous to Γ F , at a distance of 0.2 m.
Second-order tetrahedra have been adopted to discretise the two detached octants, that is, the upper and lower part of the sub-domain.More specifically, 23,660 tetrahedra have been considered for the former while 56,160 for the latter.The analysis converges in four iterations being the stopping criteria a tolerance of 10 −4 per cent for the CG solver and 10 −2 per cent for the GMRES.Figure 4 shows the contours of the electrical potential in the xz symmetry plane and on Γ F around the active GS.The floating potential at the centre of the cross of the passive GS is about V F = 4.32 V.In other words, the TP to the passive GS is the 4.32% of the voltage appearing at the ground surface of the active GS.The active GS of Figure 3 is considered while the passive system is a massive steel rail with a length of 4 m aligned along the y-axis with cross section of 10 cm × 20 cm.The bar is buried in the soil being the upper side aligned with the soil surface.The middle point of the bar is far away 15 m from the center of the GS in the x direction.The soil is modelled as an upper layer of depth h = 20 cm and σ 1 = 3.33•10 −4 S/m and a below unbounded medium of conductivity σ 0 = 0.05 S/m.The two detached subdomains are discretized by means of tetrahedral finite elements of the second order, for a total 52301 and 81154 nodes.In Figure 5, the contours of the potential are drawn on the soil surface around the bar.The mean value of the electric potential of the bar is about 5.7 V. Once again, this means the TP at the bar is 5.7% of the active GS.

Conclusions
In this paper, the hybrid method FEM-DBCI has been applied to evaluate the transferred potentials from an active grounding system.Firstly, the advantages that derive from the use of such hybrid numerical methods for grounding system have been recalled and then have been shown in case of known grounding resistance.Then, two numerical examples have been solved and the results have highlighted the accuracy of the solution provided by the method.
Due to the hybrid nature of the FEM-DBCI, the possibility of decomposing the analysis domain in separated subdomains has been exploited, avoiding the discretization of the large volume of soil between the active and the floating conductors.Moreover, the implementation of 1D elements to model the grounding grid and the rods, that had a very small section compared with the analysis

Conclusions
In this paper, the hybrid method FEM-DBCI has been applied to evaluate the transferred potentials from an active grounding system.Firstly, the advantages that derive from the use of such hybrid numerical methods for grounding system have been recalled and then have been shown in case of known grounding resistance.Then, two numerical examples have been solved and the results have highlighted the accuracy of the solution provided by the method.
Due to the hybrid nature of the FEM-DBCI, the possibility of decomposing the analysis domain in separated subdomains has been exploited, avoiding the discretization of the large volume of soil between the active and the floating conductors.Moreover, the implementation of 1D elements to model the grounding grid and the rods, that had a very small section compared with the analysis domain, hugely decreased the number of tetrahedra in the mesh.In this way, applying FEM-DBCI dramatically reduced the computing time.

Figure 1 .
Figure 1.An active GS in the vicinity of another passive GS embedded in a two-layer soil.

Figure 1 .
Figure 1.An active GS in the vicinity of another passive GS embedded in a two-layer soil.

Figure 2 .
Figure 2. Three-dimensional geometry of the active GS.Figure 2. Three-dimensional geometry of the active GS.

Figure 2 .
Figure 2. Three-dimensional geometry of the active GS.Figure 2. Three-dimensional geometry of the active GS.

Figure 3 .
Figure 3. Three-dimensional geometry of the passive GS.

Figure 3 .
Figure 3. Three-dimensional geometry of the passive GS.

Energies 2018 , 11 Figure 4 .
Figure 4. Contours of the potential around the active GS.

Figure 4 .
Figure 4. Contours of the potential around the active GS.

Figure 5 .
Figure 5. Contours of the potential around the metal rail.

Figure 5 .
Figure 5. Contours of the potential around the metal rail.

Table 1 .
Grounding resistance R g and computing effort in terms of CPU times and unknowns vs. different values of parameter s.

Table 1 .
Grounding resistance Rg and computing effort in terms of CPU times and unknowns vs. different values of parameter s.