Elektrotehniški vestnik 84(1-2): 39-47, 2017 Original scientific paper
A combined boundary element and an analytical approach to grounding mesh modeling in a multi-layer soil
Adnan Mujezinovic1, Ajdin Mulaosmanovic2, Alija Muharemovic1, Irfan Turkovic1, Zijad Bajramovic1
'Faculty of Electrical Engineering, University of Sarajevo, Department of Power Engineering, Bosnia and Herzegovina 2Public Electric Utility Elektroprivreda of Bosnia and Herzegovina, Hydropower Plants on Neretva E-mail: adnan.mujezinovic@etf.unsa.ba
Abstract. This paper deals with a combined computer-aided numerical - analytical approach to calculation of the grounding. The electromagnetic theory on which the presented mathematical model is based is described. The model which is based on boundary element method, Aitken's S2 algorithm and grounding potential non-uniformity correction factors is explained in detail. A special attention is paid to the selection of an appropriate Green's function in the calculation of parameters of a large and complex grounding mesh for a homogeneous and a two-layer soil. Finally, the model is used to calculate the grounding mesh parameters of a real power substation with complex geometry.
Keywords: grounding mesh, potential distribution, boundary element method (BEM), Aitken's S2 algorithm, non-uniformity correction factor.
Združeni numerično-analitični pristop za modeliranje
ozemljitvenih mrež
V prispevku obravnavamo združen numerično-analitičen pristop za izračun parametrov ozemljitvenih mrež z uporabo računalnika. Najprej je predstavljena elektromagnetna teorija, na kateri temelji predlagani matematični model. V nadaljevanju je podrobneje opisan matematični model za izračun parametrov ozemljitvenih mrež, ki temelji na Aitkenovem 52 algoritmu z metodo končnih elementov in upoštevanju korekcijskih faktorjev ozemljitvenega potenciala zaradi nehomogenosti ozemljitve. Posebna pozornost je namenjena izboru ustrezne Greenove funkcije pri izračunu parametrov velikih in kompleksnih ozemljitvenih mrež pri homogeni in dvoslojni ozemljitvi. Predlagani matematični model smo uporabili za izračun parametrov ozemljitvenih mrež obstoječe transformatorske postaje s kompleksno geometrijo.
1 Introduction
Thanking into account safety, the most important element of power substations are the grounding systems. The primary role of a grounding system is to provide safety of personnel and integrity of equipment during faults [1]. This important function of the grounding systems is performed by conducting the fault current into the surrounding soil in which they are placed. In order to properly perform their function, the grounding systems should have a low resistance, thus limiting the potential values at the ground surface during the highest values of the fault currents and to keep a large enough
value of fault current for the protection devices in power substations to react. To satisfy this, the grounding systems are designed in complex geometries consisting of a large number of horizontal, vertical and inclined uninsulated galvanically coupled conductors [2,3]. For large power substations, complex grounding meshes are mostly used as a grounding system [4].
To model grounding meshes, many authors have developed different analytical and numerical approaches. Nowadays analytical expressions are poorly represented in the literature because they are not able to give accurate results for most of the real cases. Inaccuracies of analytical models are usually caused by stratification of the soil in which a grounding system is placed as well as geometric factors. Therefore, numerical techniques are used to calculate the potential distribution around the grounding system. Numerical techniques that can be used for modeling the grounding system are the finite difference method (FDM), finite element method (FEM) [5-9], charge simulation method (CSM) [10-12], boundary element method (BEM) [1318] and hybrid combinations of these methods, like the FEM/BEM method [19]. The first two methods, FDM and FEM, are rarely used to calculate the current field generated by a grounding mesh because they require discretization of the entire domain (grounding mesh conductor and surrounding soil). Also, large differences in the size between subdomains that need to be discretized result in a large number of segments. All this leads to large matrix systems that need to be solved.
Received 11 April 2015 Accepted 23 January 2017
40
MUJEZINOVIC, MULAOSMANOVIC, MUHAREMOVIC, TURKOVIC, BAJRAMOVIC
The next two mentioned methods, CSM and BEM, unlike FDM and FEM, do not require discretization of the entire domain, but only the boundary surfaces. Also, these two methods do not need the infinite boundaries to be discretized. By applying an appropriate Green's function, the need for discretization of the boundary surface soil-air (earth surface) and boundaries between two soil layers is avoided, which ultimately leads to the need for discretization of only boundary area of the grounding system. The result is a significantly smaller matrix system compared to FDM and FEM. But because CSM and BEM take into account the mutual impact between segments, these matrices are completely filled, while in FDM and FEM the matrices are rarely filled. The major difference between the CSM and BEM method is the way of solving the integral equation that describes the distribution of the current field in the soil. In the CSM method, this integral is solved analytically, setting the point current source in the center of the grounding segment. This approach is very easy to implement, but its use is limited by the requirement for the discretization of the grounding system into a large number of elements [20]. This problem can be overcome by applying the BEM method in which the solution of the field integral equations is evaluated numerically.
The last method, FEM/BEM, is based on hybridization of the FEM and BEM methods. This method is organized so that the grounding system and the soil near the grounding system is solved with the FEM method, while the problem of the infinite boundaries and layered soil is solved with the BEM method. The problem of this method is the need to determinate the dimension of the FEM domain which can significantly affect the precision of calculation of the relevant parameter grounding [20].
In this paper, an approach based on the BEM method and non-uniformity potential correction factor is used to calculate the grounding mesh parameters. For the solution of the infinite series that appear in the case of modeling the grounding mesh in a multi-layer soil, the Aitken's 52 algorithm is used.
2 Grounding Mesh Current Field
The Maxwell's electromagnetic theory can be used to model the phenomenon of the fault current dissipation in the soil around a grounding mesh [13]. Calculation of the current and potential distribution of the grounding mesh comes down to the solution of the Laplace's partial differential equation with the use of appropriate boundary conditions [21]. The Laplace's partial differential equation can be obtained from Maxwell's equation:
curl E = 0	(1)
In stationary-current fields, the first Kirchhoffs law must be satisfied [22], which is in a differential form given as:
div a = 0	(2)
where a is the vector of current density.
The analytical relationship between electric potential y at some point and electric field vector E is given as:
E = —V p = —grad p	(3)
In the case of a stationary-current field, calculation of distribution of the current density and potential is based on the assumption that the Ohm's law can be applied on the flow of the current in the soil. The potential and current are distributed so that the total Ohmic losses are reduced and distribution adjusted to the boundary conditions. From the above it follows that the vector density in a linear environment must meet the generalized version of the Ohm's law:
a = y • E = —y- gradp	(4)
where y is the soil conductivity.
Substituting equation (4) in relation to the first
Kirchhoff's law in a differential form, the following is
obtained:
diva = —div(y • gradp)	(5)
Assuming that the soil conductivity is a scalar value and by applying a standard vector identity:
div(grad p) = Ap	(6)
By applying a standard vector identity on the equation (5), the Laplace's partial differential equation is obtained, which has the form:
Ap = 0	(7)
The potential at any point of the soil can be obtained by solving the Laplace's partial differential equation (7). In order to obtain a unique solution, it is necessary to use appropriate boundary conditions in addition to the Laplace's partial differential equation.
3 Mathematical Model
As mentioned above, distribution of the grounding mesh potential is described by the Laplace's partial differential equation. To solve the analyzed problem, it is necessary to introduce the Green's function as follows:
where E is the vector of the electric field.
A COMBINED BOUNDARY ELEMENT AND AN ANALYTICAL APPROACH TO GROUNDING MESH MODELING IN
41
G(p, q) =
4n
(8)
p - q
where |p - q\| is the Euclidian's distance between source pointp and observation point q. By applying the Green's symmetrical identity on both the Laplace's partial differential equation (7) and Green's function (8), the potential of any point of the domain can be calculated by using the following integral equation:
cpk ) = ^ G(p, q)-Ç,(p)^
dn
dn
dT (9)
¿q) = ^ Gp, q)dT
(10)
The following mathematical shift can be used:
dv(p) = m( p)
dn y
(11)
where a(p) is the unknown current density of source point p.
Now, integral equation (10) can be written in the following form:
1 f
v(q)=- Mp)G(p, q)dT
y J
(12)
section/length is very small, so that the grounding conductors can be discretized with the 1D boundary elements. Also, after discretization of the grounding mesh conductors, transformation from the global to the local coordinate system is performed [23].
In this paper, geometry and current density on one 1D boundary element is approximated with the second-order polynomial as follows:
^^ ())=Z xe-w, ())
i=1 3
r ())=Z ye-w, ())
i=1 3
where dq/dn is a derivative of the potential in the direction of the outward normal vector, and dG/dn is a derivative of the Green's function in the direction of the outward normal vector, r is the surface of the grounding mesh conductors.
By applying an adequate Green's function, the whole analysed system is characterized with only the Dirichlet's boundary conditions. Therefore, one of the two planar integrals can be eliminated. Integral equation (9) now takes the following form [23]:
z
'())=£ ze w ())
i=1 3
M ()) = ^M-w, ())
(13 a) (13 b) (13.c) (13 d)
The used 1D boundary elements for discretization of the grounding mesh in the global and local coordinate system are shown in Fig. 1. Shape functions yt(£) of the second-order polynomial approximation are shown in Fig. 2.
£
4 = -1
4 = 0
4 = 1
t
//2
X
//2
3
Figure 1. 1D boundary element used for discretization of the grounding mesh.
The above equation is the Fredholm's integral equation of the first kind. The Fredholm's integral equation and adequate boundary conditions are the basis for solving the stationary-current field of the grounding mesh. To solve this equation, the indirect boundary element method was applied.
3.1 Boundary element method
To solve integral field equation (12), the grounding mesh was discretized on boundary elements. As the length of the grounding mesh conductor is significantly larger than the cross-section, i.e. the ratio of the cross-
			-w) -w) w)	/
				
				
				
				
				
Figure 2. Shape function for the second-order polynomial approximation.
1
1
T
i=1
T
z
t
0.8
0.6
0.4
0.2
0
-0.5
0.5
42
MUJEZINOVIC, MULAOSMANOVIC, MUHAREMOVIC, TURKOVIC, BAJRAMOVIC
After geometry discretization and coordinate system transformation, the collocation method at the point is applied on the integral field equation. After applying the collocation method at the point, integral equation (12) can be written as:
<Pk) = 1 £ } £ Vf -Vt (^G (f, q)det J (14)
y e=1 -1 t=1
where aet is the current density on the t-th collocation point of the e-th boundary element and detJ(^) is the determinate of the Jacobean matrix.
By applying the Gauss- Legendre's quadrature rule, the above equation can be written as an algebraic equation:
1 "e ng 3 I	\
p(q )=- £££ Ve - v> & )pe &, q)det jy
y
=1 m=1 t=1
(15)
where ng is the number of the Gauss-Legendre's integration points and wm is the m-th weighting coefficient.
Finally, the above equation can be written in the following matrix form [16]:
[r]-m=w
(16)
"e	e 3
^ =£ie =£ £ v - r
(17)
=1 t=1
where If is the fault current that enters into the grounding mesh, Ie is part of the fault current on the e-th boundary element and let is the calculated length in the i-th collocation point of the e-th boundary element. Then, matrix equation (16) can be written in the following form [24]:
R11	R12 '	" R1" L"cp	-1
R21	R22 '	Ri"	-1
R"cp 1	R"cp 2 '	■■ R" " cp cp	-1
l1	l2 '	l"cp	0
v

v„
Pg
(18)
where is the electric potential of the grounding mesh. Two most important parameters to determine when designing any grounding mesh are the grounding mesh resistance and potential distribution on the earth surface [25]. When the current densities of each collocation point and potential of grounding mesh are known, the grounding resistance of the analyzed grounding mesh can be easily determined using the following relation:
R = -
Pg
"e 3
(19)
££ve - ie
e=1 t=1
where [R] is the square matrix whose dimensions are "cp x "cp, while {a} and are vectors of unknown current densities and potentials with dimensions "cp x 1, respectively. Therefore, the number of the unknown variables is 2"p, where "cp is the number of collocation points on whole segments of the grounding mesh and is equal to 3me. Assuming that all conductors of the grounding mesh are on the same potential the number of the unknown variables is significantly reduced. With this assumption, the number of the unknown variables is reduced to "cp + 1. Therefore, it is necessary to add an additional equation to provide a unique solution. As the resistance of the grounding conductors is significantly lower than the resistance of the surrounding soil, it is permissible to assume that the entire fault current that enters the grounding mesh dissipates into the surrounding soil. This can mathematically be written in the following form:
When the current densities of all boundary elements are known, the potential distribution on the earth surface can be calculated. As the potential on the earth surface varies from one point to another matrix equation (18) needs to be modified. The matrix equation to calculate the potential distribution on the earth surface can be written in the following form:
R11	R12 '	" R1 "cp
R21	R22 '	' R2"cp
R31	R3 2 '	' R3"cp
R'n„ 1 R
"„2
R'„
Pp1
Pp2
Pp3
P
p"p
(20)
where is the matrix whose elements are calculated with the Green's function that is suitable to calculate the potential on the earth surface. The dimensions of this matrix are "p x "cp where "p is number of the points on the earth surface in which the potential needs to be determined. Vector {^p} is the vector of unknown potentials on the earth surface.
From the calculated values of the potential on the earth surface, the touch voltage and step voltage can be determined according to the definitions of the touch and step voltage given in IEEE Std. 80-2000 [26], by using the following equations [27]:
Ut = max
pm y 0) - Pp ta ± ^ y 0 ^ pm y 0) - Pp (xo, y 0 ±1(m))
(21)
cp
I
F
e
a
a
e
"
cp
p cp
e
e
A COMBINED BOUNDARY ELEMENT AND AN ANALYTICAL APPROACH TO GROUNDING MESH MODELING IN
43
Us = max
\<Pp ta , yo )~Vp (xo ± 1(m) yo ^ \pp (tao, yo )~Pp (xo , y o ± 1(m))
(22)
where U is the touch voltage, U is the step voltage and ^m is the maximum value of the potential of the metal parts of power substations during a fault.
3.1.1 Multi-layer soil model
To model a grounding mesh, an adequate soil model should be assumed. The Green's function given by relation (8) is valid only for the case when the grounding mesh is placed in an infinite medium. Though not possible in practice, the Green's function can be used for a deeply buried grounding mesh in a homogeneous soil. As the grounding mesh is very seldom buried to a depth such that the soil can be computed as an infinite medium and the fact that the soil is almost always layered, the Green's function given by relation (8) is rarely used. Therefore, to model a grounding mesh buried low in a homogeneous soil, the following Green's function should be used:
1
G( p, q) = — ■
1
1
4n U|p-q 11 jipj-q 11
(23)
Gii( p,q) =
4n
1
1
iip-qii iipi-qii
œ	5	!
tiP t?iipi-qii
(24)
g22( p, q) = ~r~
(
4n
1
P
\
1 -P2 4n
iip'-qii ii p3 - qii 1
A
(25)
iip; - qii
-tPn
n=1
ii p5 - qii
where n is number of summands to consider in the evaluation of the series of the Green's function until convergence is achieved [3o], and i is the order of the source image. In relations (24, 25), coefficient ft is the reflection factor given by:
P =
n - y 2 Ï1 +Ï2
(26)
where \\p'i - q\\ is the Euclidian's distance between the image of source point p'i and observation point q. In notation p'i, index 1 represents the order of the image. In practical situations, a homogenous soil is very rarely encountered. It is often composed of multiple layers of different electrical conductivities. In such situations, it is convenient to apply a multi-layer soil model with horizontal change of the soil conductivity. In the case of a two-layer soil, the Green's functions have the shape of an infinite sum. The grounding mesh can be positioned in any layer of a horizontally layered soil, so it is necessary to consider two possible scenarios. If source point p and observation point q are both in the upper layer of a two-layer soil, then the Green's function takes the form [28-30]:
where yi is the soil conductivity of the upper layer and Y2 is the soil conductivity of the lower layer.
From equation (14), it can be noted that to calculate the elements of both expanded matrices in equation (18) and in equation (20), the soil conductivity is needed. In a two-layer soil model, to calculate the elements of these matrices, the electrical soil conductivity of the layer in which the source point (grounding mesh) is placed is used.
3.2 Grounding mesh potential non-uniformity correction factor
It can be noted from the Green's functions for the two-layer soil model that infinite series appear. These infinite series can affect the speed of calculation of the grounding mesh parameters. For the case when the grounding mesh is placed in the upper layer of a two-layer soil, the infinite series can be avoided by using non-uniformity correction factors given in [31, 32]. In this case, the grounding mesh potential is calculated by using only the Green's function given in relation (23) and soil conductivity of the upper layer. On the calculated value of the grounding mesh potential, the following relation is applied:
Pgc (n1,y2) = M 2/1 ■PG (7\)
(27)
where $gc (yi, 72) is a corrected value of the potential of the grounding mesh placed in a two-layer soil and M1/2 is equal to:
M 2/1 =
(28)
1
1
X
For the case when both a source point p and observation point q are placed in a lower layer of a two-layer soil, the Green's function takes the form [28,29]:
and x is equal to:
44
MUJEZINOVIC, MULAOSMANOVIC, MUHAREMOVIC, TURKOVIC, BAJRAMOVIC
0.075 • 0
L VP
log
V
2
h
h <4Â
otherwise
(29)
Voc (ri,r2,rs) = M3/2 Vo Oi^
(30)
where:
M 3/2 =
( \x Ï3
(3i)
and in this case, x is equal to:
0.075 • 0
< n ^
log
h1 + h2
h + h < -Ja
h + h2 otherwise
(32)
where A is the surface of the grounding mesh and h is the thickness of the upper layer of the soil. In some situations, a two-layer soil model can be sufficient and a three-layer soil model must be used. In case of a three-layer soil model, more complex Green's functions must be adopted [33]. If the grounding mesh is placed in the first layer or in the second layer of a three-layer soil, these complex procedures can be avoided in a similar way as previously by using the non-uniformity correction factors for a three-layer soil model. For the case when the grounding mesh is placed in the first layer, the previously given procedure can be used, while if grounding mesh is placed in the second layer, the Green's functions (25) must be used to calculate the grounding mesh potentials. Then, for both cases, the corrected value of the potential of the grounding mesh can be calculated as follows:
Definition of the input parameters (fault curent, soil model, soil conductivity, geometry)	
	
Geometry discretization	
	
Analysis of segment position and choice of Green's function(s)	
	
Calculation of grounding system potential	
	
Calculation of collocation point current densities	
where h1 is the thickness of the first layer of the soil and h2 is the thickness of the second layer of a three-layer soil.
After calculating the corrected value of the grounding mesh potential, the corrected value of the current density for all collocation points must be calculated by using relation (18). From these corrected values of the current densities both the grounding resistance and the potential on the earth surface can be calculated by using relations (19) and (20) and finally the touch and step voltage from relations (21) and (22).
A flow chart of the grounding mesh parameters calculation procedure is shown in Fig. 3.
Correction of the grounding system potential
Calculation of the grounding resistance, touch and step voltages
Figure 3. Flow chart of the calculation procedure. 3.3 Aitken's S2 Algorithm
It can be noted form the Green's functions for the two-layer soil model that infinite series appear, thus significantly affecting the speed of calculation of the grounding mesh parameters. To accelerate the convergence of the infinite series, different acceleration algorithms are proposed in the literature such as Euler -Maclaurin's summation formula [30], Wynn's e method [34] and Aitken's 52 algorithm [15, 34]. In this paper, the Aitken's 52 algorithm is used for convergence acceleration of the infinite series that appear in the Green's functions for a two-layer soil model. The iterative procedure of the Aitken's 52 algorithm used in this paper can be written in the following form:
f(n) - f(n) _ Jk+i = Jk
J
(n+i)
_ fkn )
j(n+2) _ 2f (n+i)
+ f(n)
(33)
where k is the number of iterations, n is the number of element in the series and f is the element of the series. For example, if calculations are carried out with the Green's function (25), the n-th element of the series has the form:
f (n ) = pn
1
IIp5_ qll
(34)
x = <
2
X = <
2
A COMBINED BOUNDARY ELEMENT AND AN ANALYTICAL APPROACH TO GROUNDING MESH MODELING IN
45
The procedure iterates until a sufficient accuracy is achieved. This procedure was done individually for every interaction for the case when the Green's function with infinite series was used.
4 Model Application
In order to show a practical application of the presented mathematical model, calculations of the potential distribution on the earth surface, touch and step voltage were performed on a real power substation grounding mesh. Calculations were performed on the grounding mesh of the Sidi Rached (60/30 kV) power substation. Geometry of the analyzed grounding mesh with the dimension (red lines) and the profile lines (blue lines) for the plot of step and touch voltage is given in Fig. 4. Other relevant input parameters required for the calculations are given in Table 1.
m
Figure 4. Geometry of the analyzed grounding mesh. Table 1. Input parameters for calculation.
Parameter	Value
Expected fault current:	18.9 (kA)
First layer soil conductivity:	0.014 (S/m)
Second layer soil conductivity:	1.58 (S/m)
Third layer soil conductivity:	1.67 (S/m)
Thickness of the soil first layer:	0.36 (m)
Thickness of the soil second layer:	9.46 (m)
The depth of the grounding mesh:	0 . 99 (m)
Duration of the fault current:	0,35 (s)
Conductivity of the crush rock on the surface:	3 . 3 10"4 (S/m)
Thickness of the crush rock on the surface:	0,12 (m)
Human body mass:	70 (kg)
Frequency:	50 (Hz)
Results of calculating of the potential distribution on the earth surface are given in Fig. 5.
0 0
Figure 5. Geometry of the analyzed grounding mesh.
As seen from the calculated potential distribution on the earth surface, the highest values of the potential are at areas around the center of the grounding mesh, the maximum value of the potential is below 800 (V). The highest potential gradients are at the edges and outside the grounding mesh. This is an indication that in these areas, are high values of the touch and step voltage. In order to see whether these values are above the allowable limits, distribution of touch and step voltage should be calculated.
Distribution of the touch and step voltage is shown in Figs. 6 and 7, respectively.
Figure 6. Touch voltage distribution on the earth surface.
Figure 7. Step voltage distribution on the earth surface.
In Figs. 8 and 9, a comparison is given of the touch and step voltage values for the profile lines outlined in the Fig. 5 and permissible values for the touch and step voltage are given.
46
MUJEZINOVIC, MULAOSMANOVIC, MUHAREMOVIC, TURKOVIC, BAJRAMOVIC
Figure 8. Distribution of the touch voltage along the profile lines.
Figure 9. Distribution of the step voltage along the profile lines.
From the results given in the Figs. 8 and 9, it is clearly visible that the values of the touch and step voltage are significantly below the limit values. Therefore, it can be concluded that the given geometry of the grounding mesh satisfies the safety criteria in the analyzed power substation.
5 Conclusion
To calculate the grounding mesh parameters, because of their important role for the safety of the staff and equipment of power substations, very accurate and precise methods should be used. In this paper, mathematical model is presented based on a combination of the boundary element method and non-uniformity correction factor to be used in calculation of the parameters of a grounding mesh placed in a homogeneous or stratified soil (two and three-layered soil). The model can be used for most cases that can be encountered in practice.
The model is organized so as to use that best characteristics of each of the two methods. The boundary element method is used to calculate the grounding mesh in a homogeneous soil and the non-uniformity correction factor is then used to correct the grounding mesh potential and current density of all colocation points on the grounding mesh. Because of the lack of an adequate non-uniformity correction factor
for the case when the grounding mesh is placed in a second layer of a two or three-layer soil, the boundary element method is used to calculate the parameters of a grounding mesh in a two-layer soil and the grounding mesh potential can be corrected if necessary (i.e. in case of a three-layer soil). In this case, in the mathematical model, infinite series exists. To minimize the time needed to calculate the grounding mesh parameters, the Aitken's 52 algorithm is used.
Acknowledgment
This research has been supported by the Ministry of Education and Science of the Federation of Bosnia and Herzegovina.
References
[1]	A. Mujezinovic, Analysis of discontinuity of electrical conductivity of soil on grounding resistance, MSc. Thesis, Faculty of Electrical Engineering, University of Sarajevo, July 2011 (in Bosnian).
[2]	A. Mujezinovic, A. Muharemovic, A. Muharemovic, I. Turkovic, Z. Bajramovic, Application of Finite Element Method in Calculation of Large and Complex Grounding Systems, International Conference and Exposition on Electrical and Power Engineering (EPE), pp. 688 - 692., Romania, Iasi, October 2012.
[3]	A. Mulaosmanovic, Calculation of grounding system parameters of large electric power facilities, MSc. Thesis, Faculty of Electrical Engineering, University of Sarajevo, September 2015 (in Bosnian).
[4]	J. Nahman, S. Skuletich, Resistance to ground and mesh voltages of ground grids, PROC IEE 126 (1), pp. 57 - 61, January 1979.
[5]	M. V. F. da Luz, P. Dular, Analytical and finite element modeling of grounding system, IX International Symposium on Lightning Protection SIPDA, Brazil, Foz do Iguacu, November 2007.
[6]	J. A. Guemes, L. Azpiazu, F. E: Hernando, Earth network resistance calculation using 3-D solid models, 17th International Conference on Electricity Distribution, CIRED, Barcelona, May
2003.
[7]	I. Paunovic, J. Nahman, Complex grounding systems analysis using finite element method, The Journal of the Union of Yugoslav Electric Power Industry 3 (3), pp. 92-100, Belgrade
2004.
[8]	M. I. Lorentzou, A. G. Klads, N. D. Hatziargyriou, Finite Element Modelling of Grounding System Considering Electrode Geometry Effect, IEEE Trans. on Magnetics, Vol. 35., No. 3., pp. 1757 -1760., May 1999.
[9]	M. Trlep, A. Hamler, B. Hribernik, The FEM analysis of complex grounding system by the program package FIELD_GF, Trans. on Engineering Sciences 22., 1999 © WIT Press, www.witpress.com, ISSN 1743 - 3533, pp. 243 - 252.
[10]	E. Bendito, A. Carmona, A. M. Encians, M.J. Jimenez, The Extremal Charges Method in Grounding Grid Design, IEEE Trans. on Power Delivery 19 (1), pp. 118 - 123, January 2004.
[11]	F. Freschi, M. Mitolo, M. Tartaglia, An Effective Semianalytical Method for Simulating Grounding Grid, IEEE Trans. on Industry Applications 49 (1), pp. 256 - 263, Jan/Feb. 2013.
[12]	F. Freschi, M. Mitolo, M. Tartaglia, Interferences Phenomena between Separete Grounding Systems, IEEE Trans. on Industry Applications 50 (4), pp. 2853 - 2860, July/Aug. 2014.
[13]	I. Colominsa, F. Navarrina, M. Casteleiro, A Boundary Element Numerical Approach for Earthing Grid Computation, Computer methods in applied mechanics and engineering 174 (1), pp. 73 -90, April 1999.
[14]	F. Navarrina, I. Colominas, M. Casteleiro, Why Do Computer Methods For Grounding Analysis Produce Anomalous Results?,
A COMBINED BOUNDARY ELEMENT AND AN ANALYTICAL APPROACH TO GROUNDING MESH MODELING IN
47
IEEE Trans. on Power Delivery 18 (4), pp. 1192 - 1202, October 2003.
[15]	I. Colominas, J. París, F. Navarrina, M. Casteleiro, Improvement of the computer methods for grounding analysis in layered soils by using high - efficient convergence acceleration techniques, Advances in Engineering Software 44(1), pp. 80 - 91, 2012.
[16]	H. Zildzo, A. Muharemovic, I. Turkovic, H. Matoruga, Numerical Calculation of Floating Potentials for Large Earthing Systems,
XXII	International Symposium on Information, Communication and Automation Technologies, Sarajevo, October 2009.
[17]	A. Muharemovic, V. Madžarevic, I. Džananovic, A. Muharemovic, M. Ramic, Grounding System - Design, Measurement and Regulations, Faculty of Electrical Engineering Tuzla, Tuzla 2011.
[18]	Z. Haznadar, S. Berberovic, Numerical procedure for calculation and design of grounding systems, Elektrotehnika 27(1-2), pp. 9 -17, 1984.
[19]	M.Trlep, A. Hamler, M. Jesenik, B. Štumberger, The FEM-BEM Analysis of Complex Grounding Systems, IEEE Trans. on Magnetics 39(3), pp. 1155 - 1158, May 2003.
[20]	S. S. M. Ghoneim, Charge and Current Simulation Method with Boundary Element Method for Grounding System Calculations in Case of Mult-Layer Soil, IOSR Journal of Engineering (IOSRJEN) 3 (6), pp. 14 - 22, June 2013.
[21]	A. Mujezinovic, A. Muharemovic, A. Muharemovic, Current Field Distribution of the Grounding System in Case Instalation of Underground Steel Pipeline in the Vicinity, BH K Elektortehnika, 6 (1), pp. 35 - 39, December 2012.
[22]	E. Hot, I. Hot, I. Hot, A. Muharemovic, A. Dautovic, Theory of Electromagnetic fields, Faculty of Electrical Engineering Sarajevo, Sarajevo 2007.
[23]	C. A. Brebbia, J. C. F. Telles, L. C. Wrobel, Boundary Element Method Theory and Application in Engineering, Berlin/Heidelberg, Springer-Verlag, 1984.
[24]	H. Zildžo, Application of Computer in Power Engineering, Faculty of Electrical Engineering Sarajevo, Sarajevo 2004.
[25]	I. Turkovic, Contribution to tracking the neutral point of 20(10) kV networks for the purpose of eliminating excessive touch voltage in the substations 20(10)/0.4 kV, MSc. Thesis, Faculty of Electrical Engineering, University of Sarajevo, December 2002 (in Bosnian).
[26]	IEEE Std. 80 - 2000, IEEE Guide for Safety in AC Substation Grounding, IEEE New York, 2000.
[27]	D. Poljak, The theory of electromagnetic fields with applications in engineering, Školska knjiga, Zagreb 2014.
[28]	Z. Haznadar, S. Berberovic, D. Vujevic, D. Markoninovic, N. Žanic. Z. Žic, Touch Voltage and potential distribution of grounding system in power systems, Journal of Energy 7 (8), Zagreb 1981.
[29]	H. Zildžo, A. Muharemovic, H. Matoruga, A. Mujezinovic, Application of coupled 1D and 2D isoperimetric discontinuous boundary elements in the calculation of large grounding systems,
XXIII	International Symposium on Information, Communication and Automation Technologies, Sarajevo, October 2011.
[30]	S. Berberovic, Z. Haznadar, Ž. Štih, Method of moments in analysis of grounding systems, Engineering Analysis with Boundary Elements 27 (3), pp. 351 - 360, 2002.
[31]	J. Nahman, I.Paunovic, Mesh voltages at earthing grids buried in multi-layer soil, Electric Power Systems Research 80 (5), pp. 556 - 561, 2010.
[32]	J. Nahman, V. B. Djordjevic, Nonuniformity Correction Factors for Maximum Mesh- and Step-Voltages of Grounding Grids and Combined Ground Electrodes, IEEE Trans. on Power Delivery 10 (3), pp. 1263 -1269, July 1995.
[33]	J. He, R. Zeng, B. Zhang, Methodology and Technology for Power System Grounding, John Wiley & Sons, Singapore 2013.
[34]	E. J. Weniger, Nonlinear Sequence Transformations for the Acceleration of Convergence and the Summation of Divergent Series, Computer Physics Reports 10 (5), pp. 189 - 371, 1989.
Adnan Mujezinovic received his B.Sc. and M.Sc. degrees in electrical engineering from University of Sarajevo (Bosnia and Herzegovina), in 2009 and 2011, respectively. From 2012 he has been working at Faculty of Electrical Engineering, University of Sarajevo, as an assistant. His topics of interest are numerical calculations of electromagnetic fields, cathodic protection and grounding systems.
Ajdin Mulaosmanovic received his B.Sc. and M.Sc. degrees in electrical engineering from the Faculty of Electrical Engineering, University of Sarajevo, Bosnia and Herzegovina, in 2013 and 2015, respectively. Currently, he is employed with the Public Electric Utility Elektroprivreda of Bosnia and Herzegovina, HPP Neretva. His areas of interest are maintenance of electrical equipment, power generation, electrical protection and grounding systems.
Alija Muharemovic received his B.Sc., M.Sc. and Ph.D. degrees in electrical engineering from University of Sarajevo (Bosnia and Herzegovina), in 1974, 1981 and 1989, respectively. Currently, he is a full Professor at the Faculty of Electrical Engineering, University of Sarajevo. His topics of interest are cathodic protection, electrical measurements, electrical installations and safety measures, and EMC in power systems.
Irfan Turkovic received his M.Sc. and Ph.D. degrees in electrical engineering from the University of Sarajevo, Bosnia and Herzegovina, in 2003 and 2010, respectively. Currently, he is an associate Professor at the University of Sarajevo -Faculty of Electrical Engineering. His topics of interest include cathodic protection, electrical measurements, safety in low voltage network and EMC in power systems.
Zijad Bajramovic received his M.Sc. and Ph.D. degrees in electrical engineering from the University of Sarajevo, Bosnia and Herzegovina, in 1986 and 2009, respectively. Currently, he is an associate Professor at the Faculty of Electrical Engineering, University of Sarajevo. His areas of interest are testing of electric power equipment, electromagnetic compatibility research, digital simulation of electromagnetic transients, and insulation coordination in power system. He is a member of the CIGRE and IEEE.