On the Electromagnetic Field of an Overhead Line Current Source

Abstract

This work presents an analytical series-form solution for the time-harmonic electromagnetic (EM) field components produced by an overhead current line source. The solution arises from casting the integral term of the complete representation for the generated axial electric field into a form where the non-analytic part of the integrand is expanded into a power series of the vertical propagation coefficient in the air space. This makes it possible to express the electric field as a sum of derivatives of the Sommerfeld integral describing the primary field, whose explicit form is known. As a result, the electric field is given as a sum of cylindrical Hankel functions, with coefficients depending on the position of the field point relative to the line source and its ideal image. Analogous explicit expressions for the magnetic field components are obtained by applying Faraday’s law. The results from numerical simulations show that the derived analytical solution offers advantages in terms of time cost with respect to conventional numerical schemes used for computing Sommerfeld-type integrals.

1. Introduction

The computation of the electromagnetic (EM) fields from overhead electric lines located above dissipative terrestrial areas is a classical problem which is still of interest today because of the public concern about the biological effects of field exposure. In fact, it is well known that the strong EM fields associated with the high-intensity currents and voltages of power lines can generate undesirable effects on humans, animals and other forms of life [1,2,3]. The problem, in all of its variants, has spawned a number of contributions to the scientific literature since a solution was first developed at the beginning of the 20th century by Carson [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. Yet, in spite of the variety of proposed approaches, to Author’s knowledge there is scarcity of purely analytical techniques in literature, and most of the published methods attempt to solve the problem through the derivation of analytical formulations amenable to numerical treatment. Excellent illustrations of these approaches are the well-known method of moments (MoM), which is used to solve surface integral formulations, and the method of auxiliary sources [10], consisting of expressing the fields as weighted superpositions of the contributions from a finite number of fictitious currents flowing on mathematical surfaces separated from the air-ground interface. Even if these methods may exhibit good performances in terms of accuracy and efficiency, as numerical procedures they suffer from the intrinsic disadvantage of involving memory requirements and computational costs significantly greater than those implied by an analytical solution. Moreover, they provide less insight in the physics of the problem and are also less suitable for sensitivity analysis. On the other hand, one important contribution in the context of purely analytical techniques is the work by Prof. J. R. Wait [5], which presents field expressions derived through usage of the complex image theory. The only drawback of the obtained formulas resides in that they are valid in the quasi-static regime only, that is when the effects of the displacement currents in the air space are negligible. The aim of the present paper is to derive rigorous series-form expressions for the time-harmonic EM field components produced in the air-space by a uniform-current line source located above a homogeneous dissipative ground. The field expressions are obtained starting from decomposing the integral representation for the axial component of the generated electric field into three parts, that is the direct field induced by the source current, the ideal reflected field induced by a negative image current and a correction term due to the imperfect conductivity of the ground. Next, after recognizing that the first two terms can be straightforwardly expressed in explicit form, the non-analytic part of the integrand of the correction term is expanded into a power series of the vertical propagation coefficient in the air space. This permits to express the electric field as a sum of derivatives of the Sommerfeld Integral describing the direct field. As a result, the field is given as a sum of cylindrical Hankel functions, with coefficients depending on the position of the field point relative to the original current source and its ideal image. Finally, analogous explicit expressions for the magnetic field components are obtained by applying Faraday’s law. The obtained solution is not subject to simplifying assumptions, and hence is valid even when the high-frequency effects due to the displacement currents in both the air and the soil are not negligible. As a consequence, the solution is also of practical use as an analytical benchmark for simulation tools employed to solve EM boundary value problems. Numerical tests are performed to show the validity of the developed method and its advantages in terms of computation time with respect to standard numerical algorithms used to evaluate Sommerfeld-type integrals.

2. Formulation of the Problem

Consider an infinite line source of current located in proximity of a flat homogeneous ground, as shown in Figure 1.

The line is situated above the ground at height h, and the material medium is characterized by the dielectric permittivity ${ϵ}_1$ and electrical conductivity ${\sigma }_1$, while the magnetic permeability is taken to be everywhere that of free-space ${\mu }_0$. Under the hypothesis that the observation point lies in the near-field region of the line source, it is reasonable to assume that the line supports a nearly uniform current [1], i.e., $I(x,t)\cong I{e}^{j\omega t}$. Due to the symmetry of the problem, the total electric field generated in the air space has only one component in the axial direction, whose integral representation is well known and given by [10]

$$\begin{array}{c}\hfill E_x=-\frac{j\omega {\mu }_{0}I}{2\pi }\left[S_d\left(|z-h|\right)+S_r\left(z+h\right)\right],\end{array}$$

(1)

where

$$\begin{array}{cc}\hfill S_d\left(\zeta \right)& ={\int }_0^{\infty }\frac{1}{u_0}{e}^{-u_0\zeta }cos\left(\lambda y\right)d\lambda =-\frac{j\pi }{2}{H}_0^{\left(2\right)}\left(k_0\sqrt{y^2+{\zeta }^2}\right)\hfill \end{array}$$

(2)

and

$$\begin{array}{cc}\hfill S_r\left(\zeta \right)& ={\int }_0^{\infty }\frac{1}{u_0}\frac{u_0-u_1}{u_0+u_1}{e}^{-u_0\zeta }cos\left(\lambda y\right)d\lambda \hfill \end{array}$$

(3)

are, respectively, the Sommerfeld Integrals describing the direct and reflected fields, being $H_0^{\left(2\right)}$ the zeroth-order Hankel function of the second kind, and $u_n=\sqrt{{\lambda }^2-k_n^2}$. The aim of this paper is to derive exact series representations for $E_x$ and the magnetic field components, the latter arising from applying Faraday’s law. To this goal, it is first convenient to use the identity

$$\begin{array}{c}\hfill \frac{u_0-u_1}{u_0+u_1}=-1+\frac{2u_0}{u_0+u_1}\end{array}$$

(4)

and split $S_r$ into two terms, that is the ideal reflected field contribution (induced by a negative image current) plus a term due to the imperfect conductivity of the ground. It reads

$$\begin{array}{c}\hfill S_r\left(z+h\right)=-S_d\left(z+h\right)+S_c\left(z+h\right),\end{array}$$

(5)

with

$$\begin{array}{c}\hfill S_c\left(\zeta \right)=2{\int }_0^{\infty }\frac{e^{-u_0\zeta }}{u_0+u_1}cos\left(\lambda y\right)d\lambda ,\end{array}$$

(6)

and where account has been taken of (2). The integral $S_c$ may be evaluated by proceeding as follows. After setting $\alpha =\sqrt{k_1^2/k_0^2-1}$ and

$$\begin{array}{c}\hfill \frac{1}{u_0+u_1}=\frac{u_0}{k_0^2{\alpha }^2}+jP\left(u_0\right),\end{array}$$

(7)

we replace the quantity $e^{u_0\beta }P\left(u_0\right)$, where $\beta$ is a non-negative real constant to be determined, with its Taylor expansion about $u_0$ = 0, in a similar fashion as in [5]. It yields

$$\begin{array}{c}\hfill f\left(u_0\right)=e^{u_0\beta }P\left(u_0\right)=\sum _{i=0}^{\infty }\frac{f^{\left(i\right)}\left(0\right)}{i!}{u}_0^i,\end{array}$$

(8)

and, as a consequence, (6) becomes

$$\begin{array}{cc}\hfill S_c\left(\zeta \right)=& \frac{2}{{\left(k_0\alpha \right)}^2}{\int }_0^{\infty }{u}_0{e}^{-u_0\zeta }cos\left(\lambda y\right)d\lambda +2j\sum _{i=0}^{\infty }{a}_i\left(\beta \right){\int }_0^{\infty }{u}_0^i{e}^{-u_0(\zeta +\beta )}cos\left(\lambda y\right)d\lambda ,\hfill \end{array}$$

(9)

being $a_i\left(\beta \right)$ = $f^{\left(i\right)}\left(0\right)/i!$. Under the assumption that $\zeta$ > 0, the second integral on the right-hand side of (9) converges regardless of the value chosen for the non-negative arbitrary constant $\beta$. This happens because, as $\lambda$ becomes sufficiently large, the exponential factor $e^{-u_0(\zeta +\beta )}$ approaches $e^{-\lambda (\zeta +\beta )}$, which rapidly decays with increasing $\lambda$. Since $\beta$ may be assumed to be as small as desired, we now take the limit of (9) as $\beta \to 0^{+}$. It is not difficult to prove that, subject to this condition, the odd coefficients $a_{2i+1}$ ($i=0,1,2,\dots$) become null. On the other hand, the limit of the even coefficients reads

$$\begin{array}{c}\hfill \underset{\beta \to 0^{+}}{lim}{a}_{2i}\left(\beta \right)=\frac{c_{2i}}{{\left(k_0\alpha \right)}^{2i+1}},\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}{c}_i=\frac{(i+1)!!}{(i^2-1)i!!},\end{array}$$

(10)

and, as a consequence, (9) may be expressed as

$$\begin{array}{cc}\hfill S_c\left(\zeta \right)=& \frac{2}{{\left(k_0\alpha \right)}^2}\frac{{\partial }^2{S}_d}{\partial {\zeta }^2}-2j\sum _{i=0}^{\infty }\frac{c_{2i}}{{\left(k_0\alpha \right)}^{2i+1}}\frac{{\partial }^{2i+1}{S}_d}{\partial {\zeta }^{2i+1}},\hfill \end{array}$$

(11)

where the $\zeta$-derivatives of $S_d$, that is the derivatives of $H_0^{\left(2\right)}$, are to be made explicit. After letting let r = $\sqrt{y^2+{\zeta }^2}$, from the differential properties of the Bessel functions [26] it follows that

$$\begin{array}{cc}\hfill \frac{\partial H_0^{\left(2\right)}}{\partial \zeta }& =-\frac{k_0\zeta }{r}{H}_1^{\left(2\right)}\hfill \\ \hfill \frac{{\partial }^2{H}_0^{\left(2\right)}}{\partial {\zeta }^2}& =-k_0\left[1+{\zeta }^2\left(\frac{1}{r}\frac{d}{dr}\right)\right]\frac{H_1^{\left(2\right)}}{r}\hfill \\ \hfill \frac{{\partial }^3{H}_0^{\left(2\right)}}{\partial {\zeta }^3}& =-k_0\left[3\zeta \left(\frac{1}{r}\frac{d}{dr}\right)+{\zeta }^3{\left(\frac{1}{r}\frac{d}{dr}\right)}^2\right]\frac{H_1^{\left(2\right)}}{r}\hfill \\ \hfill \frac{{\partial }^4{H}_0^{\left(2\right)}}{\partial {\zeta }^4}& =-k_0\left[3\left(\frac{1}{r}\frac{d}{dr}\right)+6{\zeta }^2{\left(\frac{1}{r}\frac{d}{dr}\right)}^2+{\zeta }^4{\left(\frac{1}{r}\frac{d}{dr}\right)}^3\right]\frac{H_1^{\left(2\right)}}{r}\hfill \\ \hfill & ⋮\hfill \\ \hfill \frac{{\partial }^l{H}_0^{\left(2\right)}}{\partial {\zeta }^l}& =-k_0\sum _{m=0}^{\lfloor l/2\rfloor }{(-1)}^m{d}_{lm}{\zeta }^{l-2m}{\left(\frac{1}{r}\frac{d}{dr}\right)}^{l-m-1}\frac{H_1^{\left(2\right)}}{r},\hfill \end{array}$$

(12)

where the argument of the Hankel functions, that is $k_{0}r$, has been omitted for notational simplicity, and with

$$\begin{array}{c}\hfill d_{lm}={(-1)}^m\left(\genfrac{}{}{0pt}{}{l}{l-2m}\right)\left(2m-1\right)!!=\frac{{(-1)}^{m}l!}{(l-2m)!\left(2m\right)!!}.\end{array}$$

(13)

It should be noted that the number of terms on the right-hand sides of (12) is $(l/2+1)$ for even l, and $\left[\right(l-1)/2+1]$ for odd l. Hence, it is equal to unity plus the integer part of $(l/2)$, denoted by $\lfloor l/2\rfloor$. Next, use of the tabulated result [26]

$$\begin{array}{c}\hfill {\left(\frac{1}{r}\frac{d}{dr}\right)}^{\nu }\frac{H_1^{\left(2\right)}}{r}={(-k_0)}^{\nu }\frac{H_{\nu +1}^{\left(2\right)}}{r^{\nu +1}}\end{array}$$

(14)

into (12) leads to the expression

$$\begin{array}{c}\hfill \frac{{\partial }^l{\varphi }_0}{\partial {\zeta }^l}={(-k_0)}^l\sum _{m=0}^{\lfloor l/2\rfloor }{d}_{lm}\frac{{\varphi }_{l-m}}{{\left(k_0\zeta \right)}^m},\end{array}$$

(15)

and

$$\begin{array}{c}\hfill {\varphi }_n={\left(\frac{\zeta }{r}\right)}^n{H}_n^{\left(2\right)}\left(k_{0}r\right).\end{array}$$

(16)

Finally, substituting (2) and (15) into (11) provides

$$\begin{array}{cc}\hfill S_c\left(\zeta \right)& =j\pi {\Phi }_{2i+1,1},\hfill \end{array}$$

(17)

with

$$\begin{array}{cc}\hfill {\Phi }_{ln}& =\frac{{\varphi }_n}{{\alpha }^2{k}_0\zeta }-\frac{{\varphi }_{n+1}}{{\alpha }^2}-j\sum _{i=0}^{\infty }\frac{c_{2i}}{{\alpha }^{2i+1}}\sum _{m=0}^{\lfloor l/2\rfloor }{d}_{lm}\frac{{\varphi }_{2i-m+n}}{{\left(k_0\zeta \right)}^m},\hfill \end{array}$$

(18)

and hence, the total axial electric field component $E_x$ (1) may be rewritten as

$$\begin{array}{cc}\hfill E_x=& -\frac{\omega {\mu }_{0}I}{4}\left[{\left.{\varphi }_0\right|}_{\zeta =|z-h|}-{\left({\varphi }_0+2{\Phi }_{2i+1,1}\right)}_{\zeta =z+h}\right].\hfill \end{array}$$

(19)

The non-null magnetic field components may be obtained from $E_x$ by applying Faraday’s law, as follows

$$\begin{array}{cc}\hfill H_y& =-\frac{1}{j\omega {\mu }_0}\frac{\partial E_x}{\partial z}=\frac{I}{2\pi }\left\{{\left.\frac{\partial S_d\left(\zeta \right)}{\partial \zeta }\right|}_{\zeta =z-h}-{\left[\frac{\partial S_d\left(\zeta \right)}{\partial \zeta }-\frac{\partial S_c\left(\zeta \right)}{\partial \zeta }\right]}_{\zeta =z+h}\right\},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill H_z& =\frac{1}{j\omega {\mu }_0}\frac{\partial E_x}{\partial y}=-\frac{I}{2\pi }\left\{{\left.\frac{\partial S_d\left(\zeta \right)}{\partial y}\right|}_{\zeta =z-h}-{\left[\frac{\partial S_d\left(\zeta \right)}{\partial y}-\frac{\partial S_c\left(\zeta \right)}{\partial y}\right]}_{\zeta =z+h}\right\}.\hfill \end{array}$$

(21)

In fact, differentiating (17) with respect to $\zeta$ leads to write

$$\begin{array}{c}\hfill \frac{\partial S_c\left(\zeta \right)}{\partial \zeta }=-j\pi k_0\left(\frac{2{\varphi }_2}{{\alpha }^2{k}_0\zeta }+{\Phi }_{2i+2,2}\right)\end{array}$$

(22)

and, as a consequence, (20) is turned into

$$\begin{array}{c}\hfill H_y=\frac{j{k}_{0}I}{4}\left[{\left.{\varphi }_1\right|}_{\zeta =z-h}-{\left({\varphi }_1+\frac{4{\varphi }_2}{{\alpha }^2{k}_0\zeta }+2{\Phi }_{2i+2,2}\right)}_{\zeta =z+h}\right].\end{array}$$

(23)

On the other hand, since it holds

$$\begin{array}{cc}\hfill \frac{{\partial }^{l+1}{\varphi }_0}{\partial y\partial {\zeta }^l}& =-k_0\sum _{m=0}^{\lfloor l/2\rfloor }{(-1)}^m{d}_{lm}y{\zeta }^{l-2m}{\left(\frac{1}{r}\frac{d}{dr}\right)}^{l-m}\frac{H_1^{\left(2\right)}}{r}\hfill \\ \hfill & ={(-k_0)}^{l+1}\frac{y}{\zeta }\sum _{m=0}^{\lfloor l/2\rfloor }{d}_{lm}\frac{{\varphi }_{l-m+1}}{{\left(k_0\zeta \right)}^m}\hfill \end{array}$$

(24)

and, therefore,

$$\begin{array}{c}\hfill \frac{\partial S_c\left(\zeta \right)}{\partial y}=-\frac{j\pi k_{0}y}{\zeta }{\Phi }_{2i+1,2}\left(\zeta \right),\end{array}$$

(25)

the explicit expression for the $H_z$-field reads

$$\begin{array}{c}\hfill H_z=-\frac{j\pi k_{0}yI}{4\pi \zeta }\left[{\left.{\varphi }_1\right|}_{\zeta =z-h}-{\left({\varphi }_1+2{\Phi }_{2i+1,2}\right)}_{\zeta =z+h}\right].\end{array}$$

(26)

3. Results

As validation of the theoretical development, expressions (23) and (26) were used to calculate the magnetic field components that an infinite line source carrying 1 A of current produced at a plane 1 m above the interface between air and a clay soil. The source was situated 4 m above the material medium, whose electrical conductivity and dielectric permittivity were taken to be equal to ${\sigma }_1$ = 0.1 mS/m and ${ϵ}_1$ = 40 ${ϵ}_0$, respectively [27]. At first, the fields were computed against the horizontal distance y from the line axis, assuming that the operating frequency was equal to 1 MHz. Four y-profiles were generated, each one corresponding to a different value for the truncation index L of the outer infinite sum in (18). The results of the computations, shown in Figure 2 and Figure 3, were compared with those arising from numerically evaluating the integral representations for the magnetic field components, namely [28]

$$\begin{array}{cc}\hfill H_y=& \frac{I}{2\pi }{\int }_0^{\infty }\left[e^{-u_0\left|z-h\right|}-\frac{u_0-u_1}{u_0+u_1}{e}^{-u_0(z+h)}\right]cos\left(\lambda y\right)d\lambda ,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill H_z=& \frac{I}{2\pi }{\int }_0^{\infty }\frac{\lambda }{u_0}\left[e^{-u_0|z-h|}+\frac{u_0-u_1}{u_0+u_1}{e}^{-u_0(z+h)}\right]sin\left(\lambda y\right)d\lambda .\hfill \end{array}$$

(28)

Numerical integration was carried out by applying a Gauss–Kronrod G7-K15 scheme, originating from the combination of a seven-point Gauss rule with a 15-point Kronrod rule. From the analysis of the plotted curves it emerged that increasing the truncation index L improved the accuracy of the result of the computation. In fact, if L grew the curves provided by (23) and (26) approached the outcomes from numerical quadrature, and close agreement was achieved when L = 9 in both the situations. Thus, the proposed series-form solution converged to the exact solution. This was confirmed by the curves plotted in Figure 4, which show the relative error of the outcomes from (23) as compared to numerical integration data. As is seen, the error monotonically decreased as L increased. On the other hand, for a fixed value of L the error was not substantially affected by a variation of the distance y as long as the latter was smaller than 50 m. Thereafter, the relative error grew and exhibited an oscillating behavior, pretty similar for all the considered values of the truncation index of the outer sum in (18). However, even with the rapid increase as the distance y increased, for $L=9$ the error did not exceed ${10}^{-3}$ in the considered interval.

The accuracy of the result of the computation also depended on the values of the electromagnetic parameters of the lossy ground and, for instance, better accuracy was observed for larger values of the electrical conductivity ${\sigma }_1$. This aspect is illustrated by Figure 5 and Figure 6, which show, respectively, profiles of $|H_y|$ versus the ratio ${\sigma }_1/\left(\omega {ϵ}_1\right)$, and the relative error resulting from using (23) instead of G7-K15 scheme. Here, it is assumed that the source carried 1 A and operated at 100 kHz, and that ${ϵ}_1$ = 20 ${ϵ}_0$, h = 5 m, y = 10 m, and z = 1 m. As is seen, the sequence of profiles of $|H_y|$ generated by the partial sums in (23) converged to the curve provided by the G7-K15 scheme, and perfect agreement between the analytical and numerical data was again obtained for L = 9. Furthermore, convergence was faster in the good-conductor limit (${\sigma }_1\gg \omega {ϵ}_1$), where it sufficed to use a sum constituted by seven terms to produce sufficiently accurate results. This does not necessarily imply that the relative error generated by each partial sum in (23) always decreased as the conductivity grew. In fact, as pointed out by Figure 6, the error fluctuated around a decreasing mean value rather than being monotonically decreasing.

The proposed method made it possible to achieve time savings compared to G7-K15 numerical scheme, while maintaining the same accuracy as the latter. In fact, on a single-core 2.2 GHz processor, the average CPU time taken by (23) to generate the profile of $|H_y|$ associated with L = 9 in Figure 5 was equal to 157 ms, while about 9 s were taken by numerical integration of (27) to produce the same outcome. This means that, limited to the considered example, the proposed approach was about 57 times faster than Gaussian integration, and, as a consequence, the numerical complexity of the former was significantly smaller than that exhibited by the latter. For the sake of completeness, the time costs of (23) corresponding to the remaining values of L depicted in Figure 5 are indicated in

Table 1. CPU time comparisons for the calculation of $H_y$.

Approach	Average CPU Time [s]	Speed-Up
G7-K15	$8.97$	-
(23) with $L=4$	$3.13\times {10}^{-4}$	$2.87\times {10}^4$
(23) with $L=5$	$1.26\times {10}^{-3}$	$7.12\times {10}^3$
(23) with $L=7$	$3.84\times {10}^{-2}$	$2.34\times {10}^2$
(23) with $L=9$	$1.57\times {10}^{-1}$	$57.1$

.

Finally, one would ask whether the proposed solution may still be used when both h and z approach zero, since it has been obtained starting from expression (9), whose integrals on the right-hand side, strictly speaking, converged only for $\zeta$ = $z+h>0$. This point is clarified by Figure 7, which shows the amplitude–frequency spectrum of the $H_y$-field that a unit-current line source generated at a point placed 30 m apart from it, assuming that both the source and the observation point lay on a homogeneous ground with ${\sigma }_1$ = 20 mS/m and ${ϵ}_1$ = 10 ${ϵ}_0$. As is evident from the analysis of the plotted curves, the trends arising from the sequence of partial sums in (23) still converged to the curve resulting from numerical integration of (27), even if the number of terms of (23) required to approach the numerical data grew as frequency decreased.

Hence, the rate of convergence of (23) depended on frequency, and convergence was faster at higher frequencies. This aspect is further investigated in Figure 8, which depicts the relative error of the outcomes of (23) plotted in Figure 7, as compared to the data arising from the G7-K15 rule. The curves of Figure 8 point out that, for a fixed L, the accuracy of the result of the computation significantly worsened when entering the low-frequency range. At the same time, they confirmed that the accuracy could always be enhanced by increasing L. This implies that, all over the considered frequency range, it holds [29]

$$\begin{array}{c}\hfill \underset{L\to \infty }{lim}\frac{|H_y^{(L+1)}-H_y|}{|H_y^{\left(L\right)}-H_y{|}^p}=q,\end{array}$$

(29)

where the symbol $\left\{H_y^{\left(L\right)}\right\}$ denotes the sequence of partial sums that originates from (23), as a result of truncating the infinite sum in (18), and where p and q are the order of convergence (OC) of the sequence and the asymptotic error constant (AEC), respectively. It is easily understood that the parameters p and q depend on the operating frequency, and that the knowledge of their values permits to acquire further information on the rate of convergence. As an example,

Table 2. Estimated OC and AEC for the sequence $\left\{H_y^{\left(L\right)}\right\}$.

L	${\mathit{p}}_{\mathit{L}}$	${\mathit{q}}_{\mathit{L}}$
4	0.857	0.674
5	0.871	0.591
7	0.944	0.563
10	0.989	0.495
13	0.996	0.442

illustrates Lth order estimates of the OC and the AEC for the sequence $\left\{H_y^{\left(L\right)}\right\}$ corresponding to 1 kHz in Figure 7. The estimates were calculated using the well-known expressions [29]

$$\begin{array}{cc}\hfill p_L& =\frac{log\left[|H_y^{(L+1)}-H_y^{\left(L\right)}|/|H_y^{\left(L\right)}-H_y^{(L-1)}|\right]}{log\left[|H_y^{\left(L\right)}-H_y^{(L-1)}|/|H_y^{(L-1)}-H_y^{(L-2)}|\right]},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill q_L& =\frac{|H_y^{(L+1)}-H_y^{\left(L\right)}|}{|H_y^{\left(L\right)}-H_y^{(L-1)}{|}^{p_L}}\hfill \end{array}$$

(31)

whose limits as $L\to \infty$ approach p and q [29]. As pointed out by Table 2, $p_L\to 1$ with increasing L, and this means that the sequence $\left\{H_y^{\left(L\right)}\right\}$ converges linearly. Moreover, convergence is accelerated by the small values of the $q_L$’s, which further reduce the remainder $|H_y-H_y^{\left(L\right)}|$ at any additional iteration of the sequence.

4. Conclusions

This paper has presented a rigorous analytical approach for evaluating the time-harmonic EM field components generated in the air space by an infinite current line source located above a homogeneous lossy ground. The approach consists of expanding the non-analytic part of the integrand of the correction term in the integral expression for the axial electric field into a power series of the z-directed propagation coefficient in air. This leads to express the axial electric field as a sum of derivatives of the Sommerfeld Integral describing the direct field, which may be analytically evaluated. As a result, the electric field is given as a sum of cylindrical Hankel functions, with coefficients depending on the position of the field point relative to the line source and its ideal image. Then, explicit expressions for the magnetic field components are also derived by applying Faraday’s law. The obtained solution is not subject to simplifying assumptions, and hence is valid even when the high-frequency effects due to the displacement currents in both the air and the soil are no longer negligible. Numerical simulations have been carried out to show that the proposed approach exhibits good accuracy, while taking significantly less computation time than conventional numerical quadrature schemes used to evaluate Sommerfeld-type integrals.