Electromagnetic Field Analysis of an Electric Dipole Antenna Based on a Surface Integral Equation in Multilayered Dissipative Media

In this paper, a novel method based on the Poggio–Miller–Chang-Harrington–Wu–Tsai (PMCHWT) integral equation is presented to study the electromagnetic fields excited by vertical or horizontal electric dipoles in the presence of a layered region which consists of K-layered dissipative media and the air above. To transform the continuous integral equation into a block tridiagonal matrix with the feature of convenient solution, the Rao–Wilton–Glisson (RWG) functions are introduced as expansion and testing functions. The electromagnetic fields excited by an electric dipole are calculated and compared with the available results, where the electric dipole antenna is buried in the non-planar air–sea–seabed, air–rock–earth–mine, and multilayered sphere structures. The analysis and computations demonstrate that the method exhibits high accuracy and solving performance in the near field propagation region.


Introduction
The electromagnetic fields excited by a dipole source in a layered dissipative medium structure have been investigated widely because of their applications in Through-the-Earth communication, underwater communication, ground penetrating radar technology, and antenna design [1][2][3][4][5][6]. The developments in electromagnetic field propagation in layered medium structures have been analyzed by many investigators, including Wait [7][8][9] and King [7,10,11] who proposed the asymptotic methods, and the surface-impedance technique has been used for a two-layered region. Contour integrations and ranch cuts have been also used for the electromagnetic field in a layered region [8,9]. King et al. have obtained the complete formulas for electromagnetic fields excited by horizontal and vertical electric dipoles in two-and three-layered regions by using the proposed method in the recent studies. However, it is difficult to solve the arbitrary layered region problem. To overcome this drawback, Michalski et al. [12] presented a compact formulation of the electric-type and magnetic-type dyadic Green's functions for a plane-stratified, laterally unbounded multilayered structure, which is based on the tranmission-line network. Nikita et al. [13] has presented a near-fielding radiating dipole antenna next to a three-layered lossy sphere close to a human head model. Moreover, the numerical results have been verified via a unified method of moments (MoM) model based on the electric field integral equation (EFIE) given by Khamas [14]. These algorithms are suitable for analyzing the non-planar electromagnetic problem only when the structure is a concentric sphere, which limits their application in an arbitrary shaped non-planar layered region. Quintana et al. [15] studied the electromagnetic field propagation characteristic of shallow water in detail, and the influence of seabed layer is also considered.
In most studies [16][17][18][19], surface integral equation (SIE)-based methods have been adopted in the numerical analysis of the scattering of arbitrary shaped objects and the scattering of objects buried in layered structure, but few of them focus on the electromagnetic fields excited by dipoles in arbitrary layered dissipative media. In this paper, an SIE-based method to study the electromagnetic fields excited by electric dipole submerged in arbitrary layered homogeneous dissipative media is proposed. The total electromagnetic fields in each layer with three source components including the excitation source itself, equivalent electric/magnetic surface currents at the top of the layer, and equivalent electric/magnetic surface currents at the bottom of the layer are derived in detail. Then, the formulation proposed for the interface of each layer is discretized by using the Galerkin's method. After that, a matrix equation for a layered dissipative structure with a block tridiagonal matrix feature of a convenient solution is derived [20]. Three application scenarios are constructed and computed to analyze the electromagnetic characteristics with the vertical and horizontal dipole buried in wide, practically important media including sea, wet and other dissipative materials. Finally, the numerical results are presented to discuss the difference of the proposed method in comparison with the Computer Simulation Technology (CST) Studio Suite.

SIEs for Multilayered Dissipative Medium Structures
A non-planar K-layered dissipative medium structure is illustrated in Figure 1. Let the frequency of excitation sources in homogeneous medium be ω. The electromagnetic constants at the ith-layer are denoted by conductivity σ i , permittivity ε i , and permeability µ i , with complex wave number k i and wave resistance η i . Let S i−1 and S i denote the top and bottom surfaces of ith-layer medium of the region D i , respectively, where i = 0, 1, 2, . . . , K. In addition, there is no top surface on the half-space region D 0 and no bottom surface beneath the half-space region D K .  Hence, the total electric and magnetic fields at arbitrary point on surface S i−1 or S i can be expressed as The electric and magnetic surface currents on surface S i in region D i are defined as The linear operator L i and K i can be defined as where ∇ s is the surface divergence of a vector field. The Green's function in homogeneous isotropic infinite space is given by (15) Similarly, we can also calculate the electromagnetic fields on surface S i in region D i+1 , Combining Equations (13) and (16), electric field on surface S i can be written as, The magnetic field on surface S i is obtained by combining Equations (15) and (17),

Discretization
For numerical solutions of equivalent electric and magnetic currents on the media interfaces, the Galerkin's method and the Rao-Wilton-Glisson (RWG) basis functions [16,26] are utilized to discretize the continuous formulations Equations (18) and (19) [27]. The media interface of each layer is discretized into planar triangular elements, with T i,n = T + i,n ∪ T − i,n which is defined as the nth triangular element on surface S i , i = 0, 1, 2, ..., K−1 and n = 0, 1, 2, ..., N i .
The RWG basis function f i,n assigned to a pair of triangular elements on surface S i is defined as: where p + i,n and p − i,n are the vertex of triangular elements T + i,n and T − i,n , which share a common edge of length L i,n [18]. A + i,n and A − i,n are the areas of triangular elements T + i,n and T − i,n , respectively. Spatial distribution of RWG functions are shown in Figure 2.
where α i,n and β i,n (i = 0, 1, 2, ..., K−1; n = 0, 1, 2, ..., N i ) are the nth complex expansion coefficients of J i and M i . By substituting (21) into (18) and (19), and testing them with RWG basis, the matrix equation for media interface S i can be written as where, The submatrix in the matrix equation can be expressed as: The basis f j,n (r ), (j = 0, 1, 2, ..., K−1) in medium D i is operated by the linear operator L and K: For r ∈ S i , the electric and magnetic excitation source matrix in D i can be expressed as: According to Equation (22), the matrix equation for all surfaces with block tridiagonal matrix equation is derived as follows: . .
Equation (39) is the relation between the surface current, impedance matrix and incidence fields in K-layered dissipative medium structures. The impedance matrix is a square matrix with the number of elements equal to

Numerical Examples
In this section, three experiments are constructed to analyze the performance of the proposed method with the vertical or horizontal dipole buried in media including sea, wet and other dissipative materials. The algorithm is implemented by MATLAB 2011b on a computer with CPU i3-4030U working at 1.9 GHz with 4 GB RAM and a Windows 10 operating system. The Low Frequency (LF) Domain Solver of CST Studio Suite in the same computer is used to simulate the models for comparison.
We first consider a non-planar air-sea-seabed structure with a hemispherical depression at the seabed shown in Figure 3 for the application that the naval vessel receives signals from the transmitter burried in the seabed. The air permittivity and permeability are ε 0 µ 0 , respectively, which are the same as vacuum. Generally, the conductivity of the sea is σ 1 = 4 S/m, and the permittivity and permeability of the sea are ε 1 = 74ε 0 and µ 1 = µ 0 , respectively. Conductivity of the seabed is σ 2 = 0.01 S/m, with permittivity ε 2 = 13ε 0 and permeability µ 2 = µ 0 . Media interface S 0 and S 1 are both circular with radius r = 50 m, located at plane z = 0 m and z = −20 m, respectively. The hemispherical depression is set at point (0, 0, −20) with radius r = 10 m. A vertical electric dipole excitation source is located at point (0, 0, −25), with dipole moment 1 A · m and the frequency is 100 Hz. In order to reduce the mesh amount, we refine the mesh in the source areas and other part is sparse. The number of triangles on the air-sea interface S 0 and the sea-seabed interface S 1 are 856 and 1528. As a result, 7088 unknowns are generated. In numerical calculations, the processor takes 49.4 s and 0.58 GB to solve the matrix equation. However, the calculation of the same model with 1, 046, 366 unknowns using CST takes 156 s and 2.3 GB of memory.  Figure 4 shows the x-direction electric field and the total electric field at point (x, 0, −60), 0 ≤ x ≤ 25. The results calculated by CST are also plotted in Figure 4 for comparison. These results show clearly that the numerical results calculated by the proposed method are close to the results computed by CST, with most errors within a range of ±2%. However, from the results listed in Table 1, it can be observed that there are relatively large errors in the first and last column, which are larger than 2%. As the electric field in E x close to z-axis is a minimum value, non-uniform discretized mesh can introduce a relatively large error to the first column in numerical calculations. The sparse mesh far from source also causes large error in the last column. It can be also observed from the results in Table 2 that most of the results are within the error range of ±2%. However, there is relatively large error when field point is in the sparse mesh area. The comparisons indicate that the proposed method provides good performance in solving time and high accuracy results for the electric dipole in the two-layered dissipative structure.   For the second example, we consider the air-rock-earth-mine layered structure shown in Figure 5. It can be used in Through-The-Earth communication for mines, where communication devices communicate with other nodes by using horizontal electric dipoles located in the earth layer. The frequency of excitation source is 10 kHz. The conductivity of the rock layer is σ 1 = 0.001 S/m, and the permittivity and permeability of the rock are ε 1 = 3ε 0 and µ 1 = µ 0 , respectively. Conductivity of the earth layer is σ 2 = 0.01 S/m, with permittivity ε 2 = 4ε 0 and permeability µ 2 = µ 0 . The conductivity of the mine layer is σ 3 = 0.001 S/m, and the permittivity and permeability of the mine are ε 3 = 3ε 0 and µ 3 = µ 0 , respectively. Media interfaces S 0 , S 1 and S 2 are circular of radius r = 1000 m. Media interfaces S 0 , S 1 and S 2 are located at plane z = 0 m, z = −150 m and z = −350 m. The x-directional horizontal electric dipole with dipole moment 1A · m is located at (0, 0, −250). The x-component electric field propagation at point (0, y, −500), 0 ≤ y ≤ 1000 is plotted in Figure 6. The z-component electric field propagation and the y-component magnetic field propagation at point (x, 0, −500), 0 ≤ x ≤ 1000 are plotted in Figures 7 and 8, respectively. We can notice that the numerical results calculated by the proposed method are close to the results computed by CST. It can be observed from Tables 3-5 that most of the numerical results are within the error range ±3%. The mesh is sparse at the boundary area, which results in relatively large deviation when the test point gets closer to the boundary, and the magnetic field results are more accurate than electric field results.
Each of the interfaces S 0 , S 1 and S 2 is meshed into 710 triangle patches, resulting in 6288 unknowns. It takes 37.9 s and 0.36 GB to solve the matrix equation by the use of the proposed method. However, the calculation of the same model with 461, 970 unknowns using CST takes 55 s and 1.7 GB memory.
The proposed method is used to simulate the electromagnetic fields excited by horizontal electric dipole in a dissipative sphere, which is a closed structure coated by layered materials and is shown in Figure 9. The external sphere is located at origin and has a radius of 6 m with conductivity σ 1 = 0.001 S/m, permittivity ε 1 = 3.1ε 0 and permeability µ 1 = µ 0 . The radius of the middle sphere is 3 m with conductivity σ 2 = 0.01 S/m, permittivity ε 2 = 2ε 0 and permeability µ 2 = µ 0 , which is located at point (−1, 0, 0). The internal sphere is located at origin and has a radius of 1 m, where the conductivity is σ 3 = 0.1 S/m, permittivity is ε 3 = 2ε 0 , and permeability is µ 3 = µ 0 . A horizontal electric dipole excitation source is located at origin with dipole moment of 1 A · m along x-axis, and the frequency is 1 MHz. After discretization, the numbers of triangles on the external sphere S 0 , the middle sphere S 1 and the internal sphere S 2 are 1204, 604 and 626, respectively. As a result, 7302 unknowns are generated. It takes 54.3 s and 0.66 GB of memory to solve the matrix equation in numerical calculations. Compared with the proposed method, the calculation of the same model with 923, 446 unknowns using CST takes 150 s and 1.9 GB of memory, which shows that the proposed method has a good performance in solving time and storage space.  Figure 10 shows the near magnetic field component H z as a function of the observation angle in plane z = 0, and Figure 11 illustrates the near electric field component E y , which are compared with CST. We can observe that the numerical results obtained from the proposed method are in good agreement with the results obtained by CST. Figure 12 gives the root-mean-square (RMS) error of different observation radius. The Root-mean-square (RMS) error is given by [28]: where x SIE is the calculated quantity and x CST refers to the corresponding CST result. We can see that the H z RMS error is smaller than that of E y of each observation radius. It indicates that the proposed method has a better performance in solving magnetic field than solving electric field. In the CST model, the excitation source is a quasi-ideal horizontal electric dipole. As a result, the comparison model is different from our simulation model, and the RMS error seems to be higher when the observation point is close to the source. When the observation radius are larger than 0.6 m, the numerical results are still in good agreement with the comparison model with RMS error of less than 2% for magnetic fields and 4% for electric fields.   Magnitude of E y [dBµV/m] SIE,r =0.5 CST,r =0.5 SIE,r =0.6 CST,r =0.6 SIE,r =0.7 CST,r =0.7 SIE,r =0.8 CST,r =0.8

Conclusions
In this paper, a novel method based on PMCHWT integral equations has been proposed to study the electromagnetic fields excited by the vertical or horizontal electric dipole in a layered dissipative medium region. The electromagnetic fields in each layer are excited only by the equivalent surface current at the top/bottom of the layer and the dipole in it. A block tridiagonal matrix system is derived, which results in a much more straightforward treatment to handle non-planar layered structure, multilayered planar structure and closed structure coated with arbitrary layered dissipative materials. To analyze the performance of the method, non-planar air-sea-seabed, air-rock-earth-mine and multilayered sphere structures are investigated. The numerical results show that accurate near fields in the layered dissipative medium region are obtained by the method with deviation of 2%. The proposed method can give a better performance with respect to the solving time and storage space in comparison with the CST in handling near-filled problems. Other areas for future work include filling and solving the block tridiagonal matrix using fast/parallel techniques.