**3.2. Computation of** [**Zb**] **and** [**Zs**] **matrices**

4 Computational and Numerical Simulationse

whose elements we refer to [28]–[31].

which will be introduced later.

*<sup>E</sup>*¯(*r*¯*j*) = −

*Nl* ∑ *i*=1 *jω* 

*li*

calculated by

**3.1. Mathematical model of the grounding system in the frequency domain**

[**Y**]=[**K**]

*ϕ*(*r*¯*j*) =

*A*¯(*r*¯*j*) =

The electrical field intensity (EFI) at any point can be calculated by

*<sup>G</sup>*¯¯ *<sup>A</sup>*(*r*¯*j*,*r*¯*i*) • ¯*Ibi*

*Nl* ∑ *i*=1

The vector magnetic potential (VMP) *A* at any point can be calculated by

*Nl* ∑ *i*=1 *li*

*li*

*<sup>t</sup>* · [**Zs**]

With the above considerations, according to [28]–[31], the electric circuit may be studied using the conventional nodal analysis method [35], resulting in the following equations:

<sup>−</sup><sup>1</sup> · [**K**]+[**A**] · [**Zb**]

where [**F**] is an *Np* × 1 vector of external current sources; [**Z***b*] is the *Nl* × *Nl* branch mutual induction matrix of the circuit including resistive and inductive effects, which gives a matrix relationship between branch currents [**Ib**]; [**Zs**] is an *Nl* × *Nl* mutual impedance matrix, which gives a matrix relationship between the average SEP [**U**] and leakage currents [**Is**] through the rapid Galerkin moment method [38]. Both [**A**] and [**K**] are incidence matrices, which are used to relate branches and nodes. There are rectangular matrices of order *Nl* × *Np*, for

The vector of nodal SEP [**Vn**] may be obtained by solving 2. The average SEP [**U**], leakage current [**Is**], branch voltage [**Ul**], and branch current [**Ib**] can also be calculated [28]–[31]. Once the branch currents and leakage currents are known, the SEP at any point can be

where *Gϕ*(*r*¯*j*,*r*¯*i*) is the scalar Green's function of a monopole in the multilayered earth model.

Here, *<sup>G</sup>*¯¯ *<sup>A</sup>*(*r*¯*j*,*r*¯*i*) is the dyadic Green's function of a dipole in the multilayered earth model,

(*r*¯*i*)*dti* −

*Nl* ∑ *i*=1 ∇ 

*li*

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

*li*

*dli*. (6)

*<sup>G</sup>*¯¯ *<sup>A</sup>*(*r*¯*j*,*r*¯*i*) • ¯*Ibi*

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

*li*

[**F**]=[**Y**] · [**Vn**] (2)

*<sup>t</sup>* (3)

*dli* (4)

(*r*¯*i*)*dti*. (5)

<sup>−</sup><sup>1</sup> · [**A**]

From [28]–[31], we know that each segment is modeled as a lumped resistance and self-inductance. Mutual inductances or impedances between branch segments' branch currents or leakage currents are also included in the model: 

$$[\overline{\mathbf{X\_{q}}}]\_{N\_{l}\times N\_{l}} = \begin{bmatrix} \mathbf{X\_{1,1}} & \dots & \mathbf{X\_{1,i}} & \dots & \mathbf{X\_{1,j}} & \dots & \mathbf{X\_{1,N\_{l}}} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{X\_{i,1}} & \dots & \mathbf{X\_{i,i}} & \dots & \mathbf{X\_{i,j}} & \dots & \mathbf{X\_{i,N\_{l}}} \\ \vdots & \dots & \vdots & \dots & \vdots & \vdots & \vdots \\ \mathbf{X\_{j,1}} & \dots & \mathbf{X\_{j,i}} & \dots & \mathbf{X\_{j,j}} & \dots & \mathbf{X\_{j,N\_{l}}} \\ \vdots & \dots & \vdots & \dots & \vdots & \vdots & \vdots \\ \mathbf{X\_{N\_{l}}} 1 & \dots & \mathbf{X\_{N\_{l}}} & \dots & \mathbf{X\_{N\_{l}}} & \dots & \mathbf{X\_{N\_{l}}} \end{bmatrix} \tag{8}$$

$$\text{1. The case of } [\overline{\overline{\mathbf{X\_{q}}}}] = [\overline{\overline{\mathbf{Z\_{b}}}}] \text{: } X\_{\mathbf{i},\mathbf{i}} = \mathcal{R}\_{\mathbf{i}} + \text{j}\omega L\_{\mathbf{i}\prime} \text{ } X\_{\mathbf{i},\mathbf{j}} = \mathcal{R}\_{\mathbf{i},\mathbf{j}} + \text{j}\omega M\_{\mathbf{i},\mathbf{j}}.$$

The diagonal elements consists of self impedance and self induction, the other elements belong to mutual induction between a pairs of conductor segments. The formula for self impedance and self induction can be found in [28]–[31], here, we give the formula for mutual induction:

$$M\_{\mathbf{i},\mathbf{j}} = \int\_{\mathbf{l}\_{\parallel}} \int\_{\mathbf{l}\_{\parallel}} \bar{\mathbf{G}}\_{A}(\bar{r}\_{\mathbf{j}}, \bar{r}\_{\mathbf{i}}) \bullet \bar{\mathbf{I}}\_{\mathbf{b}\_{\parallel}} dt\_{\mathbf{i}} \bullet \bar{\mathbf{I}}\_{\mathbf{b}\_{\parallel}\mathbf{i}} dt\_{\mathbf{j}} \tag{9}$$

For an infinite homogeneous conductivity medium, one has *<sup>G</sup>*¯¯ *<sup>A</sup>*(*r*¯*j*,*r*¯*i*) = *<sup>µ</sup>* <sup>4</sup>*<sup>π</sup>* <sup>1</sup> *Rij* ¯¯*I*0, where ¯¯*I*<sup>0</sup> is the diagonal unit matrix.

$$M\_{\mathbf{i},\mathbf{j}} = \frac{\mu}{4\pi\epsilon} \int\_{\mathbf{l}\_i} \int\_{\mathbf{l}\_j} \frac{1}{\mathcal{R}\_{i\mathbf{j}}} \mathbf{\bar{I}}\_{\mathbf{b}\mathbf{l}} d\mathbf{t}\_i \bullet \mathbf{\bar{I}}\_{\mathbf{b}\mathbf{j}} d\mathbf{t}\_j \tag{10}$$

The above integral can be analytically calculated, this formula can be found in [27].

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 earth, and its elements are the mutual impedance coefficients *Zi*,*j*.

$$Z\_{i,j} = \int\_{I\_l} \int\_{I\_j} \mathbf{G}\_{\Phi}(\vec{r}\_j, \vec{r}\_i) \, \frac{dt\_i}{l\_i} \frac{dt\_j}{l\_j} \tag{11}$$

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

$$Z\_{i,j} = \frac{1}{4\pi\sigma\_1 l\_j l\_j} \int\_{l\_l} \int\_{l\_l} \frac{1}{R\_{ij}} dt\_i dt\_j \tag{12}$$

**Figure 1.** The earth model

Poisson equation as:

where *R* = |*r*¯ − ¯

where *<sup>i</sup>*, *<sup>j</sup>* <sup>=</sup> 1, . . . , *Ns*, and *<sup>δ</sup>*(*r*,*r*′

the spherical coordinate system is.

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

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

*σ<sup>i</sup>* = *σ<sup>i</sup>* + *jωε<sup>i</sup>* is the complex conductivity of the *i*th layer medium. Supposing the monopole

If the monopole lies in infinite homogenous medium, its expression of Green's function in

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

<sup>∞</sup>

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

*r*′| is the distance between the source point and field point. While its

*ϕ* = *G*(*r*¯, ¯

4*π* · *σ*

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

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

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

)*δ*(*ij*) *σi*

) is the Dirac delta function. *δ*(*ij*) is Kronecker's symbol,

<sup>4</sup>*<sup>π</sup>* · *<sup>σ</sup>* · *<sup>R</sup>* (14)

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

Based on Quasi-Static Complex Image Method

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

399

(13)

<sup>∇</sup>2*ϕ<sup>j</sup>* <sup>=</sup> <sup>−</sup>*δ*(*r*,*r*′

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

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

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 discussed next.
