**4. The closed form of the Green's function of a point source in a horizontal multilayered earth model and the QSCIM**

The closed form of the Green's function of a scalar monopole and vector dipole buried in horizontal multilayered earth model will be respectively introduced.

## **4.1. The closed form of the Green's function of a scalar point source in a horizontal multilayered earth model and the QSCIM**

The main task of simulation grounding system is calculating the element *Zjk* of matrix [**Zs**] in Eq. (11). When the earth model is considered as horizontal multilayered conductivity media, influences from the interface will be considered, which will lead to infinite integral about Bessel function associated with Green's function of a scalar monopole. However, the element *Zjk* of matrix [**Zs**], which includes the Green's function, can be fast calculated by using the QSCIM.

To perform the calculation, the corresponding Green's function of a scalar monopole must be defined first. The Green's function can be regarded as the SEP of one point located at any place produced by a scalar monopole with unit current *δ* in a horizontal multilayered conducting medium. For low frequency (50 or 60 Hz and higher harmonic wave) and limited size of the substation, the propagating effect of electromagnetic wave can be neglected, so <sup>398</sup> Computational and Numerical Simulations Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi-static Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model Based on Quasi-Static Complex Image Method http://dx.doi.org/10.5772/57049 399

**Figure 1.** The earth model

6 Computational and Numerical Simulationse

discussed next.

QSCIM.

2. The case of [**Xq**]=[**Zs**]: *Xi*,*<sup>j</sup>* = *Zi*,*j*.

*Zi*,*<sup>j</sup>* is the mutual impedance coefficient between a pair of conductor segments in the grounding system. The matrix above includes the conductive and capacitive effects of the

*<sup>G</sup>ϕ*(*r*¯*j*,*r*¯*i*) *dti*

*li*

*dtj lj*

(11)

4*πσ*¯1 1 *Rij* , so

*dtidtj* (12)

earth, and its elements are the mutual impedance coefficients *Zi*,*j*.

*Zi*,*<sup>j</sup>* =

*lj*

*li*

For an infinite homogeneous conductivity medium, one has *Gϕ*(*r*¯*j*,*r*¯*i*) = <sup>1</sup>

4*πσ*¯1*ljlj*

**4. The closed form of the Green's function of a point source in a**

**4.1. The closed form of the Green's function of a scalar point source in a**

*lj* 1 *Rij*

*li*

Note the medium surrounding the point current source here was considered as homogeneous and infinite. However, in any practical case, the earth is represented via a multilayered earth model. The QSCIM can be used to dealt with the multilayered earth model, this will be

The closed form of the Green's function of a scalar monopole and vector dipole buried in

The main task of simulation grounding system is calculating the element *Zjk* of matrix [**Zs**] in Eq. (11). When the earth model is considered as horizontal multilayered conductivity media, influences from the interface will be considered, which will lead to infinite integral about Bessel function associated with Green's function of a scalar monopole. However, the element *Zjk* of matrix [**Zs**], which includes the Green's function, can be fast calculated by using the

To perform the calculation, the corresponding Green's function of a scalar monopole must be defined first. The Green's function can be regarded as the SEP of one point located at any place produced by a scalar monopole with unit current *δ* in a horizontal multilayered conducting medium. For low frequency (50 or 60 Hz and higher harmonic wave) and limited size of the substation, the propagating effect of electromagnetic wave can be neglected, so

*Zi*,*<sup>j</sup>* <sup>=</sup> <sup>1</sup>

where *σ*¯1 = *σ*<sup>1</sup> + *jωε*1. Eq. (12) can be solved analytically [36].

**horizontal multilayered earth model and the QSCIM**

horizontal multilayered earth model will be respectively introduced.

**horizontal multilayered earth model and the QSCIM**

the electromagnetic field here can be regarded as quasi-static field, so the SEP *ϕ* satisfies the Poisson equation as:

$$
\nabla^2 \varphi\_j = -\frac{\delta(r, r')\delta(i\vec{j})}{\overline{\sigma\_i}} \tag{13}
$$

where *<sup>i</sup>*, *<sup>j</sup>* <sup>=</sup> 1, . . . , *Ns*, and *<sup>δ</sup>*(*r*,*r*′ ) is the Dirac delta function. *δ*(*ij*) is Kronecker's symbol, *σ<sup>i</sup>* = *σ<sup>i</sup>* + *jωε<sup>i</sup>* is the complex conductivity of the *i*th layer medium. Supposing the monopole is located at origin of the coordinate system, seen from Fig. 1.

If the monopole lies in infinite homogenous medium, its expression of Green's function in the spherical coordinate system is.

$$\varphi = \mathcal{G}(\vec{r}, \vec{r'}) = \frac{1}{4\pi \cdot \overline{\sigma} \cdot \mathcal{R}} \tag{14}$$

where *R* = |*r*¯ − ¯ *r*′| is the distance between the source point and field point. While its expression in the cylindrical coordinates system (*ρ*, *z*) is as follows:

$$\mathfrak{g} = \mathcal{G}(\rho, z) = \frac{1}{4\pi \cdot \overline{\sigma}} \int\_0^\infty e^{-\lambda|z|} J\_0(\lambda \rho) d\lambda \tag{15}$$

where *J*0(*λρ*) is the Bessel function of the first kind of order zero.

For horizontal multilayered earth model, for examples, two-layer earth model, in the cylindrical coordinate system, the general form of Green's functions can be expressed by [1]:

$$\varphi\_{10} = G\_{10}(\rho\_\prime z) = \frac{1}{4\pi \cdot \overline{\sigma}\_1} \int\_0^\infty \left[ A\_{10}(\lambda) e^{-\lambda \cdot z} + B\_{10}(\lambda) e^{+\lambda \cdot z} \right] J\_0(\lambda \rho) d\lambda \tag{16}$$

where *α<sup>n</sup>* and *β<sup>n</sup>* are constants to be determined by choosing sample points of function *f*(*λ*). Considering the Laplace transform of Bessel function of the first kind of order *v Jv*(*λρ*), we

Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model

<sup>−</sup>*c<sup>λ</sup> <sup>J</sup>*0(*λρ*)*d<sup>λ</sup>* <sup>=</sup> <sup>1</sup>

By employing Eq. (20) and Lipschitz integration, the closed form of the Green's function (Eq.

where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of

source, whose location is indicated by *Rni*,*i*=1∽<sup>4</sup> and amplitude is *αn*, which can be seen in Fig. 1. However, in Eq. (23) *α<sup>n</sup>* and *Rni* are usually complex numbers, so that this approach

The procedure for the closed form of the Green's function of a scalar monopole in arbitrary horizontal multilayered earth model is first to derive the specific solution of the Green's function in each layer based on the general expressions (for example, the general expressions of a scalar monopole in three layers media is given in Eq. (15), Eq. (16) and Eq. (17)), and the interface and infinite conditions of SEP of the monopole is used to decide the unknown constants. Then, by employing the MP method exponential series development and the

Once the closed form of Green's function of a scalar monopole in the horizontal multilayered earth model has been gotten, the mutual impedance coefficient Eq. (11) can be fast analytical

Lipschitz integration, the final expression of the Green's functions can be obtained.

)2] 1 (*c*<sup>2</sup> + *ρ*2)

1 2

− *k*01*k*<sup>12</sup> *Rn*<sup>2</sup>

+ *k*<sup>12</sup> *Rn*<sup>3</sup>

1

) and the field point is at (*ρ*, *<sup>z</sup>*), so *<sup>R</sup>*<sup>0</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*′

*<sup>R</sup>*<sup>0</sup> of Eq. (23) can be regarded as a image scalar monopole

<sup>2</sup> , *Rn*<sup>1</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′ <sup>−</sup> *zn*)2]

<sup>2</sup> , *Rn*<sup>4</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*′ <sup>+</sup> *zn*)2]

− *k*01*k*<sup>12</sup> *Rn*<sup>4</sup>

1

*ρ <sup>c</sup>* <sup>+</sup> *<sup>ρ</sup>*<sup>2</sup> <sup>+</sup> *<sup>c</sup>*<sup>2</sup> *<sup>v</sup>* <sup>1</sup>

(*c*<sup>2</sup> + *ρ*2)

1 2

Based on Quasi-Static Complex Image Method

http://dx.doi.org/10.5772/57049

(21)

401

(22)

)] (23)

<sup>2</sup> , *Rn*<sup>2</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* −

<sup>2</sup> , in which *zn* = 2*h* − *βn*.

)2] 1 2 ,

have [39]

<sup>∞</sup>

<sup>−</sup>*c<sup>λ</sup> Jv*(*λρ*)*d<sup>λ</sup>* <sup>=</sup>

<sup>∞</sup>

0 *e*

[ 1 *R*0 − *k*<sup>01</sup> *R*′ 0 + *N* ∑ *n*=1 *αn*( *k*2 <sup>01</sup>*k*<sup>12</sup> *Rn*<sup>1</sup>

<sup>0</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′

1

0 *e*

Set *v* = 0, so we have Lipschitz integration:

where *Re*(*v*) > −1 and *Re*(*c*) > 0.

(19)) can be given as following:

*G*11(*r*¯, ¯

in the same way, we have *<sup>R</sup>*′

is named as the QSCIM.

1

*<sup>z</sup>*′ <sup>−</sup> *zn*)2]

*<sup>r</sup>*′) = <sup>1</sup> 4*πσ*<sup>1</sup>

the earth, the source point is at (0, *<sup>z</sup>*′

It can be seen that each term except <sup>1</sup>

*G*<sup>10</sup> and *G*<sup>12</sup> can be similarly gotten.

calculated, which is just same as the Eq. (12).

<sup>2</sup> , *Rn*<sup>3</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′ <sup>+</sup> *zn*)2]

$$\varphi\_{11} = G\_{11}(\rho, z) = \frac{1}{4\pi \cdot \overline{\sigma}\_1} \int\_0^\infty \left[ e^{-\lambda|z|} + A\_{11}(\lambda)e^{-\lambda z} + B\_{11}(\lambda)e^{+\lambda z} \right] J\_0(\lambda \rho) d\lambda \tag{17}$$

$$\varphi\_{12} = G\_{12}(\rho, z) = \frac{1}{4\pi \cdot \overline{\sigma}\_1} \int\_0^\infty \left[ A\_{12}(\lambda) e^{-\lambda \cdot z} + B\_{12}(\lambda) e^{+\lambda \cdot z} \right] I\_0(\lambda \rho) d\lambda \tag{18}$$

where *G*10, *G*<sup>11</sup> and *G*<sup>12</sup> express the Green's function for the field point in the air, the top layer and bottom layer, respectively.

The six constants *<sup>A</sup>*<sup>10</sup> ∼ *<sup>B</sup>*<sup>12</sup> can be determined by employing the interface and infinite conditions of the SEP (*ϕ*<sup>10</sup> = *ϕ*11, *σ*<sup>0</sup> *∂ϕ*<sup>10</sup> *<sup>∂</sup><sup>z</sup>* = *σ*<sup>1</sup> *∂ϕ*<sup>11</sup> *<sup>∂</sup><sup>z</sup>* , *ϕ*<sup>11</sup> = *ϕ*12, *σ*<sup>1</sup> *∂ϕ*<sup>11</sup> *<sup>∂</sup><sup>z</sup>* = *σ*<sup>2</sup> *∂ϕ*<sup>12</sup> *<sup>∂</sup><sup>z</sup>* , *<sup>ϕ</sup>*10|*z*�→−<sup>∞</sup> <sup>=</sup> <sup>0</sup> and *<sup>ϕ</sup>*12|*z*�→+<sup>∞</sup> = 0).

Consequently, the expression of *G*<sup>11</sup> can be given as follows:

$$G\_{11}(\rho, z) = \frac{1}{4\pi\overline{\sigma}\_1} \int\_0^\infty [e^{-k\_\rho|z|} - k\_{01}e^{-k\_\rho(2z'+z)} + f(k\_\rho)(k\_{01}^2 k\_{12} e^{-k\_\rho(2h+2z'+z)})$$

$$[-k\_{01}k\_{12}e^{-k\_{\rho}(2\hbar+z)} + k\_{12}e^{-k\_{\rho}(2\hbar-2z'-z)} - k\_{01}k\_{12}e^{-k\_{\rho}(2\hbar-z)})]l\_{0}(k\_{\rho}\rho)dk\_{\rho} \tag{19}$$

where *f*(*λ*) = *<sup>k</sup>*<sup>12</sup> (1+*k*01*k*12*e*−2*λ*·*<sup>h</sup>* ) ,*k*<sup>01</sup> <sup>=</sup> (*σ*0−*σ*1) (*σ*0+*σ*1), *<sup>k</sup>*<sup>12</sup> <sup>=</sup> (*σ*1−*σ*2) (*σ*1+*σ*2), h is the thickness of the top layer, *σ*0,*σ*<sup>1</sup> and *σ*<sup>2</sup> are the complex conductivity of air and two layers, respectively.

Now this *f*(*λ*) can be developed into an exponential series with finite terms by MP method [37] instead of Maclaurin's series to avoid verbose calculation [38]:

$$f(\lambda) = \sum\_{n=1}^{N} \alpha\_n \cdot e^{\mathcal{G}\_n \lambda} \tag{20}$$

<sup>400</sup> Computational and Numerical Simulations Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi-static Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model Based on Quasi-Static Complex Image Method http://dx.doi.org/10.5772/57049 401

> where *α<sup>n</sup>* and *β<sup>n</sup>* are constants to be determined by choosing sample points of function *f*(*λ*). Considering the Laplace transform of Bessel function of the first kind of order *v Jv*(*λρ*), we have [39]

$$\int\_0^\infty e^{-c\lambda} f\_\upsilon(\lambda \rho) d\lambda = \left[\frac{\rho}{c + \sqrt{\rho^2 + c^2}}\right]^\upsilon \frac{1}{(c^2 + \rho^2)^{\frac{1}{2}}} \tag{21}$$

where *Re*(*v*) > −1 and *Re*(*c*) > 0.

8 Computational and Numerical Simulationse

*<sup>ϕ</sup>*<sup>10</sup> <sup>=</sup> *<sup>G</sup>*10(*ρ*, *<sup>z</sup>*) = <sup>1</sup>

*<sup>ϕ</sup>*<sup>12</sup> <sup>=</sup> *<sup>G</sup>*12(*ρ*, *<sup>z</sup>*) = <sup>1</sup>

*<sup>ϕ</sup>*<sup>11</sup> <sup>=</sup> *<sup>G</sup>*11(*ρ*, *<sup>z</sup>*) = <sup>1</sup>

layer and bottom layer, respectively.

conditions of the SEP (*ϕ*<sup>10</sup> = *ϕ*11, *σ*<sup>0</sup>

*<sup>G</sup>*11(*ρ*, *<sup>z</sup>*) = <sup>1</sup>

− *<sup>k</sup>*01*k*12*<sup>e</sup>*

where *f*(*λ*) = *<sup>k</sup>*<sup>12</sup>

4*πσ*<sup>1</sup>

and *<sup>ϕ</sup>*12|*z*�→+<sup>∞</sup> = 0).

<sup>4</sup>*<sup>π</sup>* · *<sup>σ</sup>*<sup>1</sup>

<sup>∞</sup>

0

<sup>4</sup>*<sup>π</sup>* · *<sup>σ</sup>*<sup>1</sup>

Consequently, the expression of *G*<sup>11</sup> can be given as follows:

<sup>∞</sup>

0 [*e*

<sup>−</sup>*k<sup>ρ</sup>* (2*h*+*z*) <sup>+</sup> *<sup>k</sup>*12*<sup>e</sup>*

(1+*k*01*k*12*e*−2*λ*·*<sup>h</sup>* ) ,*k*<sup>01</sup> <sup>=</sup> (*σ*0−*σ*1)

[37] instead of Maclaurin's series to avoid verbose calculation [38]:

<sup>4</sup>*<sup>π</sup>* · *<sup>σ</sup>*<sup>1</sup>

<sup>∞</sup>

0

 *e*

<sup>∞</sup>

0

*∂ϕ*<sup>10</sup> *<sup>∂</sup><sup>z</sup>* = *σ*<sup>1</sup>

<sup>−</sup>*k<sup>ρ</sup>* <sup>|</sup>*z*<sup>|</sup> <sup>−</sup> *<sup>k</sup>*01*<sup>e</sup>*

−*k<sup>ρ</sup>* (2*h*−2*z*′

*σ*0,*σ*<sup>1</sup> and *σ*<sup>2</sup> are the complex conductivity of air and two layers, respectively.

*f*(*λ*) =

[1]:

For horizontal multilayered earth model, for examples, two-layer earth model, in the cylindrical coordinate system, the general form of Green's functions can be expressed by

*A*10(*λ*)*e*

<sup>−</sup>*λ*|*z*<sup>|</sup> <sup>+</sup> *<sup>A</sup>*11(*λ*)*<sup>e</sup>*

*A*12(*λ*)*e*

where *G*10, *G*<sup>11</sup> and *G*<sup>12</sup> express the Green's function for the field point in the air, the top

The six constants *<sup>A</sup>*<sup>10</sup> ∼ *<sup>B</sup>*<sup>12</sup> can be determined by employing the interface and infinite

*∂ϕ*<sup>11</sup>

−*k<sup>ρ</sup>* (2*z*′

(*σ*0+*σ*1), *<sup>k</sup>*<sup>12</sup> <sup>=</sup> (*σ*1−*σ*2)

Now this *f*(*λ*) can be developed into an exponential series with finite terms by MP method

*N* ∑ *n*=1 <sup>−</sup>*λ*·*<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*10(*λ*)*<sup>e</sup>*

<sup>−</sup>*λ*·*<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*12(*λ*)*<sup>e</sup>*

*<sup>∂</sup><sup>z</sup>* , *ϕ*<sup>11</sup> = *ϕ*12, *σ*<sup>1</sup>

<sup>+</sup>*z*) + *f*(*kρ*)(*k*<sup>2</sup>

<sup>−</sup>*z*) <sup>−</sup> *<sup>k</sup>*01*k*12*<sup>e</sup>*

<sup>−</sup>*λ<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*11(*λ*)*<sup>e</sup>*

+*λ*·*z* 

+*λ*·*z* 

*∂ϕ*<sup>11</sup> *<sup>∂</sup><sup>z</sup>* = *σ*<sup>2</sup>

<sup>01</sup>*k*12*e*

−*k<sup>ρ</sup>* (2*h*−*z*)

*∂ϕ*<sup>12</sup>

−*k<sup>ρ</sup>* (2*h*+2*z*′

(*σ*1+*σ*2), h is the thickness of the top layer,

*α<sup>n</sup>* · *eβn<sup>λ</sup>* (20)

+*z*)

)]*J*0(*kρρ*)*dk<sup>ρ</sup>* (19)

+*λz* 

*J*0(*λρ*)*dλ* (16)

*J*0(*λρ*)*dλ* (17)

*J*0(*λρ*)*dλ* (18)

*<sup>∂</sup><sup>z</sup>* , *<sup>ϕ</sup>*10|*z*�→−<sup>∞</sup> <sup>=</sup> <sup>0</sup>

Set *v* = 0, so we have Lipschitz integration:

$$\int\_0^\infty e^{-c\lambda} f\_0(\lambda \rho) d\lambda = \frac{1}{(c^2 + \rho^2)^{\frac{1}{2}}} \tag{22}$$

By employing Eq. (20) and Lipschitz integration, the closed form of the Green's function (Eq. (19)) can be given as following:

$$G\_{11}(\vec{r}, \vec{r'}) = \frac{1}{4\pi\tilde{\sigma}\_1} [\frac{1}{R\_0} - \frac{k\_{01}}{R\_0'} + \sum\_{n=1}^{N} a\_n (\frac{k\_{01}^2 k\_{12}}{R\_{n1}} - \frac{k\_{01} k\_{12}}{R\_{n2}} + \frac{k\_{12}}{R\_{n3}} - \frac{k\_{01} k\_{12}}{R\_{n4}})] \tag{23}$$

where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of the earth, the source point is at (0, *<sup>z</sup>*′ ) and the field point is at (*ρ*, *<sup>z</sup>*), so *<sup>R</sup>*<sup>0</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*′ )2] 1 2 , in the same way, we have *<sup>R</sup>*′ <sup>0</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′ )2] 1 <sup>2</sup> , *Rn*<sup>1</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′ <sup>−</sup> *zn*)2] 1 <sup>2</sup> , *Rn*<sup>2</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* − *<sup>z</sup>*′ <sup>−</sup> *zn*)2] 1 <sup>2</sup> , *Rn*<sup>3</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′ <sup>+</sup> *zn*)2] 1 <sup>2</sup> , *Rn*<sup>4</sup> = [*ρ*<sup>2</sup> + (*<sup>z</sup>* <sup>−</sup> *<sup>z</sup>*′ <sup>+</sup> *zn*)2] 1 <sup>2</sup> , in which *zn* = 2*h* − *βn*.

It can be seen that each term except <sup>1</sup> *<sup>R</sup>*<sup>0</sup> of Eq. (23) can be regarded as a image scalar monopole source, whose location is indicated by *Rni*,*i*=1∽<sup>4</sup> and amplitude is *αn*, which can be seen in Fig. 1. However, in Eq. (23) *α<sup>n</sup>* and *Rni* are usually complex numbers, so that this approach is named as the QSCIM.

*G*<sup>10</sup> and *G*<sup>12</sup> can be similarly gotten.

The procedure for the closed form of the Green's function of a scalar monopole in arbitrary horizontal multilayered earth model is first to derive the specific solution of the Green's function in each layer based on the general expressions (for example, the general expressions of a scalar monopole in three layers media is given in Eq. (15), Eq. (16) and Eq. (17)), and the interface and infinite conditions of SEP of the monopole is used to decide the unknown constants. Then, by employing the MP method exponential series development and the Lipschitz integration, the final expression of the Green's functions can be obtained.

Once the closed form of Green's function of a scalar monopole in the horizontal multilayered earth model has been gotten, the mutual impedance coefficient Eq. (11) can be fast analytical calculated, which is just same as the Eq. (12).

### **4.2. The closed form of the Green's function of a vector point source in a horizontal multilayered earth model and the QSCIM**

The another important task of simulation grounding system is calculating the element *Mjk* of matrix [**Zb**] in Eq. (9). When the earth model is considered as horizontal multilayered conductivity media, and influences from the interfaces on mutual induction between two conductors should be considered, which also lead to infinite integral about Bessel function in the Green's function of a vector dipole, the element *Mjk* of matrix [**Zb**] includes the Green's function, which can also be fast calculated by using the QSCIM.

Like Green's function of the scalar monopole, to perform the calculation of a vector dipole in horizontal multilayered earth model, its corresponding Green's function should also be defined first. The Green's function can be regarded as the VMP A of a point at any place produced by a vector dipole with unit current *δ* within a horizontal multilayered medium, whose electromagnetic field can also be regarded as quasi-static field, so the VMP A also satisfies the Poisson equation as:

$$
\nabla^2 A\_j = -\mu \bar{\delta}(r, r') \delta(i\jmath) \tag{24}
$$

*Ap*

*<sup>p</sup>*<sup>11</sup> <sup>=</sup> *<sup>G</sup><sup>p</sup>*

*Ap*

*Ap*<sup>10</sup> , *<sup>G</sup><sup>p</sup>*

*Gz Az*<sup>11</sup> (*r*¯, ¯

the earth. *<sup>R</sup>*0, *<sup>R</sup>*′

*Az <sup>x</sup>*<sup>10</sup> <sup>=</sup> *<sup>G</sup><sup>z</sup>*

*Az <sup>x</sup>*<sup>11</sup> <sup>=</sup> *<sup>G</sup><sup>z</sup>*

horizontal dipole case, *p* = *x*.

*<sup>p</sup>*<sup>12</sup> <sup>=</sup> *<sup>G</sup><sup>p</sup>*

integration, the final expression of *G<sup>p</sup>*

earth model can be expressed by [1]:

*<sup>r</sup>*′) = *<sup>µ</sup>*

<sup>4</sup>*<sup>π</sup>* [ <sup>1</sup> *R*0 + *k*<sup>01</sup> *R*′ 0 − *N* ∑ *n*=1 *αn*( *k*2 <sup>01</sup>*k*<sup>12</sup> *Rn*<sup>1</sup>

*Ax*<sup>10</sup> (*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

*Ax*<sup>11</sup> (*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

*Gx*

<sup>0</sup>, *Rni*,*i*=1∽<sup>4</sup> are just same as Eq. (23). For horizontal dipole, there is a Green's function for z component *G<sup>z</sup>*

> 4*π ∂ ∂x*

4*π ∂ ∂x*

*Ap*<sup>11</sup> and *<sup>G</sup><sup>p</sup>*

*Ap*

where *G<sup>p</sup>*

*<sup>p</sup>*<sup>10</sup> <sup>=</sup> *<sup>G</sup><sup>p</sup>*

*Ap*<sup>10</sup> (*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

4*π*

*Ap*<sup>12</sup> (*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

*Ap*<sup>11</sup> (*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

4*π*

<sup>∞</sup>

 *e*

<sup>∞</sup>

0

0

4*π*

<sup>∞</sup>

*A*10(*λ*)*e*

Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model

<sup>−</sup>*λ*·|*z*<sup>|</sup> <sup>+</sup> *<sup>A</sup>*11(*λ*)*<sup>e</sup>*

*A*12(*λ*)*e*

in the air, the top layer and bottom layer, respectively; for vertical dipole case, *p* = *z*, for

Just like the deduced procedure for Green's function of a scalar monopole, the six constants *<sup>A</sup>*<sup>10</sup> ∼ *<sup>B</sup>*<sup>12</sup> can also be decided by employing the interface and infinite conditions of the dipole's VMP. For vertical dipole, there is a *f*(*λ*), which can also be developed into an exponential series with finite terms by the MP method. Last by applying the Lipschitz

*Ax*<sup>11</sup> (*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of

coordinate system, the general form of the Green's function for z component in the two-layer

*A*10(*λ*)*e*

*A*11(*λ*)*e*

<sup>∞</sup>

0

<sup>∞</sup>

0

*Ap*<sup>11</sup> can be given as follows.

+

4*π R*<sup>0</sup>

*k*01*k*<sup>12</sup> *Rn*<sup>2</sup>

<sup>−</sup>*λ*·*<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*10(*λ*)*<sup>e</sup>*

<sup>−</sup>*λ<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*11(*λ*)*<sup>e</sup>*

− *k*<sup>12</sup> *Rn*<sup>3</sup> +

> +*λ*·*z*

+*λz*  *J*0(*λρ*)

*J*0(*λρ*)

*k*01*k*<sup>12</sup> *Rn*<sup>4</sup>

*Ax*<sup>11</sup> left, in the cylindrical

*dλ*

*dλ*

*<sup>λ</sup>* (32)

*<sup>λ</sup>* (33)

)] (30)

(31)

<sup>−</sup>*λ*·*<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*10(*λ*)*<sup>e</sup>*

<sup>−</sup>*λ*·*<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*12(*λ*)*<sup>e</sup>*

*Ap*<sup>12</sup> express the Green's function of the dipole for the field point

<sup>−</sup>*λ*·*<sup>z</sup>* <sup>+</sup> *<sup>B</sup>*11(*λ*)*<sup>e</sup>*

+*λ*·*z* 

+*λ*·*z* 

+*λ*·*z* 

Based on Quasi-Static Complex Image Method

*J*0(*λρ*)*dλ* (27)

http://dx.doi.org/10.5772/57049

403

*J*0(*λρ*)*dλ* (28)

*J*0(*λρ*)*dλ* (29)

0

where *i*, *j* = 1, . . . , *Ns*; *µ* is the permeability of the earth medium. Supposing the dipole is located at origin of the coordinate system, seen from Fig. 1

If the vector dipole is lying in homogenous infinite medium, its expression of Green's function in the spherical coordinate system is.

$$A = G\_A(r, r') = \frac{\mu}{4\pi \cdot R} \tag{25}$$

And its expression in the cylindrical coordinates system (*ρ*, *z*) is as follows:

$$A = G\_A(\rho, z) = \frac{\mu}{4\pi} \int\_0^\infty e^{-\lambda|z|} f\_0(\lambda \rho) d\lambda \tag{26}$$

Not like scalar monopole source, for horizontal multilayered earth model, vector dipole source must be considered into two cases, which are vertical and horizontal dipoles, respectively. This is because any placed dipole in horizontal multilayered earth model can be decomposed into horizontal and vertical components. For vertical dipole case, its Green's function can be defined as *G<sup>z</sup> Az* ; For horizontal dipole case, it own two components, which are horizontal x or y component and vertical z component, so its Green's function can be defined as *G<sup>x</sup> Ax* and *<sup>G</sup><sup>z</sup> Ax* , respectively.

We also take the two-layer earth model as an example, the dipole lies in the earth model. In the cylindrical coordinate system, the general form of Green's functions for the vertical dipole *G<sup>z</sup> Az* or horizontal dipole *<sup>G</sup><sup>x</sup> Ax* can be expressed by [1]:

<sup>402</sup> Computational and Numerical Simulations Numerical calculation for lightning response to grounding systems buried in horizontal multilayered earth model based on quasi-static Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model Based on Quasi-Static Complex Image Method http://dx.doi.org/10.5772/57049 403

10 Computational and Numerical Simulationse

satisfies the Poisson equation as:

function can be defined as *G<sup>z</sup>*

*Ax* and *<sup>G</sup><sup>z</sup>*

*Ax*

*Az* or horizontal dipole *<sup>G</sup><sup>x</sup>*

defined as *G<sup>x</sup>*

dipole *G<sup>z</sup>*

**4.2. The closed form of the Green's function of a vector point source in a**

The another important task of simulation grounding system is calculating the element *Mjk* of matrix [**Zb**] in Eq. (9). When the earth model is considered as horizontal multilayered conductivity media, and influences from the interfaces on mutual induction between two conductors should be considered, which also lead to infinite integral about Bessel function in the Green's function of a vector dipole, the element *Mjk* of matrix [**Zb**] includes the Green's

Like Green's function of the scalar monopole, to perform the calculation of a vector dipole in horizontal multilayered earth model, its corresponding Green's function should also be defined first. The Green's function can be regarded as the VMP A of a point at any place produced by a vector dipole with unit current *δ* within a horizontal multilayered medium, whose electromagnetic field can also be regarded as quasi-static field, so the VMP A also

*<sup>δ</sup>*(*r*,*r*′

where *i*, *j* = 1, . . . , *Ns*; *µ* is the permeability of the earth medium. Supposing the dipole is

If the vector dipole is lying in homogenous infinite medium, its expression of Green's

) = *<sup>µ</sup>*

)*δ*(*ij*) (24)

<sup>4</sup>*<sup>π</sup>* · *<sup>R</sup>* (25)

*J*0(*λρ*)*dλ* (26)

; For horizontal dipole case, it own two components, which

∇<sup>2</sup>*Aj* = −*<sup>µ</sup>* ¯

*<sup>A</sup>* <sup>=</sup> *GA*(*r*,*r*′

4*π*

<sup>∞</sup>

0 *e* −*λ*|*z*|

Not like scalar monopole source, for horizontal multilayered earth model, vector dipole source must be considered into two cases, which are vertical and horizontal dipoles, respectively. This is because any placed dipole in horizontal multilayered earth model can be decomposed into horizontal and vertical components. For vertical dipole case, its Green's

are horizontal x or y component and vertical z component, so its Green's function can be

We also take the two-layer earth model as an example, the dipole lies in the earth model. In the cylindrical coordinate system, the general form of Green's functions for the vertical

*Ax* can be expressed by [1]:

And its expression in the cylindrical coordinates system (*ρ*, *z*) is as follows:

*<sup>A</sup>* <sup>=</sup> *GA*(*ρ*, *<sup>z</sup>*) = *<sup>µ</sup>*

*Az*

, respectively.

**horizontal multilayered earth model and the QSCIM**

function, which can also be fast calculated by using the QSCIM.

located at origin of the coordinate system, seen from Fig. 1

function in the spherical coordinate system is.

$$A\_{p10}^p = G\_{A\_{p10}}^p(\rho, z) = \frac{\mu}{4\pi} \int\_0^\infty \left[ A\_{10}(\lambda) e^{-\lambda \cdot z} + B\_{10}(\lambda) e^{+\lambda \cdot z} \right] f\_0(\lambda \rho) d\lambda \tag{27}$$

$$A\_{p11}^p = G\_{A\_{p1}}^p(\rho, z) = \frac{\mu}{4\pi} \int\_0^\infty \left[ e^{-\lambda \cdot |z|} + A\_{11}(\lambda) e^{-\lambda \cdot z} + B\_{11}(\lambda) e^{+\lambda \cdot z} \right] f\_0(\lambda \rho) d\lambda \tag{28}$$

$$A\_{p12}^p = G\_{A\_{p12}}^p(\rho, z) = \frac{\mu}{4\pi} \int\_0^\infty \left[ A\_{12}(\lambda) e^{-\lambda \cdot z} + B\_{12}(\lambda) e^{+\lambda \cdot z} \right] f\_0(\lambda \rho) d\lambda \tag{29}$$

where *G<sup>p</sup> Ap*<sup>10</sup> , *<sup>G</sup><sup>p</sup> Ap*<sup>11</sup> and *<sup>G</sup><sup>p</sup> Ap*<sup>12</sup> express the Green's function of the dipole for the field point in the air, the top layer and bottom layer, respectively; for vertical dipole case, *p* = *z*, for horizontal dipole case, *p* = *x*.

Just like the deduced procedure for Green's function of a scalar monopole, the six constants *<sup>A</sup>*<sup>10</sup> ∼ *<sup>B</sup>*<sup>12</sup> can also be decided by employing the interface and infinite conditions of the dipole's VMP. For vertical dipole, there is a *f*(*λ*), which can also be developed into an exponential series with finite terms by the MP method. Last by applying the Lipschitz integration, the final expression of *G<sup>p</sup> Ap*<sup>11</sup> can be given as follows.

$$\mathcal{G}\_{A\_{211}}^{\mathcal{Z}}(\bar{r}, \bar{r'}) = \frac{\mu}{4\pi} [\frac{1}{R\_0} + \frac{k\_{01}}{R\_0'} - \sum\_{n=1}^{N} a\_n (\frac{k\_{01}^2 k\_{12}}{R\_{n1}} + \frac{k\_{01} k\_{12}}{R\_{n2}} - \frac{k\_{12}}{R\_{n3}} + \frac{k\_{01} k\_{12}}{R\_{n4}})] \tag{30}$$

$$G\_{A\_{\mathcal{X}1}}^{\chi}(\rho, z) = \frac{\mu}{4\pi \ R\_0} \tag{31}$$

where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of the earth. *<sup>R</sup>*0, *<sup>R</sup>*′ <sup>0</sup>, *Rni*,*i*=1∽<sup>4</sup> are just same as Eq. (23).

For horizontal dipole, there is a Green's function for z component *G<sup>z</sup> Ax*<sup>11</sup> left, in the cylindrical coordinate system, the general form of the Green's function for z component in the two-layer earth model can be expressed by [1]:

$$A\_{\mathbf{x}10}^{\mathbf{z}} = G\_{A\_{\mathbf{x}10}}^{\mathbf{z}}(\rho, \mathbf{z}) = \frac{\mu}{4\pi} \frac{\partial}{\partial \mathbf{x}} \int\_{0}^{\infty} \left[ A\_{10}(\lambda) e^{-\lambda \cdot \mathbf{z}} + B\_{10}(\lambda) e^{+\lambda \cdot \mathbf{z}} \right] J\_{0}(\lambda \rho) \frac{d\lambda}{\lambda} \tag{32}$$

$$A\_{\rm x11}^{\rm z} = G\_{A\_{\rm x11}}^{\rm z}(\rho, z) = \frac{\mu}{4\pi} \frac{\partial}{\partial \mathbf{x}} \int\_0^\infty \left[ A\_{11}(\lambda) e^{-\lambda z} + B\_{11}(\lambda) e^{+\lambda z} \right] J\_0(\lambda \rho) \frac{d\lambda}{\lambda} \tag{33}$$

$$A\_{\rm x12}^{z} = G\_{A\_{\rm x12}}^{z}(\rho, z) = \frac{\mu}{4\pi} \frac{\partial}{\partial \mathbf{x}} \int\_{0}^{\infty} \left[ A\_{12}(\lambda) e^{-\lambda \cdot z} + B\_{12}(\lambda) e^{+\lambda \cdot z} \right] J\_{0}(\lambda \rho) \frac{d\lambda}{\lambda} \tag{34}$$

**5. Time-domain solutions**

function *<sup>i</sup>*(*t*) = *Im*(*e*−*α<sup>t</sup>* <sup>−</sup> *<sup>e</sup>*−*β<sup>t</sup>*

the excitation function is defined by integral [40]:

Integral (37) can be evaluated analytically

the frequency response of the grounding system:

is defined by the integral [40]:

responses.

Once the transfer functions *W*(*jω*) have been determined for each calculated quantity, for example the electric field or current at specified points, the time-domain solutions can be obtained by direct application of (1). The calculation of the inverse Fourier transform is carried out by an FFT algorithm which is well-suited for the evaluation of time-domain

Numerical Calculation for Lightning Response to Grounding Systems Buried in Horizontal Multilayered Earth Model

The transient impedance, an essential parameter in grounding system design, is defined as a

*<sup>Z</sup>*(*t*) = *<sup>v</sup>*(*t*)

where *i*(*t*) represents the injected current at a grounding grid. This injected current represents the lightning channel current usually expressed by the double exponential

*α* and *β*, while *Im* denotes the amplitude of the current waveform. The Fourier transform of

*<sup>I</sup>*(*f*) = <sup>∞</sup>

*<sup>I</sup>*(*f*) = *Im*( <sup>1</sup>

−∞

*<sup>α</sup>* <sup>+</sup> *<sup>j</sup>*2*<sup>π</sup> <sup>f</sup>* <sup>−</sup> <sup>1</sup>

The frequency components up to few MHz are meaningfully present in the lightning current Fourier spectrum with strong decreasing importance from very low to highest frequencies. Multiplying the excitation function *I*(*f*) with the input impedance spectrum *Zin*(*f*) provides

Applying the IFFT, a time domain voltage counterpart is obtained. IFFT of the function *U*(*f*)

−∞

*<sup>U</sup>*(*t*) = <sup>∞</sup>

*β* + *j*2*π f*

*U*(*f*) = *I*(*f*)*Zin*(*f*) (39)

*U*(*f*)*ej*2*<sup>π</sup> f tdω* (40)

*<sup>i</sup>*(*t*) (36)

Based on Quasi-Static Complex Image Method

http://dx.doi.org/10.5772/57049

405

*i*(*t*)*ej*2*<sup>π</sup> f tdt* (37)

) (38)

), *t* ≥ 0, where pulse rise time is determined by constants

ratio of time varying voltage and current at injection point [19]:

Also like the deduced procedure for Green's function of a scalar monopole, the six constants *<sup>A</sup>*<sup>10</sup> ∼ *<sup>B</sup>*<sup>12</sup> can be decided by employing the interface and infinite conditions of the horizontal dipole's VMP. The *f*(*λ*) in the Green's function can also be developed into an exponential series with finite terms by the MP method. Last by applying the Lipschitz integration's varied form ( <sup>∞</sup> <sup>0</sup> *<sup>e</sup>*−*c<sup>λ</sup> <sup>J</sup>*0(*λρ*) *<sup>d</sup><sup>λ</sup> <sup>λ</sup>* <sup>=</sup> ln(*<sup>c</sup>* <sup>+</sup> *c*<sup>2</sup> <sup>+</sup> *<sup>ρ</sup>*2)), the final expression of *<sup>G</sup><sup>z</sup> Ax*<sup>11</sup> can be given as follows.

$$G\_{A\_{11}}^{z}(\rho, z) = \frac{\mu}{4\pi} \frac{\partial}{\partial x} [k\_{01} \ln(z\_0^{'} + R\_0^{'}) - \sum\_{n=1}^{N} a\_n (k\_{01}^2 k\_{12} \ln(z\_{n1} + R\_{n1}))]$$

$$\left[ +k\_{01}k\_{12}\ln(z\_{n2} + R\_{n2}) - k\_{12}\ln(z\_{n3} + R\_{n3}) + k\_{01}k\_{12}\ln(z\_{n4} + R\_{n4})) \right] \tag{35}$$

where the origin of the coordinate system shown in Fig. 1 has been moved to the surface of the earth. The *<sup>R</sup>*′ <sup>0</sup> and *Rni*,*i*=1∽<sup>4</sup> are same as Eq. (23). Other parameters are *z* ′ <sup>0</sup> <sup>=</sup> *<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′ , *zn*(1−4) = (*signa* · *<sup>z</sup>* <sup>+</sup> *signb* · *<sup>z</sup>*′ ) + *zn*, in which *signa* = 1 for *zn*<sup>1</sup> and *zn*3, *signb* = 1 for *zn*<sup>3</sup> and *zn*4, *signa* = *signb* = −1 for others.

It can also be seen that each term except <sup>1</sup> *<sup>R</sup>*<sup>0</sup> of Eqs. (30) and (35) can be regarded as a image vector dipole source, whose location is indicated by *Rni*,*i*=1∽<sup>4</sup> and amplitude is *αn*, which can be seen in Fig. 1. However, in Eqs. (30) and (35), *α<sup>n</sup>* and *Rni* are usually complex numbers, so that this approach is also named as the QSCIM.

The other Green's function for the dipole *G<sup>z</sup> Az*<sup>10</sup> , *<sup>G</sup><sup>z</sup> Az*<sup>12</sup> , *<sup>G</sup><sup>x</sup> Ax*<sup>10</sup> , *<sup>G</sup><sup>x</sup> Ax*<sup>12</sup> , *<sup>G</sup><sup>z</sup> Ax*<sup>10</sup> and *<sup>G</sup><sup>z</sup> Ax*<sup>12</sup> can be given in the similarly way.

Like closed form of Green's function of a scalar monopole, the procedure for the closed form of the Green's function of a vector dipole in arbitrary horizontal multilayered earth model is also first to derive the specific solution of the Green's function in each layer based on the general expressions (for example, the general expressions of Green's function for a vertical dipole or x component of horizontal dipole in three horizontal layers medium are given in Eq. (27), Eq. (28) and Eq. (29); for z component of horizontal dipole, are given in Eq. (32), Eq. (33) and Eq. (34)), and the interface and infinite conditions of VMP of the dipole are also used to decide unknown constants. Then, for z component of horizontal dipole and vertical dipole, by employing the MP method exponential series can be developed. Last by utilizing the Lipschitz integration and its varied form, the final expression of the Green's function of the dipole can be achieved.

Once the closed form of Green's function of a vector dipole in horizontal multilayered earth model has been given, the mutual impedance coefficient Eq. (9) can also be fast analytical calculated just like Eq. (10).
