# Improved Generalized Single-Source Tangential Equivalence Principle Algorithm with Contact-Region Modeling Method for Array Structures

Authors: Yao Han; Hanru Shao; Jianfeng Dong
Journal: International Journal of Antennas and Propagation
 (2018)

---

## 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.


---

## Body

## 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ε1,μ1, as illustrated in Figure 1(a). The exterior region Ω0 is the free space. The dielectric structure is illuminated by an incident plane wave Einc,Hinc. 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 Ω=Ω1∪Ω2. The surfaces of the subdomains Ω1 and Ω2 with outward-directed unit normal n^1 and n^2 are denoted by S1 and S2, respectively. As the distance δ→0, 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.Figure 1
(a) The scattering problem of a three-dimensional dielectric object. (b) The object decomposed into two subdomains using the contact-region modeling method.By using the extinction theorem, the total fields on the inner boundary of subdomainΩ1 are equal to zero in the exterior equivalence situation. J1 and M1 are the equivalent electric and magnetic currents defined on the outer side of the surface S1. The EFIE and MFIE for subdomain Ω1 can be expressed as
(1)η0ℒ110J1−K~110−M1+η0ℒ120J2−K~120+M2=−E1inc,r∈S1−,(2)K~110−J1+1η0ℒ110M1+K~120+J2+1η0ℒ120M2=−H1inc,r∈S1−,where η0 is the wave impedance in free space. S1− denotes the inner side of S1. The operators ℒ and K~ are defined as
(3)ℒjkir,r′·Xr′=iki∫SkXr′+1ki2∇∇′·Xr′·Gir,r′ds′,r∈Sj,(4)K~jki±r,r′·Xr′=±12Xr×n^k+PV∫Sk∇Gir,r′×Xr′ds′,r∈Sj,where X is J or M, i denotes the media domain, and Sj and Sk are the field surface and source surface, respectively. ki stands for the wave number in domain i. Gi is Green’s function defined as
(5)Gir,r′=eikir−r′4πr−r′.PV stands for the Cauchy principle value integration. For the residue term±1/2Xr×n^k, if the observation point r resides on the outer side of Sk, then “+” is chosen. Otherwise, “−” is taken if r resides on the inner side of Sk, while if r is not on Sk, the residue term should be equal to zero. Similarly, the EFIE and MFIE for subdomain Ω2 in the exterior equivalence situation can be written as
(6)η0ℒ210J1−K~210+M1+η0ℒ220J2−K~220−M2=−E2inc,r∈S2−,(7)K~210+J1+1η0ℒ210M1+K~220−J2+1η0ℒ220M2=−H2inc,r∈S2−.For the interior equivalence situation, the total fields are zero on the outer boundary ofΩ1. So the EFIE and MFIE for the Ω1 are written as
(8)η1ℒ111−J1−K~111+−M1=0,r∈S1+,(9)K~111+−J1+1η1ℒ111−M1=0,r∈S1+,where
(10)η1=μ1ε1.S1+ denotes the outer side of S1. In the domain Ω2, the EFIE and MFIE are
(11)η1ℒ221−J2−K~221+−M2=0,r∈S2+,(12)K~221+−J2+1η1ℒ221−M2=0,r∈S2+.Combining (1) and (8), the PMCHWT E-field equation for surface S1 can be derived as
(13)η0ℒ110+η1ℒ111J1−K~110−+K~111+M1+η0ℒ120J2−K~120+M2=−E1inc,r∈S1.The H-field equation can be obtained by combining (2) and (9) as
(14)K~110−+K~111+J1+1η0ℒ110+1η1ℒ111M1+K~120+J2+1η0ℒ120M2=−H1inc,r∈S1.Similarly, the PMCHWT E-field and H-field equations for surfaceS2 can be derived by combining (6) and (11) and (7) and (12), as follows:
(15)η0ℒ210J1−K~210+M1+η0ℒ220+η1ℒ221J2−K~220−+K~221+M2=−E2inc,r∈S2,K~210+J1+1η0ℒ210M1+K~220−+K~221+J2+1η0ℒ220+1η1ℒ221M2=−H2inc,r∈S2.Finally, the matrix form of these four PMCHWT equations is written as(16)η0ℒ110+η1ℒ111−η0K~110−+K~111+η0ℒ120−η0K~120+η0K~110−+K~111+η0ℒ110+η02η1ℒ111η0K~120+η0ℒ120η0ℒ210−η0K~210+η0ℒ220+η1ℒ221−η0K~220−+K~221+η0K~210+η0ℒ210η0K~220−+K~221+η0ℒ220+η02η1ℒ221·J11η0M1J21η0M2=−E1inc−η0H1inc−E2inc−η0H2inc,where η0 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 Figure2. The surfaces of domains Ω1 and Ω2 are divided into four parts as S1=Sd1∪Sc1 and S2=Sd2∪Sc2. In the exterior equivalence situation, the MFIE on the contact surfaces Sc1 and Sc2 are
(17)KSc1Sd10Jd1+1η0ℒSc1Sd10Md1+KSc1Sd20Jd2+1η0ℒSc1Sd20Md2+K~Sc1Sc10−Jc1+1η0ℒSc1Sc10Mc1+K~Sc1Sc20+Jc2+1η0ℒSc1Sc20Mc2=−HSc1inc,r∈Sc1−,(18)KSc2Sd10Jd1+1η0ℒSc2Sd10Md1+KSc2Sd20Jd2+1η0ℒSc2Sd20Md2+K~Sc2Sc10+Jc1+1η0ℒSc2Sc10Mc1+K~Sc2Sc20−Jc2+1η0ℒSc2Sc20Mc2=−HSc2inc,r∈Sc2−,where
(19)Kjkir,r′·Xr′=PV∫Sk∇Gir,r′×Xr′ds′,r∈Sj.Figure 2
The model for the continuity of the currents at the contact surface.Because the surfacesSc1 and Sc2 touch each other, the outer side of Sc1 is the same as that of Sc2. Therefore, the terms 1/η0ℒSc1Sd10Md1, 1/η0ℒSc1Sd20Md2, 1/η0ℒSc1Sc10Mc1, 1/η0ℒSc1Sc20Mc2, and HSc1inc in (17) are equal to the terms 1/η0ℒSc2Sd10Md1, 1/η0ℒSc2Sd20Md2, 1/η0ℒSc2Sc10Mc1, 1/η0ℒSc2Sc20Mc2, and HSc2inc in (18). Since the K operator does not have the residue term, the terms KSc1Sd10Jd1 and KSc1Sd20Jd2 are the same as KSc2Sd10Jd1 and KSc2Sd20Jd2. By subtracting (18) from (17), the following can be obtained:
(20)K~Sc1Sc10−Jc1−K~Sc2Sc10+Jc1+K~Sc1Sc20+Jc2−K~Sc2Sc20−Jc2=0.The definition of theK~ operator can be found in (4). It can be seen that the Cauchy principle value integration in K~Sc1Sc10− is the same as that in K~Sc2Sc10+, while the residue terms are different. The terms K~Sc1Sc20+ and K~Sc2Sc20− have similar relations. Therefore, (20) can be rewritten as
(21)−Jc1×n^c1+Jc2×n^c2=0.Sincen^c1=−n^c2, then (21) can be derived as
(22)Jc1=−Jc2.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(23)Mc1=−Mc2can 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 Figure3, the object is decomposed into N-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.Figure 3
The dielectric object decomposed into subdomains and each subdomain enclosed by an equivalence surface.The scattering operator can be derived from a dielectric object enclosed by an equivalence surface with three steps as shown in Figure4. 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
(24)Jls=η0ℒll+η0K~ll+ℒll−1K~ll−−1·η0ℒld−η0Kld·η0ℒdd0+η1ℒdd1−η0K~dd0−+K~dd1+η0K~dd0−+K~dd1+η0ℒdd0+η02η1ℒdd1−1·−η0ℒdlη0Kdl−η0Kdl−η0ℒdl·1−ℒll−1K~ll+·Jli=S~llTEH·Jli,where subscripts l and d denote the equivalence surface and dielectric object separately. TEH means one type of the expressions between the incident current Jli and the scattering current Jls, 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.Figure 4
A dielectric object enclosed by an equivalence surface with three steps.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 Figure5, the subdomains l and k are enclosed by equivalence surfaces ES l and ES k separately. The coupling from ES k to ES l is the fields produced by the equivalent scattered current Jks on ES k, which can be expressed via a translation operator T~lk, l≠k, as
(25)T~lkTEH=η0ℒll+η0K~ll−ℒll−1K~ll+−1·η0ℒlk−η0Klk·1−ℒll−1K~ll−.Figure 5
The translation operator of two nonadjacent subdomains.Therefore, the equivalence scattered currentJks on ES l is given by
(26)Jls=S~llTEH·Jli+S~llTEH·T~lkTEH·Jks.The expression of the translation operatorT~kl is similar as T~lk 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 δ→0. The coupling from ES k to ES l can be derived with several steps. Firstly, suppose the Jks on ES k is already known, then the equivalent electric current Jdk and magnetic current Mdk on the surface of dielectric k can be obtained by using the SRM as
(27)Jdk1η0Mdk=η0ℒkdk−η0Kkdk†·η0ℒkk+η0K~kk+ℒkk−1K~kk−·Jks,where ·† is the pseudoinverse. Secondly, once the equivalent currents on the surface of dielectric k is obtained, the coupling from the dielectric k to the dielectric l can be solved using the CRM method as
(28)−Edldki−1η0Hdldki=−η0ℒdldk0η0K~dldk0+−η0K~dldk0+−η0ℒdldk0·Jdk1η0Mdk.Figure 6
The translation operator of two adjacent subdomains.Thirdly, the equivalent current produced by the incident fieldsEdldki and Hdldki using the equivalence principle can be written as
(29)J~ls=η0ℒll+η0K~ll+ℒll−1K~ll−−1·η0ℒldl−η0Kldl·η0ℒdldl0+η1ℒdldl1−η0K~dldl0−+K~dldl1+η0K~dldl0−+K~dldl1+η0ℒdldl0+η02η1ℒdldl1−1·−Edldki−1η0Hdldki.Finally, by combing (27), (28), and (29), the equivalence scattered current Jls is equal to
(30)Jls=S~llTEH·Jli+J~ls=S~llTEH·Jli+T~lkSRM·Jks,where the translation operator T~lkSRM is defined as the multiplication of several operators as
(31)T~lkSRM=η0ℒll+η0K~ll+ℒll−1K~ll−−1·η0ℒldl−η0Kldl·η0ℒdldl0+η1ℒdldl1−η0K~dldl0−+K~dldl1+η0K~dldl0−+K~dldl1+η0ℒdldl0+η02η1ℒdldl1−1.−η0ℒdldk0η0K~dldk0+−η0K~dldk0+−η0ℒdldk0·η0ℒkdk−η0Kkdk†·η0ℒkk+η0K~kk+ℒkk−1K~kk−.Since the scattering operator and translation operator have been derived, the GSST-EPA equations for the object in Figure3 can be written as
(32)ℐ−Z12⋯−Z1N−Z21ℐ⋯−Z2N⋮⋮⋱⋮−ZN1−ZN2⋯ℐ·J1sJ2s⋮JNs=S~11TEH·J1iS~22TEH·J2i⋮S~NNTEH·JNi,where
(33)Zij=S~iiTEHT~ijTEH,ifiandjarenonadjacent,T~ijSRM,ifiandjareadjacent.ℐ 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 operatorS~iiTEH 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. TheT~ijTEH operator denotes the coupling between nonadjacent elements which can be calculated by MLFMA, while the T~ijSRM 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 T~ijSRM can be solved and stored only once, which can be used to reduce the calculation of any adjacent interactions.Figure 7
Eight kinds of near-interaction mechanism between nine elements.


## 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εr=2.0 and μr=1.0 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 x^-polarized plane wave propagating into the negative z^ 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.Figure 8
(a) The model of a dielectric sphere. (b) Subdomain partitions of the dielectric sphere.
(a)
(b)Figure 9
The bistatic RCS of a dielectric sphere at 300 MHz, HH polarization.The second example is used to investigate the accuracy and convergence of the GSST-EPA with CRM. As shown in Figure10, the dimension of the dielectric slab is 1.5×1.5×0.3m3. The relative permittivity and permeability of the slab are εr=2.0 and μr=1.0. The excitation is an x^-polarized plane wave propagating into the negative z^ direction at 300 MHz. The slab is decomposed into 3×3, 4×4, and 5×5 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 0.7×0.7×0.5m3, 0.6×0.6×0.5m3, and 0.5×0.5×0.5m3 for the three different decompositions separately, which means the distance between the subdomain and the equivalence surface is 0.1λ. 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.Figure 10
The model of a dielectric slab.Figure 11
The dielectric slab decomposed into subdomains and each subdomain enclosed by an equivalence surface: (a)3×3 subdomains, (b) 4×4 subdomains, and (c) 5×5 subdomains.
(a)
(b)
(c)Figure 12
The bistatic RCS of a dielectric slab at 300 MHz, HH polarization.Figure 13
The convergence of GSST-EPA for the dielectric slab.The third example is a3×3 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 εr=2.2 and μr=1.0. 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 x^-polarized plane wave propagating into the negative z^ 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.Figure 14
(a) The dimension of a bowtie antenna element and (b)3×3 bowtie antenna array.
(a)
(b)Figure 15
The bistatic RCS of a3×3 bowtie antenna array at 1.0 GHz, HH polarization.Table 1
The computational results of two methods.
Method
Unknowns
Memory (MB)
Iteration number (0.001)
Total time (s)
EFIE-PMCHWT
13,404
1747.9
1177
894.5
GSST-EPA
1674
53.3
10
442.6Finally, to investigate the capability of GSST-EPA with CRM to solve large array structures, the scattering of the15×15 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.Figure 16
The bistatic RCS of a15×15 bowtie antenna array at 1.0 GHz, HH polarization.


## 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.



---

*Source: 9875041-2018-06-27.xml*