Progress In Electromagnetics Research, Vol.

117, 209–236, 2011

EVALUATION OF ELECTROMAGNETIC FIELDS ASSO-

CIATED WITH INCLINED LIGHTNING CHANNEL US-
ING SECOND ORDER FDTD-HYBRID METHODS

M. Izadi* , M. Z. A. A. Kadir, and C. Gomes

Centre of Excellence on Lightning Protection (CELP), Faculty of
Engineering, University Putra Malaysia, UPM, Serdang, Selangor
43400, Malaysia

Abstract—Evaluation of electromagnetic fields caused by the

lightning channel is an appealing topic in order to consider the indirect
effects of lightning on the power lines. A common assumption for
the calculation of electromagnetic fields at the observation point is a
vertical lightning channel, but the fact is that in reality the lightning
channel is seldom vertical on the ground surface. In this study, the
electromagnetic fields due to inclined lightning channel at various
observation points with different angles and with respect to the image
of lightning channel on the ground surface were explored. This study
also proposes general equations that can estimate the electric fields
due to inclined lightning channel through the 2nd FDTD method.
The proposed method supports the notion of vertical lightning channel
while the channel angle with respect to z-axis is assumed to be zero.
This method was validated through the data gathered from five fields:
three at a close distance from inclined lightning channel and two at
intermediate distances from vertical lightning channel. Similarly, due
to inclined lightning channel, the effects of geometrical and current
parameters on the electromagnetic fields are considered. This study
substantiates different coupling models with FDTD structure directly
at the time domain without a need for extra converters.

1. INTRODUCTION

One of the appealing areas under discussion is the electromagnetic

fields due to lightning channel in order to appoint proper protection
level on the power lines and equipments while the electromagnetic
Received 21 April 2011, Accepted 17 May 2011, Scheduled 6 June 2011
* Corresponding author: Mahdi Izadi ([email protected]).
210 Izadi, Kadir, and Gomes

fields is one of the most influential factors in creating induction

voltage on the power lines [1]. Several studies have been done on
the evaluation of electromagnetic fields due to lightning channel with
the assumption that the lightning channel strikes the ground at a
vertical angle [2]. However, it is known to be tortuous on scales
ranging from less than 1 m to over 1 km [3–5]. Moreover, the lightning
channel typically strikes the ground at an angle less than ninety
degrees (with respect to the ground) [6]. Several studies on the
radiated electromagnetic fields associated with inclined and tortuosity
lightning channels have been done while generally the tortuosity
lightning channel is modeled with some individual linear inclined
channel segments [3, 5, 7, 8]. Accordingly, this paper focuses on the
computation of the electromagnetic fields due to the inclined lightning
channel at the observation point in the time domain using second
order FDTD (later referred as 2nd FDTD) Hybrid method. The
method is then validated by the data collected from three fields stricken
by lightning at a close distance from the inclined lightning channel
whereas the channel angle (between lightning channel and z-axis) can
be determined by a high speed camera [6]. In addition, the effects of
the channel angel and observation point angle (between the vector from
channel base to image of observation point on the ground and the image
of lightning channel on the ground surface) on the electromagnetic field
values are evaluated. The results due to inclined channel are compared
with the field values associated with the vertical channel. Also, the
proposed method is validated for vertical lightning channel where the
channel angle is equal to zero. The three fundamental assumptions in
this study are:
(i) The ground conductivity is perfect
(ii) The lightning channel is without any branches
(iii) The current velocity along lightning channel is constant
For the simulation of electromagnetic fields due to lightning channel, it
is necessary to consider both the return stroke current at the channel
base and the lightning channel in a special function and a current
model, respectively. In this study, the channel base current is simulated
by the DU function [9] which is an improved function of Heidler’s
function [10] and is in agreement with the measured current at the
channel base. The engineering models are used in this simulator, since
the return stroke currents at different heights along with lightning
channel are generally expressed by a special function based on the
channel base current and a height dependent attenuation factor in
these models.
Progress In Electromagnetics Research, Vol. 117, 2011 211

2. RETURN STROKE CURRENT

As mentioned earlier, the channel base current in this study

is simulated by the DU function which can be expressed by
Equation (1) [9]
 ³ ´n1 ³ ´n2 
t µ ¶ t µ ¶
i 01 Γ −t i02 Γ −t
i (0, t)=
11 21
³ ´ exp + ³ ´ exp 
η1 1+ t n1 Γ12 η2 1+ t n2 Γ22
Γ11 Γ21
(1)
where
i01 , i02 are the amplitudes of the channel base current,
Γ11 , Γ12 are the front time constants,
Γ21 , Γ22 are the decay-time constants,
n1 , n2 are· the exponents (2 ∼ 10), ¸
³ ´1
n1
η1 = exp − (Γ11 /Γ12 ) n1 × ΓΓ1211
,
· ³ ´1 ¸
n2
η2 = exp − (Γ21 Γ22 ) n2 × ΓΓ2221

For instance, the measured channel bases current from a triggered

lightning experiment in Florida is illustrated in Figure 1. The channel
base current parameters using DU function with respect to measured
current are estimated while the numerical fitting method is applied;
the results are tabulated in Table 1.
A simulation of the channel base current using initial current
parameters from Table 1 is shown in Figure 2.
In this study, the current wave shape along with the lightning
channel is predicted by the engineering return stroke current models,

Figure 1. The measured return stroke current at channel bases from

triggered lightning experiment [6].
212 Izadi, Kadir, and Gomes

Table 1. Typical values for Diendorfer and Uman channel base current
function.
i01 (A) i02 (A) Γ11 (s) Γ12 (s) Γ21 (s) Γ22 (s) n1 n2
8 × 103 3.6671 × 103 9.5 × 10−8 8.5 × 10−7 1.4 × 10−6 1.8 × 10−5 2 2

7
Channel Base Current(kA)

0
0 5 10 15
Time (µs)

Figure 2. The simulated channel base current using DU function

based on current parameters from Table 1.

Table 2. P (z 0 ) and v for five engineering models based on

Equation (2) [3].
Model P (z 0 ) v
BG 1 ∞
TCS 1 −c
TL 1 vf
z0
MTLL (1− H ) vf
z0
MTLE exp(− λ ) vf

while the focus is on the current by a special function that depends on

channel base current and a height dependent attenuation factor. The
engineering models can usually be presented by a general equation
such as Equation (2). The most common applied engineering models
are namely the Bruce-Golde model (BG) [11], Traveling Current
Source model (TCS) [12], Transmission Line model (TL) [2], Modified
Transmission Line Linear model with Linear decay with height
(MTLL) [2], Modified Transmission Line model with Exponential
decay with height (MTL) or (MTLE) [2]. Theses models are tabulated
Progress In Electromagnetics Research, Vol. 117, 2011 213

in Table 2 based on the current parameters in Equation (2) [3, 13]

where λ is the decay constant (∼ 2 km) and H is the cloud’s height.
µ ¶ µ ¶
¡ ¢ z0 ¡ ¢ z0
I z 0 , t = I 0, t− ×P z 0 × u t − (2)
v vf
where:
z 0 is the temporary charge height along the lightning channel,
I(z 0 , t) is the current distribution along the lightning channel at any
height z 0 and any time t,
I(0, t) is the channel base current;
P (z 0 ) is the height dependent attenuation factor,
v is the velocity of charge(return stroke current velocity),
vf is the velocity of return stroke front,
u is the Heaviside function defended as
µ ¶ ( 0
z0 1 for t ≥ vzf
u t− = 0
vf 0 for t < vzf
There are some drawbacks regarding these models. For instance, the
BG and TCS models are not convenient, due to the discontinuity
at current [2]. Whilst the TL model is not a complete model
because it is not considered on the removed charge along the lightning
channel [2, 14]. This model was later modified by the MTLL and
MTLE models. In the MTLL and MTLE models, some charges are
removed from channel with linear cloud height dependent factor and
exponentially rate, respectively [2, 14] This study opted for MTLE
model because of the channel height value is considered as unmeasured
parameter while the value of λ is recommended around 2 km. In
addition, It should be highly remarked that the return stroke velocity
is less than light speed in the free space (c) and is typically assumed to
be between (c/3) to (2c/3) [15] the return stroke current wave shape
at two different heights with respect to the ground surface along the
channel using TL, MTLE and MTLL models are shown in Figure 3
where the channel base current parameters are applied from Table 2.
Note that, the λ and H factors are set on 2000 m and 7000 m for MTLE
and MTLL models respectively while the return stroke velocity is fixed
at 1.3×108 m/s.

3. ELECTROMAGNETIC FIELDS DUE TO LIGHTNING

CHANNEL

Dipole method which is given in Equations (3) to (5) can be used

to show the electromagnetic fields due to vertical lightning channel;
214 Izadi, Kadir, and Gomes

9000
TL model,z'=1000 m
MTLE model,z'=1000 m
8000
MTLL model,z'=1000m
TL model,z'=3000 m
7000 MTLE model,z'=3000 m
Return stroke current(A)

MTLL model,z'=3000m
6000

5000

4000

3000

2000

1000

0
0 5 10 15 20 25 30
Time(µs)

Figure 3. The return stroke current wave shape at different heights

along lightning channel using TL, MTLE and MTLL models.

Figure 4. The geometry of the vertical lightning channel.

Progress In Electromagnetics Research, Vol. 117, 2011 215

while the geometry of problem is illustrated in Figure 4 (the ground is

assumed to be perfect) [2, 16, 17]
µ ¶Z H1 µ Z µ ¶
~ 1 3r(z−z 0) t 0 R
Er (r, z, t) = i z , τ− dτ
4πε0 H2 R5 0 c
µ ¶ ¡ ¢!
3r(z−z 0) 0 R r(z−z 0) ∂i z 0 , t− Rc
+ i z , t− + 2 3 dz 0 (3)
cR4 c c R ∂t
µ ¶Z H1 µ Z µ ¶
~ 1 2(z−z 0)2 −r2 t 0 R
Ez (r, z, t) = i z , τ− dτ
4πε0 H2 R5 0 c
µ ¶ ¡ ¢!
2(z−z 0)2 −r2 0 R r2 ∂i z 0 , t− Rc
+ 4 i z , t− − 2 3 dz 0 (4)
cR c c R ∂t
Z Ã µ ¶ ¡0 ¢!
³µ ´ H1 r R r ∂i z , t− R
0
B~ϕ (r, z, t) = i z 0 , t− + 2 c
dz 0 (5)
4π H2 R3 c cR ∂t
where q
R = r2 +(z−z 0 )2 ,
q
R = r2 +(z + z 0 )2 ,
0

E~r (r, z, t) is the horizontal electric field due to vertical lightning

channel,
E~z (r, z, t) is the vertical electric field due to vertical lightning channel,
B~ϕ (r, z, t) is the magnetic flux density due to vertical lightning channel,
z is height of observation point,
z 0 is the vertical space variable,
ε0 is the permittivity of free space,
µ0 is the permeability of free space,
β = vc , ½ ¾
q p
0 β 2 2 2
H1 =z = 1−β 2 − (βz − ct) − (βct − z) +(r 1−β ) ,
½ q ¾
β 2
p
00
H2 = z = 1−β 2 − (βz + ct) + (βct + z) +(r 1−β ) . 2 2

The different components of electric field associated with the lightning

channel at a remote observation point can be estimated by the 2nd
FDTD-Hybrid method as presented by Equations (6) to (8) and
the magnetic flux density can be calculated by Equation (5) [16].
It should be noted that, the observation point is held constant
at (m∆y, n∆x, p∆z), the lightning channel is along z-axis and the
216 Izadi, Kadir, and Gomes

processing time is equal to k∆t

(
−ε E~ k−2 (m, n, p) −4E ~ k−1 (m, n, p)
~ y (m, n, p) =
k 0 y y
E 3ε0
σ+ 2∆t 2∆t
à !)
cosγ1 B ~ k (m, n, p + 1) −B ~ k (m, n, p − 1)
ϕ ϕ
+ (6)
ε0 µ0 2∆z
(
−ε 0 E~ xk−2 (m, n, p) −4E~ xk−1 (m, n, p)
~ xk (m, n, p) =
E 3ε0 2∆t
σ+ 2∆t
à !)
sinγ1 B ~ ϕk (m, n, p + 1) −B ~ ϕk (m, n, p − 1)
− (7)
ε0 µ0 2∆z
(
−ε 0 E~ k−2 (m, n, p) −4E ~ k−1 (m, n, p)
~ zk (m, n, p) =
E z z
3ε0 2∆t
σ+ 2∆t

1  cosγ3 B ~ kϕ (m + 1, n, p) −cosγ4 B ~ kϕ (m − 1, n, p)
−
ε0 µ0 2∆y

k k
− sinγ5 B ϕ (m, n + 1, p) +sinγ6 B ϕ (m, n − 1, p) 
~ ~
−  (8)
2∆x 
where σ is the medium conductivity,
m∆y
cosγ1 = q
(m∆y)2 +(n∆x)2
n∆x
sinγ2 = q
(m∆y)2 +(n∆x)2
(m + 1)∆y
cosγ3 = q
((m + 1)∆y)2 +(n∆x)2
(m − 1)∆y
cosγ4 = q
((m − 1)∆y)2 +(n∆x)2
(n + 1)∆x
sinγ5 = q
(m∆y)2 +((n + 1)∆x)2
(n − 1)∆x
sinγ6 = q
(m∆y)2 +((n − 1)∆x)2
Progress In Electromagnetics Research, Vol. 117, 2011 217

The geometry of problem for estimation of electromagnetic fields due

to inclined lightning channel is shown in Figure 6, with the assumption
that the lightning channel usually strikes the ground surface with an
inclined shape as depicted in Figure 5. The angle between lightning
channel and z-axis is illustrated by θ and is called channel angle (see
Figure 6) Besides the observation point is located above the ground
surface with special angle between the vector from channel base to
image of observation point on the ground (~r) and the image of lightning
channel on the ground surface (y-axis).
Having applied Maxwell’s equations [19] and Lorentz’s gauge [20],
different components of magnetic flux density in the Cartesian
coordinate can be estimated by Equations (9) to (11) which is shown
in Figure 6. Note that, the equations of (9) to (11) could easily be
solved using Gauss-Lobatto quadrature method [21].
Z H1
cos θ
B~ x (x, y, z, θ, t) = 10 ×
−7
H2

  ³ cos θ ´ 
 0 ∂i r 0 , t− R(r )
0
0
µ 0
¶
y−r sin θ c y−r sin θ R (r ) 
− cos θ  2 + 0
3 i r , t−
 cR (r ) 0 ∂t 0
R (r ) c
 ³ ´ 
0 ∂i r 0 , t− R(r )
0
0
µ 0
¶
z−r cos θ c z−r cos θ R (r ) 
+ sin θ 2 + 3 i r0 , t− dr0(9)
0
cR (r ) ∂t 0
R (r ) c 

Z H1
cos θ
~ y (x, y, z, θ, t) = 10
B −7
× cos θ
H2
cos θ
 ³ ´ 
 x ∂i r0 , t− R(r
0)
µ 0
¶
c x 0 R (r )
+ 3 i r , t− dr0 (10)
cR (r0)2 ∂t 0
R (r ) c 

Figure 5. The video frame for the return stroke lightning channel
during the triggered lightning experiment [18].
218 Izadi, Kadir, and Gomes

Figure 6. The geometry of the inclined lightning channel with an

observation point above the ground surface.

Z H1
cos θ
~ z (x, y, z, θ, t) = −10−7 × sin θ
B
H2
cos θ
 ³ ´ 

0
∂i r0 , t− R(rc ) µ 0
¶
x x 0 R (r ) 0
+ 3 i r , t− c dr (11)
cR (r0 )2 ∂t 0
R (r )

where:
r is radial distance frompchannel base to image of observation point on
the ground surface (r = x2 +y 2 ),
r0 is the temporary channel length along lightning channel,
θ is the angle between lightning channel and z axis (the channel angle),
φ is the angle between y-axis and ~r while y-axis is fixed on the image of
lightning channel on the ground surface (the observation point angle),
x is the position of observation point at x axis (x = r sin φ),
Progress In Electromagnetics Research, Vol. 117, 2011 219

y is the position of observation point at y axis (y = r cos φ),

z is observation point height from ground surface,
B~ x (x, y, z, θ, t) is the magnetic flux density at x-direction due to
inclined lightning channel,
B~ y (x, y, z, θ, t) is the magnetic flux density at y-direction due to
inclined lightning channel,
B~ z (x, y, z, θ, t) is the magnetic flux density at z-direction due to
inclined lightning channel.
q
0 ~00 ~
R(r ) = |r −r | = x2 + (y−r0 sin θ)2 + (z−r0 cos θ)2 .
0

In addition, the magnetic flux density components in the cylindrical

domain can be estimated by Equations (12) to (14)
³ ´
~ r (r, z, φ, θ, t) = cos π − φ B
B ~ x (x, y, z, θ, t)
2³ ´
π ~ y (x, y, z, θ, t)
+ sin −φ B (12)
³π2 ´
B ~ ϕ (r, z, φ, θ, t) = − sin −φ B ~ x (x, y, z, θ, t)
³π2 ´
+ cos −φ B ~ y (x, y, z, θ, t) (13)
2
B~ z (r, z, φ, θ, t) = B
~ z (x, y, z, θ, t) (14)
where:
B~ r (r, z, φ, θ, t) is the magnetic flux density at r-direction due to
inclined lightning channel,
B~ ϕ (r, z, φ, θ, t) is the magnetic flux density at ϕ-direction due to
inclined lightning channel,
B~ z (r, z, φ, θ, t) is the magnetic flux density at z-direction due to
inclined lightning channel,
φ = arccos ( √ 2y 2 )
x +y
The magnetic flux density could also be expressed by Equation (15)
regarding Maxwell’s equations [16].
à ! à !
∂ ~
D ∂ ~
E
∇×B=µ~ 0 j+
~ 0
=µ0 σ E+ε (15)
∂t ∂t
where:
J~ is the current density,
D~ is the electric flux density.
The solutions to curl function in Equation (15) in the Cartesian
220 Izadi, Kadir, and Gomes

coordinates can be presented by Equations (16) to (18) [22]

à !
∂E~x 1 ∂B~ z ∂B~y
= − ~x
−µ0 σ E (16)
∂t ε0 µ0 ∂y ∂z
à !
∂E~y 1 ∂B~ x ∂B~z
= − ~y
−µ0 σ E (17)
∂t ε0 µ0 ∂z ∂x
à !
∂E~z 1 ∂B~ y ∂B~x
= − ~z
−µ0 σ E (18)
∂t ε0 µ0 ∂x ∂y
Therefore, by solving Equations (16) to (18) based on 2nd
FDTD method [22–24] and applying perfect ground conductivity
conditions, different components of electric field can be expressed by
Equations (19) to (21) At the same time, the components of magnetic
flux density are obtained from Equations (9) to (11) and the electric
fields components in the cylindrical coordinates from Equations (22)
to (24)
Ã
2∆t B~ z (x, y+∆y, z, θ,tk ) −B ~ z (x, y−∆y, z, θ,tk )
E~ x (x, y, z, θ,tk ) =
3ε0 µ0 2∆y
!
B~ y (x, y, z + ∆z, θ,tk ) −B ~ y (x, y, z − ∆z, θ,tk )
−
2∆z
1 ~h i
+ 4E x (x, y, z, θ,tk−1 ) −E~ x (x, y, z, θ,tk−2 ) (19)
3 Ã
~ ~
E~ y (x, y, z, θ,tk ) = 2∆t Bx (x, y, z+∆z, θ,tk ) −Bx (x, y, z−∆z, θ,tk )
3ε0 µ0 2∆z
!
B~ z (x + ∆x, y, z, θ,tk ) −B ~ z (x − ∆x, y, z, θ,tk )
−
2∆x
1 ~h i
+ 4E ~
y (x, y, z, θ,tk−1 ) −Ey (x, y, z, θ,tk−2 ) (20)
3 Ã
~ ~
E~ z (x, y, z, θ,tk ) = 2∆t By (x+∆x, y, z, θ,tk ) −By (x−∆x, y, z, θ,tk )
3ε0 µ0 2∆x
!
B~ x (x, y + ∆y, z, θ,tk ) −B ~ x (x, y − ∆y, z, θ,tk )
−
2∆y
1 ~h i
+ 4E ~
z (x, y, z, θ,tk−1 ) −Ez (x, y, z, θ,tk−2 ) (21)
3
Progress In Electromagnetics Research, Vol. 117, 2011 221

³ ´
~ r (r, z, φ, θ,tk ) = cos π − φ E
E ~ x (x, y, z, θ,tk )
2³ ´
π ~ y (x, y, z, θ,tk )
+ sin −φ E (22)
³ 2 ´
~ ϕ (r, z, φ, θ,tk ) = − sin π − φ E
E ~ x (x, y, z, θ,tk )
2
³π ´
+ cos −φ E ~ y (x, y, z, θ,tk ) (23)
2
~ z (r, z, φ, θ,tk ) = E
E ~ z (x, y, z, θ,tk ) (24)
where:
E~ x (x, y, z, θ,tk ) is the electric field at x-direction associated with the
inclined lightning channel,
E~ y (x, y, z, θ,tk ) is the electric field at y-direction associated with the
inclined lightning channel,
E~ z (x, y, z, θ,tk ) is the electric field at z-direction associated with the
inclined lightning channel(vertical electric field),
E~ r (r, z, φ, θ,tk ) is the electric field at r-direction(horizontal) due to in-
clined lightning channel,
E~ ϕ (r, z, φ, θ,tk ) is the electric field at ϕ-direction due to inclined light-
ning channel,
E~ z (r, z, φ, θ,tk ) is electric field at z-direction (vertical) due to inclined
lightning channel,
tk = k∆t.
Computational stability requires the condition specified in Equa-
tion (25). Note that, ∆x = ∆y = ∆z < λe /10 where λe is the wave-
length [16].
1
∆t ≤ q (25)
1 1 1
c× (∆x)2 + (∆y) 2+
(∆z)2

According to proposed algorithm, the current model (in this study

MTLE model) has been entered into magnetic flux density calculations
0)
by i(r0 , t− R(r
c ) parameter in Equations (9) to (11). Therefore,
the current model effect on the electric field calculations has
been established indirectly through magnetic flux density values in
Equations (19) to (21).

4. RESULT AND DISCUSSION

In this section, the proposed equation reflects on the return stroke

current parameters that were obtained from Table 1. Likewise, the
proposed method is validated via the data from three measured fields.
222 Izadi, Kadir, and Gomes

By the same token, new equations are corroborated by the data

gathered from the measured fields from vertical lightning channel when
the channel angle (θ) is assumed to be equal to zero. The magnetic
flux density and the vertical electric fields due to an inclined lightning
channel were estimated using TL, MTLE and MTLL models (Figures 7
and 8, respectively) by considering the current parameters obtained
from Table 1 and a fixed return stroke velocity of v = 1.3×108 m/s.
Figures 7 and 8 illustrate that the magnetic flux density values for
three models are almost the same at close distance from the lightning
channel. However, the simulated vertical electric field due to MTLE
model has higher values than the other two models. On the other
hand, by applying MTLE model as the current model, the effect of
return stroke velocity observation point height on magnetic flux density
and the vertical electric filed are taken into account as illustrated in
Figures 9 and 10, respectively. It should be noted that in order to
see the effect of observation point height, the velocity should be set at
v = 1.3×108 m/s.
As demonstrated in Figures 9 and 10, an increase in velocity
results in an increase in the magnetic flux density values but a decrease
in the vertical electric field values. Similarly, as the observation point
height with respect to the ground increases, the magnetic flux density
values boost up and the vertical electric fields decline. The magnetic
flux density and the vertical electric filed behaviors versus the channel
angle and observation point angle changes are shown in Figures 11 and
12, respectively.

-4
x 10
1.2
TL model
MTLL model
1 MTLE model
Magnetic flux density(Wb/m 2)

0.8

0.6

0.4

0.2

0
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Time(µs)

Figure 7. The Magnetic flux density for different important

engineering models (r = 15 m, θ = 20◦ , φ = 60◦ , z = 0, λ = 2000 m,
H = 7000 m, v = 1.3 × 108 m/s).
Progress In Electromagnetics Research, Vol. 117, 2011 223

4
x 10
7
TL model
MTLL model
6 MTLE model

5
Vertical electric field(v/m)

0
0 5 10 15
Time(µs)

Figure 8. The vertical electric field for different important engineering

models (r = 15 m, θ = 20◦ , φ = 60◦ , z = 0, λ = 2000 m, H = 7000 m,
v = 1.3 × 108 m/s).
-4
x 10
1.2 v=1.5E8 m/s,z=0
v=2E8 m/s,z=0 v=1.3E8 m/s,z=5 m

v=1.3E8 m/s,z=10 m

1
Magnetic flux density(Wb/m2)

v=1.3E8 m/s,z=0
0.8

v=1E8 m/s,z=0
0.6

0.4

0.2

0
0 0.5 1 1.5 2 2.5 3 3.5 4
Time(µs)

Figure 9. The Magnetic flux density behavior versus different return

stroke velocities and observation point heights (r = 15 m, θ = 20◦ ,
φ = 60◦ , λ = 2000 m).

Figure 11 exhibits that by increasing the channel angle (θ) with

respect to z-axis, the values of magnetic flux density fall. On the other
hand, Figure 12 shows that the vertical electric field is directly depend
on θ parameter and by increasing θ, the vertical electric field values are
increased. In other words, whenever the observation point angle (φ)
goes up to 90 degrees (at a constant value of θ), the values of vertical
224 Izadi, Kadir, and Gomes

4
x 10
8

7
v=1.3E8 m/s,z=0
v=1.3E8 m/s,z=5 m v=1E8 m/s,z=0
6
Vertical electric field(v/m)

5
v=1.5E8 m/s,z=0 v=1.3E8 m/s,z=10 m

4
v=2E8 m/s,z=0

0
0 1 2 3 4 5 6 7 8 9 10
Time(µs)

Figure 10. The vertical electric field behavior versus different return
stroke velocities and observation point heights (r = 15 m, θ = 20◦ ,
φ = 60◦ , λ = 2000 m).
-4
x 10
1.2
teta=0,phi=90 degrees

1
teta=20 degrees,phi=0
Magnetic flux density(Wb/m2)

0.8

0.6

teta=30 degrees,phi=90 degrees

0.4 teta=20 degrees,phi=45 degrees

teta=20 degrees,phi=90 degrees

0.2

0
0 1 2 3 4 5 6 7 8 9 10
Time(µs)

Figure 11. The magnetic flux density behavior versus different

channel angles and observation point angles (r = 15 m, z = 0,
λ = 2000 m, v = 1.3 × 108 m/s).

electric field and magnetic flux density drop, as shown in Figures 11 and
12, respectively. In addition, the proposed algorithm is validated using
the three measured fields at different observation point angles (φ) with
respect to y-axis as shown in Figures 13, 15 and 18 for dE dt , Ez and Bϕ
z

respectively as the current parameters are obtained from Table 1 and

the velocity is set at v = 1.3×108 m/s. Additionally, Figure 14 shows
Progress In Electromagnetics Research, Vol. 117, 2011 225

4
x 10
7 teta=20 degrees,phi=0
teta=30 degrees,phi=90 degrees

teta=0,phi=90 degrees
5
Vertical electric field(v/m)

teta=20 degrees,phi=45 degrees

4
teta=20 degrees,phi=90 degrees

0
0 1 2 3 4 5 6 7 8 9 10
Time(µs)

Figure 12. The vertical electric field behavior versus different channel
angles and observation point angles (r = 15 m, z = 0, λ = 2000 m,
v = 1.3 × 108 m/s).

dEz
Figure 13. The measured dt (r = 15 m, z = 0, φ = 135◦ ,
θ = 20◦ ) [6].

a comparison between the measured derivative of vertical electric field

and the time taken from Figure 13 with simulated results from inclined
and vertical channels.
According to Figure 14, an increase is observed in the simulated
values from inclined lightning channel with respect to similar field at
vertical channel. These values are in line and close with each other
comparing to measured fields. In addition, the simulated field from
inclined channel is in agreement with the increase and decrease in time
226 Izadi, Kadir, and Gomes

based on ◦measured field. Moreover, the measured vertical electric field

at φ =60 (Figure 15) illustrates two parts for vertical electric field
due to leader and return stroke (part (a)). So, the part of vertical
electric field associated with return stroke indicated in part (b) is cut
off from measured field The comparison between measured vertical
electric field (Figure 15, part (b)) and simulated vertical electric field

11
x 10
3
inclined channel

2.5

2
vertical channel
dEz/dt(v/m/s)

1.5
measured field

0.5

-0.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time(µs)

Figure 14. Comparison between simulated dE dt due to vertical and

z

inclined lightning channels with measured values (r = 15 m, z = 0,

φ = 135◦ , θ = 20◦ , v = 1.3 × 108 m/s).

(a) (b)

Figure 15. The measured vertical electric field ((a) associated with
leader and return stroke, (b) associated with return stroke, r = 15 m,
z = 0, φ = 60◦ , θ = 20◦ ) [6].
Progress In Electromagnetics Research, Vol. 117, 2011 227

as a result of inclined and vertical lightning channel is exemplified in

Figure 16, whilst the current parameters are obtained from Table 1.
Figure 16 demonstrates that the simulated vertical electric field is
compatible with the inclined lightning channel in contrast to measured
field. The vertical electric field values due to inclined channel have
increased compared to similar fields associated with vertical channel.
In other words, the vertical electric field is more effective in the coupling
models and induced voltage values. Thus, evaluation of accurate
vertical electric field values can be influential in predicting the induced
voltage on the lines and setting appropriate protection levels on the
power lines The simulated magnetic flux densities due to vertical and
inclined channel are compared with those of the measured fields in
Figure 18 when applying channel base current from Table 1 and are
depicted in Figure 17. The figure reveals that the magnetic flux density
due to inclined channel is lesser in value than in similar fields due to
vertical channel; however, the difference between the values is negligible
(See Figure 17). Likewise, the simulated results show an agreement
between the two channels in contrast to measured field.
An algorithm is proposed that can support vertical lightning
channel case when θ is set on zero. the simulated magnetic flux
density, vertical and horizontal electric fields, respectively are shown
in Figures 19, 20 and 21 using the proposed algorithm, while the
observation point is located at 10 m height above ground surface and
the radial distance from vertical lightning channel is equal 9 km as

4
x 10
7

measured field
6

5
Vertical electric field(v/m)

inclined channel

4
vertical channel

0
0 1 2 3 4 5 6 7 8 9 10
Time(µs)

Figure 16. Comparison between simulated vertical electric fields

due to vertical and inclined lightning channels with measured values
(r = 15 m, z = 0, φ = 135◦ , θ = 20◦ , λ = 2000 m, v = 1.3 × 108 m/s).
228 Izadi, Kadir, and Gomes

-4
x 10
1.4

1.2
Magnetic flux density(Wb/m2)

0.8 vertical channel

0.6

inclined channel
0.4
measured field

0.2

0
0 5 10 15
Time(µs)

Figure 17. Comparison between simulated magnetic flux density

due to vertical and inclined lightning channels with measured values
(r = 15 m, z = 0, φ = 45◦ , θ = 20◦ , λ = 2000 m, v = 1.3 × 108 m/s).

Figure 18. The measured magnetic flux density (r = 15 m, z = 0,

φ = 45◦ , θ = 20◦ ) [6].

expressed in the reference [16]. Note that, the current parameters

are Γ11 = 1 µs, Γ12 = 2 µs, Γ21 = 8 µs, Γ22 = 30 µs, i01 = 19.5 kA,
i02 = 12 kA, n1 = 2, n2 = 2, v = 1 × 108 ms , λ = 1500 m, when the
DU channel base current function and the MTLE current model are
applied.
Therefore, the comparison between simulation results using the
proposed algorithm and measured fields and other simulation results
from the method in the reference [16] show that the proposed
algorithms are in agreement. On the other hand, the proposed
Progress In Electromagnetics Research, Vol. 117, 2011 229

-7
x 10
1.4

1.2
Magnetic flux density(Wb/m2)

0.8

0.6

0.4

0.2

0
30 35 40 45 50 55 60 65 70
Time(µs)

Figure 19. The simulated magnetic flux density (r = 9000 m, z = 10,

θ = 0◦ , λ = 1500 m, v = 1 × 108 m/s).
40

35

30
Vertical electric field(v/m)

25

20

15

10

0
30 35 40 45 50 55 60 65 70
Time(µs)

Figure 20. The simulated vertical electric field (r = 9000 m, z = 10,

θ = 0◦ , λ = 1500 m, v = 1 × 108 m/s).

algorithm does not rely on observation point angle (φ) when it is

used in vertical lightning channel basis. The proposed algorithm
is more compatible with different coupling models while different
electromagnetic fields components can be estimated in the time
domain.
Figures 22, 23 and 24 show the radial distance (r), the observation
point angle (φ) and the channel angle (θ) effects on the peak
of electromagnetic field components generated by inclined lightning
channel, respectively while the basic current assumptions are obtained
230 Izadi, Kadir, and Gomes

0.05

0.045

0.04
Horizontal electric field(v/m)

0.035

0.03

0.025

0.02

0.015

0.01

0.005

0
30 30.5 31 31.5 32 32.5 33 33.5 34 34.5 35
Time(µs)

Figure 21. The simulated horizontal electric field (r = 9000 m,

z = 10, θ = 0◦ , λ = 1500 m, v = 1 × 108 m/s).
1.8

1.6 2
dB/dtp*(2*10 ) Wb/m2/s
1.4 11
-dEz/dtp*(2*10 ) v/m/s

1.2
-4
Bp*10 Wb/m2
Field peak

0.8

0.6

0.4

0.2

0
0 100 200 300 400 500 600 700 800 900 1000
r(m)

Figure 22. Radial distance effect on the electromagnetic fields peaks

associated with inclined lightning channel (z = 0, φ = 45◦ , θ = 20◦ ,
λ = 2000 m, v = 1.3 × 108 m/s).

from Table 1.
Figure 22 illustrates that by increasing the radial distance from
lightning channel, peak of electromagnetic field components drop
and this reduction in the near distance from lightning channel
has a greater rate as compared to farther distances. Whilst Figure 23
shows that by increasing the observation point angle (φ) to 90 degrees,
the peaks of electromagnetic fields have a downward trend until about
90 degrees when they are increased almost symmetrically with the first
half. It is important to mention that the value of observation point
Progress In Electromagnetics Research, Vol. 117, 2011 231

3.5

2.5 2
dB/dtp*( 2*10 ) Wb/m2/s
11
Field peak

-dEz/dtp*( 2*10 ) v/m/s

2
-4
Bp *10 Wb/m2
1.5

0.5
0 20 40 60 80 100 120 140 160 180
phi(degree)

Figure 23. Observation point angle effect on the electromagnetic

fields peaks associated with inclined lightning channel (r = 15 m, z = 0,
θ = 20◦ , λ = 2000 m, v = 1.3 × 108 m/s).
3.5

3
2
dB/dtp*(2*10 ) Wb/m2/s
2.5
11
Field peak

-dEz/dtp*(2*10 ) v/m/s
2
-4
Bp*10 Wb/m2

1.5

0.5
0 5 10 15 20 25 30 35 40
Teta(degree)

Figure 24. Channel angle effect on the electromagnetic fields peaks

associated with inclined lightning channel (r = 15 m, z = 0, φ = 45◦ ,
λ = 2000 m, v = 1.3 × 108 m/s).

angle is more effective on the value of x, y and R(r0 ) in Equations (9)

to (11).
Figure 24 illustrates that by increasing the channel angle (θ), the
peak values of magnetic flux density and dB/dt are reduced while
the peak values of dEz/dt have a rising trend with opposite behavior
compared to magnetic flux density.
In order to demonstrate the behavior of electromagnetic fields in
the presence of a channel base current with higher decay-constant
232 Izadi, Kadir, and Gomes

-Ez 800 -9
B*10
(v/m) (Wb/m2)
700
teta=30 degree

600 teta=0
teta=10 degree teta=20 degree
500 teta=20 degree
teta=30 degree teta=10 degree

400

300
teta=0

200

100

0
0.5 1 1.5 2 2.5 3 3.5
Time(µs) x 10

Figure 25. The simulated magnetic flux density (solid line) and
vertical electric field (dot line) at different channel angles (r = 2000 m,
z = 10, φ = 45◦ , λ = 1500 m, v = 1.5 × 108 m/s).

using proposed method, Figure 25 shows behavior of both magnetic

flux density and vertical electric field at different channel angles. Note
that, the current parameters are Γ11 = 2 µs, Γ12 = 4.8 µs, Γ21 = 20 µs,
Γ22 = 26 µs, i01 = 10.5 kA, i02 = 9 kA, n1 = 2, n2 = 2, v = 1.5 × 108 ms ,
λ = 1500 m [16]. The results illustrate that by increasing the channel
angle, the magnetic flux density is decreased while the vertical electric
field is increased. Noted that the results at θ = 0 are also validated
with case 3 in reference [16].
In this study, a method was proposed that posse the following
unique features:
1. The electromagnetic field components can be calculated directly
at the time domain needless of employing Fourier transform, when
a realistic channel base current function is applied (such as Heidler
current function).
2. The proposed method can support different observation point
angles (φ), while a number of previous methods just considered
φ = 0 which causes complications when simultaneous field values
at several points are desired [24].
3. The accuracy of calculated fields can be increased by changing the
values of ∆t, ∆x, ∆y, ∆z.
4. The proposed algorithm is compatible with other coupling models
for the estimation of lightning induced voltage. Since the proposed
method has FDTD structure, while the coupling models are
usually introduced by two partial differential equations the FDTD
method is just a common solution to coupling equations [1, 25–
27]. In addition, the proposed method provides different
Progress In Electromagnetics Research, Vol. 117, 2011 233

electromagnetic field components to be fed to various coupling

models [1].
5. The proposed method is applicable for close and intermediate
distances from the lightning channel.
6. The proposed method can be expanded to the lighting channels
with the shape of a zigzag.
Further, the inclined lightning channel subject can be very beneficial
in increasing the accuracy of the calculated induced voltage in the
coupling models because of these two reasons:
1. Due to the inclined lightning channel, the vertical electric field
values increased compared to the values of the vertical channel.
Thus, it can be more effective on the shape of the induced
voltage wave. Whereas, the lightning induced voltage is directly
dependent on the vertical electric field values in some coupling
models [1, 25, 28–30].
2. Due to inclined lightning channel, the observation point angle
(φ) is an effective factor for the evaluation of electromagnetic
fields. The radial distance between the charge along the lightning
channel and various observation points at different observation
point angles (φ) are different from each other; hence, the inclined
lightning channel can be more effective on the calculation of
lightning induced voltage on the power lines, while the coupling
models rely more on the electromagnetic field components at
different points along the power line [1, 25, 27, 31].

5. CONCLUSION

In this study, the electromagnetic fields due to inclined lightning

channel were estimated using the 2nd FDTD method. An algorithm
was proposed and then validated by three measured fields at a close
distance from lightning channel. The results were compared with
ones in simulated fields from vertical lightning channel. Likewise, the
proposed equations were validated by measured fields from vertical
lightning channel while the channel angle was assumed to be zero.
Moreover, the effects of different geometrical and current parameters
on the electromagnetic fields due to inclined lightning channel were
discussed. It was capable of estimating the electromagnetic fields
associated with the vertical lightning channel in the time domain.
Besides, in order to bond to different coupling models, different
electromagnetic fields components were available in the proposed
method.
234 Izadi, Kadir, and Gomes

REFERENCES

1. Nucci, C. A., “Lightning-induced voltages on overhead power

lines. Part II: Coupling models for the evaluation of the induced
voltages,” Electra, Vol. 162, 121–145, 1995.
2. Nucci, C. A., “Lightning-induced voltages on overhead power
lines. Part I: Return stroke current models with specified channel-
base current for the evaluation of the return stroke electromagnetic
fields,” Electra, Vol. 161, 75–102, 1995.
3. Rakov, V. and M. Uman, “Review and evaluation of lightning
return stroke models including some aspects of their application,”
IEEE Transactions on Electromagnetic Compatibility, Vol. 40,
403–426, 1998.
4. Hill, R., “Analysis of irregular paths of lightning channels,”
Journal of Geophysical Research, Vol. 73, 1897–1906, 1968.
5. Hill, R., “Electromagnetic radiation from erratic paths of lightning
strokes,” Journal of Geophysical Research, Vol. 74, 1922–1929,
1969.
6. Uman, M., J. Schoene, V. A. Rakov, K. J. Rambo, and
G. H. Schnetzer, “Correlated time derivatives of current, electric
field intensity, and magnetic flux density for triggered lightning
at 15 m,” Journal of Geophysical Research, Vol. 107, 4160–4172,
2002.
7. Le Vine, D. and R. Meneghini, “Electromagnetic fields radiated
from a lightning return stroke: Application of an exact solution to
Maxwell’s equations,” Journal of Geophysical Research, Vol. 83,
2377–2384, 1978.
8. Le Vine, D. and R. Meneghini, “Simulation of radiation from
lightning return strokes: The effects of tortuosity,” Radio Science,
Vol. 13, 801–809, 1978.
9. Nucci, C. A., G. Diendorfer, M. A. Uman, F. Rachidi, M. Ianoz,
and C. Mazzetti, “Lightning return-stroke models with channel-
base specified current: A review and comparison,” Journal of
Geophysical Research, Vol. 95, 20395–20408, 1990.
10. Heidler, F., “Analytische blitzstromfunktion zur LEMP-
berechnung,” 18th ICLP, Munich, Germany, 1985.
11. Bruce, C. E. R. and R. H. Golde, “The lightning discharge,”
Journal of the Institute of Electrical Engineering, Vol. 88, 487–
505, 1941.
12. Heidler, F., “Travelling current source model for LEMP
calculation,” Proc. of the 6th Symposium and Technical Exhibition
on Electromagnetic Compability, Zurich, 1985.
Progress In Electromagnetics Research, Vol. 117, 2011 235

13. Cooray, V., The Lightning Flash, IET Press, 2003.

14. Thottappillil, R., V. Rakov, and M. Uman, “Distribution of
charge along the lightning channel: Relation to remote electric
and magnetic fields and to return-stroke models,” Journal of
Geophysical Research, Vol. 102, 6987–7006, 1997.
15. Rakov, V., “lightning return stroke speed,” Journal of lightning
Research, Vol. 1, 2007.
16. Izadi, M., M. Z. A. Ab Kadir, C. Gomes, and W. F. Wan Ahmad,
“An analytical second-FDTD method for evaluation of electric and
magnetic fields at intermediate distances from lightning channel,”
Progress In Electromagnetic Research, Vol. 110, 329–352, 2010.
17. Izadi, M. and M. Z. A. Kadir, “New algorithm for evaluation of
electric fields due to indirect lightning strike,” CMES: Computer
Modeling in Engineering & Sciences, Vol. 67, 1–12, 2010.
18. Schoene, J., “Analysis of parameters of rocket triggered lightning
measured during 1999 and 2000 camp blanding experiment and
modeling of electric and magnetic field derivatives using the
transmission line model,” MSC, University of Florida, 2002.
19. Zhou, X., “On independence completeness of Maxwell’s equations
and uniqueness theorems in electromagnetics,” Progress In
Electromagnetic Research, Vol. 64, 117–134, 2006.
20. Nevels, R. and C. S. Shin, “Lorenz, Lorentz, and the gauge,” IEEE
Antennas & Propagation Magazine, Vol. 43, 70–72, 2001.
21. Eslahchi, M. M. and E. Babolian, “On numerical improvement
of Gauss-Lobatto quadrature rules,” Applied Mathematics and
Computation, Vol. 164, 707–717, 2005.
22. Kreyszig, E., Advanced Engineering Mathematics, Wiley-India,
2007.
23. Sadiku, M. N. O., Numerical Technique in Electromagnetics, CRC
Press, LLC, 2001.
24. Song, T. X., Y. H. Liu, and J. M. Xiong, “Computations of
electromagnetic fields radiated from complex lightning channels,”
Progress In Electromagnetics Research, Vol. 73, 93–105, 2007.
25. Paolone, M., C. A. Nucci, and F. Rachidi, “A new finite difference
time domain scheme for the evaluation of lightning induced
overvoltage on multiconductor overhead lines,” International
Conference on Power System Transients (IPST), 596–602, 2001.
26. Paolone, M., C. A. Nucci, E. Petrache, and F. Rachidi,
“Mitigation of lightning-induced overvoltages in medium voltage
distribution lines by means of periodical grounding of shielding
wires and of surge arresters: Modeling and experimental
236 Izadi, Kadir, and Gomes

validation,” IEEE Transactions on Power Delivery, Vol. 19, 423–

431, 2004.
27. Rachidi, F., “Formulation of the field-to-transmission line
coupling equations interms of magnetic excitation field,” IEEE
Transactions on Electromagnetic Compatibility, Vol. 35, 404–407,
1993.
28. Borghetti, A., A. Morched, F. Napolitano, C. A. Nucci, and
M. Paolone, “Lightning-induced overvoltages transferred through
distribution power transformers,” IEEE Transactions on Power
Delivery, Vol. 24, 360–372, 2008.
29. Borghetti, A., J. Gutierrez, C. A., Nucci, M. Paolone,
E. Petrache, and F. Rachidi, “Lightning-induced voltages on
complex distribution systems: Models, advanced software tools
and experimental validation,” Journal of Electrostatics, Vol. 60,
163–174, 2004.
30. Gomes, C. and M. Z. A. Ab Kadir, “Protection of naval systems
against electromagnetic effects due to lightning,” Progress In
Electromagnetics Research, Vol. 113, 333–349, 2011.
31. Rachidi, F., M. Rubinstein, S. Guerrieri, and C. A. Nucci, “Volt-
ages induced on overhead lines by dart leaders and subsequen-
treturn strokes in natural and rocket-triggered lightning,” IEEE
Transactions on Electromagnetic Compatibility, Vol. 39, 160–166,
1997.