Vector Single-Source SIE for TE Scattering From Objects Embedded in Multilayers
11 Vector Single-Source SIE for TE Scattering FromObjects Embedded in Multilayers
Zekun Zhu, Xiaochao Zhou, and Shunchuan YangPaper submitted to the IEEE Transactions on Antennas and Propagation
Abstract —A single-source (SS) surface integral equation (SIE)for transverse electric (TE) scattering from objects embeddedin multilayers is proposed in this paper. Through recursivelyapplying the surface equivalence theorem from innermost tooutermost boundaries, an equivalent model with surface currentsources enforced on the outermost boundary of the original objectcan be obtained. By further incorporated with the differentialsurface admittance operator (DSAO), only the electric currentdensity is required. In addition, an efficient singularity can-cellation approach is proposed to accurately evaluate singularand nearly singular integrals in the new formulation. Comparedwith other SIEs, like the Poggio-Miller-Chan-Harrington-Wu-Tsai (PMCHWT) formulation and other dual sources (DSs)formulations with both electric and magnetic current densitiesrequired, only single surface equivalent electric current densityis required to be enforced on the outermost boundary in theproposed formulation. It is found that the overall count ofunknowns and memory consumption can be significantly reducedand the conditioning of the final system is much better comparedwith the PMCHWT formulation. As numerical results shown,only a very small fraction amount of unknowns, memory andCPU time, ∼ of unknowns, ∼ of memory and ∼ of CPU time of the PMCHWT formulation in oursimulations, is possibly needed in the proposed formulation. Index Terms —Multilayers, method of moment, single-sourceformulation, surface equivalence theorem, singularity cancella-tion, TE polarization
I. I
NTRODUCTION T HE two-dimensional transverse electric (TE) electro-magnetic analysis exist in many practical engineeringapplications, especially for objects embedded in multilayers,like power cables [1], infinite long transmission line [2]. Itis more involved than the scalar transverse magnetic (TM)problems since vector field equations are required [3].The method of moment (MOM) is widely used to solve theTE scattering problems due to its unknows only residing insource regions and much less count of unknowns is neededcompared with the volumetric mesh based method, like thefinite-difference time-domain (FDTD) method [4], the finiteelement method (FEM) [5]. One of the most commonly
Manuscript received xxx; revised xxx.This work was supported in part by the National Natural Science Foundationof China through Grant 61801010, in part by Beijing Natural ScienceFoundation through Grant 4194082, and in part by Fundamental ResearchFunds for the Central Universities. (Corresponding author: Shunchuan Yang)
Z. Zhu and X. Zhou are with the School of Electronic and InformationEngineering, Beihang University, Beijing, 100083, China (e-mail: [email protected], [email protected].)S. Yang is with the Research Institute for Frontier Science and the School ofElectronic and Information Engineering, Beihang University, Beijing, 100083,China (e-mail: [email protected].) used surface integral equation (SIE) is the Poggio-Miller-Chan-Harrington-Wu-Tsai (PMCHWT) formulation [6]. Inthis formulation, dual sources (DSs) including both electricand magnetic current densities are required on each interfaceof different homogeneous media. Then, several variations, likethe combined tangential formulation (CTF) [7], the Mullerformulation [8], MT-PMCHWT [9] and PMCHWT-CBFM[10], are proposed to improve matrix conditioning and ef-ficiency. However, a large system equation and many two-region problems need to be solved simultaneously when a largenumber of interfaces is involved.Single-source(SS) formulations are proposed to model com-plex composite objects to improve the efficiency. For example,a single-source surface-volume-surface (SS-SVS) formulationis proposed to solve the TE scattering problems of penetrableobjects [11]. Through mapping the surface to volume integraloperator and then back to the surface integral operator, itcan significantly improve the efficiency compared with thevolume integral equation (VIE). Other various SS-SIEs areproposed to model penetrable objects [12]–[15]. Since inthose formulations only SS rather than DSs is required, theyare more efficient than their DSs counterparts. In [16], ageneralized impedance boundary condition (GIBC) is used tomodel interconnects through eliminating the magnetic currentdensity. In addition, the differential surface admittance oper-ator (DSAO) [17] with only single electric current density isproposed to model two-dimensional rectangular interconnects.Then, it is extended to model various shaped cables and inter-connects [18]–[21], and to solve three-dimensional scatteringproblems [22], model antennas [23] [24], and composite lossyconductors [25] [26].Besides those SS formulations, several equivalence principlealgorithms (EPAs) [27] [28] with both electric and magneticcurrent densities are used to solve three-dimensional scatteringproblems. In [29], the EPA combined with the connectionscheme is proposed to model periodic perfectly electric con-ductor (PEC) embedded in planar multilayer media for TMscattering. Through recursively applying the boundary condi-tion at the interfaces between two different planar media, thescattering problem by the multilayer media can be efficientlysolved. In those EPAs, both the electric and magnetic currentdensities are required on a closed surface due to Love’sequivalence theorem [30]. However, when objects are fullyembedded in multilayers, which are interested in this paper,all these formulations can not solve the problem or still sufferfrom the efficiency issue.In this paper, a new vector SS-SIE is proposed to solve the a r X i v : . [ c s . C E ] J un TE scattering problems of objects embedded in multilayersas shown in Fig. 1(a). In the proposed formulation, thesurface equivalence theorem [31] is recursively applied fromthe innermost to outermost boundaries, original objects arereplaced by the background medium as shown in Fig. 1(b).Through incorporating with the differential surface admittanceoperator (DSAO) [17], only single electric current densityis required to be enforced on the outermost boundary. Inaddition, if multiple repetitions within structures, like identicalelements in a large array [23], exist, the equivalent surfacecurrent density can only be calculated once for each typeof scatters. Therefore, it can further improve the efficiency.In addition, an efficient approach through combining thenumerical integration and the analytical integration is proposedto accurately calculate singular and nearly singular integralsin the proposed formulation.It should be noted that the proposed approach is differentfrom those in [32], in which a recursive Green’s function isconstructed through recursively uniting fictitious cells. How-ever, the proposed approach directly uses the surface equiva-lence theorem on each interface and derives the equivalent cur-rent density. Compared with the recursive blockwise inversionof the full matrix method [16], the proposed formulation hasstrictly physical interpretation rather than purely mathematicalmanipulations. This paper explores the possibility which thesurface equivalence theorem solves vector scattering problemsby objects embedded in multilayers and greatly extends ourprevious work [33] [34], which only solve the scalar TMscattering problems.Compared with our previous work and existing literature,its contributions in this paper are twofold.1) A new vector SS-SIE is proposed to solve the TEscattering problems by penetrable or PEC objects em-bedded in multilayers. In the proposed formulation,only single electric current density is required to beenforced on the outermost boundary of objects, whichretains fields the same as those in the original problem.The proposed approach decomposes the original largetwo-region problem into several small scattering ones.Therefore, significant efficiency improvements in termsof CPU time and memory usage can be obtained.2) An efficient singularity cancellation approach is pro-posed to handle various types of singular and nearlysingular integrals in the proposed formulation. Throughcombining the numerical integration and the analyticalintegration, all the matrix elements can be accuratelyand efficiently calculated. Our numerical results showthat the proposed approach can significantly improve theconvergence rate and only uses less than ten integrationpoints to reach the level of machine precision.This paper is organized as follows. In Section II, theconfigurations and preliminaries are first demonstrated. Then,the surface equivalence theorem is briefly summarized todemonstrate the derivation of the SS formulation. The pro-posed SS-SIE formulation is detailed demonstrated to modela single penetrable object and then a penetrable or PECobject embedded in another medium. At last, extension of σ , ε ,μ σ ,ε ,μ i i i σ ,ε ,μ n n n σ ,ε ,μ γ γ i γ n γ in E k σ ,ε ,μ (a) n γ in E k σ ,ε ,μ σ ,ε ,μ in-1 r in r in+1 r n l n+1 l n ( ) ρ r n J (b) Fig. 1. (a) The TE scattering on objects embedded in multilayers, (b) theequivalent model with the electric current density enforced on the outmostboundary γ n . the proposed approach to model objects embedded into anymultilayers is illustrated. In Section III, through combiningthe surface current density and the electric field integralequation (EFIE) on the outermost boundary, the proposed SS-SIE to solve the exterior scattering problem is presented. InSection IV, detailed formulations of the proposed singularitycancellation approach are demonstrated. In Section V, theaccuracy of the proposed singularity cancellation approach andthe performance of the proposed SS-SIE are comprehensivelyinvestigated through several numerical examples. At last, wedraw some conclusions in Section VI.II. M ETHODOLOGY
A. The Configurations and Preliminaries
The TE scattering on penetrable objects with arbitraryshapes embedded in multilayers is considered as shown in Fig.1(a). The penetrable object with constant material parameters ε i , σ i and µ i is enclosed by two adjacent enclosed boundaries γ i . The background medium is air with constant materialparameters ε , σ and µ . Each boundary is discretizedwith N , N , ... , N n line segments. A TE-polarized planewave incidents from the exterior region. According to thesurface equivalence theorem [31], the scattering fields canbe regarded to be generated by surface electric and magneticcurrent densities on the boundary enclosing the objects. Sincefields outside objects are interested in our study, fields insideobjects in the euqivalent problem can be arbitrary . Therefore,an equivalent model as shown in Fig. 1(b) with only theelectric current density J n enforced on the outermost boundarycan be obtained through carefully selecting arbitrary fieldsinside objects after the penetrable objects are replaced by thebackground medium. Before we move to detailed formulationderivation, some symbol notations used in this paper areintroduced.A bold character with both a superscript and a subscript,like U k ( i,j ) , denotes a matrix, where the superscript representsthat this matrix is related to the surface equivalence theoremapplied for the k th time and the subscript ( i, j ) denotesthe i th and j th boundaries are the source and testing ones,respectively. A bold character with only a subscript, like E i ,denotes a vector defined on γ i . A quantity with (cid:98) , like (cid:98) E , denotes that it is for the equivalent problem. The inner productof two vectors f and g is defined as (cid:104) f ( r ) , g ( r ) (cid:105) = (cid:90) γ f ( r ) · g ( r ) d r . (1) B. The Surface Equivalence Theorem
According to the surface equivalence theorem [30], anyenclosed boundary with surface electric and magnetic cur-rent densities can reproduce exactly the same fields in theequivalent configurations. The equivalent electric and magneticcurrent densities are expressed as (cid:98) J i ( (cid:126)r ) = H t i ( (cid:126)r ) − (cid:98) H t i ( (cid:126)r ) , (2) (cid:99) M i ( (cid:126)r ) = E t i ( (cid:126)r ) − (cid:98) E t i ( (cid:126)r ) , (3)where (cid:126)r ∈ γ i , E t i ( (cid:126)r ) , (cid:98) E t i ( (cid:126)r ) , H t i ( (cid:126)r ) , (cid:98) H t i ( (cid:126)r ) are the surfacetangential electric and magnetic fields in the original andequivalent model, respectively. When (cid:98) E t i ( (cid:126)r ) = (cid:98) H t i ( (cid:126)r ) = , itcorresponds to the Love’s equivalence theorem. In this paper,we use other options to derive a SS formulation [17] [23].One option to obtain the single electric current densityformulation is to enforce E t i ( (cid:126)r ) = (cid:98) E t i ( (cid:126)r ) . Then, (2) and (3)are rewritten as (cid:98) J i ( (cid:126)r ) (cid:54) = , (cid:99) M i ( (cid:126)r ) = . (4)When E t i ( (cid:126)r ) = , (4) corresponds to PEC objects and E t i ( (cid:126)r ) (cid:54) = for penetrable objects. Therefore, it is applicablefor both the PEC and penetrable objects.Another option to obtain the single magnetic current densityformulation is to enforce that H t i ( (cid:126)r ) = (cid:98) H t i ( (cid:126)r ) . Therefore, (2)and (3) are expressed as (cid:99) M i ( (cid:126)r ) (cid:54) = , (cid:98) J i ( (cid:126)r ) = . (5)When H t i ( (cid:126)r ) = , (5) corresponds to perfect magneticconductor (PMC) objects and H t i ( (cid:126)r ) (cid:54) = for dielectricobjects. Therefore, it is applicable for both the PMC andpenetrable objects.Both (4) and (5) allow us to derive a SS formulation. In thispaper, (4) is used to derive the SS formulation by enforcingthat all tangential electric fields inside γ i in the original andequivalent problems equal to each other. It should be noted thatfictitious fields are obtained inside the penetrable objects sincethe equivalent therom is applied. If those fields are anywayrequired, they can be recovered through the surface electriccurrent density. C. The SS-SIE for A Penetrable Object Embedded in theBackground Medium
A single penetrable object with constant parameters ε , σ ,and µ and the boundary γ as shown in Fig. 2(a) is firstconsidered. We reported the preliminary idea in [35]. Detailedformulations are presented in this subsection.According to the Stratton-Chu formulation [31], the electricfield E inside γ can be expressed as T ˆt ( E ) = ˆt ( L [ n (cid:48) × H ]) + ˆt ( K [ n (cid:48) × E ]) , (6) σ , ε , μ γ in E k σ ,ε ,μ (a) in E k σ ,ε ,μ σ ,ε ,μ γ J (b) Fig. 2. (a) The TE scattering problem on a penetrable object, (b) the equivalentmodel with the electric current density enforced on γ . where ˆt ( A ) = n × n × A , (7) L ( A ) = − jωµ (cid:73) γ (1 + 1 k ∇ (cid:48) ∇ (cid:48) · ) G ( r , r (cid:48) )( A ) d r (cid:48) , (8) K ( A ) = − (cid:90) γ ( A ) × ∇ (cid:48) G ( r , r (cid:48) ) d r (cid:48) , (9) K is Cauchy principal value without including the residualterm. E and H are electric and magnetic fields inside γ , n (cid:48) de-notes a unit normal vector pointing to the interior region of γ ,the Green’s function is G ( r , r (cid:48) ) = − ( j/ H (2)0 ( k | r − r (cid:48) | ) inside the penetrable object, and H ( · ) is the zero th orderHankel function of the second kind, T = 1 / when r and r (cid:48) are located on the same boundary, otherwise, T = 1 .The vector basis function f n ( r (cid:48) ) [3] and the dual vector basisfunction n (cid:48) × f n ( r (cid:48) ) are used to expand n (cid:48) × H and n (cid:48) × E . f n ( r (cid:48) ) is defined as f n ( r (cid:48) ) = r (cid:48) − r in − l n , r (cid:48) ∈ [ r in − , r in ] r in +1 − r (cid:48) l n +1 , r (cid:48) ∈ [ r in , r in +1 ] , (10)where f n ( r (cid:48) ) is the n th basis function, l n and l n +1 are thelength of the n th and ( n + 1)th segments and r in − , r in , r in +1 are the endpoints of the n th and ( n + 1)th segments asshown in Fig. 1(b), respectively. f n ( r (cid:48) ) mimics the RWG basisfunction in three-dimensional space and is also divergence free[3]. Therefore, n (cid:48) × E and n (cid:48) × H can be expanded as n (cid:48) × E = (cid:88) e n n (cid:48) × f n ( r (cid:48) ) , (11) n (cid:48) × H = (cid:88) h n f n ( r (cid:48) ) . (12) n (cid:48) × E is expanded through the dual vector basis function n (cid:48) × f n ( r (cid:48) ) rather than f n ( r (cid:48) ) to well test the residual term onthe left hand side (LHS) of (6) [36].Then, f n ( r ) is used to test (6) and we obtain T U , E = L , H + K , E , (13)where U , , L , and K , are square matrixes withdimension of N × N . E and H are two column vectorsincluding the expansion coefficients defined as E = [ e , e , e , · · · , e N ] T , (14) H = [ h , h , h , · · · , h N ] T . (15)The entities of matrixes U , , L , and K , are givenby [ U , ] mn = − (cid:10) f m , ˆ t ( f n ) (cid:11) , (16) [ L , ] mn = − (cid:10) f m , ˆ t ( L ( f n )) (cid:11) , (17) [ K , ] mn = − (cid:10) f m , ˆ t ( K [ n (cid:48) × f n ]) (cid:11) . (18)Accurate evaluation of the singular and nearly singularintegrals in (17) and (18) is required to obtain reliable resultswhen the source and testing segments are overlapped and nearto each other. Our numerical experiments show that traditionalsingularity approaches in [37] [38] are not enough to obtainaccurate results. To mitigate the problem, we proposed an-other singularity cancellation approach through combining thenumerical and radial integration approaches to accurately andefficiently calculate singular and nearly singular integrals in(17) and (18). Detailed formulations are presented in SectionIV.By moving the second term on the right hand side (RHS)of (13) to its LHS and inversing the square coefficient matrix,we get H = [ L , ] − [ T U , − K , ] (cid:124) (cid:123)(cid:122) (cid:125) Y E , (19)where Y is the surface admittance operator (SAO) [17] forthe original problem.After the surface equivalence theorem [31] is applied, thepenetrable object is replaced by its surrounding medium andsurface equivalent current density is introduced to ensure fieldsin the exterior region unchanged as shown in Fig. 2(b). Withsimilar procedure above, we obtain (cid:98) H = [ (cid:98) L , ] − [ T (cid:98) U , − (cid:98) K , ] (cid:124) (cid:123)(cid:122) (cid:125) (cid:98) Y E , (20)where (cid:98) Y is the SAO [17] for the equivalent problem andall the material parameters are replaced by those of thesurrounding medium.Therefore, according to (4), the surface equivalent electriccurrent density can be expressed as J = H − (cid:98) H = ( Y − (cid:98) Y ) (cid:124) (cid:123)(cid:122) (cid:125) Y γ E , (21)where J is the expansion coefficients stored in a columnvector when the equivalent current density on γ is expandedby f n ( r (cid:48) ) and Y γ is the DSAO [17]. Magnetic current densityis not required since we keep electric fields on the inner sideof γ in the original and equivalent problems unchanged asshown in (4). D. The SS-SIE for A Two-layered Object Embedded in theBackground Medium
Two scenarios, in which a dielectric or PEC object isembedded in the innermost region as shown in Fig. 3(a) and σ , ε ,μ σ , ε ,μ γ γ in E k σ ,ε ,μ (a) σ , ε ,μ γ γ in E k σ ,ε ,μ PEC (b)
Fig. 3. (a) A two-layered penetrable object is embedded in the backgroundregion, (b) a penetrable object with a PEC object inside is embedded in thebackground region. σ , ε ,μ γ γ in E k σ ,ε ,μ σ , ε ,μ J (a) in E k σ ,ε ,μ σ ,ε ,μ γ γ σ ,ε ,μ J (b) Fig. 4. (a) The TE scattering problem upon object embedded in a two-layered medium and the innermost penetrable object replaced by an equivalentcurrent density, (b) the equivalent model in which the object is replaced by thebackground medium and only an equivalent electric current density enforcedon γ . (b), are considered in this subsection. Our goal is to derivea surface equivalent electric current density J enforced onthe outermost boundary γ through recursively applying thesurface equivalence theorem from γ to γ .
1) A Dielectric Object in the Innermost Region:
Let’sconsider a dielectric object in the innermost region embeddedinto another medium in this subsection as shown in Fig. 3(a).After the inner penetrable object inside γ is modeled throughthe approach presented in Section III-C, a surface equivalentcurrent density J is obtained and enforced on γ as shownin Fig. 4(a). Compared the equivalent model in Fig. 4(a) withthat in Fig. 3(a), an equivalent current density J is enforcedon γ .Then, the electric field inside γ can be expressed throughthe Stratton-Chu formulation [31] as T ˆt ( E ) = ˆt ( L [ n (cid:48) × H ]) + ˆt ( K [ n (cid:48) × E ]) + ˆt [ L ( J )] . (22)The tangential magnetic field n (cid:48) × H and the electric currentdensity J are expanded by f n ( r (cid:48) ) and the tangential electricfield n (cid:48) × E is expanded by n (cid:48) × f n ( r (cid:48) ) . f n ( r ) is used to test(22) on γ . The following matrix equation is obtained T U , E = L , H + K , E + L , J , (23)where U , , L , , K , , and L , are matrixes withdimension of N × N , N × N , N × N , and N × N . J , E , H denote the coefficient column vectors of the currentdensity on γ and tangential electric and magnetic fields on γ , respectively. The entities of the above four matrixes are calculated by [ U , ] mn = − (cid:10) f m , ˆ t ( f n ) (cid:11) , (24) [ L , ] mn = − (cid:10) f m , ˆ t ( L ( f n )) (cid:11) , (25) [ K , ] mn = − (cid:10) f m , ˆ t ( K [ n (cid:48) × f n ]) (cid:11) , (26) [ L , ] mn = − (cid:10) f m , ˆ t ( L ( f n )) (cid:11) . (27)By substituting (21) into (23), moving the last term onthe RHS of (23) to its LSH and then inversing the squarecoefficient matrix, we obtain E = (cid:104) T U , − L , Y γ (cid:105) − (cid:124) (cid:123)(cid:122) (cid:125) C γ (cid:104) L , H + K , E (cid:105) . (28)By using f n ( r (cid:48) ) to test (22) on γ , we have T U , E = L , H + K , E + L , J , (29)where U , , L , , K , , and L , are matrixes withdimension of N × N , N × N , N × N , and N × N .Their entities are similar to U , , L , , K , , and L , and can be calculated by (24)-(27).By substituting (21) and (28) into (29) and rearrangingeach term, the relationship between the tangential electric andmagnetic field on γ can be expressed as H = Y E , (30)where Y = (cid:104) L , + L , Y γ C γ L , (cid:105) − · (cid:104) T U , − K , − L , Y γ C γ K , (cid:105) . (31)
2) A PEC Object in the Innermost Region:
When theinnermost embedded object is PEC, J physically exists andtangential electric fields vanish on γ . Therefore, E equals tozero on the LHS of (23) and we get = L , H + K , E + L , J . (32)By inversing square matrix L , , we get J = − (cid:104) L , (cid:105) − (cid:104) L , H + K , E (cid:105) . (33)Then, after substituting (33) into (29) and rearranging eachterm, we have H = Y PEC E , (34)where Y PEC = (cid:104) L , − L , [ L , ] − L , (cid:105) − (cid:104) T U , − K , − L , [ L , ] − K , (cid:105) . (35)
3) Equivalent Problem:
After the surface equivalence theo-rem is applied on γ , the object is replaced by the backgroundmedium and there is no current source inside γ as shownin Fig. 4(b). With the similar procedure in the previous subsection, we obtain (cid:98) H = [ (cid:98) L , ] − [ T (cid:98) U , − (cid:98) K , ] (cid:124) (cid:123)(cid:122) (cid:125) (cid:98) Y E , (36)where the entitites of matrix (cid:98) L , , (cid:98) U , and (cid:98) K , aregiven by [ (cid:98) U , ] mn = − (cid:10) f m , ˆ t ( f n ) (cid:11) , (37) [ (cid:98) L , ] mn = − (cid:10) f m , ˆ t ( L ( f n )) (cid:11) , (38) [ (cid:98) K , ] mn = − (cid:10) f m , ˆ t ( K [ n (cid:48) × f n ]) (cid:11) . (39)Then, (34) for penetrable objects or (30) for PEC objectscan be substituted into (4). Therefore, the surface equivalentelectric current density J on γ can be expressed as J = H − (cid:98) H = ( Y / PEC − (cid:98) Y ) (cid:124) (cid:123)(cid:122) (cid:125) Y γ E , (40) Y γ is the DSAO on γ , which relates the tangential electricfield to surface equivalence current density.Compared with the micromodeling approach proposed in[34] for antenna array modeling, which is only applicablefor the PEC object embedded in the background medium, theapproach in this paper is applicable for much more general sce-narios including both penetrable and PEC objects. In addition,as shown in the later subsection, the proposed approach is stillsuitable for objects embedded in multilayers. Therefore, themicromodeling approach in [34] can be regarded as a specialcase proposed in this paper. E. Extension of the Proposed Formulation for Objects Embed-ded in Multilayers
For objects embedded in multilayers, the TE scattering fieldscaused by objects can be still equivalent to the scattering fieldsgenerated by a single electric current density on the outermostboundary γ n . To generalize the proposed approach for suchobjects, the surface equivalence theorem is recursively appliedfrom γ to γ n . Eventually, the equivalent current density J n on γ n can be obtained. The procedure is similar to that in theprevious subsection. Then, we get J n = Y γ n E n , (41)where Y γ n is the DASO on γ n , which relates the tangentialelectric field to surface equivalence current density on γ n .When multiple scatters embedded in another medium, ablock-diagonal operator Y γ i is obtained from assembling eachblock DSAO [34].III. S CATTERING M ODELING OF THE P ROPOSED
SS-SIEOnce the equivalent current density J n on γ n is derivedthrough (41), the scattering problem shown in Fig. 1 canbe solved by combining the EFIE. The electric field E inthe exterior region of γ n is the superposition of the incidentfield E inc and the scattered field. Since the surface equivalent current density J n exists on γ n , the total electric field isexpressed as E ( (cid:126)r ) = − jωµ (cid:90) γ n J n ( (cid:126)r (cid:48) ) G ( (cid:126)r, (cid:126)r (cid:48) ) ds (cid:48) + E inc , (42)where (cid:126)r (cid:48) on the outermost boundary γ n of the object, G is the Green’s function expressed as G = − jH (2)0 ( k ρ ) / ,where k is the wavenumber in the background medium. J n is the equivalent current density on γ n . Through taking theobservation points on γ n and expanding the tangential electricfields with the dual vector basis function and the electriccurrent density J n through the vector basis function, we have U n ( n,n ) E n = L n ( n,n ) Y γ n E n + E i . (43)The entities of in E i , U n ( n,n ) and L n ( n,n ) are expressed as [ E i ] m = − (cid:10) f m , ˆ t (cid:0) E inc (cid:1)(cid:11) , (44) [ U n ( n,n ) ] mn = − (cid:104) f m , f n (cid:105) , (45) [ L n ( n,n ) ] mn = − (cid:10) f m , ˆ t ( L ( f n )) (cid:11) . (46)Therefore, the scattering problem can be solved and E n on γ n is expressed as E n = (cid:104) U n ( n,n ) − L n ( n,n ) Y γ n (cid:105) − E i . (47)IV. T HE P ROPOSED A PPROACH FOR A CCURATE C ALCULATION OF S INGULAR AND N EARLY S INGULAR I NTEGRALS
In this paper, we proposed an accurate approach to calculatethe nearly singular and singular integrals through combing thenumerical and analytical integration approaches. The Green’sfunction in the two dimensional space is the zero th orderHankel’s function of the second kind. When the source point r (cid:48) is close to the observation point r , the Green’s function hasa ln R singularity. For the gradient of the Green’s function, ∇ (cid:48) G ( r , r (cid:48) ) , / R singularity exists. Therefore, when r isclose to r (cid:48) , the Green’s function and its gradient show strongsingular and nearly singular behaviors.Two integrals shown in (8) and (9) are required in theproposed formulation. Two different scenarios are consideredaccording to the augments as follows.1) When | k R | = | k ( r − r (cid:48) ) | ≥ . , direct numericalintegration is used to evaluate (8) and (9) in the proposedformulation.2) When | k R | = | k ( r − r (cid:48) ) | ≤ . , small-argumentapproximation is used and the analytical approach tocalculate Green’s function and its gradient to ensure theaccuracy.With the small-augment approximation, the Green’s func-tion and its gradient can be expressed as, G ( | k R | ) ≈ − j (cid:20) − j π ln (cid:18) γ | k R | (cid:19)(cid:21) | R | → , (48) ∇ (cid:48) G ( | k R | ) ≈ − j (cid:32) k jπ | R | (cid:33) R | R | → , (49)where γ = 1 . . x y o r r r r P P l R R τ Fig. 5. The configuration of a test point and a source segment, where R − = r − − r , R + = r + − r , R = r (cid:48) − r , (cid:126)τ = r + − r − | r + − r − | , | P | = | P − r | , l − = R − · (cid:126)τ , l + = R + · (cid:126)τ , (cid:126)l = (cid:126)τ ( R · (cid:126)τ ) . r r r r P R x y o R k P P
Fig. 6. No small argument approximation is used.
The outer integral of the double line integral in (7) and (8)is calculated through the five point Gaussian numerical inte-gration. Careful attention should be paid to accurate evaluationof their inner integrals. Here, we consider a testing point r andthe inner integral segment l i as shown in Fig. 5. p − and p + are the endpoints of the source segment, r (cid:48) is a point on l i and p is the projecting point from the observation point r to l i . Several other quantities are defined as shown in Fig. 5.There are three scenarios for different distance between r and l i . First, the small argument is not applicable toevaluate the Green’s function and its gradient as shownin Fig. 6. Second, the small argument should be usedin the whole l i as shown in Fig. 6. Third, only partial l i is applicable for analytical evaluation of the Green’sfunction and its gradient evaluation as shown in Fig.8. Let | P min | = min ( | R − | , | R + | , | P | ) , | P max | = max ( | R − | , | R + | , | P | ) . Then, we presented the proposed ef-ficiently approach for the singular and nearly singular integralevaluation. A. No small argument approximation is used.
When | k P min | ≥ . as shown in Fig. 6, the Green’sfunction and its gradient can be accurately calculated throughthe recursive relationship of Hankel function [3]. Therefore,Gaussian numerical integration is used to calculate (8) and (9),five points Gassian numerical integration is enough to obtainaccurate results. B. The small argument approximation is used in the whole l i When | k P max | ≤ . as shown in Fig. 7, small agumentsare used to approximate the Hankel function and its gradient.We proposed the following analytical approach to evaluate theinner integral of (8) and (9). r r r r P R x y o max k R P P Fig. 7. The small argument approximation is used in the whole l i . According to (48), when the vector and dual vector basisfunctions are used to expand the tangential electric and mag-netic fields, the analytical results of the operator L ( A ) can beevaluated with small-argument approximation as follows, f m ( r ) · (cid:90) l n L [ f n ( r (cid:48) )] d r (cid:48) = ωµ (cid:26) f m ( r ) · (cid:90) l n f n ( r (cid:48) ) G ( k | P | ) d r (cid:48) −∇ · f m ( r ) · (cid:90) l n ∇ (cid:48) · f n ( r (cid:48) ) G ( k | P | ) d r (cid:48) (cid:27) = c (cid:8) f m ( r ) · (cid:2) I + c I + ( p − r in − ) · ( I + I ) (cid:3) −∇ · f m · ( r ) ( I + c I ) } . (50)Similarly, according to (49), the analytical results of theoperator K ( A ) is caculated as follows, f m ( r ) · (cid:90) l n K [ n (cid:48) × f n ( r (cid:48) )] d r (cid:48) = j f m ( r ) · (cid:90) l n [ n (cid:48) × f n ( r (cid:48) )] × ∇ (cid:48) G ( k | P | ) d r (cid:48) = c f m ( r ) · { v × ( c I + c I ) + v ( c I + c I ) − n (cid:48) ( c I + c I ) + | p | ( c I + c I ) } (51)and c , c , c , c , c are all constants and v , v are constantvectors defined as, c = ωµ l n ( n +1) , c = − jπ ln( γk/ ,c = j l n ( n +1) , c = k , c = 2 j/π, (52) v = n (cid:48) × ( p − r in − n +1) ) , v = [ n (cid:48) × ( p − r in − n +1) )] × ( p − r ) , (53)where l n ( n +1) and r in − n ) denote that when the first part ofthe basis function f n ( r (cid:48) ) is selected, l n and r in − are used,otherwise, l n +1 and r in +1 should be used. r r r p p x o r r p p r y r r p p r p p p P P P x y o (a) r r r p p x o r r p p r y r r p p r p p p P P P x y o (b) r r r p p x o r r p p r y r r p p r p p p P P P x y o (c)Fig. 8. The minimum value for | R | is | p | , | R − | or | R + | in the case of(a), (b) or (c), respectively. C. The small argument approximation is partially used on l i . When | k P min | ≤ . and | k P max | ≥ . , there are threecases shown in Fig. 8. There may be two intersection points, r − and r + in Fig. 8(a) and only one intersection point in Fig.8(b) and Fig. 8(c). Therefore, it is found that (8) and (9) aredivided into three cases as shown in Fig. 8.1) When | k P | ≤ . , | k R − | ≥ . and | k R + | ≥ . as shown in Fig. 8(a), the integral result of the wholesegment may be divided into three parts, where thenumerical integral is used for the intervals [ r − , p − ] , [ p + , r + ] and the analytic method should be applied forthe interval [ p − , p + ] .2) When | k R − | ≤ . and | k R + | ≥ . as shownin Fig. 8(b), the source segment is divided into twoparts, [ r − , p + ] and [ p + , r + ] . Therefore, the inner lineintegration of (8) and (9) includes analytical integrationin [ r − , p + ] and numerical integration in [ p + , r + ] .3) Fig. 8(c) is similar to Fig. 8(b), where | k R − | ≥ . and | k R + | ≤ . .In the proposed approach, the inner line integral can beaccurately calculated and the outer line integral is calculatedthrough the numerical integration. Through carefully con-sidering the position of the source and testing segments asstated above, extremely fast convergence rate can be obtained.As shown in our numerical results, less than ten Gaussianintegration points can reach the level of machine precision.Furthermore, the proposed approach can be easy to be used toevaluate different singular and nearly singular integrals fromformulations, like in [39] [40].V. N UMERICAL R ESULTS AND D ISCUSSION
Two numerical examples are first carried out to verify theconvergence of the proposed singularity cancellation approachproposed in Section IV. Then, three TE-polarization scatteringproblems are solved to demonstrate the performance of theproposed SS-SIE.The workstation with Intel i7-9750H 2.6 GHz CPU and 16Gmemory is used for all the experiments in this section. Our in-house codes are implemented in Matlab and full vectorization x y p p p O Fig. 9. Relative error obtained from the proposed approach and the fullynumerical integration approach with θ = 0 . π . The end points on sourcesegment p , p are ( − . , , [m] and (0 . , , [m], respectively.The testing point p is ( − . , . , [m] and . m away from p .The angle θ between the testing segment and source segments is . π . is explored to speed up the performance. Only a single threadis used for fair comparison. A. Validation of the Proposed Singularity Cancellation Ap-proach
Assuming the source and testing segments are located inthe free space and then we calculated one element of matrixes L and K with a testing point and a source segment withtwo end points of ( − . , , [m] and (0 . , , [m].The frequency used in our simulations is 300 MHz. Twochallenging scenarios are considered to validate the effective-ness of the proposed singularity cancellation approach. First,the testing segment is connected to the source segment at ( − . , , [m]. One point, which is . m away from ( − . , , [m] on the testing segment, is selected as thetesting point. Then, the integration is calculated with rotatingthe testing segment in the counterclockwise direction at thecommon vertice from θ = 0 . π to π . In this configuration, thesource segment is divided into two parts. One part satisfies thesmall argument approximation condition | k R | = | k ( r − r (cid:48) ) | ≤ . . Therefore, the integration on this part are computedthrough the analytical formulations in the proposed approach.The other part does not meet the condition and then theintegration is calculated by the Gaussian numerical integration.Fig. 9 ∼
12 present the relative error corresponding to thiscase. Second, the testing segment is parallel to the sourcesegment. We fixed the testing point at (0 , . , [m], whichis . m away from the vertice of the testing segment. Thesource segment is then divided into three parts. The middlepart is calculated by the small argument approximation andthe analytical formulations in the proposed approach, and theintegration on other two parts is calculated by the Gaussiannumerical integration. Fig. 13 shows results corresponding tothis case. Two approaches including the proposed approachand the fully numerical integration are considered in ourcomparisons. x y p p p O Fig. 10. Relative error obtained from the proposed approach and the fullynumerical integration approach with θ = 0 . π . The end points on sourcesegment p , p are ( − . , , [m] and (0 . , , [m], respectively.The testing point p is ( − . , . , [m] and . m away from p .The angle θ between the testing and source segments is . π . Two observations are made before we move to detailednumerical results:1) After the small variable approximation is applied for theintegral operator L , its imaginary part is a linear functionof | r − r (cid:48) | according to (48), which corresponds to thereal part of L as shown in (8). The integration calculatedby the Gassian integration approach is expected to beaccurate enough. Therefore, it is of little significance tocompare the real parts of results obtained from the twoapproaches. Similarly, the imaginary part of the operator K is not considered in our comparisons.2) The fully numerical integration approach suffers fromthe convergence issue when θ is small. For example, let’sus consider θ = 0 . As the count of integration points in-creases, the testing and source points are approaching toeach other and even coincide. The numerical integrationshows extremely slow convergence and even totally doesnot convergent. Therefore, this scenario is not consideredin our comparison. However, it should be noted that theproposed approach does not suffer from this issue andcan still obtain accurate integration results.Fig. 9 ∼
13 show the relative error obtained from the fullynumerical integration approach and the proposed approachfor two testing cases including the first one with θ =0 . π, . π, . π, π and the second parallel one. The perfor-mance of the fully numerical integration approach changeda lot with different θ as shown in Fig. 9 ∼
12. It requires alarge number of integration points, especially for small θ .As shown in Fig. 9, it needs over 100 integration points toobtain accurate enough results. However, as θ grows larger,the performance of the fully numerical approach in terms ofthe count of integration points becomes better. The reason isthat L and K shows stronger nearly singular behaviors forsmaller θ . However, for the proposed approach, the count ofintegration points is less than 10 points to reach the machineprecision for all scenarios as shown in Fig. 9 ∼
13. Therefore, x y p p p O Fig. 11. Relative error obtained from the proposed approach and the fullynumerical integration approach with θ = 0 . π . The end points on sourcesegment p , p are ( − . , , [m] and (0 . , , [m], respectively.The testing point p is ( − . , . , [m] and . m away from p .The angle θ between the testing and source segments is . π . x y p p p O π Fig. 12. Relative error obtained from the proposed approach and the fullynumerical integration approach with θ = π . The end points on source segment p , p are ( − . , , [m] and (0 . , , [m], respectively. The testingpoint p is ( − . , , [m] and . m away from p . The angle θ between the testing and source segments is π . the proposed approach is quite accurate, efficient and robustto calculate the nearly singular and singular integrals in theproposed SS-SIE.Another scattering example is carried out to further validatethe accuracy of the proposed singularity cancellation approach.The radar cross section (RCS) of a homogeneous dielectricobject with ε r = 4 is considered. Its cross section is anequilateral triangle with the side length of √ m. The sur-rounding medium is air. A TE-polarized plane wave with thefrequency of MHz incidents from the x -axis. The lengthof all segments used to discretize the contour of the object is λ/ , where λ is the wavelength inside the penetrable object.In our comparisons, three approaches including the proposedsingularity cancellation approach, the analytical approach [37][38] and the fully numerical integration approach are used. Forthe analytical approach, only the self segments are treated to Count of integration points -15 -10 -5 R e l a t i v e err o r Real(K)-NumericalReal(K)-ProposedReal(L)-NumericalReal(L)-Proposed x y p p p O p Fig. 13. Relative error obtained from the proposed approach and the fullynumerical integration approach with θ = 0 . π . The end points on sourcesegment p , p are ( − . , , [m] and (0 . , , [m], respectively.The testing point p is (0 , . , [m] and . m away from the endpoint p ( − . , . , [m]. The testing segment is parallel to the sourcesegment. be singular. Then, the Green’s function and its gradient areapproximated by (48) and (49) and analytical integration isapplied. Other interactions are calculated through the Gaussiannumerical integration. 14-points and 5-points Gassian numer-ical integration for the source and testing line integration areused.The reference RCS is calculated from the COMSOL withextremely fine mesh. As shown in Fig. 14, results obtainedfrom the analytical approach can only get accurate RCS from o to o and show large errors in other regions. Thereason is that numerical integration is not accurate enough toevaluate nearly singular integration. When the fully numericalapproach is applied, it is interestingly found that results aremuch more accurate than those of the analytical approach andshow good agreement in the whole range except from o to o . It maybe account for the reason that the numericalintegration has almost uniform level of accuracy. When theproposed singularity cancellation approach is used, resultsshow excellent agreement with the reference RCS. Throughsplitting the nearly singular and singular segments into theanalytical and numerical integration parts and combining thenumerical integration approach and integration cancellationapproach, the singular and nearly singular integrals can beaccurately calculated. B. A Single Pentrable Object
In this subsection, two single penetrable objects, the equilat-eral triangle including non-smoothing boundaries and a dielec-tric cylinder with only smoothing boundaries, are considered.The radius of the dielectric cylinder is m, the relativedielectric constant ε is 4. The background medium is air. Theincident wave is a TE-wave from the x -axis with the frequencyof 300 MHz. The length of all segments used to discretize thecontour of the object is λ/ , where λ is the wavelength insidethe penetrable object. a r ε 4 Fig. 14. RCS obtained from the proposed approach, the analytical approachand the COMSOL for a triangular dielectric object.
The reference solution is obtained from the COMSOL withextremely fine mesh. The RCS of the triangle and cylinderdielectric objects are calculated as shown in Fig. 14, in whichonly results obtained from the COMSOL and the proposedapproach are presented to demonstrate the effectiveness ofthe proposed singular cancellation approach, and Fig. 15 bythe proposed approach, the COMSOL and the PMCHWTformulation, respectively. Results obtained from the threeapproaches are all in excellent agreement.The overall count of unknows for the PMCHWT formula-tion is for the triangle object and for the cylinder.However, only and is used for the proposed approachsince only single current density is required. The conditionnumbers of the PMCHWT formulation and the proposed ap-proach are , , . for the cylinder object and , , for the triangle object, respectively. Therefore, the proposedSS-SIE shows better performance in terms of overall count ofunknows and condition number.Fig. 16 and Fig. 17 show that the relative error with respectto averaged count of segments per wavelength (ACSPW) inthe free space. The relative error in this subsection is definedas RE = (cid:80) i (cid:13)(cid:13) RCS cal ( φ i ) − RCS ref ( φ i ) (cid:13)(cid:13) (cid:80) i (cid:13)(cid:13) RCS ref ( φ i ) (cid:13)(cid:13) , (54)where RCS cal ( φ i ) denotes results calculated from the pro-posed approach and RCS ref ( φ i ) is the reference solutionobtained from the Comsol. As shown in Fig. 16 and Fig.17, it can be found that the relative error fast exponentiallydecreases to the level of − ∼ − as the ACSPW increasesto 10. It is expected that more accurate results can be obtainedfor smoothing than non-smoothing objects. Then, its accuracydoes not improve significantly when the mesh is furtherrefined. This is because the linear basis function (10) is usedin our current implementations. r ε 4 a Fig. 15. RCS caculated by the proposed approach, the PMCHWT formulationand the Comsol for the dielectric cylinder.
ACSPW -4 -3 R e l a t i v e E rr o r Fig. 16. Relative error with ACSPW for the dielectric cylinder.
ACSPW R e l a t i v e E rr o r -3 Fig. 17. Relative error with ACSPW for the dielectric triangular object. " r = 4
C. A Complex Structure
A complex structure with three dielectric 120-degree sector-shaped objects with r = 0 . m and ε r = 9 surrounded bythree concentric cylinders is considered as shown in Fig. 18.The parameters of dielectric cylinders are r = 1 m, r =1 . m, r = 2 m and the relative permittivity ε r = 6 . , ε r =4 , ε r = 2 . . The background medium is air. A plane wavewith f =
300 MHz incidents from the x -axis. The averagedlength of segments used to discretize all the boundaries is λ / , where λ is the wavelength in the free space. It shouldbe noted that this discretization corresponds to λ/ . insidethe sector-shaped object with ε r = 9 .Fig. 18 shows that the RCS obtained from three approaches.It is easy to find that results obtained from the proposedapproach are in excellent agreement with those from theCOMSOL and the PMCHWT. Fig. 20 shows that the relativeerror with respect to ACSPW in the free space. It can be foundthat the relative error of the proposed approach decreases as thecount of segments per wavelength increases. When the countis less than 16, the relative error exponentially converges.Then, after the count is larger than 20, the improvement ofaccuracy is less significant. When 20 segments per wavelengthin the free space is used, which corresponds to λ/ . insidethe sector-shaped object, the relative error almost reaches − . Therefore, the proposed approach shows quite goodconvergence property in terms of mesh size.Table I shows the computational cost for the PMCHWTformulation and the proposed approach. For the PMCHWTformulation, it requires , unknowns to solve this problem.However, only unknowns are required for the proposedapproach. It can be found that the proposed approach onlyuses 5% memory of the PMCHWT formulation and showssignificant memory saving. It should be noted that the vec-torization is fully explored to speed up the in-house matlabcodes. Therefore, more memory consumption are expectedin our simulations. Since only the electric current densityis enforced on the outermost boundary rather than both theelectric and magnetic current densities on all the interfacesin the PMCHWT formulation, significant less unknowns arerequired. In this simulation, the proposed approach only needs unknowns of the PMCHWT formulation. The overallcomputational time to solve this problem for the PMCHWTformulation is . seconds. Although overhead is required to (°) -10-50510152025 RC S ( d B s m ) COMSOLPMCHWTProposed
Fig. 19. RCS obtained from the proposed approach, the PMCHWT formu-lation, and the COMSOL. derive the equivalent current density on the outermost bound-ary, only . seconds is needed for the proposed approach.In addition, the condition number of the final linear systemfor the PMCHWT formulation and the proposed approach is , and , respectively. The proposed approach showsmuch better conditioning since fine geometry details are onlyimplicitly embedded in final system unlike in the PMCHWTformulation, where all these are directly coupled into the finalsystem. Therefore, like formulations [23] [27], the proposedapproach can significantly improve the conditioning and thenthe convergence properties. TABLE IC
OMPARISON OF COMPUTATIONAL COST FOR THE
PMCHWT
FORMULATION AND THE PROPOSED APPROACH
PMCHWT ProposedTotal Time [s]
Memory [MB]
Condition number
Overall Counts of Un-knows
D. An Array with × Single-layered Coated DielectricCylinders
In this subsection, an array with × single-layered coateddielectric cylinder is simulated as shown in Fig. 21. The innerobjects are dielectric cylinders with r = 0 . m and ε r = 4 surrounded by the coating dielectric medium with r = 0 . mand ε r = 2 . . The central distance between two adjacentcylinder elements is d = 1 . m and the background mediumis air. The TE-wave with f = 300 MHz incidents from the x -axis.The RCS is computed from the proposed approach, thePMCHWT formulation and the COMSOL as shown in Fig. ACSPW -3 -2 -1 R e l a t i v e E rr o r Fig. 20. Relative error with respect to ACSPW of the proposed approach inthe free space. TABLE IIC
OMPARISON OF COMPUTATIONAL COST FOR THE
PMCHWT
FORMULATION AND THE PROPOSED APPROACH
PMCHWT ProposedTotal Time [s]
Memory [MB]
Condition Number
Overall Counts of Un-knows λ/ is selected for boundary discretization, where λ is the wavelength in the background medium. It can befound that results obtained from the proposed approach arein excellent agreement with the results from the PMCHWTformulation and the COMSOL.Table II shows the computational cost of the PMCHWTformulation and the proposed approach. The count of unknowsfrom the PMCHWT formulation is , . However, the pro-posed approach only requires only , to solve this problem.Since only the electric current density is required to enforce onthe outermost boundaries, significant less count of unknowns isused in the proposed approach. For the memory consumption, , MB memory is used in the PMCHWT formulation.However, only
MB memory, . of the PMCHWTformulation, is used for the proposed method. The CPUtime to fill coefficient matrix is , . s for the PMCHWTformulation and . for the proposed approach. Therefore,the CPU time is reduced by . The computational timedoes not reduce linearly with the count of unknowns sinceoverhead is required to construct the equivalent current densityon the outermost boundaries. However, the saving in terms ofmemory and CPU time is still quite significant. It is found thatthe matrix condition number of the PMCHWT formulationis very large, which is , while the matrix conditionnumber of the proposed method is only , . It means r ε = 2.25 r ε = 4 Fig. 21. An array with 5 × (°) -20-100102030 RC S ( d B s m ) COMSOLPMCHWTProposed
Fig. 22. RCS obtained from the proposed approach, the PMCHWT formu-lation, and the COMSOL. that the proposed approach shows much better conditioningcompared with the PMCHWT formulation.VI. C
ONCLUSION
The surface equivalence theorem states that any enclosedsurface with appropriate surface equivalent currents sourcescan generate an equivalent configuration with fields exactlythe same as those in the original problems. In this paper,we fully explore this possibility in the method of momentimplementations. A novel vector SS-SIE is proposed to solveTE scattering problems from objects embedded in multilayers.The results reported in this paper shows that the fully vectorelectromagnetic waves are also valid and greatly extended ourprevious researches in [34], which only works for the scalarTM-polarized electromagnetic problems.The proposed approach in this paper only needs a singleelectric current density on the outermost boundary, which canbe derived by recursively applying the surface equivalencetheorem on each boundary from inner to exterior regions.Compared with the PMCHWT formulation, the proposedapproach can significantly improve the efficiency in terms ofoverall count of unknowns, memory consumption, CPU timeand condition number of the linear system. In addition, toaccurately handling the nearly singular and singular integrals in the proposed approach, an efficient singularity cancellationapproach is proposed. Numerical results validate its accuracy,efficiency and convergence property of the singularity cancel-lation approach and the SS-SIE.Extension of current work into three dimensional case tosolve more practical engineering problems is in progress andwe will report more results upon this topic in the future.VII. A PPENDIX
According to the definitions in Section IV, the followingeight identities [41] [42] are used. I = (cid:90) l i ln | R | d r (cid:48) = l + ln | R + | − l − ln | R − | + | p | (cid:18) tan − l + | p | − tan − l − | p | (cid:19) − ( l + − l − ) , (55) I = (cid:90) l i (cid:126)l · ln | R | d r (cid:48) = (cid:126)τ · {−
12 [( l + ) − ( l − ) ]+ | R + | ln | R + | − | R − | ln | R − |} , (56) I = (cid:90) l i | R | d r (cid:48) = 1 | p | (tan − l + | p | − tan − l − | p | ) , (57) I = (cid:90) l i (cid:126)l | R | d r (cid:48) = (cid:126)τ · ( ln | R + | − ln | R − | ) , (58) I = (cid:90) l i (cid:126)l · (cid:126)l | R | d r (cid:48) =( l + − l − ) − | p | · I , (59) I = (cid:90) l i d r (cid:48) = l + − l − , (60) I = (cid:90) l i (cid:126)ld r (cid:48) = (cid:126)τ l + ) − ( l − ) ] , (61) I = (cid:90) l i (cid:126)l · (cid:126)ld r (cid:48) = 13 [( l + ) − ( l − ) ] . (62)Based on the above eight integrals, we can accurately calculatethe nearly singular and singular integrals in (8) and (9) in theproposed SS-SIE. R EFERENCES[1] U. R. Patel, B. Gustavsen, and P. Triverio, “An equivalent surface currentapproach for the computation of the series impedance of power cableswith inclusion of skin and proximity effects,”
IEEE Trans. Power Del. ,vol. 28, pp. 2474-2482, 2013.[2] A. Siripuram and V. Jandhyala, “Surface interaction matrices for bound-ary integral analysis of lossy transmission lines,”
IEEE Trans. Microw.Theory Tech. , vol. 57, no. 2, pp. 347-382, Feb. 2009.[3] W. Gibson,
The Method of Moments in Electromagnetics . CRC press,2015.[4] A. Taflove and S. Hagness,
Computational Electrodynamics: the Finite-Difference Time-Domain Method . Artech house, 2005.[5] J. Jin,
The finite element method in electromagnetics . New York, NY,USA:Wiley, 2015.[6] A. J. Poggio and E. K. Miller, “Integral equation solutions of three-dimensional scattering problems,” in
Computer Techniques for Electro-magnetics . R. Mittra, Ed. Oxford, U.K.: Pergamon, 1973.[7] P. Yl¨a-Oijala, M. Taskinen, and S. J¨arvenp¨a¨a, “Surface integral equa-tion formulations for solving electromagnetic scattering problems withiterative methods,”
Radio Sci. , vol. 40, no. 6, pp. 119, Dec. 2005.[8] C. M¨uller,
Foundations of the Mathematical Theory of ElectromagneticWaves . New York: Springer-Verlag, 1969. vol. 61, no. 1, pp. 341350,Jan. 2013. [9] J. Hu, R. Zhao, Y. Zhang, M. Jiang and Z. Nie, “Fast solving EMscattering from penetrable objects with non-conformal multiple-tracesPMCHWT equations,” in
IEEE Int. Conf. Ubiquitous Wirel. Broadband(ICUWB) , Nanjing, 2016, pp. 1-3.[10] T. Tanaka, Y. Nishioka, Y. Inasawa and H. Miyashita, “Verification of thePMCHWT-CBFM for scattering analysis of a microstrip array antenna,”in
Eur. Conf. Antennas Propag. (EuCAP) , The Hague, 2014, pp. 3232-3236.[11] A. Menshov, V. Okhmatovski, “New single-source surface integral equa-tions for scattering on penetrable cylinders and current flow modelingin 2-D conductors,”
IEEE Trans. Microw. Theory Tech. , vol. 61, no. 1,pp. 341-350, Jan. 2013.[12] Y. Shi, C. H. Liang, “An efficient single-source integral equation solutionto EM scattering from a coated conductor,”
IEEE Antennas Wirel.Propag. Lett. , vol. 14, pp. 547-550, Nov. 2014.[13] F. S. H. Lori, M. S. Hosen, and V. Okhmatovski, “Higher order methodof moments solution of the new vector single-source surface integralequation for 2D TE scattering by dielectric objects,” in
IEEE MTT-SInt. Conf. Numer. Electromagn. Multiphysics Model. Optim. RF, Microw.,Terahertz Appl. (NEMO) , Seville, 2017, pp. 161-163.[14] L. Hosseini, A. Menshov, R. Gholami, J. Mojolagbe, and V. Okhma-tovski, “Novel single-source surface integral equation for scatteringproblems by 3-D dielectric objects,”
IEEE Trans. Antennas Propag. ,vol. 66, no. 2, pp. 797-807, Dec. 2017.[15] D. R. Swatek, I. R. Ciric, “Single source integral equation for wave scat-tering by multiply-connected dielectric cylinders,”
IEEE Trans. Magn. ,vol. 32, no. 3, pp. 878-881, May. 1996.[16] Z. G. Qian, W. C. Chew, and R. Suaya, “Generalized impedanceboundary condition for conductor modeling in surface integral equation,”
IEEE Trans. Microw. Theory Tech. , vol. 55, pp. 2354-2364, Nov. 2007.[17] D. De Zutter, L. Knockaert, “Skin effect modeling based on a differentialsurface admittance operator,”
IEEE Trans. Microw. Theory Tech. , vol. 53,no. 8, pp. 2526-2538, Aug. 2005.[18] U. R. Patel and P. Triverio, “Skin effect modeling in conductors ofarbitrary shape through a surface admittance operator and the contourintegral method,”
IEEE Trans. Microw. Theory Tech. , vol. 64, no. 9, pp.2708-2717, Aug. 2016.[19] U. R. Patel, B. Gustavsen, and P. Triverio, “Proximity-aware calculationof cable series impedance for systems of solid and hollow conductors,”
IEEE Trans. Power Del. , vol. 29, no. 5, pp. 2101-2109, Oct. 2014.[20] P. Utkarsh, and P. Triverio, “MoM-SO: a complete method for computingthe impedance of cable systems including skin, proximity, and groundreturn effects,”
IEEE Trans. Power Del. , vol. 30, no. 5, pp. 2110-2118,Oct. 2015.[21] H. Martijn, D. D. Zutter, D. V. Ginste, “A fully 3-D BIE evaluation ofthe resistance and inductance of on-board and on-chip interconnects,” in
IEEE Workshop Signal Power Integrity (SPI - Proc.) , Brest, 2018, pp.1-4.[22] U. R. Patel, P. Triverio, and S. V. Hum, “A novel single-source surfaceintegral method to compute scattering from dielectric objects,”
IEEEAntennas Wirel. Propag. Lett. , vol. 16, pp. 1715-1718, Feb. 2017.[23] U. R. Patel, P. Triverio, and S. V. Hum, “A macromodeling approach toefficiently compute scattering from large arrays of complex scatterers,”
IEEE Trans. Antennas Propag. , vol. 66, no. 11, pp. 6158-6169, Nov.2018.[24] H. Martijn, M. Gossye, D. D. Zutter, and D. V. Ginste, “A 3-D differen-tial surface admittance operator for lossy dipole antenna analysis,”
IEEEAntennas Wirel. Propag. Lett. , vol. 16, pp. 1052-1055, 2017.[25] H. Martijn, D. D. Zutter, and D. V. Ginste, “Broadband 3-D boundaryintegral equation characterization of composite conductors,” in
IEEEConf. Electr. Perform. Electron. Packag. Syst.( EPEPS) , Montreal, 2019,pp. 1-3.[26] H. Martijn, K. Y. Kapusuz, X. Sun, G. V. D. Plas, E. Beyne, D. D.Zutter, and D. V. Ginste, “Entire domain basis function expansion of thedifferential surface admittance for efficient broadband characterizationof lossy interconnects,”
IEEE Trans. Microw. Theory Tech. , vol. 68, pp.1217-1233, Apr. 2020.[27] M. K. Li, W. C. Chew, “Wave-field interaction with complex structuresusing equivalence principle algorithm,”
IEEE Trans. Antennas Propag. ,vol. 55, no. 1, pp. 130-138, Jan. 2007.[28] M. K. Li and W. C. Chew, “Multiscale simulation of complex struc-tures using equivalence principle algorithm with high-order field pointsampling scheme,”
IEEE Trans. Antennas Propag. , vol. 56, no. 8, pp.2389-2397, Aug. 2008.[29] F. G. Hu and J. M. Song, “Integral equation analysis of scattering frommultilayered periodic array using equivalence principle and connection scheme,” IEEE Trans. Antennas Propag. , vol. 58, no. 3, pp. 848-856,Mar. 2010.[30] C. Balanis,
Antenna Theory: Analysis and Design . 3rd ed. Wiley, 2005.[31] G. W. Hanson, A. B. Yakovlev,
Operator Theory for Electromagnetics .Springer, 2002.[32] M. A. Jensen, J. D. Freeze, “A recursive green’s function method forboundary integral analysis of inhomogeneous domains,”
IEEE Trans.Antennas Propag. , vol. 46, pp. 1810-1816, Dec. 1998.[33] X. C. Zhou, Z. K. Zhu, S. C. Yang, and D. L. Su, “A novel single sourcesurface integral equation for electromagnetic analysis by multilayerembedded objects,” in
Int. Appl. Comput. Electromagn. Soc. Symp.-China (ACES) , Nanjing, 2019.[34] X. C. Zhou, Z. K. Zhu, S. C. Yang, and D. L. Su, “Towards a unifiedapproach to electromagnetic analysis by multilayer embedded objects,”arXiv preprint, arXiv:1909.01854.[35] Z. K. Zhu, X. C. Zhou, J. Wang, and S. C. Yang, “Single-source SIE forTE electromagnetic analysis of penetrable objects,” in
IEEE Int. Symp.Antennas Propag. USNC-URSI Radio Sci. Meet. (APSURSI) , Montrealin,2020, in press.[36] . Ergl, L. Grel, “The use of curl-conforming basis functions for themagnetic-field integral equation,”
IEEE Trans. Antennas Propag. , vol.54, pp. 1917-1926, July. 2006.[37] R. F. Harrington,
Field Computation by Moment Methods . New York:Macmillan, 1968.[38] Y. Liu, H. Z. Yang, “A hierarchicalbutterflyLU preconditioner fortwo-dimensional electromagnetic scattering problems involving opensurfaces,”
J. Comput. Phys. , vol. 401, Jan. 2020.[39] Y. F. Jing, T. Z. Huang, Y. Duan, “A novel integration method for weaksingularity arising in two-dimensional scattering problems,”
IEEE Trans.Antennas Propag. , vol. 58, pp. 2725-2731, Aug. 2010.[40] F. S. H. Lori, A. Menshov, V. Okhmatovski, “New vector single-sourcesurface integral equation for scattering problems on dielectric objects in2D,”
IEEE Trans. Antennas Propag. , vol. 65, pp. 3794-3799, July. 2017.[41] D. R. Wilton, S. M. Rao, A. W. Glisson, “Potential integrals for uniformand linear sourcs distribution on polygonal and polyhedral domains,”
IEEE Trans. Antennas Propag. , vol. AP-32, pp. 276-281, Mar. 1984.[42] J. Stewart,