**5. Time-domain solutions**

12 Computational and Numerical Simulationse

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

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

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

4*π ∂ ∂x* [*k*<sup>01</sup> ln(*z* ′ <sup>0</sup> + *R* ′ <sup>0</sup>) <sup>−</sup>

4*π ∂ ∂x* <sup>∞</sup>

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

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

*<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>*

*N* ∑ *n*=1

+ *<sup>k</sup>*01*k*<sup>12</sup> ln(*zn*<sup>2</sup> + *Rn*2) − *<sup>k</sup>*<sup>12</sup> ln(*zn*<sup>3</sup> + *Rn*3) + *<sup>k</sup>*01*k*<sup>12</sup> ln(*zn*<sup>4</sup> + *Rn*4))] (35)

<sup>0</sup> and *Rni*,*i*=1∽<sup>4</sup> are same as Eq. (23). Other parameters are *z*

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

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,

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

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

*Az*<sup>10</sup> , *<sup>G</sup><sup>z</sup>*

*Az*<sup>12</sup> , *<sup>G</sup><sup>x</sup>*

*αn*(*k*<sup>2</sup>

) + *zn*, in which *signa* = 1 for *zn*<sup>1</sup> and *zn*3, *signb* = 1 for *zn*<sup>3</sup>

*Ax*<sup>10</sup> , *<sup>G</sup><sup>x</sup>*

*<sup>R</sup>*<sup>0</sup> of Eqs. (30) and (35) can be regarded as a image

*Ax*<sup>12</sup> , *<sup>G</sup><sup>z</sup>*

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

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

+*λ*·*z* 

<sup>01</sup>*k*<sup>12</sup> ln(*zn*<sup>1</sup> + *Rn*1)

*J*0(*λρ*)

*dλ*

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

*Ax*<sup>11</sup> can be given

′

<sup>0</sup> <sup>=</sup> *<sup>z</sup>* <sup>+</sup> *<sup>z</sup>*′

*Ax*<sup>12</sup> can be

,

0

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

varied form ( <sup>∞</sup>

*Gz*

of the earth. The *<sup>R</sup>*′

*zn*(1−4) = (*signa* · *<sup>z</sup>* <sup>+</sup> *signb* · *<sup>z</sup>*′

given in the similarly way.

the dipole can be achieved.

calculated just like Eq. (10).

and *zn*4, *signa* = *signb* = −1 for others. It can also be seen that each term except <sup>1</sup>

so that this approach is also named as the QSCIM.

The other Green's function for the dipole *G<sup>z</sup>*

as follows.

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

The transient impedance, an essential parameter in grounding system design, is defined as a ratio of time varying voltage and current at injection point [19]:

$$Z(t) = \frac{v(t)}{i(t)}\tag{36}$$

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 function *<sup>i</sup>*(*t*) = *Im*(*e*−*α<sup>t</sup>* <sup>−</sup> *<sup>e</sup>*−*β<sup>t</sup>* ), *t* ≥ 0, where pulse rise time is determined by constants *α* and *β*, while *Im* denotes the amplitude of the current waveform. The Fourier transform of the excitation function is defined by integral [40]:

$$I(f) = \int\_{-\infty}^{\infty} i(t)e^{j2\pi ft}dt\tag{37}$$

Integral (37) can be evaluated analytically

$$I(f) = I\_m(\frac{1}{\mathfrak{a} + j2\pi f} - \frac{1}{\mathfrak{F} + j2\pi f})\tag{38}$$

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 the frequency response of the grounding system:

$$\mathcal{U}(f) = I(f)Z\_{\text{in}}(f) \tag{39}$$

Applying the IFFT, a time domain voltage counterpart is obtained. IFFT of the function *U*(*f*) is defined by the integral [40]:

$$
\mathcal{U}I(t) = \int\_{-\infty}^{\infty} \mathcal{U}(f) e^{j2\pi ft} d\omega \tag{40}
$$

As the frequency response *U*(*f*) is represented by a discrete set of values the integral (41) cannot be evaluated analytically and the discrete Fourier transform, in this case the IFFT algorithm, is used, i.e.,

$$\mathcal{U}I(t) = IFFT(\mathcal{U}(f))\tag{41}$$

Earth model Mea.[42] Our model

ratio for the upper and the lower soil layers is *ρ*1/*ρ*<sup>2</sup> = 50/20, the upper layer thickness being *H* = 0.6 m. The inject lightning current parameter was set to *T*<sup>1</sup> = 3.5*µ* s, *T*<sup>2</sup> = 73*µ* s and *Im* = 12.1 A, the feed point is at the corner of the grid. For our model, the permittivities of the two layer earth model were set to *ε*<sup>1</sup> = 30*ε*<sup>0</sup> and *ε*<sup>2</sup> = 20*ε*0. The transient SEP can be seen in Fig. 3, which ultimately agreed with the measured curve in Fig. 4 (a) in [41]; meanwhile, the impulse grounding impedance was 2.12 Ω as given by [41], and it is 2.08Ω

The maximum frequency of applicability of the method is limited by the quasi-stationary approximation of the electromagnetic fields, which means the propagation effect of the

<sup>−</sup>*γeR* <sup>≈</sup> <sup>0</sup> (44)

electromagnetic field around the grounding system can be neglected, so

*e*

For most of the usual electrodes, this may be applied up to some hundreds of kHz.

**Table 1.** Comparison with published measurement result: Frequency is 80Hz

for our model.

**Figure 2.** Typical grounding grid configuration

**6.2. Validation of our method**

*<sup>e</sup>* = *jωµ* (*jωε<sup>e</sup>* + *σe*),*e* = 1, . . . , *Ne*.

where *γ*<sup>2</sup>

a 0.390 (0.387, 8.6 × <sup>10</sup>−3) b 0.658 (0.652, 1.6 × <sup>10</sup>−2)

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

407

Implementation of this algorithm inevitably causes an error due to discretization and truncation of essentially unlimited frequency spectrum. The discrete set of the time domain voltage values is defined as [40]:

$$\mathcal{U}I(t) = \mathcal{F} \sum\_{k=0}^{N} \mathcal{U}(k\Delta f) e^{jk\Delta f n \Delta t} \tag{42}$$

where *n* = 0, . . . , *N*, *F* is the highest frequency taken into account, *N* is the total number of frequency samples, ∆*f* is sampling interval and ∆*t* is the time step.

Finally, the impulse impedance, an also essential parameter in grounding system design, is defined by the following expression [41]:

$$Z\_{\mathbb{C}} = \frac{\mathcal{U}}{I} \tag{43}$$

where *U* is the voltage maximum at the discharge point and *I* is the injected current magnitude at the time instant when *U* has been reached.
