## Abstract

An algorithm is proposed in this paper for calculating the impulse earthing resistances of vertical earthing electrodes. The proposed algorithm employs the average potential method to derive the formula of the low current earthing resistance. Unlike the previous algorithm, the soil ionization effect under high impulse current is taken into account by introducing a nonlinear characteristic to represent the relationship between the electric field and current density in the ionization zone around the earthing electrode. On the basis of the nonlinear characteristic, the effective radius is evaluated for the equivalent earthing electrode. Then, the impulse earthing resistance can be calculated by substituting the effective radius into the formula of the low current earthing resistance. A comparison is also made between calculated and measured results to confirm the validity of the proposed algorithm.

**Keywords:** Earthing Resistance; Vertical Earthing Electrode; Average Potential Integral; Soil Ionization; Current Density;

## Introduction

The impulse earthing resistance of earthing electrodes is a very important factor which has to be taken into consideration in the lightning protection design of civil buildings and electric substations. An appropriate choice of protection measures against lightning overvoltage depends to a large degree on the knowledge of the values of the impulse earthing resistances of earthing electrodes. As the simplest form of earthing electrodes, vertical erathing electrodes are widely used in practical earthing systems. The impulse resistances of vertical earthing electrodes have been investigated for many years. Although a few algorithms were presented by different authors [1−3], their modeling for soil ionization effect under high impulse current is still a problem. These previous algorithms usually utilized a linear characteristic to represent the relationship between the electric field and current density in the ionization zone around the earthing electrode. In fact, this relationship has pronounced nonlinearity for typical kinds of soils in terms of the experimental investigation [4]. The nonlinear characteristic should therefore betaken into account in order to perform a more accurate calculation of the impulse earthing resistance.

The aim of this paper is to propose an efficient algorithm for calcu-lating the impulse resistances of vertical earthing electrodes. For the sake of engineering application, the inhomogeneity of soil is neglected. In the algorithm, the vertical earthing electrode is divided into a series of segments and the average potential method is employed to calculate the resistance matrix. The formula for the low current earthing resistance is derived from the resistance matrix. The soil ionization effect under high impulse current is further considered by using a nonlinear characteristic to represent the relationship between the electric field and current density in the ionization zone. With a simplified treatment made for the ionization zone, the effective radius of the equivalent earthing electrode is evaluated from the nonlinear characteristic. The impulse earthing resistance is then obtained by substituting the effective radius into the formula of the low current earthing resistance. Validity of the proposed algorithm has been verified by comparing calculated and measured results.

## Low Current Earthing resistance

Consider a vertical earthing electrode in a homogeneous soil with the resistivity ρ, as shown in Figure 1. When a low current I diffuses from it into the soil, the soil ionization cannot be caused. According to the electromagnetic field theory, the presence of the interface between

**Figure 1:** Vertical earthing

air and soil can be taken into account by installing the image electrode[5], as illustrated in Figure 2. The image electrode depicted by the dotted line is installed above the earth surface; meanwhile the actual and image electrodes are symmetrical about the earth surface. In consideration of

**Figure 2:** Division of actual and image electrodes into N segments.

the diffused current distribution along the electrode length, the actual electrode and its image are divided into N segments, respectively (see Figure 2）. At an arbitrary point on the surface of jth actual segment, as shown in Figure 3, the potential generated by the currents of jth actual segment and its image segment (j′ th segment) is calculated as [5]

${\varphi}_{\text{j}}({z}_{s})=\frac{\text{\rho}}{4\text{\pi}}\left[{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}\text{}\frac{{\tau}_{\text{j}}dz}{\sqrt{{(z-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}+{\displaystyle {\int}_{{\text{z}}_{\text{1}}}^{\text{}{z}_{2}}\frac{{{\tau}^{\prime}}_{\text{j}}\text{}d{z}^{\prime}}{\sqrt{{({z}^{\prime}-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}}}\right]\text{(1)}$

where τ_{j} and τ_{j}^{′} are the linear current densities of jth segment and its image, respectively. As the segment is rather short, the diffused current distribution on the segment is approximately considered to be uniform. Letting I_{j} and I_{j}^{′} denote the diffused currents of jth segment and its image, respectively,

**Figure 3:** Sketch for calculating self resistance of jth segment.

τ_{j} and τ_{j}^{′} become τ_{j}= I_{j}/(z_{2}-z_{1}) and τ_{j}^{′}= I_{j}^{′}/(z_{2}-z_{1}). Owing to the image symmetry, we have I_{j}= I_{j}^{′}. Thus, (2.1) is rewritten as

${\varphi}_{\text{j}}({z}_{\text{s}})=\frac{\text{\rho}{I}_{\text{j}}}{4{\text{\pi (z}}_{\text{2}}-{z}_{1})}\left[{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}\text{}\frac{dz}{\sqrt{{(z-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}+{\displaystyle {\int}_{{\text{z}}_{\text{1}}}^{\text{}{z}_{2}}\frac{d{z}^{\prime}}{\sqrt{{({z}^{\prime}-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}}}\right]\text{(2)}$

By taking the integral average of φ_{j}(z_{s}) over the length of jth segment, the average potential is obtained as

$\begin{array}{l}{\varphi}_{\text{j}}^{\text{a}}=\frac{1}{{\text{(z}}_{\text{2}}-{z}_{1})}{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}\text{}{\varphi}_{\text{j}}({z}_{\text{s}})dz}\\ \text{}=\frac{\text{\rho}{I}_{\text{j}}}{4{\text{\pi (z}}_{\text{2}}-{z}_{1}{)}^{2}}\left[{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}\text{}\frac{dzd{z}_{\text{s}}}{\sqrt{{(z-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}+{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}{\displaystyle {\int}_{{\text{z}}_{\text{1}}}^{\text{}{z}_{2}}\frac{d{z}^{\prime}d{z}_{\text{s}}}{\sqrt{{({z}^{\prime}-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}}}}}\right]\\ \text{}=\frac{\text{\rho}{r}_{\text{0}}{I}_{\text{j}}}{4{\text{\pi (z}}_{\text{2}}-{z}_{1}{)}^{2}}\left[2+2\Psi \left(\frac{{z}_{2}-{z}_{1}}{{r}_{0}}\right)-2\Psi \left(\frac{{z}_{1}+{z}_{2}}{{r}_{0}}\right)+\Psi \left(\frac{2{z}_{1}}{{r}_{0}}\right)+\Psi \left(\frac{2{z}_{2}}{{r}_{0}}\right)\right]\text{(3)}\end{array}$

where

$\Psi \left(\xi \right)=\xi {\text{sinh}}^{-\text{1}}\xi -\sqrt{1+{\xi}^{2}}$The self resistance of jth segment is given as

${R}_{\text{jj}}=\frac{{\varphi}_{\text{j}}^{\text{a}}}{{I}_{\text{j}}}\text{(4)}$

In a similar manner, the average potential generated by the currents of kth segment and its image on the surface of jth segment, as shown in Figure 4, is found as

$\begin{array}{l}{\varphi}_{\text{jk}}^{\text{a}}=\frac{\text{\rho}{I}_{\text{k}}}{4{\text{\pi (z}}_{\text{2}}-{z}_{1})\left({\text{z}}_{\text{4}}-{z}_{3}\right)}\left[{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{4}}}^{\text{}-{z}_{3}}\text{}\frac{dzd{z}_{\text{s}}}{\sqrt{{(z-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}+{\displaystyle {\int}_{\text{}-{\text{z}}_{\text{2}}}^{\text{}-{z}_{1}}{\displaystyle {\int}_{{\text{z}}_{\text{3}}}^{\text{}{z}_{4}}\frac{d{z}^{\prime}d{z}_{\text{s}}}{\sqrt{{({z}^{\prime}-{z}_{\text{s}})}^{2}+{r}_{0}^{2}}}}}}}\right]\\ \text{}=\frac{\text{\rho}{r}_{\text{0}}{I}_{\text{k}}}{4{\text{\pi (z}}_{\text{2}}-{z}_{1})\left({\text{z}}_{\text{4}}-{z}_{3}\right)}\left[\Psi \left(\frac{{z}_{3}-{z}_{2}}{{r}_{0}}\right)-\Psi \left(\frac{{z}_{3}-{z}_{1}}{{r}_{0}}\right)+\Psi \left(\frac{{z}_{4}-{z}_{1}}{{r}_{0}}\right)-\Psi \left(\frac{{z}_{3}-{z}_{1}}{{r}_{0}}\right)\right.\\ \text{}\left.\Psi \left(\frac{{z}_{2}+{z}_{3}}{{r}_{0}}\right)-\Psi \left(\frac{{z}_{1}+{z}_{3}}{{r}_{0}}\right)+\Psi \left(\frac{{z}_{1}+{z}_{4}}{{r}_{0}}\right)-\Psi \left(\frac{{z}_{4}+{z}_{2}}{{r}_{0}}\right)\right]\text{(5)}\end{array}$

where I_{k} is the current of kth segment. The mutual resistance between jth and kth segments is given as

${R}_{\text{jk}}=\frac{{\varphi}_{\text{jk}}^{\text{a}}}{{I}_{\text{k}}}\text{(6)}$

**Figure 4:** Sketch for calculating mutual resistance between jth and kth segments.

By taking the self and mutual resistances as the matrix elements, the resistance matrix of the vertical earthing electrode is formed as

$R=\left[\begin{array}{cccc}{R}_{11}& {R}_{12}& \cdot \cdot \cdot & {R}_{1N}\\ {R}_{21}& R{}_{22}& \cdot \cdot \cdot & {R}_{2N}\\ \vdots & \vdots & \vdots & \vdots \\ {R}_{N1}& {R}_{N2}& \cdot \cdot \cdot & {R}_{NN}\end{array}\right]\text{(7)}$Using **R** gives the relationship between voltage and current

$\left[\begin{array}{c}{U}_{1}\\ {U}_{2}\\ \vdots \\ {U}_{N}\end{array}\right]=\left[\begin{array}{cccc}{R}_{11}& {R}_{12}& \cdot \cdot \cdot & {R}_{1N}\\ {R}_{21}& R{}_{22}& \cdot \cdot \cdot & {R}_{2N}\\ \vdots & \vdots & \vdots & \vdots \\ {R}_{N1}& {R}_{N2}& \cdot \cdot \cdot & {R}_{NN}\end{array}\right]\text{}\left[\begin{array}{c}{I}_{1}\\ {I}_{2}\\ \vdots \\ {I}_{N}\end{array}\right]\text{(8)}$

where U_{j} (j=1, 2, ..., N) is the voltage of jth segment. As the conducing surfaces are equipotential, we have U_{1}=U_{2}= ··· =U_{N}=U [5]. Calculation of the inverse matrix R^{−1} gives the following expression

$\left[\begin{array}{c}{I}_{1}\\ {I}_{2}\\ \vdots \\ {I}_{N}\end{array}\right]=\left[\begin{array}{cccc}{G}_{11}& {G}_{12}& \cdot \cdot \cdot & {G}_{1N}\\ {G}_{21}& G{}_{22}& \cdot \cdot \cdot & {G}_{2N}\\ \vdots & \vdots & \vdots & \vdots \\ {G}_{N1}& {G}_{N2}& \cdot \cdot \cdot & {G}_{NN}\end{array}\right]\text{}\left[\begin{array}{c}U\\ U\\ \vdots \\ U\end{array}\right]=\left[\begin{array}{c}{\displaystyle \sum _{k=1}^{N}{G}_{1k}}\\ {\displaystyle \sum _{k=1}^{N}{G}_{2k}}\\ \vdots \\ {\displaystyle \sum _{k=1}^{N}{G}_{Nk}}\end{array}\right]U\text{(9)}$

where G_{jk} (j, k=1, 2, ..., N) is the element of inverse matrix R^{−1}. In terms of (9), the total current is found as

$I={\displaystyle \sum _{j=1}^{N}{I}_{\text{j}}}=U{\displaystyle \sum _{j=1}^{N}{\displaystyle \sum _{k=1}^{N}{G}_{\text{jk}}}}\text{(10)}$

As a result, the earthing resistance can be obtained as

${R}_{\text{e}}=\frac{U}{I}=\frac{1}{{\displaystyle \sum _{j=1}^{N}{\displaystyle \sum _{k=1}^{N}{G}_{\text{jk}}}}}\text{(11)}$

As shown in Figure 1, a vertical earthing electrode of r_{0}=0.015m is considered here. The earthing electrode is made of plain carbon steel (con-ductivity σ=0.5×10^{7} S•m^{-1}). The soil is kind of clay. Its resistivity is ρ=42Ω·m and water content is about 20%. The earthing resistances calculated from (11) are given in Figure 5. The corresponding measured values are given simultaneously in Figure 5, which was obtained by the fall of potential method [6]. As indicated in Figure 5, the calculated values are close to the measured ones.

**Figure 5:** Calculated and measured values of low current earthing resistance.

## Impulse Earthing Resistance

For a high impulse current with crest value I, representative of lightning, the impulse resistance R_{ie} of the earthing electrode is defined as a ratio of the crest value U_{m} of the impulse potential on the earthing electrode to I, i.e. R_{ie}=U_{m}/I [7]. The impulse current can produces great current density and high electric field intensity near the earthing electrode. When the electric field intensity on the surface of the earthing electrode exceeds the critical value E_{c} of soil ionization gradient, the breakdown will occur. This process can be illustrated by Figure 6 [8].

**Figure 6:** Impulse breakdown of soil around a vertical earthing electrode, 1—arc zone, 2—streamer zone, 3—semiconductive zone, 4—constant resistivity zone.

As the current increases, streamers are developed and in turn arcs are generated. Within the streamer and arc zones, the resistivity decreases from its original value to a limit approaching conductor [8]. In addition, there is a semiconductive zone between the streamer zone and the constant resistivity zone. For the purpose of simplifying calculation, this process can be described by an equivalent model shown in Figure 7. In the equivalent model, the semiconductive zone is neglected since it is small.

**Figure 7:** Ionization zone, (a) non-uniform along electrode length, (b) uniform along electrode length.

The streamer and arc zones are equivalently modeled as an ionization zone [8,9]. As the vertical earthing electrodes used in the actual cases usually have a shorter length (less than 5m), the non-uniformity of the ionization zone along the electrode length, as shown in Figure 7(a), is not pronounced and may be neglected. Therefore, the ionization zone is approximately considered to be uniform along the electrode length, as shown in Figure 7(b) [8−10]. The border of theionizationzoneis delimitedbythecriticalelectric field value E_{c}. A nonlinear characteris-tic is introduced to represent the relationship between the electric field andcurrent density in theionizationzone. The nonlinear characteristic is given as [4]

$E=a{J}^{b}\text{(12)}$

where a and b are constants and J is the current density

$J=\frac{I}{2\text{\pi}{r}_{}l+2\text{\pi}{r}_{}^{2}}=\frac{I}{2\text{\pi}r\left(r+l\right)}\text{(13)}$

Substituting (12) into (13), the electric field intensity is rewritten as

$E=a\frac{{I}^{b}}{{\left(2\pi \right)}^{b}{r}^{b}{\left(r+l\right)}^{b}}\text{(14)}$On the border of the ionization zone, E should satisfy the boundary condition.

${\left.E\right|}_{r={r}_{i}}={E}_{\text{c}}\text{(15)}$

The values of a, b and E_{c} can be found in [4] for typical kinds of soils. Thus, the boundary radius of the ionization zone can be derived from (14) and (15)

${r}_{i}=\frac{-l+\sqrt{{l}^{2}+4Q}}{2}\text{(16)}$

where

$Q=\frac{I}{2\text{\pi}}{\left(\frac{a}{{E}_{\text{c}}}\right)}^{\frac{1}{b}}\text{(17)}$

The soil ionization is basically equivalent to an increase in the dimension of the earthing electrode, which can be taken into account by the effective radius r_{e} of the equivalent earthing electrode in constant resistivity zone, asshowninFigure 8. The voltage between the earthing electrode and the border of the ionization zone, as shown in Figure 8(a), is expressed as

${U}_{i}={\displaystyle {\int}_{\text{}{r}_{0}}^{\text{}{r}_{i}}\frac{a{I}^{b}}{{\left(2\text{\pi}\right)}^{b}}\frac{1}{{r}^{b}{\left(r+1\right)}^{b}}dr=\frac{a{I}^{b}}{{\left(2\text{\pi}\right)}^{b}}\text{}}W\text{(18)}$

where

$W={\displaystyle {\int}_{\text{}{r}_{0}}^{\text{}{r}_{i}}\frac{1}{{r}^{b}{\left(r+1\right)}^{b}}dr\text{}}\text{(19)}$On the other hand, the voltage between the equivalent earthing electrode and the border of the ionization zone, as shown in Figure 8(b), is expressed as

${U}_{\text{e}}={\displaystyle {\int}_{\text{}r\text{e}}^{\text{}{r}_{i}}\frac{\rho I}{2\text{\pi}r\left(r+l\right)}\text{}}dr=\frac{\rho I}{2\text{\pi}l}\mathrm{ln}\frac{{r}_{i}\left({r}_{\text{e}}+l\right)}{{r}_{\text{e}}\left({r}_{i}+l\right)}\text{(20)}$

**Figure 8:** Sketch for evaluation of the effective radius, 1—ionization zone, 2—constant resistivity zone, (a) ionization zone around the actual earthing electrode, (b) constant resistivity zone around the equivalent earthing electrode.

Because U_{i} and U_{e} must be equal, the effective radius of the equivalent eathing electrode can be determined from (18) and (20)

${r}_{\text{e}}=\frac{{r}_{i}l}{{r}_{i}\left(\lambda -1\right)+\lambda l}\text{(21)}$

where λ=exp(χ) and

$\chi =\frac{al{I}^{b-1}}{\rho {(2\text{\pi})}^{b-1}}W\text{(22)}$

As the integration for (19) is difficult to be evaluated analytically, a numerical solution is given to obtain W

$W\approx \frac{\Delta r}{2}{\displaystyle \sum _{n=1}^{N}\left[\frac{1}{{r}_{\text{n}-1}^{b}{\left({r}_{\text{n}-\text{1}}+l\right)}^{b}}+\frac{1}{{r}_{n}^{b}{\left({r}_{n}+l\right)}^{b}}\right]}\text{(23)}$

where Δr is the radial step and N=(r_{i}-r_{0})/Δr.

By replacing r_{0} with r_{e} in (11), the impulse earthing resistance can be obtained. To check the validity of the proposed algorithm, two numerical examples are given here. The data for the first example are: r_{0}=0.0127m, l=3.05m, ρ=87.2Ω·m, E_{c}=127kV/m, a=3094.6 and b=0.51 [4]. At different current crest values,the calculated impulse earthing resistances are shown in Figure 9.

**Figure 9:** Calculated and measured impulse earthing resistances.

The data for the second example are: r_{0}=0.025m, l=1m, ρ=43.5Ω·m, E_{c}=350kV/m, a=219.05 and b=0.82 [4]. The calculated crest potentials on the earthing electrode (U_{m}=R_{ie}I) are shown in Figure 10. Furthermore, the corresponding measured results [11,12] are given in Figures 9 and 10, respectively, for comparison. It can be seen from Figures 9 and 10 that a better agreement appears between calculated and measured results.

**Figure 10:** Calculated and measured crest potentials.

## Conclusions

An algorithm for calculating the impulse earthing resistances of vertical earthing electrodes has been proposed in this paper. It uses the average potential method to derive the formula of the low current earthing resistances. In the ionization zone, a nonlinear characteristic is intro-duced in representingtherelationship between electric field and current density and has the capability of giving an appropriate consideration for the soil ionization effect under high impulse current. On the basis of the nonlinear characteristic, a simplified procedure has been developed to evaluate the effective radius of the equivalent earthing electrode. The impulse earthing resistance can then be obtained by substituting the effective radius into the formula of the low current earthing resistance. The calculated results are compared whith the measured ones and an agreement appears between them, which shows the applicability of the proposed algorithm in practical lightning protection design.