# Electromagnetic Field Coupling to Overhead Wire Configurations: Antenna Model versus Transmission Line Approach

Authors: Dragan Poljak; Khalil El Khamlichi Drissi
Journal: International Journal of Antennas and Propagation
 (2012)

---

## Abstract

The paper deals with two different approaches for the analysis of electromagnetic field coupling to finite length overhead wire: the wire antenna theory (AT) and the transmission line (TL) method. The analysis is carried out in the frequency and time domain, respectively. Within the frequency domain analysis the wire antenna formulation deals with the corresponding set of Pocklington integrodifferential equation, while the transmission line model uses the telegrapher's equations. The set of Pocklington equations is solved via the Galerkin-Bubnov scheme of the Indirect Boundary Element Method (GB-IBEM), while the telegrapher’s equations are treated using the chain matrix method and the modal equation to derive per-unit-length parameters. For the case of the time domain analysis AT model uses the space-time Hallen integral equation set, while TL approach deals with the time domain version of the telegrapher’s equations. Hallen equations are handled via time domain version of GB-IBEM, while time domain telegrapher’s equations are solved by using Finite Difference Time Domain (FDTD) method. Many illustrative computational examples for the frequency and time domain response, respectively, for several configurations of overhead wires, obtained via different approaches, are given in this paper.


---

## Body

## 1. Introduction

The electromagnetic field coupling to overhead wires is of great practical interest for many EMC applications [1–11], such as transient excitation of antennas, power, or communications cables. The electromagnetic field coupling to finite length overhead wires can be determined by means of the transmission line model or the thin wire antenna theory in either frequency or time domain [1]. In particular, the transient response of a wire configuration of interest can be computed directly, by solving the related time domain equations or by the indirect approach, that is, by solving their frequency domain counterpart. When the indirect approach is used the frequency spectrum has to be calculated, and then the transient response is computed by means of the Inverse Fourier Transform (IFT).Many practical engineering problems dealing with electromagnetic field coupling to thin wires can be analyzed by using the Transmission Line (TL) models [1–6]. These models include the analysis of incident electromagnetic field exciting the line and the propagation of induced currents and voltages along the line.The TL models yield valid results if the line length is significantly larger than the separation between the wires and also larger than the actual height above ground [8].On the other hand the TL approximation cannot provide a complete solution for the excitation of a given wire configuration by an incident field if the wavelength of the electromagnetic field exciting a wire structure is comparable to or less than the transverse electrical dimensions of the structure. Namely, the TL model fails to predict resonances and accounts for the presence of a lossy ground only approximately [1]. One of the serious problems with TL approach occurs due to the fact that current grows to infinity at resonant points as there are no losses and radiation resistance to limit its flow [8]. The full wave approach, based on the wire antenna theory and related integral equations, is more rigorous and should be used whenever the above-ground transmission lines of the finite length are considered. However, a serious drawback of AT approach is rather long computational time required for the calculations pertaining to long lines.This paper deals with the analysis of electromagnetic field coupling to overhead wires in either frequency or time domain by using both antenna model and transmission line approach, respectively.A number of illustrative computational examples regarding electromagnetic coupling to overhead wires are given in the paper.The aboveground wires are subjected to electromagnetic fields arriving from a distant source and inducing current to flow along the wires. The key to understanding the behaviour of induced fields is the knowledge of current distribution induced along the wires. These currents generate scattered fields propagating away from the equipment.The paper is organized as follows: Section2 deals with the frequency domain analysis followed by related numerical solution methods for overhead wires. Section 2 ends up with many illustrative examples related to the aboveground lines and PLC (power line communications) systems.Section3 outlines the time domain analysis and related method of solutions of governing equations. Some computational examples pertaining to the multiconductor aboveground lines are given. Finally, the conclusion summarizes what has been discussed throughout this work.


## 2. Frequency Domain Models and Methods

This section deals with the wire antenna theory and transmission line (TL) approximation, respectively, for the analysis of electromagnetic field coupling to overhead lines of finite length in the frequency domain. The formulation arising from the wire antenna theory is based on the set of coupled Pocklington integrodifferential equation for half-space problems. The effect of a two-media configuration is taken into account by means of the reflection coefficient approximation [12]. The resulting integro-differential expressions are numerically handled via the frequency domain Galerkin-Bubnov scheme of the Indirect Boundary Element Method (GB-IBEM) [8].Transmission line model in the frequency domain is based on the corresponding telegrapher's equations which are handled by using the chain matrix method [10].


### 2.1. Antenna Theory Approach: Set of Coupled Pocklington Equations

Modeling of arbitrarily shaped wires located at different heights above a lossy ground is an important task in both antenna and electromagnetic compatibility (EMC) studies [1].This section firstly deals with an assessment of the current induced along multiple wire configurations above a lossy ground. Once the currents along the wire array have been obtained, the radiated field components could be determined.The set of Pocklington equations for a configuration of overhead wires can be obtained as an extension of the Pocklington integro-differential equation for a single wire of arbitrary shape. The Pocklington equation for a single wire above a lossy ground can be derived by enforcing the continuity conditions for the tangential components of the electric field along the perfectly conducting (PEC) wire surface. First, a single wire of arbitrary shape, insulated in free space, as shown in Figure1 is considered.Figure 1
Single wire of arbitrary shape in free space.For the PEC wire the total field composed from the excitation fieldE⃗exc and scattered field E⃗sct vanishes [1, 8]:
(1)e⃗x⋅(E⃗exc+E⃗sct)=0on the wire surface.
Starting from Maxwell’s equations and Lorentz gauge the scattered electric field can be expressed in terms of the vector potential A⃗:
(2)E⃗sct=-jωA⃗+1jωμε∇(∇A⃗).
The vector potential is defined by the particular integral over a given path C (considered conductive wire structure):
(3)A⃗(s)=μ4π∫CI(s′)g0(s,s′,s*)s′⃗ds′,
where I(s′) is the induced current along the line and g0(s,s′) denotes the lossless medium Green function:(4)g0(s,s′)=e-jkRR,
and R is the distance from the source point to the observation point, respectively, while the propagation constant of the homogeneous medium is given by
(5)k2=ω2μ0ε0.
Inserting (3) into (2) gives the relation for the scattered electric field:
(6)E⃗sct=1j4πωε0∫C′I(s′)⋅s′⃗⋅[k2+∇∇]g0(s,s′)ds′.
Combining (6) and (1) results in the Pocklington integral equation for the unknown current distribution along the wire of arbitrary shape insulated in free space:
(7)Etanexc(s)=-1j4πωε0∫C′I(s′)⋅s⃗⋅s′⃗⋅[k2+∇∇]g0(s,s′)ds′,
where Etanexc denotes the tangential component of the electric field illuminating the wire.Now the case of curved wire located above an imperfectly conducting ground can be analyzed by extending integro-differential equation (7) using the reflection coefficient approach [12]. The geometry of an arbitrary wire and its image, respectively, is shown in Figure 2.Figure 2
The wire of arbitrary shape and its image.The excitation functionEexc is now composed from incident and reflected field, respectively,(8)Eexc=Einc+Eref.
Performing certain mathematical manipulations the Pocklington integro-differential equation for a curved wire above a lossy ground becomes [12](9)Esexc(s)=j4πωε0∫0L{[k2e⃗se⃗s′-∂2∂s∂s′]g0(s,s′)+RTM[k2e⃗se⃗s*-∂2∂s∂s*]gi(s,s*)+(RTE-RTM)e⃗se⃗p⋅[k2e⃗pe⃗s*-∂2∂p∂s*]gi(s,s*)}I(s′)ds′,
where  e⃗p is the unit vector normal to the incident plane, while gi(s,s*)  arises from the image theory and is given by
(10)gi(s,s*)=e-jkR*R*,
and R* is the distance from the image source point to the observation point, respectively.An extension to the case of multiple curved wires is straightforward, that is, it follows [12](11)Esmexc(s)=j4πωε0∑n=1Nw∫0Ln{[ke⃗sme⃗sn′-∂2∂sm∂sn′]g0n(sm,sn′)+RTM[k2e⃗sme⃗sn*-∂2∂sm∂sn*]gin(sm,sn*)+(RTE-RTM)e⃗sme⃗p⋅[k2e⃗pe⃗s*-∂2∂p∂s*]gi(sm,sn*)}I(sn′)ds′,
where Nw  is the total number of wires and In(sn′) is the unknown current distribution induced on the nth wire. Furthermore, g0mn(x,x′) and gimn(s,s′)are the Green functions of the form (12)g0mn(sm,sn′)=e-jkR1mnR1mn,gimn(sm,sn′)=e-jkR2mnR2mn,
where  R1mn and R2mn are distances from the source point and from the corresponding image, respectively, to the observation point of interest.The influence of a lossy half-space is taken into account via the Fresnel plane wave reflection coefficient (RC) for TM and TE polarization, respectively [12],
(13)RTM′=n̲cos⁡θ′-n̲-sin2θ′n̲cos⁡θ′+n̲-sin2θ′,(14)RTE=cos⁡θ′-n̲-sin2θ′cos⁡θ′+n̲-sin2θ′,
where θ′ is the angle of incidence and n̲ is given by,
(15)n̲=εeffε0,εeff=εrε0-jσω,
and εeff is the complex permittivity of the ground.For the special case of single horizontal straight wire above a lossy half-space, Figure3 integro-differential equation (9) simplifies into
(16)Exexc=jωμ4π∫0LI(x′)g(x,x′)dx′-1j4πωε∂∂x∫0L∂I(x′)∂x′g(x,x′)dx′,
where I(x′) is the induced current along the horizontal wire and g(x,x′) denotes the Green’s function given by
(17)g(x,x′)=g0(x,x′)-RTMgi(x,x′).Figure 3
Horizontal wire above a lossy ground.Furthermore, if an array of multiple horizontal wires is considered (Figure4), system of (11) becomes
(18)Exmexc=-1j4πωε0∑n=1Nw∫0Ln[∂2∂x2+k12]gmn(x,x′)In(x′)dx′,m=1,2,…,M,
where In(x′) is the unknown current distribution induced along nth wire, Exmexc is the known excitation field tangential to the jth wire surface, and gmn is the corresponding Green function:(19)gmn(x,x′)=g0mn(x,x′)-RTM′gimn(x,x′)m=1,2,…,M.
It is worth noting that a trade-off between the rigorous Sommerfeld integral approach and approximate RC approach is presented in [12]. Although reflection coefficient approximation causes certain error (up to 10%) it takes a significantly less computational effort then a rigorous Sommerfeld approach [8].Figure 4
Horizontal wires above a lossy half-space at different heights.The total electric field irradiated by configuration of multiple wires of arbitrary shape is given by [13, 14]
(20)E⃗=∑n=1Nw[E⃗0n+RTME⃗in+(RTE-RTM)(E⃗in⋅e⃗p)e⃗p],
where
(21)E⃗0n=1j4πωε0[k12∫0Lne⃗sn′I(sn′)g0n(r⃗,r′⃗)dsn′+∫0L∂I(sn′)∂sn′∇g0n(r⃗,r′⃗)dsn′],E⃗in=1j4πωε0[k12∫0Lne⃗sn*I(sn′)gin(r⃗,r⃗*)dw′-∫0Ln∂I(s*)∂sn*∇gin(r⃗,r⃗*)ds′].
Note that index 0 and i are related to the source and image wire, respectively.For the special case of single horizontal straight wire above a lossy half-space (Figure3), it follows [15](22)Ex=1j4πωε0[-∫0L∂I(x′)∂x′∂g(x,x′)∂x′dx′+k2∫0Lg(x,x′)I(x′)dx′],(23)Ey=1j4πωε0∫0L∂I(x′)∂x′∂g(x′,y)∂ydx′,(24)Ez=1j4πωε0∫0L∂I(x′)∂x′∂g(x′,z)∂zdx′.
For the case of multiple horizontal wires the expressions for electric field are given by [15](25)Ex=1j4πωε0∑n=1Nw[-∫0Ln∂In(x′)∂x′∂gnm(x,x′)∂x′dx′+k2∫0Lngnm(x,x′)In(x′)dx′],(26)Ey=1j4πωε0∑n=1Nw∫0Ln∂In(x′)∂x′∂gnm(x,x′)∂ydx′,(27)Ez=1j4πωε0∑n=1Nw∫0Ln∂In(x′)∂x′∂gnm(x,x′)∂zdx′.
The radiated magnetic field of the curved wire system can be written as follows [13, 14]:(28)H⃗=∑n=1Nw[H⃗Sn+RTEH⃗In+(RTM-RTE)(H⃗In⋅e⃗p)e⃗p],
where(29)H⃗Sn=-14π∫0LnI(sn′)e⃗s′×∇g0n(r⃗,r′⃗)ds′,H⃗In=-14π∫0LnI(sn′)e⃗s*×∇gin(r⃗,r⃗*)ds′.
The reduction to the case of a single straight wire or straight wire array is straightforward, as in the case of electric field given by (22)–(27).


### 2.2. Numerical Solution

The set of Pocklington integro-differential equations (11) has been solved by using the Galerkin-Bubnov scheme of the Indirect Boundary Element Method (GB-IBEM). An outline of the method is given here, for the sake of completeness while the method has been presented in detail elsewhere, for example, in [8].Performing the Galerkin-Bubnov scheme of (GB-IBEM) in the frequency domain the set of coupled integro-differential equations (11) is transformed into the following matrix equation [13]
(30)∑n=1M∑i=1Nn[Z]jie{I}ie={V}je,
where the mutual impedance matrix is given by [13]: (31)[Z]ije=-∬-11{D}j{D′}iTg0nm(sn,sm′)dsm′dξ′dξ′dsndξdξ+k12e⃗sne⃗sm∬-11{f}j{f′}iTg0nm(sn,sm′)dsm′dξ′dξ′dsndξdξ-RTM∬-11{D}j{D′}iTginm(sn,sm*)dsm′dξ′dξ′dsndξdξ+RTMk12e⃗sne⃗sm*×∬-11{f}j{f′}iTginm(sn,sm*)dsm′dξ′dξ′dsndξdξ+j4πωε0∫-11ZT′{f}j{f′}jTdsndξdξ,
while the voltage vector is given by [13] (32){V}jn=-j4πωε0∫-11Esnexc(sn)fjn(sn)dsndξdξn.
Once the current distribution is obtained, the radiated field can be obtained applying the similar BEM formalism [13]. Thus, the total field is given by
(33)E⃗=∑k=1N[E⃗eSk+RTME⃗eIk+(RTE-RTM)(E⃗eIk⋅e⃗p)e⃗p],
where the field components due to a wire segment radiation are given by(34)E⃗eSk=1j4πωε0∑i=1n[k2∫-11e⃗ks′Iikefi(ξ)g0k(r⃗,r′⃗)dsk′dξdξ+∫-11Iike∂fi(ξ)∂ξ∇g0k(r⃗,r′⃗)dsk′dξdξ],E⃗eI=1j4πωε0∑i=1n[k2∫-11e⃗ks*Iikefi(ξ)gik(r⃗,r⃗*)dsk′dξdξ-∫-11Iike∂fi(ξ)∂ξ′∇gik(r⃗,r⃗*)dsk′dξdξ].
The total magnetic field is given by [13](35)H⃗=∑k=1N[H⃗eSk+RTEH⃗eIk+(RTM-RTE)(H⃗eIk⋅e⃗p)e⃗p],
while the magnetic field components are given by [13](36)H⃗eSk=-14π∑i=1n∫-11Iikfi(ξ)e⃗sk′×∇g0k(r⃗,r′⃗)dsk′dξdξ,H⃗eIk=-14π∑i=1n∫-11Iikefi(ξ)e⃗ks*×∇gik(r⃗,r⃗*)dsk′dξdξ.The reduction to the case of a single straight wire or straight wire array is straightforward and can be found elsewhere, for example, in [8].


### 2.3. Transmission Line Approximation: Telegrapher’s Equations in the Frequency Domain

Voltages and currents along the multiconductor transmission line shown in Figure4 induced by an external field excitation can be obtained using the field-to-transmission line matrix equations in the frequency domain [10]:
(37)ddx[V̂(x)]+[Ẑ]⋅[Î(x)]=-jωμ0∫0h[Ĥyexc(x,z)]dz,ddx[Î(x)]+[Ŷ]⋅[V̂(x)]=-jωμ0∫0h[Êzexc(x,z)]dz,
where the longitudinal impedance matrix is given by(38)[Ẑ]=jω[L]+[Ẑw]+[Ẑg],
and the transversal admittance matrix can be written as
(39)[Ŷ]=jω[C]+[G],
where [L] is the per-unit-length longitudinal inductance matrix for a perfect soil [C] and [G] are the per-unit length transverse capacitance and conductance matrix of the multiconductor line, respectively. Furthermore, ⌊Ẑw⌋ is the per-unit length internal impedance matrix of the conductors,  and ⌊Ẑg⌋ is the per-unit length ground impedance matrix. Finally, [Ĥyexc(x,z)] and [Êzexc(x,z)] are sources vectors expressed in terms of the incident magnetic and electric field, respectively [1, 8].


### 2.4. Computational Examples

The first computational example is related to the analysis of an overhead wire  (Figure3) of length L=20m, radius a=0.005m located at height h=1 m above PEC ground and illuminated by the plane wave. The amplitude of the electric field is E0=1V/m and it is parallel to x-axis. Figure 5 shows the frequency response at the center of the line. The results computed via GB-IBEM and TL are compared to the results obtained via NEC using RC and Sommerfeld integral approach, respectively, to account for the presence of a lossy half-space. The agreement between the results obtained via the different approaches is found to be satisfactory.Figure 5
Current induced at the center of the line above a PEC ground versus frequency.Figure6 shows the frequency response for the same line located above an imperfectly conducting half-space for various values of ground conductivity σ=1 mS/m. The results calculated via different approaches agree satisfactorily again.Figure 6
Current induced at the center of the line above a lossy ground versus frequency(σ=0.001S/m,εr=10).Next computational example is related to a simple Power Line Communications (PLCs) system. PLC technology aims to provide users with necessary communication means by using the already existing and widely distributed power line network and electrical installations in houses and buildings. However, one of the principal drawbacks of this technology is related to electromagnetic interference (EMI) problems, as overhead power lines at the PLC frequency range (1 MHz to 30 MHz) act as transmitting or receiving antennas, respectively [13].Figure7 shows the geometry of a simple PLC system consisting of two conductors placed in parallel above each other at the distance d. The conductors are suspended between two poles of equal height, thus heaving the shape of the catenary.Figure 7
Simple PLC circuit.The geometry of a catenary is fully defined by such parameters as the distance between the points of suspension,L,the sag of the conductor, s, and the height of the suspension point, h, as shown in Figure 7. The imperfectly conducting ground is characterized with electrical permeability εr and conductivity σ.The conductors are modeled as thin wire antennas excited by the voltage generatorVg at one end and terminated by the load impedance ZL at the other end.The influence of the load impedance is taken into account by modifying continuity condition for the tangential components of the electric field at the wire surface:(40)Esinc+Essct=ZL′I(s),
where ZL′ is the corresponding conductor per length impedance of the conductor.The modified Pocklington equation for the wire containing the load impedance is now given by(41)Einc(s)=-1j4πωε0∫0L{[k12e⃗se⃗s′-∂2∂s∂s′]g0(s,s′)+RTM[k12e⃗se⃗s*-∂2∂s∂s*]gi(s,s*)}×I(s′)ds′+ZL′I(s).
Set of integral equations (41) is numerically solved using via GB-IBEM.The actual example is related to the simple PLC circuit shown in Figure7. The distance between poles is L=200m, with the radii of wires a=6.35mm. The wires are suspended on the poles at heights h1=10m and h2=11m. The maximum sag of the conductor is assumed to be s=2m. Ground parameters are εr=13 and =0.005S/mσ. The power of the applied voltage generator is 2.5 μW (minimum power required for the PLC system operation) and operating frequency is chosen to be 14 MHz. The value of the terminating load ZL is 500 Ω. Figure 8 shows the current distribution along the simple PLC system for different values of sag.Figure 8
The current distribution along a simple PLC system.Radiated electric and magnetic fields at the distance of 30 m from the wires and 10 m above ground are shown at the Figures9 and 10, respectively.Figure 9
Radiated electric field.Figure 10
Radiated magnetic field.Analysis of the radiated field distributions shows that the conductor sag does not influence the far-field region significantly while the near-field distribution is mainly determined by the conductor geometry. Finally, the power of the applied voltage generator is changed to 1 mW (average power used at the actual PLC systems) and operating frequency is varied between 1 and 30 MHz. The values of the terminating loadZL are chosen to be 50 Ω, 500 Ω, 5000 Ω, thus simulating different conditions within the power grid. The maximum values of the radiated electric field at the distance of 30 m for different arrangements are shown in Table 1.Table 1
Maximum values of the radiated electric field at the 30 m distance.
FrequencyZL(Ω)|E|max⁡ (mV/m)7 MHz500.4595000.34150000.38014 MHz500.4775000.45850000.54128 MHz502.3945000.85350002.043According to the available international standards [16, 17], radiated electric fields should not excide level of 30 μV/m at the distance of 30 m. Obviously, the radiated field levels are at best case more than 10 times higher than the proposed limit. The spatial distributions of the radiated electric field have been calculated for the number of frequencies in the frequency range from 1 to 30 MHz. Maximum levels of the calculated electric fields values are shown to excide the limits defined by the standard for the disturbances caused by information technology equipment.


## 2.1. Antenna Theory Approach: Set of Coupled Pocklington Equations

Modeling of arbitrarily shaped wires located at different heights above a lossy ground is an important task in both antenna and electromagnetic compatibility (EMC) studies [1].This section firstly deals with an assessment of the current induced along multiple wire configurations above a lossy ground. Once the currents along the wire array have been obtained, the radiated field components could be determined.The set of Pocklington equations for a configuration of overhead wires can be obtained as an extension of the Pocklington integro-differential equation for a single wire of arbitrary shape. The Pocklington equation for a single wire above a lossy ground can be derived by enforcing the continuity conditions for the tangential components of the electric field along the perfectly conducting (PEC) wire surface. First, a single wire of arbitrary shape, insulated in free space, as shown in Figure1 is considered.Figure 1
Single wire of arbitrary shape in free space.For the PEC wire the total field composed from the excitation fieldE⃗exc and scattered field E⃗sct vanishes [1, 8]:
(1)e⃗x⋅(E⃗exc+E⃗sct)=0on the wire surface.
Starting from Maxwell’s equations and Lorentz gauge the scattered electric field can be expressed in terms of the vector potential A⃗:
(2)E⃗sct=-jωA⃗+1jωμε∇(∇A⃗).
The vector potential is defined by the particular integral over a given path C (considered conductive wire structure):
(3)A⃗(s)=μ4π∫CI(s′)g0(s,s′,s*)s′⃗ds′,
where I(s′) is the induced current along the line and g0(s,s′) denotes the lossless medium Green function:(4)g0(s,s′)=e-jkRR,
and R is the distance from the source point to the observation point, respectively, while the propagation constant of the homogeneous medium is given by
(5)k2=ω2μ0ε0.
Inserting (3) into (2) gives the relation for the scattered electric field:
(6)E⃗sct=1j4πωε0∫C′I(s′)⋅s′⃗⋅[k2+∇∇]g0(s,s′)ds′.
Combining (6) and (1) results in the Pocklington integral equation for the unknown current distribution along the wire of arbitrary shape insulated in free space:
(7)Etanexc(s)=-1j4πωε0∫C′I(s′)⋅s⃗⋅s′⃗⋅[k2+∇∇]g0(s,s′)ds′,
where Etanexc denotes the tangential component of the electric field illuminating the wire.Now the case of curved wire located above an imperfectly conducting ground can be analyzed by extending integro-differential equation (7) using the reflection coefficient approach [12]. The geometry of an arbitrary wire and its image, respectively, is shown in Figure 2.Figure 2
The wire of arbitrary shape and its image.The excitation functionEexc is now composed from incident and reflected field, respectively,(8)Eexc=Einc+Eref.
Performing certain mathematical manipulations the Pocklington integro-differential equation for a curved wire above a lossy ground becomes [12](9)Esexc(s)=j4πωε0∫0L{[k2e⃗se⃗s′-∂2∂s∂s′]g0(s,s′)+RTM[k2e⃗se⃗s*-∂2∂s∂s*]gi(s,s*)+(RTE-RTM)e⃗se⃗p⋅[k2e⃗pe⃗s*-∂2∂p∂s*]gi(s,s*)}I(s′)ds′,
where  e⃗p is the unit vector normal to the incident plane, while gi(s,s*)  arises from the image theory and is given by
(10)gi(s,s*)=e-jkR*R*,
and R* is the distance from the image source point to the observation point, respectively.An extension to the case of multiple curved wires is straightforward, that is, it follows [12](11)Esmexc(s)=j4πωε0∑n=1Nw∫0Ln{[ke⃗sme⃗sn′-∂2∂sm∂sn′]g0n(sm,sn′)+RTM[k2e⃗sme⃗sn*-∂2∂sm∂sn*]gin(sm,sn*)+(RTE-RTM)e⃗sme⃗p⋅[k2e⃗pe⃗s*-∂2∂p∂s*]gi(sm,sn*)}I(sn′)ds′,
where Nw  is the total number of wires and In(sn′) is the unknown current distribution induced on the nth wire. Furthermore, g0mn(x,x′) and gimn(s,s′)are the Green functions of the form (12)g0mn(sm,sn′)=e-jkR1mnR1mn,gimn(sm,sn′)=e-jkR2mnR2mn,
where  R1mn and R2mn are distances from the source point and from the corresponding image, respectively, to the observation point of interest.The influence of a lossy half-space is taken into account via the Fresnel plane wave reflection coefficient (RC) for TM and TE polarization, respectively [12],
(13)RTM′=n̲cos⁡θ′-n̲-sin2θ′n̲cos⁡θ′+n̲-sin2θ′,(14)RTE=cos⁡θ′-n̲-sin2θ′cos⁡θ′+n̲-sin2θ′,
where θ′ is the angle of incidence and n̲ is given by,
(15)n̲=εeffε0,εeff=εrε0-jσω,
and εeff is the complex permittivity of the ground.For the special case of single horizontal straight wire above a lossy half-space, Figure3 integro-differential equation (9) simplifies into
(16)Exexc=jωμ4π∫0LI(x′)g(x,x′)dx′-1j4πωε∂∂x∫0L∂I(x′)∂x′g(x,x′)dx′,
where I(x′) is the induced current along the horizontal wire and g(x,x′) denotes the Green’s function given by
(17)g(x,x′)=g0(x,x′)-RTMgi(x,x′).Figure 3
Horizontal wire above a lossy ground.Furthermore, if an array of multiple horizontal wires is considered (Figure4), system of (11) becomes
(18)Exmexc=-1j4πωε0∑n=1Nw∫0Ln[∂2∂x2+k12]gmn(x,x′)In(x′)dx′,m=1,2,…,M,
where In(x′) is the unknown current distribution induced along nth wire, Exmexc is the known excitation field tangential to the jth wire surface, and gmn is the corresponding Green function:(19)gmn(x,x′)=g0mn(x,x′)-RTM′gimn(x,x′)m=1,2,…,M.
It is worth noting that a trade-off between the rigorous Sommerfeld integral approach and approximate RC approach is presented in [12]. Although reflection coefficient approximation causes certain error (up to 10%) it takes a significantly less computational effort then a rigorous Sommerfeld approach [8].Figure 4
Horizontal wires above a lossy half-space at different heights.The total electric field irradiated by configuration of multiple wires of arbitrary shape is given by [13, 14]
(20)E⃗=∑n=1Nw[E⃗0n+RTME⃗in+(RTE-RTM)(E⃗in⋅e⃗p)e⃗p],
where
(21)E⃗0n=1j4πωε0[k12∫0Lne⃗sn′I(sn′)g0n(r⃗,r′⃗)dsn′+∫0L∂I(sn′)∂sn′∇g0n(r⃗,r′⃗)dsn′],E⃗in=1j4πωε0[k12∫0Lne⃗sn*I(sn′)gin(r⃗,r⃗*)dw′-∫0Ln∂I(s*)∂sn*∇gin(r⃗,r⃗*)ds′].
Note that index 0 and i are related to the source and image wire, respectively.For the special case of single horizontal straight wire above a lossy half-space (Figure3), it follows [15](22)Ex=1j4πωε0[-∫0L∂I(x′)∂x′∂g(x,x′)∂x′dx′+k2∫0Lg(x,x′)I(x′)dx′],(23)Ey=1j4πωε0∫0L∂I(x′)∂x′∂g(x′,y)∂ydx′,(24)Ez=1j4πωε0∫0L∂I(x′)∂x′∂g(x′,z)∂zdx′.
For the case of multiple horizontal wires the expressions for electric field are given by [15](25)Ex=1j4πωε0∑n=1Nw[-∫0Ln∂In(x′)∂x′∂gnm(x,x′)∂x′dx′+k2∫0Lngnm(x,x′)In(x′)dx′],(26)Ey=1j4πωε0∑n=1Nw∫0Ln∂In(x′)∂x′∂gnm(x,x′)∂ydx′,(27)Ez=1j4πωε0∑n=1Nw∫0Ln∂In(x′)∂x′∂gnm(x,x′)∂zdx′.
The radiated magnetic field of the curved wire system can be written as follows [13, 14]:(28)H⃗=∑n=1Nw[H⃗Sn+RTEH⃗In+(RTM-RTE)(H⃗In⋅e⃗p)e⃗p],
where(29)H⃗Sn=-14π∫0LnI(sn′)e⃗s′×∇g0n(r⃗,r′⃗)ds′,H⃗In=-14π∫0LnI(sn′)e⃗s*×∇gin(r⃗,r⃗*)ds′.
The reduction to the case of a single straight wire or straight wire array is straightforward, as in the case of electric field given by (22)–(27).


## 2.2. Numerical Solution

The set of Pocklington integro-differential equations (11) has been solved by using the Galerkin-Bubnov scheme of the Indirect Boundary Element Method (GB-IBEM). An outline of the method is given here, for the sake of completeness while the method has been presented in detail elsewhere, for example, in [8].Performing the Galerkin-Bubnov scheme of (GB-IBEM) in the frequency domain the set of coupled integro-differential equations (11) is transformed into the following matrix equation [13]
(30)∑n=1M∑i=1Nn[Z]jie{I}ie={V}je,
where the mutual impedance matrix is given by [13]: (31)[Z]ije=-∬-11{D}j{D′}iTg0nm(sn,sm′)dsm′dξ′dξ′dsndξdξ+k12e⃗sne⃗sm∬-11{f}j{f′}iTg0nm(sn,sm′)dsm′dξ′dξ′dsndξdξ-RTM∬-11{D}j{D′}iTginm(sn,sm*)dsm′dξ′dξ′dsndξdξ+RTMk12e⃗sne⃗sm*×∬-11{f}j{f′}iTginm(sn,sm*)dsm′dξ′dξ′dsndξdξ+j4πωε0∫-11ZT′{f}j{f′}jTdsndξdξ,
while the voltage vector is given by [13] (32){V}jn=-j4πωε0∫-11Esnexc(sn)fjn(sn)dsndξdξn.
Once the current distribution is obtained, the radiated field can be obtained applying the similar BEM formalism [13]. Thus, the total field is given by
(33)E⃗=∑k=1N[E⃗eSk+RTME⃗eIk+(RTE-RTM)(E⃗eIk⋅e⃗p)e⃗p],
where the field components due to a wire segment radiation are given by(34)E⃗eSk=1j4πωε0∑i=1n[k2∫-11e⃗ks′Iikefi(ξ)g0k(r⃗,r′⃗)dsk′dξdξ+∫-11Iike∂fi(ξ)∂ξ∇g0k(r⃗,r′⃗)dsk′dξdξ],E⃗eI=1j4πωε0∑i=1n[k2∫-11e⃗ks*Iikefi(ξ)gik(r⃗,r⃗*)dsk′dξdξ-∫-11Iike∂fi(ξ)∂ξ′∇gik(r⃗,r⃗*)dsk′dξdξ].
The total magnetic field is given by [13](35)H⃗=∑k=1N[H⃗eSk+RTEH⃗eIk+(RTM-RTE)(H⃗eIk⋅e⃗p)e⃗p],
while the magnetic field components are given by [13](36)H⃗eSk=-14π∑i=1n∫-11Iikfi(ξ)e⃗sk′×∇g0k(r⃗,r′⃗)dsk′dξdξ,H⃗eIk=-14π∑i=1n∫-11Iikefi(ξ)e⃗ks*×∇gik(r⃗,r⃗*)dsk′dξdξ.The reduction to the case of a single straight wire or straight wire array is straightforward and can be found elsewhere, for example, in [8].


## 2.3. Transmission Line Approximation: Telegrapher’s Equations in the Frequency Domain

Voltages and currents along the multiconductor transmission line shown in Figure4 induced by an external field excitation can be obtained using the field-to-transmission line matrix equations in the frequency domain [10]:
(37)ddx[V̂(x)]+[Ẑ]⋅[Î(x)]=-jωμ0∫0h[Ĥyexc(x,z)]dz,ddx[Î(x)]+[Ŷ]⋅[V̂(x)]=-jωμ0∫0h[Êzexc(x,z)]dz,
where the longitudinal impedance matrix is given by(38)[Ẑ]=jω[L]+[Ẑw]+[Ẑg],
and the transversal admittance matrix can be written as
(39)[Ŷ]=jω[C]+[G],
where [L] is the per-unit-length longitudinal inductance matrix for a perfect soil [C] and [G] are the per-unit length transverse capacitance and conductance matrix of the multiconductor line, respectively. Furthermore, ⌊Ẑw⌋ is the per-unit length internal impedance matrix of the conductors,  and ⌊Ẑg⌋ is the per-unit length ground impedance matrix. Finally, [Ĥyexc(x,z)] and [Êzexc(x,z)] are sources vectors expressed in terms of the incident magnetic and electric field, respectively [1, 8].


## 2.4. Computational Examples

The first computational example is related to the analysis of an overhead wire  (Figure3) of length L=20m, radius a=0.005m located at height h=1 m above PEC ground and illuminated by the plane wave. The amplitude of the electric field is E0=1V/m and it is parallel to x-axis. Figure 5 shows the frequency response at the center of the line. The results computed via GB-IBEM and TL are compared to the results obtained via NEC using RC and Sommerfeld integral approach, respectively, to account for the presence of a lossy half-space. The agreement between the results obtained via the different approaches is found to be satisfactory.Figure 5
Current induced at the center of the line above a PEC ground versus frequency.Figure6 shows the frequency response for the same line located above an imperfectly conducting half-space for various values of ground conductivity σ=1 mS/m. The results calculated via different approaches agree satisfactorily again.Figure 6
Current induced at the center of the line above a lossy ground versus frequency(σ=0.001S/m,εr=10).Next computational example is related to a simple Power Line Communications (PLCs) system. PLC technology aims to provide users with necessary communication means by using the already existing and widely distributed power line network and electrical installations in houses and buildings. However, one of the principal drawbacks of this technology is related to electromagnetic interference (EMI) problems, as overhead power lines at the PLC frequency range (1 MHz to 30 MHz) act as transmitting or receiving antennas, respectively [13].Figure7 shows the geometry of a simple PLC system consisting of two conductors placed in parallel above each other at the distance d. The conductors are suspended between two poles of equal height, thus heaving the shape of the catenary.Figure 7
Simple PLC circuit.The geometry of a catenary is fully defined by such parameters as the distance between the points of suspension,L,the sag of the conductor, s, and the height of the suspension point, h, as shown in Figure 7. The imperfectly conducting ground is characterized with electrical permeability εr and conductivity σ.The conductors are modeled as thin wire antennas excited by the voltage generatorVg at one end and terminated by the load impedance ZL at the other end.The influence of the load impedance is taken into account by modifying continuity condition for the tangential components of the electric field at the wire surface:(40)Esinc+Essct=ZL′I(s),
where ZL′ is the corresponding conductor per length impedance of the conductor.The modified Pocklington equation for the wire containing the load impedance is now given by(41)Einc(s)=-1j4πωε0∫0L{[k12e⃗se⃗s′-∂2∂s∂s′]g0(s,s′)+RTM[k12e⃗se⃗s*-∂2∂s∂s*]gi(s,s*)}×I(s′)ds′+ZL′I(s).
Set of integral equations (41) is numerically solved using via GB-IBEM.The actual example is related to the simple PLC circuit shown in Figure7. The distance between poles is L=200m, with the radii of wires a=6.35mm. The wires are suspended on the poles at heights h1=10m and h2=11m. The maximum sag of the conductor is assumed to be s=2m. Ground parameters are εr=13 and =0.005S/mσ. The power of the applied voltage generator is 2.5 μW (minimum power required for the PLC system operation) and operating frequency is chosen to be 14 MHz. The value of the terminating load ZL is 500 Ω. Figure 8 shows the current distribution along the simple PLC system for different values of sag.Figure 8
The current distribution along a simple PLC system.Radiated electric and magnetic fields at the distance of 30 m from the wires and 10 m above ground are shown at the Figures9 and 10, respectively.Figure 9
Radiated electric field.Figure 10
Radiated magnetic field.Analysis of the radiated field distributions shows that the conductor sag does not influence the far-field region significantly while the near-field distribution is mainly determined by the conductor geometry. Finally, the power of the applied voltage generator is changed to 1 mW (average power used at the actual PLC systems) and operating frequency is varied between 1 and 30 MHz. The values of the terminating loadZL are chosen to be 50 Ω, 500 Ω, 5000 Ω, thus simulating different conditions within the power grid. The maximum values of the radiated electric field at the distance of 30 m for different arrangements are shown in Table 1.Table 1
Maximum values of the radiated electric field at the 30 m distance.
FrequencyZL(Ω)|E|max⁡ (mV/m)7 MHz500.4595000.34150000.38014 MHz500.4775000.45850000.54128 MHz502.3945000.85350002.043According to the available international standards [16, 17], radiated electric fields should not excide level of 30 μV/m at the distance of 30 m. Obviously, the radiated field levels are at best case more than 10 times higher than the proposed limit. The spatial distributions of the radiated electric field have been calculated for the number of frequencies in the frequency range from 1 to 30 MHz. Maximum levels of the calculated electric fields values are shown to excide the limits defined by the standard for the disturbances caused by information technology equipment.


## 3. Time Domain Models and Methods

This section deals with direct time domain analysis of transient electromagnetic field coupling to straight overhead wires using the wire antenna theory and the transmission line method, respectively. The time domain antenna theory formulation is based on a set of the space-time Hallen integral equations. The transmission line approximation is based on the corresponding time domain Telegrapher's equations. The space-time integral equations arising from the wire antenna theory are handled by the time domain scheme of GB-IBEM. The time domain Telegrapher's equations are solved using the Finite Difference Time Domain (FDTD) method. Time domain numerical results obtained with both approaches are compared to the results computed via NEC 2 code combined with Inverse Fourier Transform procedure. Some illustrative comparisons of results obtained by means of antenna theory and transmission line approach are presented in this section.It is worth mentioning that, for the sake of simplicity, only straight wires are analyzed in this paper.


### 3.1. The Antenna Theory Model

Generally, a direct time-domain analysis of thin wire in the presence of a lossy half-space can be carried out via the appropriate space-time integral equations of either Pocklington or Hallen type [1, 8]. When applied to the solution of the Hallen integral equation the Galerkin-Bubnov Indirect Boundary Element Method (GB-IBEM) [8] results in relatively complex procedures compared to various procedures for the solution of Pocklington equations, but, at the same time, it is proven to be highly efficient, accurate, and unconditionally stable  [8, 18, 19]. On the other hand, the implementation of GB-IBEM to the solution of the Pocklington-type equation is relatively simple, but suffers from serious numerical instabilities. The origin of these instabilities is the discretization of space-time differential operator [8]. The GB-IBEM solution of the Pocklington equation in free space for certain values of time domain integration parameters has been presented elsewhere, for example, in [19], while the Hallen integral equation solution by means of GB-IBEM has been obtained for thin wire structures in the presence of a dielectric half-space, for example, in [11]. In both cases, the influence of imperfect ground has been taken into account via the corresponding reflection coefficient. The numerical solution was mostly limited to scenarios in which the finite conductivity of the ground could be ignored. This approximation involves cases where the wires are sufficiently far from the two-media interface or where the ground conductivity is appreciably low or very high, that is, where the approximation of pure dielectric medium or perfect ground is applied. Through these approximations the time-dependent part of the reflection coefficient function vanishes, and the resulting matrix equation simplifies significantly.However, for the cases where these approximations are not valid, modifications to the original methods are required in order to include the ground conductivity [8]. Namely, the related reflection coefficient is space-time dependent, and the resulting convolution integrals have to be included in the matrix system and numerically computed. This leads to a significant increase in the overall computational cost of the method and consequently requires several modifications.This section deals with the transient analysis of multiple horizontal wires above a lossy ground using the Hallen integral equation approach.The set of space-time Hallen's integral equations can be derived as an extension of the single wire case. First, a single wire insulated in free space is considered.Thin wire antenna or scatterer of lengthL and radius a, oriented along the x-axis, is considered. The wire is assumed to be perfectly conducting and excited by a plane wave electric field. For the sake of simplicity, the analysis is restricted to the case of a normally incident electric field.The tangential component of the total field vanishes on the PEC wire surface, that is,(42)Exinc+Exsct=0,
where Exinc is the incident and Exsct scattered field on the metallic wire surface. Starting from Maxwell equations and obeying the Lorentz gauge one obtains a time domain version of  (2):(43)(∂2A⃗∂t2-1με∇(∇A⃗))|tan=∂E⃗inc∂t|tan,
where A⃗ is the space-time-dependent vector potential.According to the thin wire approximation, only the axial component of the vector potential exists, that is, (43) becomes  (44)∂2Ax∂x2-1c2∂2Ax∂t2=-1c2∂Exinc∂t,
where c denotes the velocity of light.The corresponding solution of (44) can be expressed in terms of a sum of the general solution of the homogeneous equation and the particular solution of the inhomogeneous equation:(45)Ax(x,t)=Axh(x,t)+Axp(x,t).
The solution of the homogeneous wave equation is given as a superposition of incident and reflected wave [8]:(46)Axh(x,t)=F1(t-xc)+F2(t+xc),
while the particular solution is given by the integral [8]:(47)Axp(x,t)=12Z0∫0LExinc(x′,t-|x-x′|c)dx′,
where L denotes the total antenna length.On the other hand, the magnetic vector potential on the PEC wire surface is given by the particular integral:(48)Ax(x,t)=μ4π∫SI(x′,t-R/c)Rdx′.Combining (45)–(48) yields the space-time Hallen equation:(49)∫0LI(x′,t-R/c)4πRdx′=F0(t-xc)+FL(t-L-xc)+12Z0∫0LExinc(x′,t-|x-x′|c)dx′,
where I(x′) is the equivalent axial current to be determined, Exinc is the known tangential incident field, R=[(x-x′)2+a2]1/2 is the distance from the source point (the equivalent current in the antenna axis) to the observation point, and Z0 is the wave impedance of a free space.The unknown functionsF0(t) and FL(t) account for the multiple reflections of the current at the free ends of the wire.A direct time formulation for a straight thin wire above a dissipative half-space can be obtained as the extension of the free-space Hallen equation (49).The free space Hallen equation (49) is first transferred into the Laplace frequency domain:(50)∫0LI(x′,s)e-sR/c4πRdx′=F0(s)e-sx/c+FL(s)e-s((L-x)/c)+12Z0∫0LExinc(x′,s)e-s|x-x′|/cdx′,
where s=jω is the Laplace variable.According to the theory of images the free space integral equation (50) is extended by an additional term multiplying the Green function of the image source by space-frequency-dependent reflection coefficient RTM(θ′,s) for TM polarization. The integral equation in the frequency domain is given by(51)∫0LI(x,s)esR/c4πRdx′-∫0LRTM(θ,s)I(x,s)esR*/c4πR*dx′=F0(s)e-sx/c+FL(s)e-s((L-x)/c)+12Z0∫0LExexc(x′,s)e-s(|x-x′|/c)dx′,
where R*=(x-x′)2+4h2 and  RTM(θ′,s) is determined by the expression [1](52)RTM(θ′,s)=εr(1+σ/εs)cos⁡⁡θ′-εr(1+σ/εs)-sin⁡2⁡θ′εr(1+σ/εs)cos⁡⁡θ′+εr(1+σ/εs)-sin⁡2⁡θ′,
where σ and  ε are the lossy medium conductivity and permittivity, respectively, andθ′=arctg(|x-x′|/2h).The reflection coefficient (RC) approach is a satisfactory approximation in half-space calculations, as long as the field is calculated far away from the source, and the imperfect ground, respectively, to ensureθ′<π/2 [8].Performing the convolution, the time domain counterpart of (52) is obtained in the form(53)∫0LI(x′,t-R/c)4πRdx′-∫-∞t∫0Lr(θ,τ)I(x′,t-R*/c-τ)4πR*dx′dτ=12Z0∫0LExexc(x′,t-|x-x′|c)dx′+F0(t-xc)+FL(t-L-xc),
where r(θ,t) is the space-time reflection coefficient which, for convenience, can be written in the form [18] (54)r(θ,τ)=r′(θ,τ)+r′′(θ,τ),
where(55)r′(θ,t)=Kδ(t),r′′(θ,t)=4β1-β2e-αtt∑n=1∞(-1)n+1nKnIn(αt),τ=σε0εr,β=εr-sin2θεrcos⁡θ,γ=τ1-(sin2θ/εr),θ=arctg|x-x′|2h,K=1-β1+β,α=τ2.Note thatIn is the modified Bessel function of the first order and nth degree.For the case of normal incidence, which is considered for the sake of simplicity, the excitation term is given by(56)Exexc(t)=Exinc(t)-Exref(t*),
where t*=t-R*/c.The transient ground-reflected field is obtained as the convolution of the incident field and the space-time reflection coefficient for the angle of incidenceθ=0 (in accordance with the parallel incidence of the electric field), as is proposed in [20](57)Exref(t)=∫-∞tExinc(t-τ)r(θ=0,τ)dτ
and the integral equation (53) becomes
(58)∫0LI(x′,t-R/c)4πRdx′-∫-∞t∫0Lr(θ,τ)I(x′,t-R*/c-τ)4πR*dx′dτ=F0(t-xc)+FL(t-L-xc)+12Z0∫0LExinc(x′,t-|x-x′|c)dx′-12Z0∫-∞t∫0LExinc(x′,t-|x-x′|c-τ)×r(θ=0,τ)dx′dτ.The unknown time functionsF0(t),FL(t),F0(t-(L/c)), and FL(t-(L/c))  can be obtained in the same manner, as in the case of free space, in terms of as auxilliary functions K0(t) and KL(t) [8]:(59)F0(t)=∑n=0∞K0(t-2nLc)-∑n=0∞KL(t-(2n+1)Lc),FL(t)=∑n=0∞KL(t-2nLc)-∑n=0∞K0(t-(2n+1)Lc),
where(60)KL(t)=∫0LI(x′,t-R0/c)4πR0dx′-∫-∞t∫0Lr(θ′,τ)I(x′,t-R0*/c-τ)4πR0*dx′dτ-12Z0∫0LExinc(x′,t-x′c)dx′,KL(t)=∫0LI(x′,t-RL/c)4πRLdx′-∫-∞t∫0Lr(θ′,τ)I(x′,t-RL*/c-τ)4πRL*dx′dτ-12Z0∫0LExexc(x′,t-L-x′c)dx′,
while R0 and RL are the distances from the wire ends to the source point and R0*,RL* are the distances from the image wire ends to the image source point.If the case of a perfect (ideal) dielectric half-space is considered, the Hallen equation (53) simplifies into(61)∫0LI(x′,t-R/c)4πRdx′-∫0Lr(θ)I(x′,t-R*/c)4πR*dx′=F0(t-xc)+FL(t-L-xc)+12Z0∫0LExexc(x′,t-|x-x′|c)dx′.
Space-time integral equation (53) or (61), respectively, can be solved assuming the zero current at the free ends of the wire and with the initial conditions requiring the wire not to be excited before the certain instant t=t0  [8].The transient behavior ofM straight horizontal thin wires located at different heights above an infinite ground plane is determined by a set of the coupled space-time integral equations of the Hallen type [11]:(62)∑s=1M∫x0sxLsIs(x′,t-Rvs/c)4πRvsdx′-∑s=1M∫-∞t∫x0sxLsrvs(θ,τ)Is(x′,t-Rvs*/c-τ)4πRvs*dx′dτ=F0v(t-x-x0vc)+FLv(t-xLv-xc)+12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
where v,s=1,2,…,M denote the index of the observed and source wire, respectively. Furthermore, Ls and Lv are the lengths of the sth and vth wire, and x′,x are the x-coordinates of the source and observation points on respective wires. The distances between observation point (x,y,z) on the wire v and source point (x′,y′,z′) on the wire s are given by(63)Rvs={(x-x′)2+(y-y′)2+(z-z′)2:v≠s(x-x′)2+a2:v=s,Rvs*=(x-x′)2+(y-y′)2+(z+z′)2,
where asterisk is related to source points are located on the image wire.Unknown time signalsF0v(t) and FLv(t) account for the multiple reflections of transient currents at the wire open ends and can be written in the form(64)F0v(t)=∑n=0∞K0v(t-2nLvc)-∑n=0∞KLv(t-(2n+1)Lvc),(65)FLv(t)=∑n=0∞KLv(t-2nLvc)-∑n=0∞K0v(t-(2n+1)Lvc),
while the auxiliary functions K are defined, as follows:(66)K0v(t)=∑s=1M∫x0sxLsIs(x′,t-Rvs(0)/c)4πRvs(0)dx′-∑s=1M∫-∞t∫x0sxLsrvs(θ,τ)Is(x′,t-Rvs*(0)/c-τ)4πRvs*(0)dx′dτ-12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,(67)KLv(t)=∑s=1M∫x0sxLsIs(x′,t-Rvs(L)/c)4πRvs(L)dx′-∑s=1M∫-∞t∫x0sxLsrvs(θ,τ)Is(x′,t-Rvs*(L)/c-τ)4πRvs*(L)dx′dτ-12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
where Rvs(0) and Rvs(L) are distances from considered source point on each wire s to a corresponding observation point at the ends of the wire v:(68)Rvs(0)=Rvs|x=x0v,Rvs(L)=Rvs|x=xLv,
while Rvs*(0) and Rvs*(L) are distances between the source point at the image of the wire sand observation point located at the ends of the wire v:(69)Rvs*(0)=Rvs*|x=x0v,Rvs*(L)=Rvs*|x=xLv.The space-time reflection coefficientrvs(θ,t) accounts for the influence of the interface and is given by [11] (70)rvs(θvs′,t)=Aδ(t),
where(71)A=1-β1+β,β=εr-sin2θ′εrcos⁡θ′,θvs′=Arctg(x′-x)2+(y′-y)2z′+z.The angleθvs′is the angle between the source point on the image of the wire s(x′,y′,-z′) and the observation point (x,y,z) on wire v.Substituting (70) into (62) yields(72)∑s=1M∫x0sxLsIs(x′,t-Rvs/c)4πRvsdx′-∑s=1M∫x0sxLsrvs(θ)Is(x′,t-Rvs*/c)4πRvs*dx′=F0v(t-x-x0vc)+FLv(t-xLv-xc)+12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
where Exvexc is the space-time-dependent tangential electric field on the vth wire.For the case of normal incidence the total excitation fieldExvexc(x′,t) is given as the sum of the incident field Exvinc(x′,t) and the field reflected from the interface Exvref(x′,t) [11],(73)Exvexc(x′,z,t)=Exvinc(x′,t-T)+Exvref(x′,t-T).The time shiftT represents the time required for the wave to travel from the highest wire to the height z of the observed vth wire. Assigning the highest wire with index U, it can be written(74)T=zU-zc,zU=max⁡(z1,z2,…,z,…,zM).The field reflected from the interface for the case of normal incidence is given by(75)Exvref(x′,t-T)=r(θ=0)⋅Exvinc(x′,t-T-2zc),
where t-T-2z/c is the time needed for the wave to travel from observed vth wire to the interface.For the case of PEC ground plane, the space-time reflection coefficient (54) simply becomes(76)rvs(θvs′,t)=1.Thus, the set of (72) simplifies into (77)∑s=1M∫x0sxLsIs(x′,t-Rvs/c)4πRvsdx′-∑s=1M∫x0sxLsIs(x′,t-Rvs*/c)4πRvs*dx′=F0v(t-x-x0vc)+FLv(t-xLv-xc)+12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
and the field reflected from PEC ground is simply given by(78)Exvref(x′,t-T)=Exvinc(x′,t-T-2z/c).Given the dielectric constant of the medium and the known excitationExvinc(x′,t), a system of M coupled Hallen integral equations can be solved using time domain version of GB-IBEM and by taking into account appropriate boundary and initial conditions. Boundary conditions assume zero currents at the end of each wire, while initial conditions assume all the currents to be zero for t≤0.


### 3.2. The Numerical Solution

First, numerical procedure for single wire Hallen equation is outlined. Applying the weighted residual approach in the spatial domain and GB-IBEM procedure [8], the following local matrix system is obtained:
(79)[A]{I}i|t-R/c-[A*]{I}i|t-R*/c-{Â}|t-R*/c=[B]{E}|t-|x-x′|/c+[C]{∑n=0∞In}i|t-(R0/c)-(2nL/c)-(x/c)-[C*]{∑n=0∞In}i|t-(R0*/c)-(2nL/c)-(x/c)-{∑n=0∞Ĉn}|t-(R0*/c)-(2nL/c)-(x/c)-[B]{∑n=0∞En}|t-(x′/c)-(2nL/c)-(x/c)-[D]{∑n=0∞In}i|t-(RL/c)-((2n+1)L/c)-(x/c)+[D*]{∑n=0∞In}i|t-(RL*/c)-((2n+1)L/c)-(x/c)+[B]{∑n=0∞En}|t-((L-x′)/c)-((2n+1)L/c)-(x/c)+{∑n=0∞D̂n}|t-(RL*/c)-((2n+1)L/c)-(x/c)+[D]{∑n=0∞In}i|t-(RL/c)-(2nL/c)-(L-x/c)-[D*]{∑n=0∞In}i|t-(RL*/c)-(2nL/c)-(L-x/c)+{∑n=0∞D̂n}|t-(RL*/c)-(2nL/c)-(L-x/c)-[B]{∑n=0∞En}|t-((L-x′)/c)-(2nL/c)-(L-x/c)-[C]{∑n=0∞In}i|t-(R0/c)-((2n+1)L/c)-(L-x/c)+[C*]{∑n=0∞In}i|t-(R0*/c)-((2n+1)L/c)-(L-x/c)+[B]{∑n=0∞En}|t-(x′/c)-((2n+1)L/c)-(L-x/c)-{∑n=0∞Ĉn}|t-(R0*/c)-((2n+1)L/c)-(L-x/c).The space-dependent local matrices representing the interaction betweenith source and jth observation element are defined as follows:(80)[A]=∫Δlj∫Δli{f}j{f}iT14πRdx′dx,[B]=12Z0∫Δlj∫Δli{f}j{f}iTdx′dx,[C]=∫Δlj∫Δli{f}j{f}iT14πR0dx′dx,[D]=∫Δlj∫Δli{f}j{f}iT14πRLdx′dx,[A*]=∫Δlj∫Δli{f}j{f}iTr(θ)4πR*dx′dx,[C*]=∫Δlj∫Δli{f}j{f}iTr(θ)4πR0*dx′dx,[D*]=∫Δlj∫Δli{f}j{f}iTr(θ)4πRL*dx′dx,
where {f} stands for the shape functions, while additional time dependent vectors are given by
(81){Â}=∫0t-R*/c∫Δlj∫Δli{f}j{f}iTH1dx′dx{I(τ)}idτ{Ĉn}=∫0t-(R0*/c)-(2nL/c)-(x/c)×∫Δlj∫Δli{f}j{f}iTH2dx′dx{I(τ)}idτ{D̂n}=∫0t-(RL*/c)-((2n+1)L/c)-(x/c)×∫Δlj∫Δli{f}j{f}iTH3dx′dx{I(τ)}idτ,
where(82)H1=r′′(θ,t-(R*/c)-τ)4πR*,H2=r′′(θ,t-(R0*/c)-(2nL/c)-(x/c)-τ)4πR0*,H3=r′′(θ,t-(RL*/c)-((2n+1)L/c)-(x/c)-τ)4πR0*.
Assembling the local matrices and vectors into the global ones the following global matrix system is formed:(83)[A]{I}|t-R/c={g}|previous timeinstants+{ĝ}|previous timeinstants,
where
(84){g}=[A*]{I}|t-R*/c+[B]{E}|t-|x-x′|/c+[C]{∑n=0∞In}|t-(R0/c)-(2nL/c)-(x/c)-[C*]{∑n=0∞In}|t-(R0*/c)-(2nL/c)-(x/c)-[B]{∑n=0∞En}|t-(x′/c)-(2nL/c)-(x/c)-[D]{∑n=0∞In}|t-(RL/c)-((2n+1)L/c)-(x/c)+[D*]{∑n=0∞In}|t-(RL*/c)-((2n+1)L/c)-(x/c)+[B]{∑n=0∞En}|t-((L-x′)/c)-((2n+1)L/c)-(x/c)+[D]{∑n=0∞In}|t-(RL/c)-(2nL/c)-(L-x/c)-[D*]{∑n=0∞In}|t-(RL*/c)-(2nL/c)-(L-x/c)-[B]{∑n=0∞En}|t-(L-x′/c)-(2nL/c)-(L-x/c)-[C]{∑n=0∞In}|t-(R0/c)-((2n+1)L/c)-(L-x/c)+[C*]{∑n=0∞In}|t-(R0*/c)-((2n+1)L/c)-(L-x/c)+[B]{∑n=0∞En}|t-(x′/c)-((2n+1)L/c)-(L-x/c),(85){ĝ}={Â}|t-(R*/c)-{∑n=0∞Ĉn}|t-(R0*/c)-(2nL/c)-(x/c)+{∑n=0∞D̂n}|t-(RL*/c)-((2n+1)L/c)-(x/c)+{∑n=0∞D̂n}|t-(RL*/c)-(2nL/c)-(L-x′/c)-{∑n=0∞Ĉn}|t-(R0*/c)-((2n+1)L/c)-(L-x′/c).Applying the weighted residual approach in the time domain and using the Dirac impulses as weight functions, the time sampling is provided, and the following recurrent formula is obtained:(86)Ij|tk=∑i=1NsajiIj|tk-R/c-gj|previous timeinstants-ĝj|previous timeinstantsajj,
where Ij|tk is current for the jth space node at kth time instant, N is total number of space segments, while the overbar indicates the absence of diagonal members.It is worth noting that the numerical calculation of convolution integrals is rather tedious task leading to tremendously large computational time of the overall method. The main advantage of the method, on the other hand, is its unconditional stability.Time domain GB-BEM procedure for the set of Hallen equations is undertaken in a similar manner as in the case of a single wire.The solution of (72) and (77), respectively, is also carried out using the GB-IBEM technique.Applying the boundary element discretisation to (72) and (77), respectively, leads to a local system of linear equations for the vth observed wire:(87)∑s=1M[∫Δli∫Δlj14πRvs{f}j{f}iTdx′dx{Is}|t-Rvs/c-∫Δli∫Δljrvs(θ)4πRvs*{f}j{f}iTdx′dx{Is}|t-Rvs*/c]=12Z0∫Δli∫ΔljExvexc(x′,t-|x-x′|c){f}jdx′dx+∫ΔljF0(t-x-x0vc){f}jdx+∫ΔljFL(t-xLv-xc){f}jdx,
where i,j=1,2,…,Ndenotes the index of the elements located on the sth source wire and the vth observed wire, respectively, with N as the total number of space segments, while M is the actual number of wires.Finally, substituting (64)–(67) into (87), the following local matrix system is obtained:(88)∑s=1M[Avs]{Is}|t-Rvs/c-∑s=1M[Avs*]{Is}|t-Rvs*/c=[Bv]{Ev}|t-|x-x′|/c+∑s=1M[Cvs]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n/c)Lv)-(Rvs(0)/c)-∑s=1M[Cvs*]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n/c)Lv)-(Rvs*(0)/c)-[Dv]{∑n=0∞Evn}|t-((x-x0v)/c)-((2n/c)Lv)-(|x′-x0v|/c)-∑s=1M[Evs]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n+1/c)Lv)-Rvs(L)/c+∑s=1M[Evs*]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n+1/c)Lv)-Rvs*(L)/c+[Dv]{∑n=0∞Evn}|t-((x-x0v)/c)-((2n+1/c)Lv)-(|xLv-x′|/c)+∑s=1M[Evs]{∑n=0∞Isn}|t-((xLv-x)/c)-((2n/c)Lv)-(Rvs(L)/c)-∑s=1M[Evs*]{∑n=0∞Isn}|t-((xLv-x)/c)-((2n/c)Lv)-(Rvs*(L)/c)-[Dv]{∑n=0∞Evn}|t-((xLv-x)/c)-((2n/c)Lv)-(|xLv-x′|/c)-∑s=1M[Cvs]{∑n=0∞Isn}|t-((xLv-x)/c)-((2n+1/c)Lv)-(Rvs(0)/c)+∑s=1M[Cos*]{∑n=0∞Isn}|t-((xLo-x)/c)-((2n+1/c)Lo)-(Ros*(0)/c)+[Do]{∑n=0∞Eon}|t-((xLo-x)/c)-((2n+1/c)Lo)-(|x′-x0o|/c),
where {E} denotes excitation vector and space-dependent matrices are of the form(89)[Avs]=∫Δlj∫Δli14πRvs{f}j{f}iTdx′dx,[Avs*]=∫Δlj∫Δlirvs(θ)4πRvs*{f}j{f}iTdx′dx,[Bv]=12Z0∫Δlj∫Δli{f}j{f}iTdx′dx,[Cvs]=∫Δlj∫Δli14πRvs(0){f}j{f}iTdx′dx,[Cvs*]=∫Δlj∫Δlirvs(θ)4πRvs*(0){f}j{f}iTdx′dx,[Dv]=12Z0∫Δlj∫Δli{f}j{f}iTdx′dx,[Evs]=∫Δlj∫Δli14πRvs(L){f}j{f}iTdx′dx,[Evs*]=∫Δlj∫Δlirvs(θ)4πRvs*(L){f}j{f}iTdx′dx.Relations containing summation fromn=0 to infinity pertain to the reflections of transient current from the wire ends. Note as the observed time interval is always finite, only a finite number of reflections occurs within a given observation interval. A shorter observed interval requires smaller number of summands and vice versa.According to GB-IBEM, a global matrix system is assembled from the local matrix systems for all wiresv=1,2,…,M. Finally, the resulting global matrix system can be written as follows:(90)[A]{I}|t-Rvs/c-[A*]{I}|t-Rvs*/c={g}.The time-domain solution on theith boundary element is given by(91)Ii(t)=∑k=1NtIikTk(t′),
where Iik are unknown coefficients, Tk are the linear time-domain shape functions, and Nt is the total number of time samples.Applying the weighted residual approach to (90) leads to the expression(92)∫tktk+Δt([A]{I}|t-Rvs/c-[A*]{I}|t-Rvs*/c-{g}θk)dt=0,k=1,2,…,Nt,
where θk  denotes the set of time-domain weights.Using the set of Dirac impulses for the test functions, time sampling is ensured and (92) becomes(93)[A]{I}|tk-Rvs/c-[A*]{I}|tk-Rvs*/c=g|all previous discrete  instants.If the space-time discretization is performed by satisfying the Courant condition,Δx≥cΔt, the transient current for a jth space node and kth time node can be obtained from a recurrence formula. Separating the terms relating to the current induced at the instant tk in (93) yields(94)AjjIj|tk+[A̅]{I}|tk-Rvs/c-[A*]{I}|tk-Rvs*/c={g}|allpreviousdiscreteinstants,
where overbar indicates the absence of diagonal terms.The first term in (94) pertains to the current at the jth space node and kth time node, that is, the present instant. Other terms are related to all previous instants. Finally, the recurrence formula for the transient current at jth space node and kth time node is obtained in the forms(95)Ij|tk=-∑i=1N(Aji¯Ii|tk-Rvs/c+Aji*Ii|tk-Rvs*/c)+gj|allpreviousdiscreteinstantsAjj,whereN is total number of space elements, k=1,2,…,Nt is the index of the kth time instant.


### 3.3. The Transmission Line Model

The time-domain field-to-transmission line coupling equations can be written in the matrix form [11] (96)∂∂x[V(x,t)]+[R][I(x,t)]+[L]⋅∂∂t[I(x,t)]=[EF(x,t)]-[z′(t)]*[I(x,t)],(97)∂∂x[I(x,t)]+[G]⋅[V(x,t)]+[C]⋅∂∂t[V(x,t)]=[HF(x,t)],
where “*” stands for the convolution product, [z′(t)] is the transient inverse Fourier transform of the ground, conductors matrix [Ẑw(s)+Ẑg(s)], and s=jω is the Laplace variable.[EF(x,t)]and [HF(x,t)] are the excitation terms, given by [11](98)[EF(x,t)]=-∂∂x[VT(x,t)]+[EL(x,t)],[HF(x,t)]=-[G][VT(x,t)]-[C]∂∂t[VT(x,t)],
where [VT(x,t)] is the transverse voltage derived from the transverse incident field excitation [11], and [EL(x,t)] represents the longitudinal electrical field excitation.The classical per-unit-length inductances of[L] matrix have been used:(99)Lii=μ02πln⁡(2hiai),Lij=μ02πln⁡(Dijdij),
where hi,ai are, respectively, the height and the radius of the ith conductor.  Dij is the distance between the ith conductor and the image of jth conductor and dij is the distance between the ith conductor and the jth conductor.The capacitance matrix⌈C⌉ is derived from the inductance matrix ⌈L⌉ by [C]=ε0μ0[L]-1.In case of PEC ground, the transient ground matrix[z′(t)] is equal to [ZW(t)] while if a dielectric half-space is of interest, the ground matrix [Ẑg(s)] is given by Carson integral [11].It is worth emphasizing that, for both cases, the conductivity matrix is neglected.


### 3.4. The Solution of Transmission Lineequations via FDTD

Discretizing each conductor of the multiconductor transmission line (MTL) intoNx  sections each of length Δx and discretizing the entire time interval into increments of duration Δt, the FDTD method is applied to (88) and (97).The solutions of (91) and (99) for lossy MTL by FDTD are given by (100)[Z+][Ikn+1/2]=-[Z-][Ikn-1/2]+[Vk-1/2n]+[VT,k-1/2n]-[Vk+1/2n]-[VT,k+1/2n]+([EL,kn]-[Skn])Δxfor0≤k≤Nx-1,n≥1,
where (101)[Z±]=[R]2±[L]+[Z0(1)]ΔtΔt,(102)[Y+]([Vk+1/2n+1]+[VT,k+1/2n+1])=-[Y-]([Vk+1/2n]+[VT,k+1/2n])-[Ik+1n+1/2]+[Ikn+1/2]for0≤k≤Nx-1,n≥1,(103)[Y±]=[G]2±[C]Δt.The expressions for[V],[I],[VT],[EL],[S] are given for time n≥1 and space k≥0:(104)[Vk+1/2n]=[V((k+1/2)Δx,nΔt)],0≤k≤Nx-1.Voltage and adjacent current are interlaced in time and space, respectively, byΔt/2 and Δx/2, then(105)[Ikn-1/2]=[I(kΔx,(n-12)Δt)],[VT,k+1/2n]=[VT((k+12)Δx,nΔt)],[EL,kn]=[EL(kΔx,nΔt)],[Skn]=[S(kΔx,nΔt)],0≤k≤Nx.The convolution product[Skn] appearing in (100) can be written as follows:(106)[Skn]=[Z0(2)][Ikn-1/2]+∑l=1n-1[Z0(l+1)-Z0(l)][Ikn-l+1/2]-[Z0(n)][Ik1/2],
where (107)[Z0(l)]=[Z0w(l)]+[Z0g(l)]=∫l-1l[Zw+Zg](uΔt)du,1≤l≤n,t=nΔt, Nx is the number of space steps, and n is the number of time steps.The corresponding components of theM×1 vectors [ET,kn], [EL,kn] are as follows:(108)[VT,kn]i=yiEyinc(xi=kΔx,yi,zi,nΔt)+ziEzinc(xi=kΔx,yi,zi,nΔt),(109)[EL,kn]i=Exinc(xi=kΔx,yi,zi,nΔt)-Ezinc(xi=kΔx,0,0,nΔt).Equations (107) and (108) are valid for 1≤i≤M, where M is the number of conductors, and  yi,zi are the positions of ith conductor. Exinc,Eyinc,Ezinc are the components of the incident electromagnetic field evaluated in the absence of conductors. [R],[L],[G],and[C] are, respectively, the per-unit-length resistance, inductance, conductance, and capacitance matrices of dimension  M×M.[Zg(t)] is the transient ground resistance and is equal to the inverse Fourier of [Zg(s)/s]:(110)Zgij(t)=F-1(Zgij(s)s),
where the ground impedance in frequency domain is given by Carson formula:(111)Zgij(s)s=μ0π∫0∞e-(hi+hj)λλ2+γg2+λ⋅cos⁡(dijλ)dλ.[Zw(t)] is a diagonal matrix and corresponds to the transient conductor resistance, where each element is given by(112)Zwi(t)=1πσwiai2∑m=1∞e-xm2(t/τwi),whereτwi=μ0σwiai2,
and hi,hj,dij are the corresponding heights of the two conductors i,j and the distance between the two conductors in the horizontal plane.Termsxm stand for the zeros of J1, Bessel function of first kind.Finally,γg is the propagation constant defined as(113)γg2=sμ0(σg+sε0εrg),
where σg and εrg are, respectively, ground conductivity and permittivity.For the cases considered in this paper, all conductors are in open circuit at both ends, so the currents vanishes at near and far ends.By using the boundary formulation [11], voltages at both ends are simply expressed as follows:(114)[V0n]=[V(1/2)n],[VNxn]=[VNx-(1/2)n],(115)[I0n]=[0],[INxn]=[0].


## 3.1. The Antenna Theory Model

Generally, a direct time-domain analysis of thin wire in the presence of a lossy half-space can be carried out via the appropriate space-time integral equations of either Pocklington or Hallen type [1, 8]. When applied to the solution of the Hallen integral equation the Galerkin-Bubnov Indirect Boundary Element Method (GB-IBEM) [8] results in relatively complex procedures compared to various procedures for the solution of Pocklington equations, but, at the same time, it is proven to be highly efficient, accurate, and unconditionally stable  [8, 18, 19]. On the other hand, the implementation of GB-IBEM to the solution of the Pocklington-type equation is relatively simple, but suffers from serious numerical instabilities. The origin of these instabilities is the discretization of space-time differential operator [8]. The GB-IBEM solution of the Pocklington equation in free space for certain values of time domain integration parameters has been presented elsewhere, for example, in [19], while the Hallen integral equation solution by means of GB-IBEM has been obtained for thin wire structures in the presence of a dielectric half-space, for example, in [11]. In both cases, the influence of imperfect ground has been taken into account via the corresponding reflection coefficient. The numerical solution was mostly limited to scenarios in which the finite conductivity of the ground could be ignored. This approximation involves cases where the wires are sufficiently far from the two-media interface or where the ground conductivity is appreciably low or very high, that is, where the approximation of pure dielectric medium or perfect ground is applied. Through these approximations the time-dependent part of the reflection coefficient function vanishes, and the resulting matrix equation simplifies significantly.However, for the cases where these approximations are not valid, modifications to the original methods are required in order to include the ground conductivity [8]. Namely, the related reflection coefficient is space-time dependent, and the resulting convolution integrals have to be included in the matrix system and numerically computed. This leads to a significant increase in the overall computational cost of the method and consequently requires several modifications.This section deals with the transient analysis of multiple horizontal wires above a lossy ground using the Hallen integral equation approach.The set of space-time Hallen's integral equations can be derived as an extension of the single wire case. First, a single wire insulated in free space is considered.Thin wire antenna or scatterer of lengthL and radius a, oriented along the x-axis, is considered. The wire is assumed to be perfectly conducting and excited by a plane wave electric field. For the sake of simplicity, the analysis is restricted to the case of a normally incident electric field.The tangential component of the total field vanishes on the PEC wire surface, that is,(42)Exinc+Exsct=0,
where Exinc is the incident and Exsct scattered field on the metallic wire surface. Starting from Maxwell equations and obeying the Lorentz gauge one obtains a time domain version of  (2):(43)(∂2A⃗∂t2-1με∇(∇A⃗))|tan=∂E⃗inc∂t|tan,
where A⃗ is the space-time-dependent vector potential.According to the thin wire approximation, only the axial component of the vector potential exists, that is, (43) becomes  (44)∂2Ax∂x2-1c2∂2Ax∂t2=-1c2∂Exinc∂t,
where c denotes the velocity of light.The corresponding solution of (44) can be expressed in terms of a sum of the general solution of the homogeneous equation and the particular solution of the inhomogeneous equation:(45)Ax(x,t)=Axh(x,t)+Axp(x,t).
The solution of the homogeneous wave equation is given as a superposition of incident and reflected wave [8]:(46)Axh(x,t)=F1(t-xc)+F2(t+xc),
while the particular solution is given by the integral [8]:(47)Axp(x,t)=12Z0∫0LExinc(x′,t-|x-x′|c)dx′,
where L denotes the total antenna length.On the other hand, the magnetic vector potential on the PEC wire surface is given by the particular integral:(48)Ax(x,t)=μ4π∫SI(x′,t-R/c)Rdx′.Combining (45)–(48) yields the space-time Hallen equation:(49)∫0LI(x′,t-R/c)4πRdx′=F0(t-xc)+FL(t-L-xc)+12Z0∫0LExinc(x′,t-|x-x′|c)dx′,
where I(x′) is the equivalent axial current to be determined, Exinc is the known tangential incident field, R=[(x-x′)2+a2]1/2 is the distance from the source point (the equivalent current in the antenna axis) to the observation point, and Z0 is the wave impedance of a free space.The unknown functionsF0(t) and FL(t) account for the multiple reflections of the current at the free ends of the wire.A direct time formulation for a straight thin wire above a dissipative half-space can be obtained as the extension of the free-space Hallen equation (49).The free space Hallen equation (49) is first transferred into the Laplace frequency domain:(50)∫0LI(x′,s)e-sR/c4πRdx′=F0(s)e-sx/c+FL(s)e-s((L-x)/c)+12Z0∫0LExinc(x′,s)e-s|x-x′|/cdx′,
where s=jω is the Laplace variable.According to the theory of images the free space integral equation (50) is extended by an additional term multiplying the Green function of the image source by space-frequency-dependent reflection coefficient RTM(θ′,s) for TM polarization. The integral equation in the frequency domain is given by(51)∫0LI(x,s)esR/c4πRdx′-∫0LRTM(θ,s)I(x,s)esR*/c4πR*dx′=F0(s)e-sx/c+FL(s)e-s((L-x)/c)+12Z0∫0LExexc(x′,s)e-s(|x-x′|/c)dx′,
where R*=(x-x′)2+4h2 and  RTM(θ′,s) is determined by the expression [1](52)RTM(θ′,s)=εr(1+σ/εs)cos⁡⁡θ′-εr(1+σ/εs)-sin⁡2⁡θ′εr(1+σ/εs)cos⁡⁡θ′+εr(1+σ/εs)-sin⁡2⁡θ′,
where σ and  ε are the lossy medium conductivity and permittivity, respectively, andθ′=arctg(|x-x′|/2h).The reflection coefficient (RC) approach is a satisfactory approximation in half-space calculations, as long as the field is calculated far away from the source, and the imperfect ground, respectively, to ensureθ′<π/2 [8].Performing the convolution, the time domain counterpart of (52) is obtained in the form(53)∫0LI(x′,t-R/c)4πRdx′-∫-∞t∫0Lr(θ,τ)I(x′,t-R*/c-τ)4πR*dx′dτ=12Z0∫0LExexc(x′,t-|x-x′|c)dx′+F0(t-xc)+FL(t-L-xc),
where r(θ,t) is the space-time reflection coefficient which, for convenience, can be written in the form [18] (54)r(θ,τ)=r′(θ,τ)+r′′(θ,τ),
where(55)r′(θ,t)=Kδ(t),r′′(θ,t)=4β1-β2e-αtt∑n=1∞(-1)n+1nKnIn(αt),τ=σε0εr,β=εr-sin2θεrcos⁡θ,γ=τ1-(sin2θ/εr),θ=arctg|x-x′|2h,K=1-β1+β,α=τ2.Note thatIn is the modified Bessel function of the first order and nth degree.For the case of normal incidence, which is considered for the sake of simplicity, the excitation term is given by(56)Exexc(t)=Exinc(t)-Exref(t*),
where t*=t-R*/c.The transient ground-reflected field is obtained as the convolution of the incident field and the space-time reflection coefficient for the angle of incidenceθ=0 (in accordance with the parallel incidence of the electric field), as is proposed in [20](57)Exref(t)=∫-∞tExinc(t-τ)r(θ=0,τ)dτ
and the integral equation (53) becomes
(58)∫0LI(x′,t-R/c)4πRdx′-∫-∞t∫0Lr(θ,τ)I(x′,t-R*/c-τ)4πR*dx′dτ=F0(t-xc)+FL(t-L-xc)+12Z0∫0LExinc(x′,t-|x-x′|c)dx′-12Z0∫-∞t∫0LExinc(x′,t-|x-x′|c-τ)×r(θ=0,τ)dx′dτ.The unknown time functionsF0(t),FL(t),F0(t-(L/c)), and FL(t-(L/c))  can be obtained in the same manner, as in the case of free space, in terms of as auxilliary functions K0(t) and KL(t) [8]:(59)F0(t)=∑n=0∞K0(t-2nLc)-∑n=0∞KL(t-(2n+1)Lc),FL(t)=∑n=0∞KL(t-2nLc)-∑n=0∞K0(t-(2n+1)Lc),
where(60)KL(t)=∫0LI(x′,t-R0/c)4πR0dx′-∫-∞t∫0Lr(θ′,τ)I(x′,t-R0*/c-τ)4πR0*dx′dτ-12Z0∫0LExinc(x′,t-x′c)dx′,KL(t)=∫0LI(x′,t-RL/c)4πRLdx′-∫-∞t∫0Lr(θ′,τ)I(x′,t-RL*/c-τ)4πRL*dx′dτ-12Z0∫0LExexc(x′,t-L-x′c)dx′,
while R0 and RL are the distances from the wire ends to the source point and R0*,RL* are the distances from the image wire ends to the image source point.If the case of a perfect (ideal) dielectric half-space is considered, the Hallen equation (53) simplifies into(61)∫0LI(x′,t-R/c)4πRdx′-∫0Lr(θ)I(x′,t-R*/c)4πR*dx′=F0(t-xc)+FL(t-L-xc)+12Z0∫0LExexc(x′,t-|x-x′|c)dx′.
Space-time integral equation (53) or (61), respectively, can be solved assuming the zero current at the free ends of the wire and with the initial conditions requiring the wire not to be excited before the certain instant t=t0  [8].The transient behavior ofM straight horizontal thin wires located at different heights above an infinite ground plane is determined by a set of the coupled space-time integral equations of the Hallen type [11]:(62)∑s=1M∫x0sxLsIs(x′,t-Rvs/c)4πRvsdx′-∑s=1M∫-∞t∫x0sxLsrvs(θ,τ)Is(x′,t-Rvs*/c-τ)4πRvs*dx′dτ=F0v(t-x-x0vc)+FLv(t-xLv-xc)+12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
where v,s=1,2,…,M denote the index of the observed and source wire, respectively. Furthermore, Ls and Lv are the lengths of the sth and vth wire, and x′,x are the x-coordinates of the source and observation points on respective wires. The distances between observation point (x,y,z) on the wire v and source point (x′,y′,z′) on the wire s are given by(63)Rvs={(x-x′)2+(y-y′)2+(z-z′)2:v≠s(x-x′)2+a2:v=s,Rvs*=(x-x′)2+(y-y′)2+(z+z′)2,
where asterisk is related to source points are located on the image wire.Unknown time signalsF0v(t) and FLv(t) account for the multiple reflections of transient currents at the wire open ends and can be written in the form(64)F0v(t)=∑n=0∞K0v(t-2nLvc)-∑n=0∞KLv(t-(2n+1)Lvc),(65)FLv(t)=∑n=0∞KLv(t-2nLvc)-∑n=0∞K0v(t-(2n+1)Lvc),
while the auxiliary functions K are defined, as follows:(66)K0v(t)=∑s=1M∫x0sxLsIs(x′,t-Rvs(0)/c)4πRvs(0)dx′-∑s=1M∫-∞t∫x0sxLsrvs(θ,τ)Is(x′,t-Rvs*(0)/c-τ)4πRvs*(0)dx′dτ-12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,(67)KLv(t)=∑s=1M∫x0sxLsIs(x′,t-Rvs(L)/c)4πRvs(L)dx′-∑s=1M∫-∞t∫x0sxLsrvs(θ,τ)Is(x′,t-Rvs*(L)/c-τ)4πRvs*(L)dx′dτ-12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
where Rvs(0) and Rvs(L) are distances from considered source point on each wire s to a corresponding observation point at the ends of the wire v:(68)Rvs(0)=Rvs|x=x0v,Rvs(L)=Rvs|x=xLv,
while Rvs*(0) and Rvs*(L) are distances between the source point at the image of the wire sand observation point located at the ends of the wire v:(69)Rvs*(0)=Rvs*|x=x0v,Rvs*(L)=Rvs*|x=xLv.The space-time reflection coefficientrvs(θ,t) accounts for the influence of the interface and is given by [11] (70)rvs(θvs′,t)=Aδ(t),
where(71)A=1-β1+β,β=εr-sin2θ′εrcos⁡θ′,θvs′=Arctg(x′-x)2+(y′-y)2z′+z.The angleθvs′is the angle between the source point on the image of the wire s(x′,y′,-z′) and the observation point (x,y,z) on wire v.Substituting (70) into (62) yields(72)∑s=1M∫x0sxLsIs(x′,t-Rvs/c)4πRvsdx′-∑s=1M∫x0sxLsrvs(θ)Is(x′,t-Rvs*/c)4πRvs*dx′=F0v(t-x-x0vc)+FLv(t-xLv-xc)+12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
where Exvexc is the space-time-dependent tangential electric field on the vth wire.For the case of normal incidence the total excitation fieldExvexc(x′,t) is given as the sum of the incident field Exvinc(x′,t) and the field reflected from the interface Exvref(x′,t) [11],(73)Exvexc(x′,z,t)=Exvinc(x′,t-T)+Exvref(x′,t-T).The time shiftT represents the time required for the wave to travel from the highest wire to the height z of the observed vth wire. Assigning the highest wire with index U, it can be written(74)T=zU-zc,zU=max⁡(z1,z2,…,z,…,zM).The field reflected from the interface for the case of normal incidence is given by(75)Exvref(x′,t-T)=r(θ=0)⋅Exvinc(x′,t-T-2zc),
where t-T-2z/c is the time needed for the wave to travel from observed vth wire to the interface.For the case of PEC ground plane, the space-time reflection coefficient (54) simply becomes(76)rvs(θvs′,t)=1.Thus, the set of (72) simplifies into (77)∑s=1M∫x0sxLsIs(x′,t-Rvs/c)4πRvsdx′-∑s=1M∫x0sxLsIs(x′,t-Rvs*/c)4πRvs*dx′=F0v(t-x-x0vc)+FLv(t-xLv-xc)+12Z0∫x0vxLvExvexc(x′,t-|x-x′|c)dx′,
and the field reflected from PEC ground is simply given by(78)Exvref(x′,t-T)=Exvinc(x′,t-T-2z/c).Given the dielectric constant of the medium and the known excitationExvinc(x′,t), a system of M coupled Hallen integral equations can be solved using time domain version of GB-IBEM and by taking into account appropriate boundary and initial conditions. Boundary conditions assume zero currents at the end of each wire, while initial conditions assume all the currents to be zero for t≤0.


## 3.2. The Numerical Solution

First, numerical procedure for single wire Hallen equation is outlined. Applying the weighted residual approach in the spatial domain and GB-IBEM procedure [8], the following local matrix system is obtained:
(79)[A]{I}i|t-R/c-[A*]{I}i|t-R*/c-{Â}|t-R*/c=[B]{E}|t-|x-x′|/c+[C]{∑n=0∞In}i|t-(R0/c)-(2nL/c)-(x/c)-[C*]{∑n=0∞In}i|t-(R0*/c)-(2nL/c)-(x/c)-{∑n=0∞Ĉn}|t-(R0*/c)-(2nL/c)-(x/c)-[B]{∑n=0∞En}|t-(x′/c)-(2nL/c)-(x/c)-[D]{∑n=0∞In}i|t-(RL/c)-((2n+1)L/c)-(x/c)+[D*]{∑n=0∞In}i|t-(RL*/c)-((2n+1)L/c)-(x/c)+[B]{∑n=0∞En}|t-((L-x′)/c)-((2n+1)L/c)-(x/c)+{∑n=0∞D̂n}|t-(RL*/c)-((2n+1)L/c)-(x/c)+[D]{∑n=0∞In}i|t-(RL/c)-(2nL/c)-(L-x/c)-[D*]{∑n=0∞In}i|t-(RL*/c)-(2nL/c)-(L-x/c)+{∑n=0∞D̂n}|t-(RL*/c)-(2nL/c)-(L-x/c)-[B]{∑n=0∞En}|t-((L-x′)/c)-(2nL/c)-(L-x/c)-[C]{∑n=0∞In}i|t-(R0/c)-((2n+1)L/c)-(L-x/c)+[C*]{∑n=0∞In}i|t-(R0*/c)-((2n+1)L/c)-(L-x/c)+[B]{∑n=0∞En}|t-(x′/c)-((2n+1)L/c)-(L-x/c)-{∑n=0∞Ĉn}|t-(R0*/c)-((2n+1)L/c)-(L-x/c).The space-dependent local matrices representing the interaction betweenith source and jth observation element are defined as follows:(80)[A]=∫Δlj∫Δli{f}j{f}iT14πRdx′dx,[B]=12Z0∫Δlj∫Δli{f}j{f}iTdx′dx,[C]=∫Δlj∫Δli{f}j{f}iT14πR0dx′dx,[D]=∫Δlj∫Δli{f}j{f}iT14πRLdx′dx,[A*]=∫Δlj∫Δli{f}j{f}iTr(θ)4πR*dx′dx,[C*]=∫Δlj∫Δli{f}j{f}iTr(θ)4πR0*dx′dx,[D*]=∫Δlj∫Δli{f}j{f}iTr(θ)4πRL*dx′dx,
where {f} stands for the shape functions, while additional time dependent vectors are given by
(81){Â}=∫0t-R*/c∫Δlj∫Δli{f}j{f}iTH1dx′dx{I(τ)}idτ{Ĉn}=∫0t-(R0*/c)-(2nL/c)-(x/c)×∫Δlj∫Δli{f}j{f}iTH2dx′dx{I(τ)}idτ{D̂n}=∫0t-(RL*/c)-((2n+1)L/c)-(x/c)×∫Δlj∫Δli{f}j{f}iTH3dx′dx{I(τ)}idτ,
where(82)H1=r′′(θ,t-(R*/c)-τ)4πR*,H2=r′′(θ,t-(R0*/c)-(2nL/c)-(x/c)-τ)4πR0*,H3=r′′(θ,t-(RL*/c)-((2n+1)L/c)-(x/c)-τ)4πR0*.
Assembling the local matrices and vectors into the global ones the following global matrix system is formed:(83)[A]{I}|t-R/c={g}|previous timeinstants+{ĝ}|previous timeinstants,
where
(84){g}=[A*]{I}|t-R*/c+[B]{E}|t-|x-x′|/c+[C]{∑n=0∞In}|t-(R0/c)-(2nL/c)-(x/c)-[C*]{∑n=0∞In}|t-(R0*/c)-(2nL/c)-(x/c)-[B]{∑n=0∞En}|t-(x′/c)-(2nL/c)-(x/c)-[D]{∑n=0∞In}|t-(RL/c)-((2n+1)L/c)-(x/c)+[D*]{∑n=0∞In}|t-(RL*/c)-((2n+1)L/c)-(x/c)+[B]{∑n=0∞En}|t-((L-x′)/c)-((2n+1)L/c)-(x/c)+[D]{∑n=0∞In}|t-(RL/c)-(2nL/c)-(L-x/c)-[D*]{∑n=0∞In}|t-(RL*/c)-(2nL/c)-(L-x/c)-[B]{∑n=0∞En}|t-(L-x′/c)-(2nL/c)-(L-x/c)-[C]{∑n=0∞In}|t-(R0/c)-((2n+1)L/c)-(L-x/c)+[C*]{∑n=0∞In}|t-(R0*/c)-((2n+1)L/c)-(L-x/c)+[B]{∑n=0∞En}|t-(x′/c)-((2n+1)L/c)-(L-x/c),(85){ĝ}={Â}|t-(R*/c)-{∑n=0∞Ĉn}|t-(R0*/c)-(2nL/c)-(x/c)+{∑n=0∞D̂n}|t-(RL*/c)-((2n+1)L/c)-(x/c)+{∑n=0∞D̂n}|t-(RL*/c)-(2nL/c)-(L-x′/c)-{∑n=0∞Ĉn}|t-(R0*/c)-((2n+1)L/c)-(L-x′/c).Applying the weighted residual approach in the time domain and using the Dirac impulses as weight functions, the time sampling is provided, and the following recurrent formula is obtained:(86)Ij|tk=∑i=1NsajiIj|tk-R/c-gj|previous timeinstants-ĝj|previous timeinstantsajj,
where Ij|tk is current for the jth space node at kth time instant, N is total number of space segments, while the overbar indicates the absence of diagonal members.It is worth noting that the numerical calculation of convolution integrals is rather tedious task leading to tremendously large computational time of the overall method. The main advantage of the method, on the other hand, is its unconditional stability.Time domain GB-BEM procedure for the set of Hallen equations is undertaken in a similar manner as in the case of a single wire.The solution of (72) and (77), respectively, is also carried out using the GB-IBEM technique.Applying the boundary element discretisation to (72) and (77), respectively, leads to a local system of linear equations for the vth observed wire:(87)∑s=1M[∫Δli∫Δlj14πRvs{f}j{f}iTdx′dx{Is}|t-Rvs/c-∫Δli∫Δljrvs(θ)4πRvs*{f}j{f}iTdx′dx{Is}|t-Rvs*/c]=12Z0∫Δli∫ΔljExvexc(x′,t-|x-x′|c){f}jdx′dx+∫ΔljF0(t-x-x0vc){f}jdx+∫ΔljFL(t-xLv-xc){f}jdx,
where i,j=1,2,…,Ndenotes the index of the elements located on the sth source wire and the vth observed wire, respectively, with N as the total number of space segments, while M is the actual number of wires.Finally, substituting (64)–(67) into (87), the following local matrix system is obtained:(88)∑s=1M[Avs]{Is}|t-Rvs/c-∑s=1M[Avs*]{Is}|t-Rvs*/c=[Bv]{Ev}|t-|x-x′|/c+∑s=1M[Cvs]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n/c)Lv)-(Rvs(0)/c)-∑s=1M[Cvs*]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n/c)Lv)-(Rvs*(0)/c)-[Dv]{∑n=0∞Evn}|t-((x-x0v)/c)-((2n/c)Lv)-(|x′-x0v|/c)-∑s=1M[Evs]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n+1/c)Lv)-Rvs(L)/c+∑s=1M[Evs*]{∑n=0∞Isn}|t-((x-x0v)/c)-((2n+1/c)Lv)-Rvs*(L)/c+[Dv]{∑n=0∞Evn}|t-((x-x0v)/c)-((2n+1/c)Lv)-(|xLv-x′|/c)+∑s=1M[Evs]{∑n=0∞Isn}|t-((xLv-x)/c)-((2n/c)Lv)-(Rvs(L)/c)-∑s=1M[Evs*]{∑n=0∞Isn}|t-((xLv-x)/c)-((2n/c)Lv)-(Rvs*(L)/c)-[Dv]{∑n=0∞Evn}|t-((xLv-x)/c)-((2n/c)Lv)-(|xLv-x′|/c)-∑s=1M[Cvs]{∑n=0∞Isn}|t-((xLv-x)/c)-((2n+1/c)Lv)-(Rvs(0)/c)+∑s=1M[Cos*]{∑n=0∞Isn}|t-((xLo-x)/c)-((2n+1/c)Lo)-(Ros*(0)/c)+[Do]{∑n=0∞Eon}|t-((xLo-x)/c)-((2n+1/c)Lo)-(|x′-x0o|/c),
where {E} denotes excitation vector and space-dependent matrices are of the form(89)[Avs]=∫Δlj∫Δli14πRvs{f}j{f}iTdx′dx,[Avs*]=∫Δlj∫Δlirvs(θ)4πRvs*{f}j{f}iTdx′dx,[Bv]=12Z0∫Δlj∫Δli{f}j{f}iTdx′dx,[Cvs]=∫Δlj∫Δli14πRvs(0){f}j{f}iTdx′dx,[Cvs*]=∫Δlj∫Δlirvs(θ)4πRvs*(0){f}j{f}iTdx′dx,[Dv]=12Z0∫Δlj∫Δli{f}j{f}iTdx′dx,[Evs]=∫Δlj∫Δli14πRvs(L){f}j{f}iTdx′dx,[Evs*]=∫Δlj∫Δlirvs(θ)4πRvs*(L){f}j{f}iTdx′dx.Relations containing summation fromn=0 to infinity pertain to the reflections of transient current from the wire ends. Note as the observed time interval is always finite, only a finite number of reflections occurs within a given observation interval. A shorter observed interval requires smaller number of summands and vice versa.According to GB-IBEM, a global matrix system is assembled from the local matrix systems for all wiresv=1,2,…,M. Finally, the resulting global matrix system can be written as follows:(90)[A]{I}|t-Rvs/c-[A*]{I}|t-Rvs*/c={g}.The time-domain solution on theith boundary element is given by(91)Ii(t)=∑k=1NtIikTk(t′),
where Iik are unknown coefficients, Tk are the linear time-domain shape functions, and Nt is the total number of time samples.Applying the weighted residual approach to (90) leads to the expression(92)∫tktk+Δt([A]{I}|t-Rvs/c-[A*]{I}|t-Rvs*/c-{g}θk)dt=0,k=1,2,…,Nt,
where θk  denotes the set of time-domain weights.Using the set of Dirac impulses for the test functions, time sampling is ensured and (92) becomes(93)[A]{I}|tk-Rvs/c-[A*]{I}|tk-Rvs*/c=g|all previous discrete  instants.If the space-time discretization is performed by satisfying the Courant condition,Δx≥cΔt, the transient current for a jth space node and kth time node can be obtained from a recurrence formula. Separating the terms relating to the current induced at the instant tk in (93) yields(94)AjjIj|tk+[A̅]{I}|tk-Rvs/c-[A*]{I}|tk-Rvs*/c={g}|allpreviousdiscreteinstants,
where overbar indicates the absence of diagonal terms.The first term in (94) pertains to the current at the jth space node and kth time node, that is, the present instant. Other terms are related to all previous instants. Finally, the recurrence formula for the transient current at jth space node and kth time node is obtained in the forms(95)Ij|tk=-∑i=1N(Aji¯Ii|tk-Rvs/c+Aji*Ii|tk-Rvs*/c)+gj|allpreviousdiscreteinstantsAjj,whereN is total number of space elements, k=1,2,…,Nt is the index of the kth time instant.


## 3.3. The Transmission Line Model

The time-domain field-to-transmission line coupling equations can be written in the matrix form [11] (96)∂∂x[V(x,t)]+[R][I(x,t)]+[L]⋅∂∂t[I(x,t)]=[EF(x,t)]-[z′(t)]*[I(x,t)],(97)∂∂x[I(x,t)]+[G]⋅[V(x,t)]+[C]⋅∂∂t[V(x,t)]=[HF(x,t)],
where “*” stands for the convolution product, [z′(t)] is the transient inverse Fourier transform of the ground, conductors matrix [Ẑw(s)+Ẑg(s)], and s=jω is the Laplace variable.[EF(x,t)]and [HF(x,t)] are the excitation terms, given by [11](98)[EF(x,t)]=-∂∂x[VT(x,t)]+[EL(x,t)],[HF(x,t)]=-[G][VT(x,t)]-[C]∂∂t[VT(x,t)],
where [VT(x,t)] is the transverse voltage derived from the transverse incident field excitation [11], and [EL(x,t)] represents the longitudinal electrical field excitation.The classical per-unit-length inductances of[L] matrix have been used:(99)Lii=μ02πln⁡(2hiai),Lij=μ02πln⁡(Dijdij),
where hi,ai are, respectively, the height and the radius of the ith conductor.  Dij is the distance between the ith conductor and the image of jth conductor and dij is the distance between the ith conductor and the jth conductor.The capacitance matrix⌈C⌉ is derived from the inductance matrix ⌈L⌉ by [C]=ε0μ0[L]-1.In case of PEC ground, the transient ground matrix[z′(t)] is equal to [ZW(t)] while if a dielectric half-space is of interest, the ground matrix [Ẑg(s)] is given by Carson integral [11].It is worth emphasizing that, for both cases, the conductivity matrix is neglected.


## 3.4. The Solution of Transmission Lineequations via FDTD

Discretizing each conductor of the multiconductor transmission line (MTL) intoNx  sections each of length Δx and discretizing the entire time interval into increments of duration Δt, the FDTD method is applied to (88) and (97).The solutions of (91) and (99) for lossy MTL by FDTD are given by (100)[Z+][Ikn+1/2]=-[Z-][Ikn-1/2]+[Vk-1/2n]+[VT,k-1/2n]-[Vk+1/2n]-[VT,k+1/2n]+([EL,kn]-[Skn])Δxfor0≤k≤Nx-1,n≥1,
where (101)[Z±]=[R]2±[L]+[Z0(1)]ΔtΔt,(102)[Y+]([Vk+1/2n+1]+[VT,k+1/2n+1])=-[Y-]([Vk+1/2n]+[VT,k+1/2n])-[Ik+1n+1/2]+[Ikn+1/2]for0≤k≤Nx-1,n≥1,(103)[Y±]=[G]2±[C]Δt.The expressions for[V],[I],[VT],[EL],[S] are given for time n≥1 and space k≥0:(104)[Vk+1/2n]=[V((k+1/2)Δx,nΔt)],0≤k≤Nx-1.Voltage and adjacent current are interlaced in time and space, respectively, byΔt/2 and Δx/2, then(105)[Ikn-1/2]=[I(kΔx,(n-12)Δt)],[VT,k+1/2n]=[VT((k+12)Δx,nΔt)],[EL,kn]=[EL(kΔx,nΔt)],[Skn]=[S(kΔx,nΔt)],0≤k≤Nx.The convolution product[Skn] appearing in (100) can be written as follows:(106)[Skn]=[Z0(2)][Ikn-1/2]+∑l=1n-1[Z0(l+1)-Z0(l)][Ikn-l+1/2]-[Z0(n)][Ik1/2],
where (107)[Z0(l)]=[Z0w(l)]+[Z0g(l)]=∫l-1l[Zw+Zg](uΔt)du,1≤l≤n,t=nΔt, Nx is the number of space steps, and n is the number of time steps.The corresponding components of theM×1 vectors [ET,kn], [EL,kn] are as follows:(108)[VT,kn]i=yiEyinc(xi=kΔx,yi,zi,nΔt)+ziEzinc(xi=kΔx,yi,zi,nΔt),(109)[EL,kn]i=Exinc(xi=kΔx,yi,zi,nΔt)-Ezinc(xi=kΔx,0,0,nΔt).Equations (107) and (108) are valid for 1≤i≤M, where M is the number of conductors, and  yi,zi are the positions of ith conductor. Exinc,Eyinc,Ezinc are the components of the incident electromagnetic field evaluated in the absence of conductors. [R],[L],[G],and[C] are, respectively, the per-unit-length resistance, inductance, conductance, and capacitance matrices of dimension  M×M.[Zg(t)] is the transient ground resistance and is equal to the inverse Fourier of [Zg(s)/s]:(110)Zgij(t)=F-1(Zgij(s)s),
where the ground impedance in frequency domain is given by Carson formula:(111)Zgij(s)s=μ0π∫0∞e-(hi+hj)λλ2+γg2+λ⋅cos⁡(dijλ)dλ.[Zw(t)] is a diagonal matrix and corresponds to the transient conductor resistance, where each element is given by(112)Zwi(t)=1πσwiai2∑m=1∞e-xm2(t/τwi),whereτwi=μ0σwiai2,
and hi,hj,dij are the corresponding heights of the two conductors i,j and the distance between the two conductors in the horizontal plane.Termsxm stand for the zeros of J1, Bessel function of first kind.Finally,γg is the propagation constant defined as(113)γg2=sμ0(σg+sε0εrg),
where σg and εrg are, respectively, ground conductivity and permittivity.For the cases considered in this paper, all conductors are in open circuit at both ends, so the currents vanishes at near and far ends.By using the boundary formulation [11], voltages at both ends are simply expressed as follows:(114)[V0n]=[V(1/2)n],[VNxn]=[VNx-(1/2)n],(115)[I0n]=[0],[INxn]=[0].


## 4. Numerical Results

Figure11 shows the transient response at the centre of the straight wire L=20 m, a=0.005m, located at height h=1m above a dielectric half-space (εr=10) excited by the electromagnetic pulse (EMP):(116)Exinc=E0(e-at-e-bt),
with E0=1.1V/m, a=7.92*104s-1,b=4*104s-1.Figure 11
The transient current induced at the center of the line above dielectric half-space (εr=10).The next example is related to a transient scattering from a straight thin wire of length isL=1m, radius a=2mm, located at heighth=0.25m above ground with permittivity εr=10, while the conductivity is varied. The wire is illuminated by the tangential electromagnetic pulse (EMP) plane wave with E0=1V/m, a=4·107s-1, b=6·108s-1.Figure12 shows the transient current induced at the wire center for different ground conductivities.Figure 12
Transient current at the wire center,L=1 m, a=2mm, h=0.25, εr=10.The influence of the ground conductivity to the transient response is particularly visible from around 0.1 S/m to 1 S/m.The time domain results obtained via different approaches are found to agree satisfactorily.Next set of examples is related to a two-wire array above a PEC ground Figure13, (Geometry No. 1) and dielectric-half space (εr=10), Figure 14 (Geometry No. 2), respectively.Figure 13
Geometry No. 1: Two-wire array above a PEC ground (a=2cm, L=10m, d=1m, h1=1m, and h2=2m).Figure 14
Geometry No. 2: Two-wire array above a dielectric half-space (εr=10, a=2cm, L=10m, d=1m, h1=1m, and h2=2m).Figures15 and 16 show the transient current induced at the center of wire 2 for the case of Geometry No. 1 and No. 2, respectively, obtained via TD GB-IBEM, TL, and NEC 2 combined with inverse fast Fourier transform (IFFT).Figure 15
Transient current induced at the center of wire 2 (Geometry No. 1)—comparison between IFFT-NEC2, TD GB-IBEM and TL results.Figure 16
Transient current induced at the center of wire 2 (Geometry No. 2)—comparison between IFFT-NEC2, TD GB-IBEM and TL results.Generally, the results calculated via different approaches are in relatively acceptable agreement. Nevertheless, some discrepancies can be noticed, in particular for the case of PEC ground. In this analysis, the applied TL model accounts not only for classical propagation effect but also for skin effects and for a correction resistance representing the radiation effect.In order to include the radiation effect in TL model, a small DC resistance (1 Ω/m) has been used to to represent the attenuation effect.It is known that TL model accounts for coupling between transverse cells only while AT takes into account mutual effects. This phenomenon is assumed to be the source of the differences in propagation velocities which are observed.


## 5. Conclusion

The paper reviews the models and methods used for the analysis of electromagnetic field coupling to overhead wires in the frequency and time domain, respectively, using the wire antenna theory and the transmission line approximation, respectively. The frequency domain wire antenna model is based on the space-time Pocklington integral equations, while the transmission line model is based on the frequency domain Telegrapher's equations. The time domain wire antenna model is based on a set of the space-time Hallen integral equations, while the transmission line model is based on the time domain Telegrapher's equations. The set of Pocklington equations is solved via the Galerkin-Bubnov variant of the Indirect Boundary Element Method (GB-IBEM), while the frequency domain transmission line equations are treated using the chain matrix method and modal equation to derive per-unit-length parameters. A number of illustrative computational examples for the frequency response of several configurations of overhead wires, obtained via different approaches, are given in this paper.The coupled space-time Hallen integral equations are numerically solved via the time domain Galerkin-Bubnov scheme of the Indirect Boundary Element Method (GB-IBEM), while the time-domain transmission line equations are solved by the finite difference time domain (FDTD) method.Some numerical results pertaining to the transient behaviour of overhead wires, obtained via different approaches, are given in this paper.



---

*Source: 730145-2012-04-12.xml*