The IEEE Model for a Ground Rod in a Two Layer Soil – A FEM Approach

The calculation of a ground electrode resistance, using a two layer soil model, has been widely presented in literature. Several methods had been used. Formulas for grid in two layers soil using the synthetic-asymptote approach have been developed in (Salama et al., 1995). Berberovic explored the Method of Moments in the calculation of ground resistance, using higher order polynomials approximation in the unknown current distribution in (Berberovic et al., 2003), and the Galerkin’s Moment Method with a variation was used in (Sharma & De Four, 2006). Another theoretical tool commonly used is the Boundary Element Method, as in (Colominas et al., 1998, 2002a, 2002b; Adriano et al., 2003). These authors transformed the differential equation that governs the physical phenomenon into an equivalent boundary integral equation. The Matrix/Integration Method for calculating the mutual resistance segment in one and two layered soil was adopted by (Coa, 2006) and an optimised method of images for multilayer soils was used in (Ma et al., 1996). Even in the study of ionization phenomena, the two layer ground model was used in (Liu et al., 2004). In general these works used the theory of images, which implies infinite series for the expanded Green function as in (Berberovic et al., 2003). Recently, a work presented the effect of low resistance materials filling in a pit surrounding a rod, working with two different soil resistivity’s (Al-Arayny et al., 2011). This type of research was also presented in (Zhenghua et al., 2011), that even considered the use of electrolytic materials. A Finite Element Method (FEM) application to grounding can also be found in (Manikandan et al., 2011) to the analysis of wind turbines grounding. In this chapter the FEM is presented in a theoretical basis for cylindrical symmetry problems, using a ground rod resistance calculation as an example. Comparison with experimental result is also made.


Introduction
The calculation of a ground electrode resistance, using a two layer soil model, has been widely presented in literature. Several methods had been used. Formulas for grid in two layers soil using the synthetic-asymptote approach have been developed in (Salama et al., 1995). Berberovic explored the Method of Moments in the calculation of ground resistance, using higher order polynomials approximation in the unknown current distribution in (Berberovic et al., 2003), and the Galerkin's Moment Method with a variation was used in (Sharma & De Four, 2006). Another theoretical tool commonly used is the Boundary Element Method, as in (Colominas et al., 1998(Colominas et al., , 2002a(Colominas et al., , 2002bAdriano et al., 2003). These authors transformed the differential equation that governs the physical phenomenon into an equivalent boundary integral equation. The Matrix/Integration Method for calculating the mutual resistance segment in one and two layered soil was adopted by (Coa, 2006) and an optimised method of images for multilayer soils was used in (Ma et al., 1996). Even in the study of ionization phenomena, the two layer ground model was used in (Liu et al., 2004). In general these works used the theory of images, which implies infinite series for the expanded Green function as in (Berberovic et al., 2003). Recently, a work presented the effect of low resistance materials filling in a pit surrounding a rod, working with two different soil resistivity's (Al-Arayny et al., 2011). This type of research was also presented in (Zhenghua et al., 2011), that even considered the use of electrolytic materials. A Finite Element Method (FEM) application to grounding can also be found in (Manikandan et al., 2011) to the analysis of wind turbines grounding. In this chapter the FEM is presented in a theoretical basis for cylindrical symmetry problems, using a ground rod resistance calculation as an example. Comparison with experimental result is also made.
This chapter presents the Finite Element Method for the calculation of a rod resistance in a two layer soil. Theoretical basis of the method are presented. The FEM mesh was tested in homogeneous soil, first a cylindrical approach to IEEE model was considered, in order to guarantee that mesh is voltage and energy adapted, and then the whole model was discretized. Resistance was calculated using Joule's law in a FEM energetic approach. The zero volt Dirichlet boundary was located within a distance of 3 cm, 15 cm and 7.6 m from the ground electrode, to analyze if the resistance relations indicated by IEEE in homogeneous soil could be also used in a two layer soil. The simulated results were compared with those obtained from Tagg formula. In homogeneous soil the errors are less than 1%. In two layer soil, resistance error to Tagg formula decreased from values of 28% with zero volt boundary at 3 cm to 22% with zero volt boundary at 15 cm; for the whole model discretized the error is near -18%. The results were unsatisfactory, thus the percentage resistances at these distances cannot be generalized in a two layer soil. However the error between simulated and field measured values is of 4.6%, turning the FEM analysis a valuable simulation tool.

Finite element method
The first order finite elements method using triangular elements may be regarded as twodimensional generalizations of piecewise-linear approximation techniques as in reference (Sylvester & Ferrari, 1990), widely used in Electrical Engineering. The method allows several choices for mesh types and an easy treatment for boundary shapes.
Several problems in electrical engineering require the solution of Laplace equation in two or three dimensions with two kinds of boundary conditions, such as prescribed potential values along the referred boundaries, Dirichlet conditions, and vanishing normal derivative along the symmetry planes, Neumann condition. As an example, in power system grounding, where capacitive and inductive effects are not considered, since industrial frequency is too low, the soil potential satisfies Laplace equation. On the other hand grounding systems dimensions are much smaller than power line wavelength, so that propagation phenomenon is not considered. Laplace equation solution is equivalent, according to the minimum potential energy principle, to the following energy (   W u ) functional minimization, that stores field energy per unit volume, as in (Sylvester & Ferrari, 1990): where u i s t h e potential and V the volume.
The integral is evaluated over all the volume defined by the problem boundary.

Discretization
In Fig. 1 the IEEE model for a ground rod is presented. The boundary conditions within a distance of 7.6 m from the electrode have zero potential. At surface, normal derivative potential is also zero, since there is no current flowing in the air.
The IEEE Model for a Ground Rod in a Two Layer Soil -A FEM Approach 145 Figure 1. IEEE model as in (IEEE Std 142, 2007) Due to isotropy, the problem was discretized in the 'rz' plane, since cylindrical coordinates were used. To obtain an approximate solution by FEM the problem region is subdivided into triangular elements. The solution mesh is presented in Fig. 2, near the rod top end.

First order triangular elements
The essence of the method lies in first approximating the potential u within each element in a standardized fashion, and thereafter interrelating the potential distribution in the various elements so as to constrain the potential to be continuous across interelement boundaries, as in reference (Sylvester & Ferrari, 1990). This standardized fashion is done using a first order polynomial for potential within each element, which is the approximate solution to the actual one.
( , ) u r z a br cz    The true potential distribution is thus replaced by a piecewise planar function and within each triangle side, potential is obtained by a linear interpolation of node potentials.
Considering a generic triangular finite element, as shown in Fig. 3, equation (2) must satisfy node potentials. Using the equation in the three nodes a system is obtained, allowing the constants 'a', 'b' and 'c' to be calculated as functions of node potentials.
where U1, U2 and U3 are node potentials. Using Cramer´s rule the referred constants are easily obtained. The formula for 'a' is: The denominator is recognized as twice the triangle area (2A) as in reference (Sylvester & Ferrari, 1990). It must be pointed that this formula was found by integration so that area become negative if the nodes are clockwise numbered. Equation (2) is useful when node potentials are already known. In order to calculate these potentials is useful to represent potential within the element as a function of nodes potentials, so equation (2) became: The factors multiplying nodes potentials, after being divided by 2A, are known as shape functions, so we can rewrite previous equation using these functions.
Or in a more elegant manner: Shape functions have two important properties; 1. Their value is one in the associated node and zero in the two others. 2. In any point inside the finite element their sum is one.
The finite element energy can now be calculated using equation (1). Node potentials, although unknown, are constants so the potential gradient is The energy can now be evaluated, using equation (8) and (1). Considering that the integrand has nine sums, due to the square of potential gradient, energy functional for a single element is: Integrating by decomposition: Defining the S variable as: As shape functions are first degrees polynomials, gradients are constants, such as their inner product, so the integral is obtained by the calculus of the revolution volume of the finite element, which is 2 c r  where c r is the centroid radius. The S variable became: The remaining terms are obtained by index cyclic rotation. The finite element energy is finally represented as: where S is a 3x3 matrix and U the nodes potential vector.

Elements assembly
The total energy model is the sum of all finite elements energy. Consider the following two elements with disjoint nodal numbering: The nodes potential vector is: The dis S matrix is formed by the elementary elements matrices.
This matrix is tridiagonal and is called the Dirichlet matrix, as in reference (Sylvester & Ferrari, 1990). The energy of the two elements is: Potentials continuity implies that the values for the same node are equal. So potential at node three is equal to potential at node six, in Fig. 4, and so on. It can be defined a global nodal numbering for element assembly illustrated in Fig. 5. The potential continuity for corresponding nodes will be guaranteed by a linear transformation that relates disjoint nodes numbering with global numbering. In this case it would be: In all cases would be: Introducing this relation in equation (15) it can be obtained an energy formulation in function of global numbering nodal potentials.

 
In order to minimize energy expression   W u , it is necessary to calculate the derivative, in equation (21), with respect to potential vector and solve the equation Also, in order to avoid trivial solution 0 U  the S matrix should be partitioned in blocks, as well as potential vector: Where n U is a vector of unknown node potentials and k U is the known potential vector. Unknown potential nodes must be the first to be numbered and after the known potentials that satisfy boundary conditions. The n S and K S dimensions must allow matrix multiplication. Final solution for the unknown potential is given by: This FEM solution is called stored energy approach.

Electric field
Knowing the node potentials, the constants in equation (8) are easy to find as well as the electric field strength, which is given by:

IEEE model in homogeneous soil
The analysis of a ground rod was carried out with first order triangular finite element only in the 3 cm near the rod, since 25 % of the rod resistance is within this region as in reference (IEEE Std 142, 2007). This avoids the discretization of the IEEE entire model with a zero volt boundary condition at 7.6 m away from the rod, being the rod resistance four times the calculated value.

The IEEE model for cylindrical region
In order to test the FEM mesh, it was considered only the cylindrical region in homogeneous soil. This problem has a theoretical solution, since it is considered as two cylinders centered in the same axis. The inner cylinder has the rod radius and the outer cylinder has a 3 cm radius, which is the zero volt boundary condition. The finite elements are shown in Fig. 6.
The voltage at any point between the cylindrical surfaces is easily obtained by: with 'b' the outer cylinder radius, 'a' the inner cylinder radius, and 0 V the potential difference applied between the cylindrical surfaces.
With 'b'=3 cm, 'a'=8 mm and a voltage difference of 220 V between the two cylinders, equation (27) where 'r' is the distance to common axis.
The numerical values for voltage were calculated by FEM. For 2 m length cylinders, the potential at 1 cm away from the common axis is presented in Fig. 7. The potential is almost constant along a parallel line 1 cm away from rod axis, as expected.
In the first point its potential is 183.2 V and the last has the value of 182.1 V. The potential variation with the distance to axis 'r' is presented in Table 1.   The obtained figures for potential at distances of 1.4 cm and 2.1 cm are similar to Fig. 7. It was concluded that this mesh is voltage adapted. For the resistance between the two cylinders was used the FEM energetic approach given by (Martins & Antunes, 1997), where v is the voltage between the cylinders,  the conductivity of the material between the two cylinders, E the electric field intensity and V the revolution volume due to axial symmetry, generated by each finite element. The obtained value for the electric resistance was 20.8 Ω considering a medium with 200 Ωm for the resistivity. That value should be compared with the theoretical value for the resistance between two cylinders given by (Purcell, 1998) with P the electric resistivity, L the length, b the outer surface radius and a the inner cylindrical radius.
The value given by equation (30) is 21.0 Ω. The simulated value of 20.8 Ω is 1% less, which allows the conclusion that the mesh is also adapted in energy.

The complete IEEE model
To achieve the discretization of the entire IEEE model the mesh was altered in the rod bottom as shown in Fig. 8. Triangular finite elements have inner angles greater than 5º avoiding triangle areas close to zero, allowing stiffness matrix to be well defined. For a rod with these dimensions the Dwight formula (Dwight, 1936)

Zero volt equipotential at 3 cm
As in previous examples, it was supposed that the Dirichlet border was 3 cm away from the rod axis, accounting for 25% of rod resistance, a supposition that needs validation in a two layer soil. The results, for a 2 m length rod, 8 mm radius, buried at ground level, with an upper soil layer of 100 Ωm resistivity, and a 500 Ωm resistivity in the lower layer, are summarized in Table 2. Theoretical resistance was obtained using Tagg formula, as in reference (Tagg, 1964). The results are unacceptable. The assumption that 25% of resistance in a two layer soil is also in the first 3 cm is probably wrong.
Changing the values of the layers resistivity for the same rod, the results are presented in Table 3.
The results are acceptable but would be better if they were greater than the theoretical ones, as a safe margin.

Zero volt equipotential at 15 cm
The zero volt Dirichlet border was moved to 15 cm, where in homogeneous soil remains 50 % of the rod resistance. The mesh was changed improving the rod bottom discretization as shown in Fig. 2. The results in Table 4 were obtained for a 2 m length rod, 8 mm radius, buried at ground level, with an upper soil layer of 100 Ωm resistivity, and a 500 Ωm resistivity in the lower layer. The analysis of the results shows that the errors are quite high. The assumption that the first 15 cm contain 50 % of rod resistance seems to be incorrect. Interchanging the values of the layers resistivity, new results were obtained, and are shown in

IEEE entire model discretization
In the last simulation the IEEE whole model was discretized. The complete solution mesh is presented in Fig. 9.
In order to validate the mesh, it was used the same conductance for all finite elements. It was chosen the value of 0.005 S/m. With this value the used rod, 2 m length and 8 mm radius, has a resistance, using Dwight formula, of 94 Ω. According to IEEE, this formula has a 13% excess to real value that should be of 83 Ω. The computer program developed returned a value of 87 Ω, which is 4.4 % higher than IEEE value, making the mesh sufficiently accurate.
The next step was to simulate a two layer soil model, with an upper soil layer of 100 Ωm resistivity, and a 500 Ωm resistivity in the lower layer. The results are presented in the

Field measurements
In order to experimentally validate the model, resistivity measurements were done in a sandy soil. The obtained values are presented in Table 8.  The distance referred in the first row is the distance between test rods in Wenner method, as in (Telford et al., 1990). The resistivity curve is presented in Fig. 10 The upper layer thickness 'h' is 1.53 m. This value was rounded to 1.5 m.
The ground rod was discretized as indicated in Fig. 9. The measured value was 108 Ω and the simulated value using FEM is 113 Ω, which is 4.6 % higher. The equipotential lines were calculated and presented in Fig. 11.
Equipotential lines in bottom layer are closer to rod, since this layer has a smaller resistivity.

Conclusions
In this chapter the Finite Element Method for the calculation of a rod resistance in a two layer soil was presented. The developed meshes were tested in homogeneous soil and in a cylindrical problem, which has theoretical solutions, and they are well adapted considering potential distribution or energy dissipation.
In homogeneous soil, 25% of rod resistance is in the 3 cm rod closest soil, according with IEEE and as validated in this work.
For an upper layer soil, with resistivity smaller than the resistivity of the lower layer soil, the assumption that 25% of rod resistance is in the nearest 3 cm and 50% in the nearest 15 cm is wrong and cannot be generalized.
For an upper layer soil, with resistivity bigger than the resistivity of the lower layer soil, the assumption that 25% of rod resistance is in the nearest 3 cm and 50% in the nearest 15 cm is acceptable, but the results are not conservative.
Discretizing the whole IEEE model allowed obtaining results with less than 20% error, but these results are not conservative. The whole mesh was tested considering equal resistivity, obtaining results similar to the ones gotten from the homogeneous soil simulation.
The comparison with the field measurement is good, since simulated value for resistance is only 4.6% higher.