Abstract

The augmented electric field integral equation (A-EFIE) with charge neutrality enforcement provides a stable formulation to conquer low-frequency breakdown characteristic of conventional EFIE. It is augmented with additional charge unknowns through current continuity equation. The A-EFIE combined with the multilevel adaptive cross-approximation (MLACA) algorithm is developed to further reduce the memory requirement and computation time for analyzing electromagnetic problems. Numerical examples are given to demonstrate the accuracy and efficiency of the proposed method.

1. Introduction

Electromagnetic integral equations are often discretized with the method of moments (MoM) [1, 2], one of the most widespread and generally accepted techniques for electromagnetic problems. The integral equation is first discretized into a matrix equation using the Galerkin-based MoM with subdomain basis functions such as rooftop functions [3] for curvilinear quad patches and Rao-Wilton-Glisson (RWG) functions [4] for triangular patches. It is convenient to model objects with arbitrary shape using triangular patches; hence, RWG functions are widely used for representing unknown current distributions. When iterative solvers are used to solve the MoM matrix equation, the fast multipole algorithm (FMA) or multilevel fast multipole algorithm (MLFMA) [58] can be used to accelerate the calculation of matrix-vector multiplies in iterative algorithm.

However, the electric field integral equation (EFIE) only works well at mid-frequencies due to the well-known low-frequency breakdown. To circumvent this disadvantage, some methods have been proposed to remedy low-frequency characteristic of electrical field integral equation (EFIE) [912]. Wilton and Glisson proposed loop-tree methods for the Helmholtz decomposition feature at low frequency regime [10]. Mautz and Harrington proposed a loop-star method and analyzed the reason of instability for integral equation in low-frequency regime [11]. Lim et al. introduced loop-star method into RWG functions [12]. Wu et al. provided the comparison between the loop-tree method and loop-star method for the EFIE [13]. However, the loop-star or loop-tree basis functions only partly remove the problem at very low frequencies. Furthermore, the resulting system matrix in MoM is still ill-conditioned leading to slowly converging iterative solutions. Zhao and Chew utilized a basis rearrangement technique [14] to speed up convergence speed of the iterative solver, but it is also time consuming.

A well-conditioned EFIE (WEFIE) to overcome the shortcomings of the EFIE has been proposed in [1523]. Due to the fact that the square of the EFIE operator does not have eigenvalues accumulating at zero or infinity, the WEFIE can give rise to a well-conditioned EFIE system independently of the discretization density [15, 1719]. An efficient Calderón preconditioner is introduced based on the Calderón identities for scattering by perfect electric conductor (PEC) [16, 2023]. The Calderón preconditioner is made purely multiplicative by using the Buffa-Christiansen (BC) basis functions. It has also been used to accelerate the convergence rate of the iterative solution of the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) equations for wave scattering by homogeneous dielectric objects [24].

Some other fast methods have also been proposed to cover the regime from low frequency to mid frequency, such as the plane wave methods based on the generalized Gaussian quadrature rules [25] and the low-frequency fast inhomogeneous plane wave algorithm (LF-FIPWA) [26]. Since evanescent waves are highly direction dependent, much memory is required for LF-FIPWA. The fast multipole algorithm is numerically unstable due to the oscillatory characteristic of the spherical Hankel function for small arguments. As an attempt for a possible remedy, the mixed-form fast multipole algorithm (MF-FMA) [27] is proposed to cover a wide band from the low frequency to mid frequency. However, it is noticed that the MF-FMA need to construct multipole expansions for each nonempty box at the low-frequency regime. This results in a source distribution of high-density for the whole object. The adaptive MF-FMA (AMF-FMA) is proposed for efficient analysis of electromagnetic scattering from objects containing fine structure [28]. The AMF-FMA applies MF-FMA for some parts of the object which belongs to the low-frequency regime and the rest of the object which belongs to the mid-frequency regime is still analyzed by MLFMA. Hence, for a cuboid with many fine structures in [28], AMF-FMA takes advantage of adaptive grouping scheme to set the finest level box size as 0.05λ for fine structure, while the finest level box size is 0.2λ for the rest of the object. Much memory can be saved for AMF-FMA compared to MLFMA.

The augmented electric field equation (A-EFIE) [16, 2933] method has been proposed to solve the low-frequency problem without the loop-tree decomposition. As a new full-wave method, A-EFIE provides a stable formulation over a wide frequency band. In A-EFIE, the contributions of the vector potential and the scalar potential are separated; thus the imbalance between them in the conventional EFIE is avoided. It is augmented with additional charge unknowns through current continuity equation. Qian and Chew emphasize the charge neutrality enforcement (CNE) to remove the rank-deficient probelm in A-EFIE [30]. The complex interconnect packaging problem with over one million unknowns is solved without the help of parallel computers. Qian and Chew also proposed a perturbation method for solving the possible low-frequency inaccuracy problem of A-EFIE [31]. Chen et al. extended the A-EFIE method into the layered medium problem [33]. Yan et al. proposed a Calderón preconditioned A-EFIE algorithm (CP-AEFIE) for analysis of scattering by electrically large objects [16]. The resonance problem of the Calderón preconditioned EFIE is removed by introducing a reasonable complex wavenumber into the Calderón preconditioner [16].

Though the Calderón preconditioning method can eliminate the low-frequency breakdown and the dense-grid breakdown, it is highly technical. For the MLFMA and its low-frequency extensions, however, a priori knowledge of Green’s function is needed. As a result, it cannot be easily applied to analyze the layered media problems. Compared to the MLFMA, the multilevel adaptive cross-approximation (MLACA) algorithm [3436] is another popular technique for analyzing scattering/radiation problems. It makes use of the well-known fact that the approximate rank of the submatrices is low when the subscatterers are sufficiently separated. The MLACA is purely algebraic and therefore does not depend on a priori knowledge of Green’s function. In this paper, the MLACA algorithm is introduced into the A-EFIE with CNE to accelerate the matrix-vector multiplies in iterative algorithm. The ACA algorithm is an adaptive and on-the-fly rank-revealing block factorization of the rank-deficient submatrices. Compared to the singular value decomposition (SVD) technique, the ACA algorithm only requires partial knowledge of the original matrix. It has been shown that for moderate electrical size problems the memory and CPU time requirements for the ACA algorithm scale as O( log N) [35]. However, the set-up time of the MLFMA is less than that of ACA due to the fact that the MLFMA reuses multipole and local expansion information across levels. Jiang et al. utilize the predetermined interaction list supported octree (PILOT) [36, 37] to reduce the set-up time and the memory consumption of the ACA algorithm.

The aim of this paper is to apply the MLACA algorithm to improve the efficiency of the A-EFIE with CNE for analyzing electromagnetic problems. This paper is organized as follows. Section 2 describes the theory and implementation of the A-EFIE combined with the ACA algorithm in more details. Section 3 gives some numerical examples to demonstrate the accuracy and efficiency of our approach. Section 4 gives some conclusions.

2. A-EFIE with MLACA Algorithm

2.1. A-EFIE Formulation

The conventional EFIE is given as where the is free space wave impedance and is the wavenumber. denotes the right-hand vector (excitation). is the unknown current. The and correspond to the magnetic vector potential and electric scalar potential: where the subscripts , denote the interaction between the th line and the th line. denotes free space scalar Green’s function. denotes RWG basis function. and represent the relative permeability and the relative permittivity, respectively. In A-EFIE, we use a normalized RWG function to expand the surface current by removing the edge length from the original RWG basis function: where the is the area of the triangle and is the free vertices of the triangle pair. Then, we use pulse function to expand the charge on each element which is defined as follows: The expansion equations of surface current and charge are given as where the is the number of inner lines and is the number of mesh elements. and are the unknown expansion coefficients for surface current and charge, respectively. According to the current continuity equation: we define a transform matrix as follow [29]: Then, the scalar potential matrix can be factorized as where the denotes a patch-to-patch scalar potential matrix as And we substitute (5) and (6); into (7) we get where the is the light speed in vacuum. Substituting the above equations into the EFIE matrix equation and enforcing the current continuity equation explicitly, we can arrive at the following A-EFIE system. So, we get the form of A-EFIE as where the is an identity matrix with dimension of .

At very low frequencies, charge neutrality causes rank deficiency in the above forms of the A-EFIE [30]. Therefore, we drop a charge unknown from each unattached objects due to the charge neutrality. Then, the vector of charge is reduced to be . is the dimension of reduced charge vector and equal to ( denotes the number of disconnected objects for a problem). We define two mapping matrixes and : where the represents the previous charge density coefficient in (6). The matrix maps the full vector forward to the reduced one, and the matrix projects the reduced vector backward to the full one. Both of them are highly sparse matrices. Then, the AEFIE in (12) is given as where the is an identity matrix with dimension of . The linear system of equations in (14) can be solved by the restart GMRES iterative method with preconditioning technique to accelerate the convergence of the iterative solver [5, 6, 31]. A simple and efficient preconditioning matrix has been given as [31] where the is the diagonal matrix of . Since is a sparse matrix, its inverse matrix can be solved by fast sparse matrix solver such as UMFPACK solver in this paper.

For the analysis of microstrip antenna and microwave integrated circuit (MIC) in layered medium, free space scalar Green’s function in (2) should be substituted by layered media Green’s function [38]. In general, spatial domain Green’s functions are expressed in terms of Sommerfeld integrals. Due to the highly oscillatory nature of the integrand, numerical integration is very time consuming. In this paper, a two-level generalized pencil of function (GPOF) method is utilized to realize DCIM [39]. Then, spatial domain Green’s functions can be obtained in closed forms from their spectral-domain counterparts via the Sommerfeld identity.

2.2. MLACA Acceleration

The MLACA employs the same octree data structure as in the MLFMA [7]. The octal-tree algorithm is used to subdivide a box that encloses an object into smaller boxes. Figure 1 shows the decomposition of the problem domain at different levels. With reference to Figure 1, far interactions exit at levels 2 and higher. Far interactions can be computed using the MLACA [40].

According to the above octree data structure, the impedance matrix in (14) can be decomposed into near-field interactions and far-field interactions . It can be simply written as where the and are computed directly and compressed by the MLACA algorithm, respectively. When the two boxes are sufficiently separated, the impedance matrix associated with them can be expressed using low-rank representations. This feature is utilized by the MLACA. In the MLACA implementation, the impedance matrix between two sufficiently distant boxes can be expressed in terms of two small matrices where the is the interaction matrix between the observation and source boxes. and denote the number of the basis functions in the observation and source boxes, respectively. The index denotes the rank of and is much smaller than and . The impedance matrix in (16) can be rewritten in multilevel operations as where the is the number of nonempty groups at level and denotes the number of far interaction groups of the th nonempty group for each observer group at level . Since the matrices and generated by the MLACA are usually not orthogonal, they may contain redundancies that can be removed by the SVD technique [41].

Finally, the linear system of equations in (14) can be solved by preconditioned restated generalized minimal residual algorithm (GMRES) using the MLACA to accelerate the matrix-vector multiplies in iterative algorithm.

3. Numerical Experiments

In this section, we show some numerical results for electromagnetic structures that illustrate the effectiveness of the proposed A-EFIE with CNE accelerated by the MLACA algorithm. All numerical experiments are performed on Intel(R) Core (TM) 2 Quad CPU at 2.83 GHz and 3.5 GB of RAM in double precision and the truncating tolerance of the ACA is (relative to the largest singular value). The restart number of the generalized minimal residual (GMRES) is set to be 30. The inner-outer flexible GMRES algorithm (FGMRES) is used as the iterative solver for microstrip circuit problems.

Additional details and comments on the implementation are given below:(i)zero vector is taken as initial approximate solution for all examples and all systems in each example;(ii)the iteration process is terminated when the normalized backward error is reduced by for all the examples.

First, electromagnetic scattering from a perfect electric conducting (PEC) sphere is analyzed in order to demonstrate the advantages of the proposed A-EFIE with the ACA algorithm (denoted as A-EFIE-ACA). The sphere has a radius of 2 m. It is discretized with 7332 triangular patches leading to 18330 unknowns. The angle of incident plane wave is , and with the frequency  KHz. Figure 2 shows the VV-polarization bistatic RCS solved by Mie series and A-EFIE-ACA. It can be found that the A-EFIE-ACA gives right results at low frequency. The condition numbers versus frequency are also demonstrated in Figure 3. It is observed from Figure 3 that the A-EFIE with charge neutrality enforcement makes the condition number almost constant when decreasing the frequency to 1 Hz. The condition number grows slowly when increasing the frequency from 1 GHz to 10 GHz. Figure 4 further shows the VV-polarization bistatic RCS when the number of unknowns for the above PEC sphere becomes larger. It can be found that the proposed A-EFIE-ACA works well when the number of unknowns increases from 19435 to 65130. However, when the number of unknowns is increased to 76480, the procedure for A-EFIE-ACA cannot work due the dense-grid breakdown problem. Table 1 shows the comparison of memory requirement and solution time of GMRES between the A-EFIE-ACA and A-EFIE without the ACA algorithm (denoted as A-EFIE-MoM) in bistatic RCS computation for the sphere. The 2-level ACA is used to reduce the solution time and memory requirement in A-EFIE-ACA. It can be found that the A-EFIE-ACA can save much time by a factor of 7.9 than that of the A-EFIE-MoM. The A-EFIE-ACA can also save much solution time by a factor of 3.0 than that of the A-EFIE-MoM.

Next, we investigate the performance of the proposed A-EFIE-ACA on electromagnetic scattering at high frequency. The mathematic description of NASA almond is defined in [42] and shown in Figure 5. The total length of the almond is 9.936 inches. The almond is discretized with 5216 triangular patches with 13040 unknowns. The angle of incident plane wave is , and . Figure 5 gives the VV-polarization bistatic RCS for an almond at 3 GHz. It can be found that the results obtained by the A-EFIE-ACA agree well with ones from traditional EFIE with ACA algorithm (denoted as EFIE-ACA). The 2-level ACA is used in the A-EFIE-ACA and EFIE-ACA. Table 2 shows the comparison of memory requirement and solution time of GMRES among different methods in bistatic RCS computation for the almond. It can be found that the A-EFIE-ACA can save much memory requirement by a factor of 2.1 than that of the A-EFIE-MoM, 6.6 than that of the EFIE- ACA. It can be also found that the A-EFIE-ACA can save much memory requirement by a factor of 2.9 than that of the A-EFIE-MoM. However, the memory requirement for the A-EFIE-ACA is more than that for the EFIE-ACA due to additional charge unknowns introduced into the A-EFIE.

The third example is a microstrip bandpass filter (BPF) based on five intercoupled split-ring resonators (SRR) proposed in [43]. The periodic SRR structure is shown in Figure 6. The distance between each SRR’s arms is 0.05 mm. The relative dielectric constant of substrate is 3.9. The structure is implemented on the substrate Corning 1059 with a relative dielectric constant of 3.9 and thickness of 0.37 mm. Microstrip line width is 0.8 mm. The width of the interdigital and the gap in each SRR structure are both 0.114 mm. The length of each SRR structure is 0.5 wavelengths at center frequency of 12 GHz. The structure is discretized with 924 triangular patches with 1077 unknowns. Frequency is from 5 GHz to 20 GHz. The sampled frequency interval is 0.1 GHz and a total of 151 frequencies need to be calculated in the proposed A-EFIE-ACA. Figure 7 gives scattering parameters S11 and S12 obtained from the A-EFIE-ACA method for the above BPF structure. It can be found that the results from A-EFIE-ACA agree well with those from Ansoft Designer software. It is observed from Figure 6 that the BPF exhibits a passband region from 9.0 to 15.2 GHz, which corresponds to a fractional bandwidth of 49.6%. The insertion losses in the passband region are around 1.2 dB. Therefore, the simulated results confirm the feasibility of utilizing the intercoupled split-ring resonator structure in microwave bandpass filter applications. The BPF is a typical electrical small structure when the wavelength is 0.015 m at 20 GHz. The condition number of impedance matrix is very bad due to the fine structure in the BPF with SRR structure. Therefore, the inner-outer Flexible GMRES algorithm (FGMRES) is used as the iterative solver. In FGMRES, the near-field impedance matrix is used as the preconditioner in the inner iteration. In the FGMRES algorithm, the stop precision is 1.E-4. Table 3 shows the comparison of memory requirement and solution time among the A-EFIE-ACA, A-EFIE-MoM, and EFIE-ACA methods for the BPF with SRR structure. It can be found that the proposed A-EFIE-ACA can save much memory requirement by a factor of 2.54 than that of the A-EFIE-MoM. The A-EFIE-ACA can also save much time for average solution time at each frequency by a factor of 1.48 than that of A-EFIE-MOM, by a factor of 39.4 than that of EFIE-ACA. Therefore, the A-EFIE-ACA can save much time for the total iterative solution time by a factor of 39.5 than that of the EFIE-ACA.

The last example is a dual-mode dual-band BPF [44] as shown in Figure 8. The dual-mode dual-band BPF is implemented on the substrate Duroid 6010 with a relative dielectric constant of 10.2, loss tangent of 0.0023, and thickness of 1.27 mm. The related dimensions as shown in Figure 6 are determined as follows:  mm,  mm,  mm,  mm,  mm,  mm,  mm,  mm, and  mm. The horizontal and vertical lengths of the slot are 9.2 mm and 6.2 mm, respectively. The size of the patch including stubs is , where is 48 mm, the guided wavelength at the lower passband. The simulated results obtained from the A-EFIE-ACA method are given in Figure 9 compared to the results from Ansoft HFSS software. The BPF structure is discretized with 1476 triangular patches with 2070 unknowns. Frequency is from 1.5 GHz to 6.0 GHz. The interval of frequency is 0.1 GHz. Good agreement is achieved between these two kinds of results. It is observed from Figure 8 that the lower and upper passbands have the fractional bandwidth of 12.8% and 6.5% centered at 2.3 and 5.05 GHz, respectively. The simulated insertion loss is 0.73 dB and 0.86 dB, respectively. The return loss is better than 10 dB. Two transmission zeros are generated near the passband edges, guaranteeing high selectivity. It also can be found that two transmission poles in both passbands demonstrate the effectiveness of the dual-mode dual-band filter. The 3-level ACA is used. The SVD is used to recompress the far interaction matrix in the A-EFIE-ACA method. Table 4 shows the comparison of memory requirement and solution time among the A-EFIE-ACA, A-EFIE-MoM, and EFIE-ACA methods for the BPF structure. It can be found that the proposed A-EFIE-ACA can save much memory requirement by a factor of 2.44 than that of the A-EFIE-MoM. The average solution time at each frequency for the A-EFIE-MOM is about 1.34 times than that of of the A-EFIE-ACA. Furthermore, the A-EFIE-ACA can save much time by a factor of 31.3 than that of the EFIE-ACA. It can be concluded that the proposed A-EFIE-ACA is very efficient compared to other methods.

4. Conclusions

In this paper, the A-EFIE with CNE accelerated by the MLACA algorithm is presented for solving electromagnetic scattering and microstrip circuit problems. Numerical experiments are performed and the comparison is made with the EFIE-ACA and A-EFIE-MoM methods. It can be found that the proposed A-EFIE-ACA method is more efficient and can significantly reduce the overall computational cost. Due to the fact that the set-up time of the ACA is much more than that of MLFMA, it is known that the ACA does not maintain its favorable complexity once the problem becomes electrically large. The PILOT algorithm can be used to reduce the set-up time and the memory consumption of ACA. The PILOT algorithm utilizes the idea that the far-field interaction lists of siblings share many common cubes to regroup a new far-field interaction list based on an octal tree data structure. Higher compression is achieved by using the PILOT algorithm when the dimension of the matrix is large. Further investigations deserve to be undertaken to study efficiency of the PILOT algorithm when the proposed A-EFIE-ACA method is used for the analysis of electromagnetic problems from electrically large objects.

Acknowledgments

The authors would like to thank Jiangsu Natural Science Foundation (BK2012034) and Natural Science Foundation (61171041) for support.