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

1. 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, 2002 a, 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.

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.

2. Finite element method

The first order finite elements method using triangular elements may be regarded as two-dimensional 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\left(u\right)$$) functional minimization, that stores field energy per unit volume, as in (Sylvester & Ferrari, 1990):

$$W(u)=\frac{1}{2}{\displaystyle {\iiint }_V{\left|\nabla u\right|}^{2}dV}$$E1

where $$u$$ is the potential and $$V$$ the volume.

The integral is evaluated over all the volume defined by the problem boundary.

2.1. 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.

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.

Figure 2.

Model discretization near rod top end

2.2. 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$$E2

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.

$$\begin{array}{l}{U}_1=a+b{r}_1+c{z}_1\\ U_2=a+b{r}_2+c{z}_2\\ U_3=a+b{r}_3+c{z}_3\end{array}$$E3

where U₁, U₂ and U₃ are node potentials.

Figure 3.

Finite element with nodes numbered

Using Cramer´s rule the referred constants are easily obtained. The formula for ‘a’ is:

$$a=\frac{\left|\begin{array}{ccc}{U}_1& r_1& z_1\\ U_2& r_2& z_2\\ U_3& r_3& z_3\end{array}\right|}{\left|\begin{array}{ccc}1& r_1& z_1\\ 1& r_2& z_2\\ 1& r_3& z_3\end{array}\right|}$$E4

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:

$$u=\left\{\begin{array}{l}\left[\left(r_2{z}_3-r_3{z}_2\right)+\left(z_2-z_3\right)r+\left(r_3-r_2\right)z\right]U_1+\\ \left[\left(r_3{z}_1-r_1{z}_3\right)+\left(z_3-z_1\right)r+\left(r_1-r_3\right)z\right]U_2+\\ \left[\left(r_1{z}_2-r_2{z}_1\right)+\left(z_1-z_2\right)r+\left(r_2-r_1\right)z\right]U_3\end{array}\right\}/2A$$E5

The factors multiplying nodes potentials, after being divided by 2A, are known as shape functions, so we can rewrite previous equation using these functions.

$$u\left(r,z\right)={\alpha }_1\left(r,z\right)U_1+{\alpha }_2\left(r,z\right)U_2+{\alpha }_3\left(r,z\right)U_3$$E6

Or in a more elegant manner:

$$u\left(r,z\right)={\displaystyle \sum _{i=1}^3{U}_i{\alpha }_i\left(r,z\right)}$$E7

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

$$\overrightarrow{{\displaystyle \nabla }}u={\displaystyle \sum _{i=1}^3{U}_i\overrightarrow{{\displaystyle \nabla }}{\alpha }_i\left(r,z\right)}$$E8

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:

$$W^e\left(u\right)=\frac{1}{2}{\displaystyle {\iiint }_V{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{j=1}^3{U}_i{U}_j\overrightarrow{{\displaystyle \nabla }}{\alpha }_i}}}.\overrightarrow{{\displaystyle \nabla }}{\alpha }_{j}dV$$E9

Integrating by decomposition:

$$W^e\left(u\right)=\frac{1}{2}{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{j=1}^3{U}_i{\displaystyle {\iiint }_V\overrightarrow{{\displaystyle \nabla }}{\alpha }_i.\overrightarrow{{\displaystyle \nabla }}{\alpha }_j}}}dV.U_j$$E10

Defining the S variable as:

$$S_{ij}={\displaystyle {\iiint }_V\overrightarrow{{\displaystyle \nabla }}{\alpha }_i.}\overrightarrow{{\displaystyle \nabla }}{\alpha }_{j}dV$$E11

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\pi r_c$$ where $$r_c$$ is the centroid radius. The S variable became:

1. For equal index

$$S_{11}^e=\frac{{\left(z_2-z_3\right)}^2+{\left(r_2-r_3\right)}^2}{2A}\pi r_c$$E12

1. For different index

$$S_{12}^e=\frac{\left(z_2-z_3\right)\left(z_3-z_1\right)+\left(r_3-r_2\right)\left(r_1-r_3\right)}{2A}\pi r_c$$E13

The remaining terms are obtained by index cyclic rotation. The finite element energy is finally represented as:

$$W^e\left(u\right)=\frac{1}{2}{U}^{T}SU$$E14

where S is a 3x3 matrix and U the nodes potential vector.

2.3. Elements assembly

The total energy model is the sum of all finite elements energy. Consider the following two elements with disjoint nodal numbering:

Figure 4.

Disjoint numbering of finite element nodes

The nodes potential vector is:

$$U_{dis}^T=\left[U_1\stackrel{}{{\displaystyle }}{U}_2\stackrel{}{{\displaystyle }}{U}_3\stackrel{}{{\displaystyle }}{U}_4\stackrel{}{{\displaystyle }}{U}_5\stackrel{}{{\displaystyle }}{U}_6\right]$$E15

The $$S_{dis}$$ matrix is formed by the elementary elements matrices.

$$S_{dis}=\left[\begin{array}{cc}{S}_1& 0\\ 0& S_2\end{array}\right]$$E16

This matrix is tridiagonal and is called the Dirichlet matrix, as in reference (Sylvester & Ferrari, 1990). The energy of the two elements is:

$$W\left(u\right)=\frac{1}{2}{U}_{dis}^T{S}_{dis}{U}_{dis}$$E17

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.

Figure 5.

Global node numbering

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:

$${\left[\begin{array}{l}{U}_1\\ U_2\\ U_3\\ U_4\\ U_5\\ U_6\end{array}\right]}_{dis}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]{\left[\begin{array}{l}{U}_1\\ U_2\\ U_3\\ U_4\end{array}\right]}_{conj}$$E18

In all cases would be:

$$U_{dis}=C{U}_{conj}$$E19

Introducing this relation in equation (15) it can be obtained an energy formulation in function of global numbering nodal potentials.

$$W\left(u\right)=\frac{1}{2}{U}_{conj}^T.S{U}_{conj}$$E20

with

$$S=C^T{S}_{dis}{C}_{}$$E21

In order to minimize energy expression$$W\left(u\right)$$, it is necessary to calculate the derivative, in equation (21), with respect to potential vector and solve the equation

$$SU=0$$E22

Also, in order to avoid trivial solution $$U=0$$ the S matrix should be partitioned in blocks, as well as potential vector:

$$\left[S_n\stackrel{}{{\displaystyle }}{S}_k\right]\left[\begin{array}{l}{U}_n\\ U_k\end{array}\right]=0$$E23

Where $$U_n$$ is a vector of unknown node potentials and $$U_k$$ 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 $$S_n$$ and $$S_K$$ dimensions must allow matrix multiplication. Final solution for the unknown potential is given by:

$$U_n=-S_n^{-1}{S}_K{U}_K$$E24

This FEM solution is called stored energy approach.

2.4. 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:

$$\overrightarrow{E}=-b\stackrel{\wedge }{{\displaystyle i}}-c\stackrel{\wedge }{{\displaystyle j}}$$E25

3. 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.

3.1. 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:

$$v(r)=\frac{V_0}{\ln \frac{b}{a}}\ln \frac{b}{r}$$E26

with ‘b’ the outer cylinder radius, ‘a’ the inner cylinder radius, and $$V_0$$ 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) becomes:

$$v(r)=166.4\ln \frac{0.03}{r}$$E27

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.

Figure 6.

Cylindrical region Discretization

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.

Distance (m)	Theoretical value (Volt)	Mean value (Volt)	Mean error (%)	Maximum error (%)
0.01	182.9	182.6	-0.2	-0.4
0.014	126.9	126.6	-0.2	-0.6
0.021	59.4	59.3	-0.2	-0.2

Table 1.

Potential variation between cylindrical surfaces

Figure 7.

Potential along a parallel line 1 cm away from rod axis

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),

$$R=\frac{v^2}{P_{JOULE}}=\frac{v^2}{\sigma {\displaystyle \underset{V}{\int }{\left|E\right|}^{2}dV}}$$E28

where $$v$$ is the voltage between the cylinders, $$\sigma$$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):

$$R=\frac{\rho }{2\pi L}\ln \frac{b}{a}$$E29

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.

3.2. 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.

Figure 8.

Bottom rod FEM discretization

For a rod with these dimensions the Dwight formula ( Dwight, 1936) gives a value for resistance of 94.0 Ω. According to the standard (IEEE Std 142, 2007), this formula has 13% excess, so the corrected theoretical value is 83.2 Ω. The FEM simulated value for resistance, with the zero volt boundary at 0.03 m, is 20.8 Ω, representing 25 % of total resistance, so the numerical simulated value is four times 20.8 which results in 83.2 Ω. This value is equal to the theoretical one.

4. IEEE model in a two layer soil

4.1. 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.

Upper layer thickness (m)	Theoretical resistance (Ω)	FEM simulated value x4 (Ω)	error (%)
0.5	134	104	-22.4
1.0	93.4	69.4	-25.7
1.5	72.0	52.0	-27.8

Table 2.

Resistance variation for positive voltage reflexion coefficient

Upper layer thickness (m)	Theoretical resistance (Ω)	FEM simulated valuex4 (Ω)	error (%)
0.5	60	51.9	-13.5
1.0	77	69.1	-10.3
1.5	114	103	-9.6

Table 3.

Resistance variation for negative voltage reflexion coefficient

4.2. 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 Table 5. In this case the errors are acceptable bur not conservative.

Upper layer thickness (m)	Theoretical resistance (Ω)	FEM simulated value x2 (Ω)	error (%)
0.6	124	103	-17
1.2	83.4	67.1	-20
1.8	63.4	49.6	-22

Table 4.

Resistance variation for positive voltage reflexion coefficient

Upper layer thickness (m)	Theoretical resistance (Ω)	FEM simulated value x2 (Ω)	error (%)
0.6	62.8	59.0	-6
1.2	89.1	85.3	-4
1.8	162	154	-5

Table 5.

Resistance variation for negative voltage reflexion coefficient

4.3. 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 Table 6.

Upper layer thickness (m)	Theoretical resistance (Ω)	FEM simulated value (Ω)	error (%)
0.5	134	118	-11.9
1.0	93.4	81.4	-12.8
1.5	72.0	62.0	-13.9

Table 6.

Resistance variation for positive voltage reflexion coefficient

Interchanging the values of resistivity layers, new results were obtained and are shown in Table 7. The results are acceptable, unfortunately by default. The last one is surprisingly high, but within 20% of theoretical value. It can be concluded than FEM meshes can provide useful results, considering the whole IEEE model in discretizing the problem.

Figure 9.

Entire IEEE model discretization

Upper layer thickness (m)	Theoretical resistance (Ω)	FEM simulated value (Ω)	error (%)
0.5	60	52.6	-12.5
1.0	77	66.7	-13.4
1.5	114	93.2	-18.3

Table 7.

Resistance variation for negative voltage reflexion coefficient

5. 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.

Distance (m)	0.5	1	2	3	4	5	6	8
Resistance (Ω)	552	207	41	7	4	2	3	2
Resistivity (Ωm)	1734	1301	515	132	101	63	113	108

Table 8.

Resistivity variation with depth

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. To choose a value for top layer resistivity, it was considered only the first resistivity measurement, since it wasn’t found an upper asymptote. For bottom layer resistivity it was considered the average resistivity value of the last four points values, which are almost in a horizontal line. The average value is 94.5 Ωm. The upper layer thickness is obtained considering the point where concavity changes. The three first points are almost in a straight line, and concavity is detected only after the third point. It was considered the point where resistivity curve crosses the 400 Ωm ordinate where abscissae seem to be 2.3 m. Using this value in Lancaster-Jones rule one obtains (Lancaster-Jones, 1930):

$$\frac{3h}{2}=2.3$$E30

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.

Figure 10.

Resistivity measurements with depth

Figure 11.

Equipotential lines with depth for a 2 m rod

6. 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.