#### Abstract

An improved generalized single-source tangential equivalence principle algorithm (GSST-EPA) is proposed for analyzing array structures with connected elements. In order to use the advantages of GSST-EPA, the connected array elements are decomposed and computed by a contact-region modeling (CRM) method, which makes that each element has the same meshes. The unknowns of elements can be transferred onto the equivalence surfaces by GSST-EPA. The scattering matrix in GSST-EPA needs to be solved and stored only once due to the same meshes for each element. The shift invariant of translation matrices is also used to reduce the computation of near-field interaction. Furthermore, the multilevel fast multipole algorithm (MLFMA) is used to accelerate the matrix-vector multiplication in the GSST-EPA. Numerical results are shown to demonstrate the accuracy and efficiency of the proposed method.

#### 1. Introduction

Surface integral equation (SIE) methods have been widely used in the simulation of complex electromagnetic phenomena in practical engineering. Application of SIE methods always leads to a dense matrix which limits its applicability. Many fast solvers have been developed to reduce the memory requirement and CPU time, such as fast multipole methods [1, 2], FFT-based methods [3, 4], low-rank-based methods [5–8], and direct solvers [9, 10]. However, for many multiscale electromagnetic problems, the dimension of the matrix is large and the matrix is ill-conditioned, which are still the challenges for the fast solvers. It is desirable for more efficient and robust SIE solutions to deal with those challenges.

In recent years, domain decomposition (DD) methods have been attracting much attention and regarded as effective techniques to solve multiscale systems. The DD methods rely on decomposing the problem into several smaller subdomains and solving each subdomain independently. In this manner, the dimension of the matrix can be reduced and the condition number of the matrix can be improved. Due to these advantages, the DD methods have been combined with finite element (FE) methods [11, 12] and finite difference (FD) methods [13, 14] successfully. To extend the DD methods to SIE methods, many studies have been proposed. Among these studies, several works are worth mentioning. The overlapping IE-DDM [15] needs to construct buffer regions. The IE-DDM in [16] uses artificial interfaces to do the decomposition. A nonoverlapping IE-DDM [17] has been proposed by using the explicit boundary condition. The discontinuous Galerkin integral equation (IEDG) [18] has been formulated by employing an interior penalty term, which shows great flexibility to geometries. The equivalence principle algorithm (EPA) [19] is a kind of DD method based on Huygens’ principle. By using the virtual equivalence surfaces that enclose the subdomains, the EPA can isolate the subdomain from others. Then, the original unknowns of subdomains are transferred to the unknowns on the equivalence surfaces, which is different from the above methods. The linear embedding via Green’s operator (LEGO) [20] and generalized transition matrix (GTM) algorithm [21] are similar with EPA. However, the connected array structures have not been dealt with in those methods.

In this paper, the improved generalized single-source tangential equivalence principle algorithm (GSST-EPA) is proposed to address the array structures with connected elements. Although the GSST-EPA has been developed to solve the connected structures [22], it is inefficient for the array structures. For the array structures with connected elements, if each element is decomposed as a subdomain and enclosed by the same equivalence surfaces, then the meshes of elements are different from each other. Therefore, the scattering matrix of each equivalence surface is different and needs to be solved and stored one by one, which leads to huge computational time and memory requirement. To overcome the above problem and make full use of the GSST-EPA, the contact-region modeling (CRM) [23] method is applied to make the decomposition for the array structures. Unlike many other DD methods, the CRM method automatically contains the boundary conditions for the contact surface between subdomains, so the additional transmission conditions are unnecessary. There are two major benefits of the proposed scheme: (a) Each element has the same meshes, so the scattering matrix of each equivalence surface is the same with each other and needs to be solved and stored only once, and (b) the shift invariant of translation matrices can be applied to reduce the computation of coupling between neighbor elements. Finally, the proposed method is accelerated by the multilevel fast multipole algorithm (MLFMA) [2].

#### 2. The Domain Decomposition Based on the CRM Method

Without loss of generality, consider the scattering problem of a three-dimensional dielectric domain with the media , as illustrated in Figure 1(a). The exterior region is the free space. The dielectric structure is illuminated by an incident plane wave . This problem can be solved by the commonly used PMCHWT formulation [24]. In this work, the CRM method is used for the domain partitioning. As shown in Figure 1(b), the dielectric domain is decomposed into two nonoverlapping subdomains such that . The surfaces of the subdomains and with outward-directed unit normal and are denoted by and , respectively. As the distance , the solution of Figure 1(b) is the same as that of Figure 1(a) due to the same physics. Next, the variational PMCHWT formulation based on the CRM method is derived.

By using the extinction theorem, the total fields on the inner boundary of subdomain are equal to zero in the exterior equivalence situation. and are the equivalent electric and magnetic currents defined on the outer side of the surface . The EFIE and MFIE for subdomain can be expressed as where is the wave impedance in free space. denotes the inner side of . The operators and are defined as where is or , denotes the media domain, and and are the field surface and source surface, respectively. stands for the wave number in domain . is Green’s function defined as

PV stands for the Cauchy principle value integration. For the residue term , if the observation point resides on the outer side of , then “+” is chosen. Otherwise, “−” is taken if resides on the inner side of , while if is not on , the residue term should be equal to zero. Similarly, the EFIE and MFIE for subdomain in the exterior equivalence situation can be written as

For the interior equivalence situation, the total fields are zero on the outer boundary of . So the EFIE and MFIE for the are written as where

denotes the outer side of . In the domain , the EFIE and MFIE are

Combining (1) and (8), the PMCHWT E-field equation for surface can be derived as

The H-field equation can be obtained by combining (2) and (9) as

Similarly, the PMCHWT E-field and H-field equations for surface can be derived by combining (6) and (11) and (7) and (12), as follows:

Finally, the matrix form of these four PMCHWT equations is written as where is used to balance the magnitude of operators and make the matrix better conditioned.

#### 3. Boundary Conditions at the Contact Surface in CRM

It is worth mentioning that the additional transmission conditions are unnecessary for the contact surface, because the continuity of the currents at the contact surface is contained in the above PMCHWT equations automatically. To prove this fact, a simple model is considered as shown in Figure 2. The surfaces of domains and are divided into four parts as and . In the exterior equivalence situation, the MFIE on the contact surfaces and are where

Because the surfaces and touch each other, the outer side of is the same as that of . Therefore, the terms , , , , and in (17) are equal to the terms , , , , and in (18). Since the operator does not have the residue term, the terms and are the same as and . By subtracting (18) from (17), the following can be obtained:

The definition of the operator can be found in (4). It can be seen that the Cauchy principle value integration in is the same as that in , while the residue terms are different. The terms and have similar relations. Therefore, (20) can be rewritten as

Since , then (21) can be derived as

It can be seen that the continuity of the electric currents at the contact surface is included in the MFIE. Similarly, the continuity of the magnetic currents at the contact surface can be derived from the EFIE equations.

#### 4. GSST-EPA with the CRM Method

Consider the scattering of a dielectric object in free space. As shown in Figure 3, the object is decomposed into -connected subdomains using the CRM method. Then, each subdomain is enclosed by an equivalence surface. The basic idea of GSST-EPA is to transform the original scattering problem into a new equivalent problem with a scattering operator and a translation operator. Next, the definitions of the scattering and translation operators are illustrated in the following.

The scattering operator can be derived from a dielectric object enclosed by an equivalence surface with three steps as shown in Figure 4. The three steps are outside-in propagation, solving for the currents on the object, and inside-out propagation. The equivalence current on the equivalence surface generates the same scattered field as that generated by original currents on the object. Therefore, the scattering of a dielectric object can be characterized by the scattering operator as
where subscripts and denote the equivalence surface and dielectric object separately. *TEH* means one type of the expressions between the incident current and the scattering current , which can be found in [25]. The dielectric object is solved by PMCHWT equations. The scattering operator can capture the full physics of the object, so the inside object can be removed.

The coupling between objects via the equivalence surfaces can be described by a translation operator. The translation operator has two different forms: one is suitable for nonadjacent cases, and the other is suitable for adjacent cases. For the two nonadjacent subdomains, the translation operator is quite simple. As shown in Figure 5, the subdomains and are enclosed by equivalence surfaces ES and ES separately. The coupling from ES to ES is the fields produced by the equivalent scattered current on ES , which can be expressed via a translation operator , , as

Therefore, the equivalence scattered current on ES is given by

The expression of the translation operator is similar as given above.

For the two adjacent subdomains, the translation operator can be derived by the source reconstruction method (SRM) [26–28]. As shown in Figure 6, the inside object has been decomposed into two subdomains with the distance . The coupling from ES to ES can be derived with several steps. Firstly, suppose the on ES is already known, then the equivalent electric current and magnetic current on the surface of dielectric can be obtained by using the SRM as where is the pseudoinverse. Secondly, once the equivalent currents on the surface of dielectric is obtained, the coupling from the dielectric to the dielectric can be solved using the CRM method as

Thirdly, the equivalent current produced by the incident fields and using the equivalence principle can be written as

Finally, by combing (27), (28), and (29), the equivalence scattered current is equal to where the translation operator is defined as the multiplication of several operators as

Since the scattering operator and translation operator have been derived, the GSST-EPA equations for the object in Figure 3 can be written as where

is the identity operator. The matrix form can be obtained by expanding the currents with the Rao-Wilton-Glisson (RWG) basis function and using Galerkin’s method to discretize the equations.

#### 5. Computational Considerations

In this paper, the periodic array structures will be solved by the GSST-EPA. If the array elements are connected with each other, then the CRM method can be used to decompose the elements. Next, each element is enclosed by an equivalence surface. Due to the same geometry and meshes of each element and equivalence surface, the scattering operator is also the same with each other. Therefore, the scattering matrix needs to be solved and stored only once, which can reduce the computational time and memory requirement.

Since the array structures have a lot of elements, the MLFMA is used to accelerate the matrix-vector multiplication in the GSST-EPA. It is notable that the translation operators have two different forms. The operator denotes the coupling between nonadjacent elements which can be calculated by MLFMA, while the operator denotes the near-field coupling between adjacent elements that should be solved directly. Due to characteristic of the periodic array structure, the translation matrices are shift invariant which can be used to reduce the computation of near-field coupling. As shown in Figure 7, it is easy to find that the eight kinds of near-interaction mechanism between nine elements is a typical representation of different near interactions between all elements. In fact, the four corner elements in the array only have three kinds of near interactions. The edge elements have five kinds of near interactions. The other elements have eight kinds of near interactions. According to the shift invariant, these eight kinds of can be solved and stored only once, which can be used to reduce the calculation of any adjacent interactions.

#### 6. Numerical Results

In this section, the performance of the proposed method is studied via several numerical experiments. All of the simulations are performed on PC with Intel Core i7-2760 2.4 GHz CPU and 32.0 GB memory.

Firstly, the accuracy of the CRM is investigated. The scattering of a dielectric sphere with and in free space is considered as shown in Figure 8(a). The sphere is partitioned into two subdomains (two hemispheres) with a fictitious contact surface by CRM as shown in Figure 8(b). The excitation is an -polarized plane wave propagating into the negative direction at 300 MHz. The curvilinear Rao-Wilton-Glisson (CRWG) [29] is used to discretize the sphere surface. The bistatic radar cross section (RCS) of HH polarization is given in Figure 9. It can be seen that the result of CRM agrees very well with Mie series solution, which also proves that the decomposition does not change the physics of the original problem.

**(a)**

**(b)**

The second example is used to investigate the accuracy and convergence of the GSST-EPA with CRM. As shown in Figure 10, the dimension of the dielectric slab is . The relative permittivity and permeability of the slab are and . The excitation is an -polarized plane wave propagating into the negative direction at 300 MHz. The slab is decomposed into , , and subdomains separately by CRM as shown in Figure 11. The subdomains have the same size as each other. Then, each subdomain is enclosed by an equivalence surface with the same size. The sizes of each equivalence surface are , , and for the three different decompositions separately, which means the distance between the subdomain and the equivalence surface is . Since each subdomain and each equivalence surface have the same size, the same meshes for each subdomain and each equivalence surface can be constructed. Therefore, the scattering matrix in GSST-EPA only needs to be calculated only once. The bistatic RCS of HH polarization is given in Figure 12, which shows good accuracy of GSST-EPA with CRM using different decomposition schemes. Figure 13 shows the comparison of the iteration number among the three different decompositions in GSST-EPA. The generalized minimal residual (GMRES) method with the tolerance 10^{−6} is used to solve the final equations. With the increase in subdomains, the iteration number only grows from 20 to 25, which shows the convergence stability of the method, while the PMCHWT equations converge to 10^{−6} with 1219 steps.

**(a)**

**(b)**

**(c)**

The third example is a bowtie antenna array located in the *x*-*y* plane. The dimension of a bowtie antenna element is shown in Figure 14(a). The yellow part is the PEC patch, while the red part is the dielectric substrate. The thickness of the antenna substrate is 0.03 m with and . The elements are connected with each other. The array is divided into 9 subdomains by CRM, which means that each element is a subdomain. Then, each element is enclosed by an equivalence surface as shown in Figure 14(b). The excitation is an -polarized plane wave propagating into the negative direction at 1.0 GHz. Because the element is a composite target, the EFIE-PMCHWT is used to solve the element. As shown in Figure 15, the bistatic RCS of GSST-EPA agrees very well with that of EFIE-PMCHWT. The comparison of computational results is given in Table 1. Because the number of unknowns on the equivalence surfaces is much less than that on the element surfaces, GSST-EPA can reduce the solution time and memory requirement significantly. It is well known that EFIE-PMCHWT is accurate but suffers from poor spectrum property. By transferring the unknowns on the element surfaces to the unknowns on the equivalence surfaces, GSST-EPA with CRM can improve the spectrum property of EFIE-PMCHWT significantly.

**(a)**

**(b)**

Finally, to investigate the capability of GSST-EPA with CRM to solve large array structures, the scattering of the bowtie antenna array is solved. The parameters of the array are the same as that in the above example, except for the number of elements. The total number of unknowns of the equivalence surfaces is 41,850. The MLFMA is applied to accelerate the calculation of the coupling between nonadjacent elements. The bistatic RCS of HH polarization is given in Figure 16. The total solution time is 1972.6 s, and memory consumption is 3.9 GB. GSST-EPA converges to 10^{−3} using 74 iteration steps. It can be seen that GSST-EPA still exhibits fast convergence as the increment of the element number.

#### 7. Conclusions

In this paper, GSST-EPA with the CRM method is introduced to solve the array structures with connected elements. Firstly, the array structures are decomposed into subdomains by CRM. Then, the unknowns on the subdomains are transferred to the unknowns on the equivalence surfaces by GSST-EPA. Moreover, the meshes of subdomains can be the same with each other, which makes the scattering matrix solved and stored only once. The shift invariant of translation matrices is also used according to the characteristic of array structures. Using this strategy, the number of unknowns can be reduced and the conditioning of the final matrix equation is significantly better than the traditional PMCHWT and EFIE-PMCHWT. Therefore, the improved GSST-EPA can be viewed as an efficient and robust solver to solve array structures.

#### Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

#### Acknowledgments

This work is supported by National Natural Science Foundation of China (Grant nos. 61501267, 61631012, 61671257, and 61425010), Natural Science Foundation of Zhejiang Province (Grant no. LQ16F010003), and Ningbo University Discipline Project (Grant nos. XKL14D2058 and XYL15008) and is partially sponsored by a K.C. Wong Magna Fund in Ningbo University.