# Full-Wave Analysis of Field-to-Line Coupling Effects Using 1D FDTD Method under Exciting Source with Different Bandwidths

Authors: Qi Zhang; Bihua Zhou; Jianbao Wang; Cheng Gao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

With the aim to analyze field-to-line coupling effects based on energy spectrum, parallel finite-difference time-domain (FDTD) method is applied to calculate the induced voltage on overhead lines under high-power electromagnetic (HPEM) environment. Firstly, the energy distribution laws of HEMP (IEC 61000-2-9), HEMP (Bell Laboratory), HEMP (Paulino et al., 2010), and LEMP (IEC61000-4-5) are given. Due to the air-earth stratified medium, both the absorbing boundary and the connecting boundary applied to scattering by finite-length objects are separately set in aerial and underground parts. Moreover, the influence of line length on induced voltage is analyzed and discussed. The results indicate that the half-peak width is wider with the increase of the line length. But the steepness of induced voltage on the overhead line is invariable. There is no further increase in the peak of induced voltage especially when the line length increases to be equivalent to the wavelength of the frequency bands with the maximum energy.


---

## Body

## 1. Introduction

In recent years, the problem of coupling effects on overhead lines caused by high-power electromagnetic (HPEM) environment has been studied extensively due to the increasing demand by the public for good reliability in power supply and communication [1–4]. In terms of a specific overhead line, although direct strikes have a high probability of producing an insulation flashover, indirect effect caused by high-power electromagnetic pulse is much more frequent [5, 6]. So, it is of significant importance to enhance the line protection against induced voltage. But before this, what should be known is that how much the induced voltage can be produced on the transmission line under HPEM environment and which factors are related to the voltage magnitude.As a matter of fact, it is difficult to get reliable data of induced voltage on overhead lines through actual observation or experiment [7]. Therefore, using a variety of numerical simulation techniques has become one of the important approaches to study HPEM protection [8–14]. In [8], the transmission line (TL) model for analyzing the field-to-line coupling has been explained in detail by Vance, in which the coupling of internal and external electromagnetic field is simplified into current and voltage on line. Then, the internal response caused by electromagnetic fields in the outer space can be solved by the coupling relationship between the transfer impedance and transfer admittance of the shield line. Based on the research of Sali, the transfer impedance model of ideal braided shielded lines has been given in [9], and it has also been verified through actual measurements. In [2], Green’s function for analyzing the coupling effects on overhead lines under the excitation of space electromagnetic fields has been presented by Tesche et al. However, TL model is an approximate method due to the fact that the antenna current, namely, the secondary radiation of the conductor, is ignored.With the rapid development of computer and computational electromagnetics, it is possible to regard field-to-line coupling as an electromagnetic scattering problem which can be solved by time-domain full-wave analysis. The finite-difference time-domain (FDTD) method, which has been extensively studied in the past few years, is an efficient tool applicable to the transient electromagnetic problems [10]. Compared with TL model, the FDTD method is much more accurate and can be used for electromagnetic scattering model with much wider frequency bands. Moreover, the induced voltage and current at any point on overhead lines can be obtained simultaneously at any time, and the physical processes of field-to-line coupling can be wholly described as well [11]. Therefore, the FDTD method is widely applied to analyze the coupling on shielded lines in 3D space. In [12], two-step FDTD method has been used to calculate the electromagnetic fields over lossy ground, and the resultant horizontal and vertical electric fields were deemed to forcing functions to compute the coupling voltage on transmission lines. In addition, a point of view has been presented in [13], that the peak of induced voltage is unchanged when the line length is 10 times longer than the distance between the lightning strike point and the line. However, most of the previous work was focused on the coupling effects under the lightning electromagnetic pulse (LEMP) environment, and less attention has been paid to the relationship between the line configuration and induced voltage under electromagnetic pulses with different frequency bands.This work deals with field-to-line coupling effects under four kinds of HPEM environments with different frequency bands by applying parallel FDTD method. Due to the fact that the computational domain involves two different mediums, both absorbing boundary condition and total fields/scattering fields (TF/SF) connecting boundary condition are set in aerial part and lossy ground separately, in which incident fields are introduced using 1D FDTD method. The influence of the line length on induced voltage is discussed, and the reasons for those changes are also analyzed based on the energy distribution of incident source.The paper is organized as follows. In Section2, the energy spectrum of exciting source is given. In Section 3, the field-to-line coupling model is set up. Besides, the computational methodology is also introduced, which consists of the parallel FDTD technique, extended grid, the setting of the TF/SF connecting boundary condition, and the absorbing boundary condition. In Section 4, the influence of the line length on coupling effects is analyzed and discussed. Finally, general conclusions are presented in Section 5.


## 2. Exciting Source

For the purpose of comparing the coupling effects on overhead lines under incident waves with different bandwidths, HEMP (IEC 61000-2-9), HEMP (Bell Laboratory), HEMP [13], and LEMP (IEC61000-4-5) are selected for calculation.The four waves are denoted by WAVE1, WAVE2, WAVE3, and LEMP, respectively, out of which the WAVE1 is fully defined in IEC 61000-2-9 and consists of three electric field pulses that are referred to as the early-time, intermediate-time, and late-time waveforms, and it is more and more widely adopted in civilian areas. WAVE2 is Bell Laboratory HEMP waveform, which is commonly used for studying system response to electromagnetic pulses. Besides, WAVE3 given in [14] is recommended in our country and LEMP is the standard 8/20 μs pulse defined in IEC61000-4-5.The time-domain expression of the exciting source isE
=
k
E
0
(
e
-
α
t
-
e
-
β
t
), which is a double-exponential shape given for the electric field. The characteristic parameters are listed in Table 1.Table 1
Characteristic parameters of the four exciting sources.
Parameters
Sources
WAVE1
WAVE2
WAVE3
LEMP
k
1.3
1.05
1.04
2.33
α (s−1)
4 × 107
4 × 106
1.5 × 106
7.714 × 104
β (s−1)
6 × 108
4.76 × 108
2.6 × 108
2.489 × 105
t
r (ns)
2.47
4.1
7.78
8 × 103
t
hw (ns)
22.98
184
483.30
20 × 103
Note. The rise time and the half-peak width are denoted by t
r and t
hw, respectively.As listed in Table1, the rise time of LEMP is the longest, and the half-peak width is the widest. On the contrary, WAVE1 is with the shortest rise time and the narrowest half-peak width.The spectrum analysis formula of exciting source is deduced from the time-domain expression; see the following equation:(1)
E
(
i
ω
)
=
k
[
(
α
α
2
+
ω
2
-
β
β
2
+
ω
2
)
+
i
(
-
ω
α
2
+
ω
2
+
ω
β
2
+
ω
2
)
]
.
According to (1), the corresponding cumulative formula of energy spectrum is
(2)
W
(
ω
)
=
∫
0
ω
|
E
(
ω
)
|
2
d
ω
W
T
=
2
π
(
β
-
α
)
[
β
arctan
ω
α
-
α
arctan
ω
β
]
,
where W
T
=
∫
0
+
∞
|
E
(
ω
)
|
2
d
ω is the total energy of incident wave and ω is frequency.The percentages of energy in different frequency bands of the four exciting sources are calculated using (2). The results are listed in Table 2.Table 2
Percentages of energy in different frequency bands of the four exciting sources.
Frequency (Hz)
WAVE1
WAVE2
WAVE3
LEMP
100-101
1.528 × 10−7
1.444 × 10−6
3.841 × 10−6
9.729 × 10−5
101-102
1.528 × 10−6
1.444 × 10−5
3.841 × 10−5
9.729 × 10−4
102-103
1.528 × 10−5
1.444 × 10−4
3.841 × 10−4
9.729 × 10−3
103-104
1.528 × 10−4
1.444 × 10−3
3.841 × 10−3
9.663 × 10−2
104-105
1.528 × 10−3
0.014
3.835 × 10−2
0.626
105-106
0.015
0.141
0.334
0.265
106-107
0.149
0.607
0.534
107-108
0.638
0.218
108-109
0.193
109-1010
Total percent
0.997
0.982
0.915
0.999As listed in Table2, about 99.7% of the energy of WAVE1 is concentrated below 1 GHz, and the corresponding wavelength of the maximum energy band (107 Hz-108 Hz) is several meters to dozens of meters. About 98.2% of the energy of WAVE2 is concentrated below 100 MHz, and the corresponding wavelength of the maximum energy band (106 Hz-107 Hz) is dozens of meters to hundreds of meters. About 91.5% of the energy of WAVE3 is concentrated below 10 MHz, and the corresponding wavelength of the maximum energy band (106 Hz-107 Hz) is dozens of meters to hundreds of meters. About 99.9% of the energy of LEMP is concentrated below 1 MHz, and the corresponding wavelength of the maximum energy band (104 Hz-105 Hz) is several kilometers to dozens of kilometers.


## 3. Computational Methodology

### 3.1. Field-to-Line Coupling Model

A schematic view of the coupling of the field due to HPEM environment to an overhead line is shown in Figure1. The length and height of the overhead line paralleling the ground are denoted as l and h, respectively. R
l is the terminating resistance which is connected with grounding body through the deflectors. To calculate conveniently, the grounding body is replaced with a rectangle structure buried underground, and the depth is 1 m.Figure 1
Schematic representation of coupling to the overhead line under HPEM environment.The wave direction of high-power electromagnetic pulse is determined by angleθ and ψ. A reference vector k
×
z is defined in the equal phase surface, where k is the propagation vector and z is the unit vector along z-axis. Then, polarization direction of the incident wave is determined by the angle between its electric field vector and the reference vector.


### 3.2. Absorbing Boundary Conditions

In order to simulate the process of electromagnetic scattering in the infinite open-domain through the finite grid space, the calculation region of FDTD method is truncated using perfectly matched layer.Since the computational domain is divided into aerial and underground parts, absorbing boundary conditions (ABC) are set separately. Modified perfectly matched layer (MPML) [15] ABC in lossless medium is applied to the aerial part, and MPML in lossy medium is applied to the underground part.


### 3.3. Connecting Boundary Conditions

As shown in Figure2, the numerical space lattice of FDTD method is zoned into distinct regions by connecting boundary. The inner part is total fields region (TF) which is not connected with the absorbing boundary. The outer part is scattering fields region (SF). The incident fields and scattering fields are both included in total fields, and the scattering object also lies in TF. Since the outer boundary of SF is the absorbing boundary, there are only outward propagating waves arriving ABC. Therefore, as long as the information of source is known, the calculation can be started to solve Maxwell equations.Figure 2
Grid division of FDTD for scattering calculation.According to the equivalence principle, the incident fields are introduced only to TF, if the tangential component of the incident electromagnetic fields is set on the connecting boundary. However, in terms of scattering in the air-earth stratified medium, it is unreasonable to consider the initial incident fields as the sole source on connecting boundary due to the influence of ground surface [16]. Therefore, in the half space, in which the incident wave arrives first, the initial incident fields plus Snell’s law reflected component should be included, while, in the subsurface, the incident fields are represented by the transmitted component.Actually, the 3D fields can be converted into a 2D problem in parallel to the incident plane, as shown in Figure3, due to the fact that the incident fields are often regarded as plane wave excitation near the ground.Figure 3
Schematic representation of 2D incident plane with oblique incidence.For the sake of realizing vertical incidence, as shown in Figure4, coordinate system of scattering object needs to be rotated by the angle φ. Then, the calculation of incident fields is performed in the new coordinate system y
′
z
′.Figure 4
Schematic representation of 2D incident plane and connecting boundary with vertical incidence.In the incident plane (y
′-z
′plane), the 2D incident fields can be further converted into 1D problem according to laws of reflection and transmission. By adding a point source in 1D FDTD, the incident fields over the ground are initial incident fields plus the reflected fields, and those underground are transmitted fields. Then, the incident fields in y
-
z plane can be given by time delay in y
′-axis and coordinate transformation. The time delay and coordinate transformation are as follows:(3)
F
2
D
(
y
,
z
,
t
)
=
F
1
D
(
t
-
sin
θ
y
c
,
z
)
,
(4)
[
F
x
F
y
F
z
]
=
[
cos
⁡
φ
sin
φ
0
-
sin
φ
cos
⁡
φ
0
0
0
1
]
[
F
x
′
F
y
′
F
z
′
]
,
where F
x, F
y, and F
z are the field components in y
-
z plane and F
x
′, F
y
′, and F
z
′ are the field components in y
′
-
z
′ plane.Likewise, the 3D incident fields could be obtained by time delay ofF
2
D in x-axis. The time delay is x
sin
θ
cos
⁡
φ
/
c. Moreover, because the delay times may locate between two time steps, linear interpolation formula is applied to calculate each field component. The 1D field equations are given in appendix.


### 3.4. Parallel Technique and Extended Grid

The coupling on overhead lines under HPEM environment is an electromagnetic scattering problem which needs to be solved in a large space. For improving the efficiency, the parallel algorithm is adopted. By dividing the calculation domain according to the number of the CPU, each CPU only needs to calculate one subdomain, which makes that the exchange of tangential field data undertaken only at the boundary of each subregion.In terms of 3D computational domain, data exchange in six surfaces is needed if the area division is performed in three dimensions. Similarly, data exchange involves four neighboring surfaces when the segmentation is conducted in two dimensions. Conceivably, data exchange is only carried out between the current domain and two neighboring areas in 1D segmentation. In analyzing the field-to-line coupling, there is only one direction with large spatial scale in 3D computational domain. Therefore, it is especially suitable to perform 1D area division along the line axis and thereby just a small amount of data needs to be exchanged among computing processes.After the area is split into several subdomains by the parallel FDTD algorithm, both connecting boundary and absorbing boundary are also divided into the subregions which require special handling in the programming. For a desired computational efficiency, the coordinate mapping is applied to reconstruct serial program.As plotted in Figure5, cycle boundaries of each subdomain are converted into i
0 and i
m. The coordinates (i
0~i
m) in each subdomain are then mapped to the coordinates of the initial serial program (A
0~A
n). The initial serial program is still applied to calculate the field component in each subregion by using this method. Only a data exchange function is added, and the modification is only performed on cycle boundaries of each subdomain.Figure 5
Schematic of coordinate mapping.In this study, the computational domain is divided into 16 subdomains according to the threads of the two computers with CPU i7 3770. Take the result under LEMP with the line lengthl
=
4000
m; for example, the required time of using serial program is about 120 min, while that of adopting parallel algorithm under 16 threads on 8 CPU cores is about 16 min. Consequently, the computational time is greatly reduced. Moreover, extended grid could be used along the line axis when the numerical dispersion and Courant stability condition are both satisfied.


## 3.1. Field-to-Line Coupling Model

A schematic view of the coupling of the field due to HPEM environment to an overhead line is shown in Figure1. The length and height of the overhead line paralleling the ground are denoted as l and h, respectively. R
l is the terminating resistance which is connected with grounding body through the deflectors. To calculate conveniently, the grounding body is replaced with a rectangle structure buried underground, and the depth is 1 m.Figure 1
Schematic representation of coupling to the overhead line under HPEM environment.The wave direction of high-power electromagnetic pulse is determined by angleθ and ψ. A reference vector k
×
z is defined in the equal phase surface, where k is the propagation vector and z is the unit vector along z-axis. Then, polarization direction of the incident wave is determined by the angle between its electric field vector and the reference vector.


## 3.2. Absorbing Boundary Conditions

In order to simulate the process of electromagnetic scattering in the infinite open-domain through the finite grid space, the calculation region of FDTD method is truncated using perfectly matched layer.Since the computational domain is divided into aerial and underground parts, absorbing boundary conditions (ABC) are set separately. Modified perfectly matched layer (MPML) [15] ABC in lossless medium is applied to the aerial part, and MPML in lossy medium is applied to the underground part.


## 3.3. Connecting Boundary Conditions

As shown in Figure2, the numerical space lattice of FDTD method is zoned into distinct regions by connecting boundary. The inner part is total fields region (TF) which is not connected with the absorbing boundary. The outer part is scattering fields region (SF). The incident fields and scattering fields are both included in total fields, and the scattering object also lies in TF. Since the outer boundary of SF is the absorbing boundary, there are only outward propagating waves arriving ABC. Therefore, as long as the information of source is known, the calculation can be started to solve Maxwell equations.Figure 2
Grid division of FDTD for scattering calculation.According to the equivalence principle, the incident fields are introduced only to TF, if the tangential component of the incident electromagnetic fields is set on the connecting boundary. However, in terms of scattering in the air-earth stratified medium, it is unreasonable to consider the initial incident fields as the sole source on connecting boundary due to the influence of ground surface [16]. Therefore, in the half space, in which the incident wave arrives first, the initial incident fields plus Snell’s law reflected component should be included, while, in the subsurface, the incident fields are represented by the transmitted component.Actually, the 3D fields can be converted into a 2D problem in parallel to the incident plane, as shown in Figure3, due to the fact that the incident fields are often regarded as plane wave excitation near the ground.Figure 3
Schematic representation of 2D incident plane with oblique incidence.For the sake of realizing vertical incidence, as shown in Figure4, coordinate system of scattering object needs to be rotated by the angle φ. Then, the calculation of incident fields is performed in the new coordinate system y
′
z
′.Figure 4
Schematic representation of 2D incident plane and connecting boundary with vertical incidence.In the incident plane (y
′-z
′plane), the 2D incident fields can be further converted into 1D problem according to laws of reflection and transmission. By adding a point source in 1D FDTD, the incident fields over the ground are initial incident fields plus the reflected fields, and those underground are transmitted fields. Then, the incident fields in y
-
z plane can be given by time delay in y
′-axis and coordinate transformation. The time delay and coordinate transformation are as follows:(3)
F
2
D
(
y
,
z
,
t
)
=
F
1
D
(
t
-
sin
θ
y
c
,
z
)
,
(4)
[
F
x
F
y
F
z
]
=
[
cos
⁡
φ
sin
φ
0
-
sin
φ
cos
⁡
φ
0
0
0
1
]
[
F
x
′
F
y
′
F
z
′
]
,
where F
x, F
y, and F
z are the field components in y
-
z plane and F
x
′, F
y
′, and F
z
′ are the field components in y
′
-
z
′ plane.Likewise, the 3D incident fields could be obtained by time delay ofF
2
D in x-axis. The time delay is x
sin
θ
cos
⁡
φ
/
c. Moreover, because the delay times may locate between two time steps, linear interpolation formula is applied to calculate each field component. The 1D field equations are given in appendix.


## 3.4. Parallel Technique and Extended Grid

The coupling on overhead lines under HPEM environment is an electromagnetic scattering problem which needs to be solved in a large space. For improving the efficiency, the parallel algorithm is adopted. By dividing the calculation domain according to the number of the CPU, each CPU only needs to calculate one subdomain, which makes that the exchange of tangential field data undertaken only at the boundary of each subregion.In terms of 3D computational domain, data exchange in six surfaces is needed if the area division is performed in three dimensions. Similarly, data exchange involves four neighboring surfaces when the segmentation is conducted in two dimensions. Conceivably, data exchange is only carried out between the current domain and two neighboring areas in 1D segmentation. In analyzing the field-to-line coupling, there is only one direction with large spatial scale in 3D computational domain. Therefore, it is especially suitable to perform 1D area division along the line axis and thereby just a small amount of data needs to be exchanged among computing processes.After the area is split into several subdomains by the parallel FDTD algorithm, both connecting boundary and absorbing boundary are also divided into the subregions which require special handling in the programming. For a desired computational efficiency, the coordinate mapping is applied to reconstruct serial program.As plotted in Figure5, cycle boundaries of each subdomain are converted into i
0 and i
m. The coordinates (i
0~i
m) in each subdomain are then mapped to the coordinates of the initial serial program (A
0~A
n). The initial serial program is still applied to calculate the field component in each subregion by using this method. Only a data exchange function is added, and the modification is only performed on cycle boundaries of each subdomain.Figure 5
Schematic of coordinate mapping.In this study, the computational domain is divided into 16 subdomains according to the threads of the two computers with CPU i7 3770. Take the result under LEMP with the line lengthl
=
4000
m; for example, the required time of using serial program is about 120 min, while that of adopting parallel algorithm under 16 threads on 8 CPU cores is about 16 min. Consequently, the computational time is greatly reduced. Moreover, extended grid could be used along the line axis when the numerical dispersion and Courant stability condition are both satisfied.


## 4. Analysis and Discussion

In this section, parallel FDTD method is applied to analyze the model shown in Figure2. By calculating the electromagnetic fields in 3D space and performing the line integral of electric field at the load, the induced voltage on the overhead line is obtained. The variation law of induced voltage, in the case of different line length, is summarized, and the reason for these changes is also given based on the energy spectrum of exciting sources. Suppose that the excitation is vertical-incident electromagnetic wave, and the direction of incident electric field is consistent with the line axis. The computational parameters are listed in Table 3. And the results yielded by WAVE1, WAVE2, WAVE3, and LEMP are shown in Figure 6.Table 3
Parameters used for calculation.
Symbol
Quantity
Value
γ
Polarization angle
0°
θ
Incident angle
30°
φ
Angle between the projection of the propagation vector on the ground surface andx-axis
0°
ε
r
Relative electric permittivity
10
σ
g
Conductivity of lossy ground
0.001 S/m
R
l
Terminating resistance
300 Ω
N
pml
PML depth
10
h
Height of the overhead line
1 mThe changes of induced voltage on the overhead line caused by line length. (a) is the results under WAVE1 with the line lengthl
=
3, 5, 10, 15, 20 m, (b) is the results under WAVE2 with the line length l
=
50, 100, 200, 300, 400 m, (c) is the results under WAVE3 with the line length l
=
100, 200, 400, 600 m, 800 m, and (d) is the results under LEMP with the line length l
=
4000, 8000, 12000, 16000, 20000 m.
(a)
(b)
(c)
(d)Figure6 illustrates that with the increase of line length, the steepness of induced voltage on the overhead line is invariable, but the half-peak width becomes wider. Besides, compared with the wavelength of incident wave, the peak of induced voltage will increase rapidly when the line length is electrically small. However, the peak will tend to be stable, if the line length increases to be equivalent to the wavelength of the frequency bands with the maximum energy. Moreover, it is noteworthy that, apart from the fact that 53.4% of energy is concentrated in (106-107) Hz, as listed in Table 2, the energy of WAVE3 concentrated in (105-106) Hz is still up to 33.4%. Therefore, the peak of induced voltage on the overhead line under WAVE3 is not stable until the line length increases to be nearly 600 m, as shown in Figure 6(c).Above-mentioned results and analysis present that the induced voltage will no longer increase, when the maximum energy of incident wave has been coupled on the overhead line. In addition, the reason why the half-peak width will be wider is that the transmission time of the reflected voltage becomes longer with the increase of the line length, namely, the oscillation period of the induced voltage increases. Therefore, the oscillation phenomenon of the induced voltage waveform under WAVE1 is the most obvious among the four exciting sources, as shown in Figure6(a), due to the shortest line length, which makes the induced voltage reflect back and forth on the overhead line.In terms of the four high-power electromagnetic pulses, the induced voltage on overhead lines caused by LEMP is the strongest. The frequency bands of LEMP, as listed in Table2, is the narrowest. Therefore, the energy in unit band and the total energy of LEMP are both the largest.


## 5. Conclusion

In this paper, parallel FDTD method is used for analyzing the induced voltage on overhead lines under four exciting sources. Due to the overground and underground parts, absorbing boundary conditions in computational domain are set separately, and incident fields on connecting boundary are calculated using 1D FDTD method. Moreover, the rationality of analyzing the field-to-line coupling from the perspective of energy distribution is demonstrated through numerical simulation. And the setting method of absorbing boundary and connecting boundary is also proven to be practicable.



---

*Source: 890181-2014-06-25.xml*