Instruction FACTA UNIVERSITATIS Series: Electronics and Energetics Vol. 29, N o 4, December 2016, pp. 675 - 688 DOI: 10.2298/FUEE1604675S NONRIGOROUS SYMMETRIC SECOND-ORDER ABC APPLIED TO LARGE-DOMAIN FINITE ELEMENT MODELING OF ELECTROMAGNETIC SCATTERERS  Slobodan V. Savić 1 , Milan M. Ilić 1,2 1 University of Belgrade, School of Electrical Engineering, Belgrade, Serbia 2 Colorado State University, Department of Electrical and Computer Engineering, Fort Collins, CO, USA Abstract. Nonrigorous symmetric second-order absorbing boundary condition (ABC) is presented as a feasible local mesh truncation in the higher-order large-domain finite element method (FEM) for electromagnetic analysis of scatterers in the frequency domain. The ABC is implemented on large generalized curvilinear hexahedral finite elements without imposing normal field continuity and without introducing new variables. As the extension of our previous work, the method is comprehensively evaluated by analyzing several benchmark targets, i.e., a metallic sphere, a dielectric cube, and NASA almond. Numerical examples show that radar cross section (RCS) of analyzed scatterers can be accurately predicted when the divergence term is included in computations nonrigorously. An influence of specific terms in the second-order ABC, which absorb transverse electric (TE) and transverse magnetic (TM) spherical modes, is also investigated. Examples show significant improvements in accuracy of the nonrigorous second-order ABC over the first- order ABC. Key words: absorbing boundary condition, electromagnetic scattering, finite element method, numerical methods 1. INTRODUCTION The finite element method (FEM) is a widely used computational tool in the frequency-domain analysis of electromagnetic (EM) problems [1-4]. To preserve the sparsity of the FEM system when analyzing open-region (radiating and scattering) problems, the necessary artificial truncation of the computational domain is often done by applying approximate local absorbing boundary conditions (ABCs) [4]. The symmetric second-order vector absorbing boundary condition (ABC) is a very popular choice among ABCs because it preserves the symmetry of the FEM system while maintaining Received September 29, 2015; received in revised form December 29, 2015 Corresponding author: Slobodan V. Savić University of Belgrade, School of Electrical Engineering, Belgrade, Serbia (email: ssavic@etf.rs) 676 S. SAVIĆ, M. ILIĆ satisfactory accuracy of the solution [5, 6]. However, this formulation requires computation of the divergence term on the faces of finite elements (FEs) belonging to the absorbing boundary surface (ABS). This, in turn, is a problem on its own because the required normal continuity of the fields is generally not enforced across the edges of adjacent elements in a standard weak-form FEM discretization where edge-based curl-conforming vector basis functions are employed. In addition, a divergence calculation of the nonconforming basis functions in such formulations cannot be done analytically for the generalized curved FEs, even across the faces of elements at the ABS (excluding the troublesome edges) where these functions are continuous and differentiable. This problem has been addressed before, however all reported conclusions pertain to evaluation of the second-order ABC in small-domain spatial discretization frameworks [7-9], where the FEM volume elements are electrically small (e.g., their edges are on the order of /10,  being the wavelength at the operating frequency of the implied time- harmonic excitation). This spatial discretization results in a rather fine mesh throughout the computational domain and at the ABS as well. It appears that in such meshes omitting the divergence term in the second-order ABC, or computing it nonrigorously without enforcing the normal continuity of the fields yields approximately the same error [8]. On the other hand, the method which rigorously implements the second-order ABC on small curved tetrahedra, while preserving the symmetry of the system, has been recently proposed in [9]. However, this method employs auxiliary variables thus mandating significant changes in the existing FEM code. Conversely however, in the open literature there appear to be no analyses of the second-order ABC performance in coarse large-domain FEM meshes, although fine meshes and small elements are really not required at the ABS, which is typically moved away from the analyzed structure and resides in a homogeneous free space. The EM field is usually not changing rapidly at the ABS, hence the advantages of large-domain modeling can be fully exploited. With the above in mind, we proposed that large-domain discretization utilizing curved elements whose edges are up to 2 long, coupled with truly higher order (e.g., up to the 10 th order) polynomial field expansion, can be efficiently used in the ABS tessellation. The number of edges shared by faces of adjacent finite elements at the ABS is thus reduced, which can, in turn, significantly reduce the error introduced by direct computation of required derivatives, because these edges are the sole locations where discontinuities of the normal field components actually arise when the second-order ABC is implemented nonrigorously. Preliminary results of the proposed method applied to a simple metallic spherical scatterer can be found in [10]. In this work we present the implementation details of the nonrigorous symmetric second-order ABC applied on large curvilinear hexahedra in higher-order FEM and evaluate its performance on a comprehensive set of benchmark targets which include: a metallic sphere, a dielectric cube (as an example of penetrable structure with sharp edges and vertices), and a metallic NASA almond as a standard nontrivial benchmark target of the Electromagnetic Code Consortium (EMCC). Nonrigorous Symmetric Second-Order ABC Applied to Large-Domain Finite Element Modeling ... 677 2. THEORY AND IMPLEMENTATION 2.1. Higher-order large-domain FEM formulation When solving three-dimensional (3-D) linear steady-state EM problems by the FEM, we first geometrically discretize the domain of interest using Lagrange-type generalized curved hexahedra of arbitrary orders, Ku, Kv, and Kw (Ku, Kv, Kw  1). These hexahedra are geometrically flexible and can be used for large-domain modeling of arbitrary shapes [11]. They are analytically described by position vector [11]      u v w wvu K i K j K k K k K j K iijk wLvLuLwvu 0 0 0 )()()(),,( rr ,       u u K il l li lK i uu uu uL 0 )( , 1,,1  wvu , (1) where ),,( kjiijk wvurr  are position vectors of interpolation nodes and u K i L represent Lagrange interpolation polynomials in the u coordinate, of the local parametric u-v-w coordinate system, with l u being the uniformly spaced interpolating nodes defined as uul KKlu /)2(  , u Kl ,...,1,0 , and similarly for )(vL v K j and )(wL w K k . We then solve the electric field vector wave equation within each of the finite elements [1, 3]. In every hexahedron we expand the electric field vector as 1 1 1 , , , , , , 0 0 0 0 0 0 0 0 0 u v w u v w u v wN N N N N N N N N u ijk u ijk v ijk v ijk w ijk w ijk i j k i j k i j k                     E f f f , (2) where f are curl-conforming (and generally div-nonconforming) hierarchical polynomial vector basis functions defined as r w k jiijkw r vk j iijkv r ukj i ijku wvPuP wPvuP wPvPu af af af )()( )( )( )()( , , ,    ,             odd ,3, even ,2,1 1,1 0,1 )( iuu iu iu iu uP i ii , 1,,1  wvu , (3) Nu, Nv, and Nw are the adopted degrees of the polynomial approximation, which are entirely independent of the element geometrical orders, Ku, Kv, and Kw, and ijku, , ijkv, and ijkw,  are unknown field-distribution coefficients (to be determined by the FEM). The reciprocal unitary vectors r u a , r v a and r w a in (3) are defined as J wv r u /)( aaa  , J uw r v /)( aaa  and J vu r w /)( aaa  , where wvu J aaa  )( is the Jacobian of the covariant transformation and u a , v a and w a are unitary vectors defined as u u  ra , v v  ra and w w  ra . By adopting higher-order polynomial field expansion [Nu, Nv, and Nw in (2) can be up to 10 th order], through the process of p-refinement, FEs could be up to 2 long in each direction [11]. Applying the standard Galerkin-type discretization yields the disconnected system of linear equations for each of the finite elements [1] 2 0 ([ ] [ ]) { } { } s A k B G   , (4) 678 S. SAVIĆ, M. ILIĆ where k0 represents the free-space wave number and {} is the column vector of electric field distribution coefficients from (2). Disconnected system of linear equations does not take into account boundary conditions which fields must satisfy on the interfaces between two adjacent FEs, but considers each finite element (FE) separately. In order to facilitate implementation (and coding), matrices [A] and [B] can be represented using submatrices as in [11] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] UUA UVA UWA A VUA VVA VWA WUA WVA WWA            , [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] UUB UVB UWB B VUB VVB VWB WUB WVB WWB            . (5) The entries in the submatrices [UVA] and [UVB] are given as 1 ˆ ˆ r ,ˆ̂ ˆ̂, , rˆ ˆ ,ˆ̂ ˆ̂, , ( ) d , d , v ijkijk ijk u ijk V v ijkijk ijk u ijk V UVA V UVB V                      f f f f ,,...,1,0,ˆ ,1,...,1,0 ,,...,1,0ˆ ,,...,1,0 ,1,...,1,0ˆ w v v u u Nkk Nj Nj Ni Ni      (6) where V stands for the volume of the FE and r and r are relative permittivity and permeability tensors [12, 13], respectively. The electric field expansion orders Nu, Nv, and Nw in (2) are selected in accordance with reduced-gradient criterion [14, 15] and by following the recipes in [16] which facilitate optimal higher-order computation. The remaining entries of matrices [A] and [B] are calculated in a similar manner. Analogously, column vector {GS} can be represented as { } { } { } { } S s S S UG G VG WG            , (7) and the entries in the column vector {UGS} are given as                     S kjiukjis SUG d 1 rˆˆ̂,ˆˆ̂, nEf , ,,...,1,0ˆ ,,...,1,0ˆ ,1,...,1,0ˆ w v u Nk Nj Ni    (8) where S stands for the boundary surface of an element, E is the electric field vector at S (generally not known in advance) and n is the unit normal on S pointing outwards of the element. The remaining entries of the column vector {Gs} are calculated in a similar manner. Connected system of linear equations [1] is then assembled from (4) and the surface integrals in {Gs} [as in (8)] are calculated only at the outer boundary of the FEM domain, and not at the boundary of each element [3]. Connected system of linear equations takes into account natural boundary conditions, i.e., tangential continuity of electric fields (explicitly) and magnetic fields (implicitly) which must be satisfied at the interfaces Nonrigorous Symmetric Second-Order ABC Applied to Large-Domain Finite Element Modeling ... 679 between finite elements. Consequently, {Gs} is calculated only at the outer FEM domain boundary, thus it represents a natural connection (interface) between the FEM domain and the surrounding space. Finally, to obtain a well-defined numerical problem, appropriate EM field boundary conditions must be imposed at the outer FEM boundary. These boundary conditions can be (i) exact and nonlocal, as in the hybrid finite element method-method of moments (FEM-MoM) [17], (ii) exact and local, when the FEM domain is surrounded by a perfect electric conductor (PEC) or a perfect magnetic conductor (PMC), or (iii) approximate and local, e.g., when EM field propagation through free space, far from EM sources and media discontinuities, is approximated by an ABC placed relatively close to the scatterer. The local boundary conditions do not reduce sparsity in the final system of linear equations, which is a highly desirable property [18, 19] and one of the strongest benefits of the FEM compared to MoM. 2.2. Symmetric second-order absorbing boundary condition Consider an EM scatterer (or generally EM field sources) occupying a finite volume, surrounded by free space and illuminated by an incident EM field (E inc and H inc ), as shown in Fig. 1. In most cases the incident EM field is a uniform plane wave, but the theory presented here applies to a general case as well. Let SABC be a fictitious spherical surface of radius rABC, centered at the origin and surrounding the scatterer. We truncate the FEM computational domain by applying ABC at SABC. Symmetric (resulting in symmetric system of linear equations) second-order ABC, obtained by approximation of the term sc ( ) r  i E utilizing the Wilcox expansion [20], given as [6]  sc sc scABC0 0 ABC scABC 0 ABC ( ) j ( ) [ ( )] 2(1 j ) ( ) , 2(1 j ) r r r r r t t r k k r r k r                i E i i E i i E E (9) will be applied at SABC, where incsc EEE  represents the scattered electric field, r i is spherical coordinate system radial unit vector, t in subscripts represents the tangential (to SABC) part of a vector or gradient operator and j is the imaginary unit. Fig. 1 With the analysis of open EM problems using ABC. 680 S. SAVIĆ, M. ILIĆ Note that for the connected system of linear equations, the surface integrals in {Gs} are calculated (only) at the entire outer FEM domain boundary SABC, and that they are zero at two finite elements junction. On the other hand, the basis and testing functions appearing in the integrals are taken locally, from a specific element, as the integration progresses. Terms in surface integrals in {Gs} [as in (8)] can be rearranged for easier implementation of the second-order ABC (9) as ABC ABC 1 ˆ ˆ ˆrˆ̂ ˆ̂ ˆ̂, , , d [ ( )] d rs ijk u ijk u ijk S S UG S S                    f E n i E f , (10) since r in  and 1 r [I]   at SABC, with ]I[ being the identity matrix. Applying (10) and imposing the second-order ABC (9), the system of linear equations (4) becomes 2 ABC 0 0 ([ ] [ ] j [ ]) { } { } s A k B k S G    . (11) Matrix [S] in (11) is the sum of three parts: the part corresponding to the first-order ABC, the part corresponding to the second-order ABC, which absorbs transverse electric (TE) spherical modes, and the part corresponding to the second-order ABC, which absorbs transverse magnetic (TM) spherical modes [6, 10]. In the matrix notation this can be written as   TE TM 1ABC 2ABC 2ABCABC 0 0 ABC [ ] [ ] [ ] [ ] , 2 ( j) r S S S S k k r     (12) where the corresponding terms are self explanatory. Analogously as in (5), matrix [S] can be represented using submatrices, namely [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] [ ] UUS UVS UWS S VUS VVS VWS WUS WVS WWS            , (13) where the entries in the submatrix [UVS], for example, are given [in accordance with (12)] as   TE TM 1ABC 2ABC 2ABCABC ˆ ˆ ˆ ˆˆ̂ ˆ̂ ˆ̂ ˆ̂, , , , 0 0 ABC , 2 ( j)ijk ijk ijk ijk ijk ijk ijk ijk r UVS UVS UVS UVS k k r     (14) and analogously for all other submatrices in (13). The entries corresponding to the first- order ABC, the TE part corresponding to the second-order ABC, and the TM part corresponding to the second-order ABC, respectively, are calculated as   ABC TE ABC TM ABC 1ABC ˆ ˆ ,ˆ̂ ˆ̂, , 2ABC ˆ ˆ ,ˆ̂ ˆ̂, , 2ABC ˆ ˆ ,ˆ̂ ˆ̂, , ( ) ( )d , [ ( )][ ( )] d , ( )( ) d , r r v ijkijk ijk u ijk S r r v ijkijk ijk u ijk S t v ijkijk ijk t u ijk S UVS S UVS S UVS S                    i f i f i f i f f f .,...,1,0,ˆ ,1,...,1,0 ,,...,1,0ˆ ,,...,1,0 ,1,...,1,0ˆ w v v u u Nkk Nj Nj Ni Ni      (15) Nonrigorous Symmetric Second-Order ABC Applied to Large-Domain Finite Element Modeling ... 681 The column vector ABC { } s G in (11) can be written in the form shown in (7), with the addition of the superscript “ABC” to distinguish the column vectors in (4) and (11). Hence, similarly as in (12), the column vector ABC { } s G can be represented as the sum of part corresponding to the first-order ABC, the TE part corresponding to the second-order ABC, and the TM part corresponding to the second-order ABC, respectively, as TE TM ABC 1ABC 2ABC 2ABC { } { } { } { } . s s s s G G G G   (16) The entries in the column vector  ABC s UG , for example, are given as   ABC TE ABC TM 1ABC inc inc ˆ ˆ ˆ0ˆ̂ ˆ̂ ˆ̂, , , 2ABC incABC ˆ ˆˆ̂ ˆ̂, , 0 ABC 2ABC ABC ˆ ˆˆ̂ ˆ̂, , 0 ABC ( ) ( ) j ( ) ( ) d , [ ( )][ ( )] d , 2(1 j ) ( ) 2(1 j ) r r rs ijk u ijk u ijk S r rs ijk u ijk S s ijk t u ijk UG k S r UG S k r r UG k r                        i f E i f i E i f i E f ABC inc ( ) d , t S S     E ,,...,1,0ˆ ,,...,1,0ˆ ,1,...,1,0ˆ w v u Nk Nj Ni    (17) and analogously for the remaining entries in ABC { } s G . 2.3. Computation of the surface integrals appearing in the symmetric second-order absorbing boundary condition applied to curvilinear elements Consider the surface integrals appearing in (11) when computing entries in [S] and ABC { } s G . The utilized basis and testing functions are curl-conforming and generally div- nonconforming, hence the divergences in the TM parts of (15) and (17), and all similar terms, cannot be expressed in the closed form. Moreover, as already discussed, these surface integrals are calculated over the entire SABC surface; in other words, they are calculated not only over the finite element surfaces belonging to SABC, but across the junctions (edges between the elements) as well. Since the basis and testing functions possess only tangential continuity, this results in appearance of squares of delta-functions ( 2 ) in the kernels of the surface-integral terms at all edges enveloping the surfaces of the finite elements belonging to ABC S [9]. In order to rigorously treat the divergence of the basis and testing functions at the edges of elements over SABC, the basis and testing functions must be adopted to enforce the normal continuity of the EM field over SABC [8] or additional auxiliary (scalar) variables need to be introduced as in [9]. Nevertheless, since the utilized higher-order polynomial basis and testing functions are continuous and differentiable over FEs faces, their divergence can be readily calculated numerically. For example, from (3) it follows that the divergence of fu,ijk is given as 1 , 1 ( ) ( ) ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) ( ) . i r r i r r u ijk j k u u j k u u ji r r i r r k u v j k u v i r r i r rk j u w j k u w iu P v P w u P v P w J J u P v u P w u P v P w J v J v P w u P v u P v P w J w J w                         f a a a a a a a a a a a a (18) 682 S. SAVIĆ, M. ILIĆ Partial derivatives in (18) are calculated numerically utilizing the symmetric finite difference. For example, ' d ' d ( ) ( ) ( ) , 2 d r r r r u v u vr r v v v v v v u v J J J v v           a a a a a a (19) where vd is a numerical-differentiation step. Since these divergences are computed only at the FEM domain-truncation boundary SABC, numerical differentiation represents minimal addition to the complexity of the overall algorithm, and computation time for the surface integrals ABC { } s G is almost negligible compared to the computation time for the FEM volume integrals appearing in matrices [A] and [B]. The procedure is similar when divergence is calculated for the functions ijkv,f and ijkw,f . 3. NUMERICAL RESULTS AND DISCUSSION 3.1. PEC spherical scatterer As the first numerical example, consider a PEC spherical scatterer of radius a = 1 m. The scatterer is situated in free space, with permittivity 0  and permeability 0  , and illuminated by a time-harmonic plane-wave of a free space wavelength m1 0  (f = 299.792 MHz), as shown in Fig. Error! Reference source not found. (a). When constructing numerical model, infinite free space surrounding the scatterer is truncated at the artificial spherical boundary SABC, of radius m5.1b , where the nonrigorous symmetric second-order ABC is imposed. The normalized thickness of the free space layer between the scatterer and SABC is 5.0)( 0  ab and it is meshed by only six cushion-like triquadratic curved hexahedral FEs. 0 1 2 3 4 5 6 7 8 9 10 10 -3 10 -2 10 -1 10 0 1 st ord. ABC 1 st ord. ABC with G s 2ABC, TE and S 2ABC, TE 1 st ord. ABC with G s 2ABC, TM and S 2ABC, TM Nonrigorous 2 nd ord. ABC Unknowns FEM-ABC L 2 n o rm ( B iR C S ) / L 2 n o rm (M ie B iR C S ) N 10 1 10 2 10 3 10 4 10 5 10 6 U n k n o w n s (a) (b) Fig. 2 (a) Large-domain FEM-ABC model of a PEC spherical scatterer. (b) Normalized L 2 error norm of the computed bistatic RCS for the PEC spherical scatterer and the number of unknowns. Nonrigorous Symmetric Second-Order ABC Applied to Large-Domain Finite Element Modeling ... 683 First, we will consider far field results. A bistatic radar cross section (RCS) of the scatterer is computed by the proposed FEM-ABC technique. The order of the polynomial expansion of the electric field for all FEs and in all directions is Nu = Nv = Nw = N. Numerical integration is performed by means of the 13 th order Gauss-Legendre quadrature. The bistatic RCS is computed in all directions uniformly (from  0start to  180 stop with the resolution of  5 , and from  0start to  360stop with the resolution of  5 ), and its error (with respect to the analytical Mie’s series solution) is calculated as a normalized 2 L norm 180 360 MieBiRCS 2 MoMBiRCS2 0 0 2 MieRCS 180 360 MieBiRCS 2MoMRCS MoMBiRCS 0 0 (FEMBiRCS( , ) ( , )) L norm( BiRCS) L norm( ) ( , )                        , (20) where FEMBiRCS stands for the numerical solution for the bistatic RCS obtained by the proposed FEM-ABC technique and MieBiRCS stands for the analytical (reference) results in the form of Mie’s series. In the following subsection, when analytical MieBiRCS solution is not available, the results obtained by MoM, denoted as MoMBiRCS, will be used as a reference, as indicated in (20). In Fig. Error! Reference source not found. (b) numerical results are compared for the first- and nonrigorous second-order ABC, along with results for the first-order ABC with only one term included from the nonrigorous second-order ABC [ , TE ABC2 s G TE ABC2 S and , TM ABC2 s G TM ABC2 S from (12) and (16)]. To validate the convergence of the method with p-refinement, the solutions are obtained for various orders N, ranging from N = 1 to N = 9. From Fig. Error! Reference source not found. (b) it can be concluded that, although not being implemented rigorously and not contributing independently to the accuracy of the solution, the TM part of the symmetric second-order ABC together with the TE part synergistically contributes to the overall solution accuracy. In addition, due to very rough mesh in this example, the FEM solution becomes sufficiently accurate for 97  N with N = 8 yielding the lowest error, which is consistent with the results reported in [16]. Moreover, the lowest errors obtained with the proposed large-domain FEM with the nonrigorous second-order ABC are of the same order of magnitude as those reported in the first example in [9], where the same scatterer was analyzed utilizing the rigorously implemented second-order ABC. In this example the nonrigorous second-order ABC performs significantly better in far field compared to the first-order ABC, and for N = 8 the solution error is 2.7 times lower compared to results obtained utilizing the first-order ABC. Note that this error difference is even greater (8.8 times in favor of the nonrigorous second- order ABC) when the ABC is set closer to the scatterer, i.e., when 1.0)( 0  ab , as reported in [10]. Noting that far fields, and related derived parameters, are less sensitive to computational errors than near fields, in order to obtain and demonstrate an even more rigorous and complete validation of the proposed FEM-ABC technique, we next analyze the accuracy of the computed near field of the presented PCE spherical scatterer. Using the mesh from Fig. Error! Reference source not found. (a) and setting N = 8 (for all elements in all direction) we compute the near electric field numerically and analytically 684 S. SAVIĆ, M. ILIĆ and show the comparison of obtained results in Fig. 3. Shown in Fig. 3 is the magnitude of the x-component of the total electric field, in the 0x plane, obtained (a) analytically (Mie’s series solution) and numerically using (b) the first-order ABC and (c) the proposed second-order ABC. The incident electric field is ]m/V[1 inc xiE  ( xi being the Cartesian unit vector in the x-direction) traveling in the z-direction, as shown in Fig. 3 (d). In Figs. 3 (e) and (f) the error of the electric field computed by the FEM (relative to the reference Mie’s series solution) for the first-order and second-order ABC models are plotted, respectively. The error is calculated as 2ImMie, Im FEM, 2Re Mie, Re FEM, )()( xxxxx EEEEE  , where Ex,FEM and Mie,xE Ex,Mie are x-components of the electric fields obtained numerically and analytically, respectively, and Re and Im stand for the real and imaginary part of the complex quantities, respectively. (a) (b) (c) (d) (e) (f) Fig. 3 Near field results for the PEC spherical scatterer from Fig. Error! Reference source not found. obtained (a) analytically and numerically using (b) the first-order and (d) the proposed second-order ABC. (d) Large-domain FEM-ABC model of a PEC spherical scatterer with illustrated incident field. Electric field error (relative to the reference Mie’s series solution) for (e) the first-order and (f) the proposed second- order ABC. From Fig. 3, it can be concluded that the proposed second-order ABC significantly outperforms the first-order ABC. The results obtained using nonrigorous second-order ABC are more accurate than those using the first-order ABC in the complete x = 0 plane, and especially for z > 0. Note that, due to symmetry, the remaining two Cartesian components of the electric field vanish in the 0x plane (Ey = 0, Ez = 0), hence they are not shown. Also, note that other field components in different planes exhibit similar Nonrigorous Symmetric Second-Order ABC Applied to Large-Domain Finite Element Modeling ... 685 errors, hence they are not shown here for brevity. In addition, the errors in the near field can be further reduced employing p-refinement. 3.2. Dielectric cubical scatterer As the second numerical example, consider a dielectric cubical scatterer with relative permittivity 25.2 r  and relative permeability 1 r  , of edge length m2a . The scatterer is situated in free space and illuminated by a time-harmonic plane-wave of a free space wavelength m2 0  (f =149.896 MHz), as shown in Fig. 4 (a). When constructing the numerical model, infinite free space surrounding the scatterer is truncated at the artificial spherical boundary SABC, of radius m2b , where the nonrigorous symmetric second-order ABC is imposed. Free space between the scatterer and the ABC S is again meshed by only six cushion-like triquadratic curved hexahedral FEs and the dielectric scatterer is meshed by only one trilinear FE. Minimal normalized distance between the scatterer and ABC S is 13.0)35.0( 0  ab and this maximal distance is (b  0.5a)/0 = 0.5. 0 1 2 3 4 5 6 7 8 9 10 10 -3 10 -2 10 -1 10 0 1 st ord. ABC 1 st ord. ABC with G s 2ABC, TE and S 2ABC, TE 1 st ord. ABC with G s 2ABC, TM and S 2ABC, TM Nonrigorous 2 nd ord. ABC Unknowns FEM-ABCL 2 n o rm ( B iR C S ) / L 2 n o rm (M o M B iR C S ) N 10 1 10 2 10 3 10 4 10 5 10 6 U n k n o w n s (a) (b) Fig. 4 (a) Large-domain FEM-ABC model of a dielectric cubical scatterer. (b) Normalized L 2 error norm of the computed bistatic RCS for the dielectric cubical scatterer and the number of unknowns. Normalized L 2 error norm of the computed bistatic RCS for the cubical scatterer is calculated as discussed in subsection 0 and shown in Fig. 4 (b). The error is calculated with respect to the fully converged MoM solutions obtained by WIPL-D software [21]. Numerical parameters regarding the field expansion and integration in the FEM model are kept the same as in the previous example. It can be concluded based on Fig. 4 (b) that the nonrigorously implemented TM part of the second-order ABC independently contributes to the quality of solutions and that, together with TE part of the second-order ABC, both parts synergistically contribute to the overall solution accuracy. In this example, the nonrigorous second-order ABC performs significantly better compared to the first-order 686 S. SAVIĆ, M. ILIĆ ABC, and for 7N the error obtained using the second-order ABC is 5.6 times smaller than that for the first-order ABC. 3.3. PEC NASA almond scatterer As the last example, consider a PEC NASA almond scatterer, which is one of the standard benchmarks of the EMCC. The NASA almond is geometrically described by the parametric equations given above Fig. 2 in [22]. The almond of length mm37.252d (parameter d from equations in [22]), situated in free space, and illuminated by horizontally and vertically (in  90 plane) polarized incident EM field at the operating frequency GHz19.1f ( mm252 0  ) will be considered, as shown in Fig. 5. Fig. 5 PEC NASA almond scatterer. Higher-order FEM-ABC model of the PEC NASA almond scatterer consists of 96 triquadratic large-domain Lagrange-type FEs. These FEs model the free space between the almond and the spherical surface ABC S , where nonrigorous symmetric second-order ABC is applied. The radius of ABC S is mm220b . Minimum and maximum distances from the almond to ABC S are 0 373.0  and 0 801.0  , respectively, and the field expansion orders are set to 6N (for all finite elements and in all directions), which results in 62220 unknown field distribution coefficients. Using the proposed nonrigorous second-order ABC coupled with the 0 30 60 90 120 150 180 -50 -45 -40 -35 -30 -25 -20 -15 -10 WIPL-D FEKO ----------------------------------- Higher Order FEM-ABC N u =N v =N w =6, 62220 Unkn. Nonrigorous 2 nd ord. ABC = 90 0 M o n o st a ti c R C S [ d B m 2 ]  0 30 60 90 120 150 180 -50 -45 -40 -35 -30 -25 -20 -15 -10 WIPL-D FEKO ----------------------------------- Higher Order FEM-ABC N u =N v =N w =6, 62220 Unkn. Nonrigorous 2 nd ord. ABC M o n o st a ti c R C S [ d B m 2 ]  = 90 0 (a) (b) Fig. 6 Computed monostatic RCS of the PEC NASA almond from Fig. 5 for the (a) horizontal and (b) vertical incident field polarization; comparison of proposed FEM-ABC and two MoM results obtained by WIPL-D [21] and FEKO [23] software. Nonrigorous Symmetric Second-Order ABC Applied to Large-Domain Finite Element Modeling ... 687 large-domain higher-order FEM technique, the monostatic RCS in the horizontal plane (  90 , )1800  is computed. The results are compared with results obtained by MoM technique [21, 23] for both horizontal and vertical incident field polarizations, and shown in Fig. 6. From Fig. 6 it can be concluded that a very good matching between the FEM-ABC and MoM results is achieved in all directions, and that scatterers of relatively complex shapes can also be accurately analyzed by the proposed FEM-ABC method. 4. CONCLUSIONS We have presented, implemented, and validated by representative numerical experiments, a nonrigorous symmetric second-order ABC in combination with large- domain higher-order FEM technique for frequency domain EM scattering analysis. In the proposed method, the ABC is implemented nonrigorously, without imposing the normal field continuity and without introducing additional variables. The required divergence of the nonconformal field components is computed numerically on the faces of elements belonging to the ABS, using simple finite differences. Numerical experiments have shown that the nonrigorous second-order ABC performs significantly better compared to the first-order ABC and that the proposed method results mach very good with referent numerical solution of high accuracy. Moreover, the examples have shown that the errors in computation of the RCS can be significantly lower if the divergence term is included in the ABC, as described, than if it is omitted. This conclusion is in contrast with results reported thus far in the literature, where examples with small-domain FEM meshes have been utilized exclusively. Finally, examples with a dielectric cubical scatterer and the NASA almond have shown that the proposed method can be successfully applied in analysis of scatterers with sharp edges and tips. Acknowledgement: This work was supported by the Serbian Ministry of Science and Technological Development under grant TR-32005. REFERENCES [1] P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers, 3 ed. New York: Cambridge University Press, 1996. [2] J. L. Volakis, A. Chatterjee, and L. C. Kempel, Finite Element Method for Electromagnetics (Antennas, Microwave Circuits, and Scattering Applications), 1 ed. New York: IEEE Press, 1998. [3] J.-M. Jin, The Finite Element Method in Electromagnetics. Hoboken, New Jersey: John Wiley & Sons, 2014. [4] J.-M. Jin and D. J. Riley, Finite Element Analysis of Antennas and Arrays, 1 ed. Hoboken, New Jersey: Wiley-IEEE Press, 2009. [5] J. P. Webb and V. N. Kanellopoulos, "Absorbing Boundary Conditions for The Finite Element Solution of The Vector Wave Equation," Microwave and Optical Technology Letters, vol. 2, no. 10, pp. 370-372, October 1989. [6] A. F. Peterson, "Accuracy of 3-D Radiation Boundary Conditions for Use with the Vector Helmholtz Equation," IEEE Transactions on Antennas and Propagation, vol. 40, no. 3, pp. 351-355, March 1992. [7] V. N. Kanellopoulos and J. P. Webb, "3D finite element analysis of a metallic sphere scatterer: comparison of first and second order vector absorbing boundary conditions," Journal de Physique III, vol. 3, no. 3, pp. 563-572, March 1993. 688 S. SAVIĆ, M. ILIĆ [8] V. N. Kanellopoulos and J. P. Webb, "The Importance of the Surface Divergence Term in the Finite Element-Vector Absorbing Boundary Condition Method," IEEE Transactions on Microwave Theory and Techniques, vol. 43, no. 9, pp. 2168-2170, September 1995. [9] M. M. Botha and D. B. Davidson, "Rigorous, Auxiliary Variable-Based Implementation of a Second- Order ABC for the Vector FEM," IEEE Transactions on Antennas and Propagation, vol. 54, no. 11, pp. 3499-3504, November 2006. [10] S. V. Savić, B. M. Notaroš, and M. M. Ilić, "Accuracy Analysis of the Nonrigorous Second-Order Absorbing Boundary Condition Applied to Large Curved Finite Elements," in 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), Turin, Italy, 2015, pp. 58-61. [11] M. M. Ilić and B. M. Notaroš, "Higher order hierarchical curved hexahedral vector finite elements for electromagnetic modeling," IEEE Transactions on Microwave Theory and Techniques, vol. 51, no. 3, pp. 1026-1033, March 2003. [12] S. V. Savić, A. B. Manić, M. M. Ilić, and B. M. Notaroš, "Efficient Higher Order Full-Wave Numerical Analysis of 3-D Cloaking Structures," Plasmonics, vol. 8, no. 2, pp. 455-463, June 1 2013. [13] S. V. Savić, B. M. Notaroš, and M. M. Ilić, "Conformal cubical 3D transformation-based metamaterial invisibility cloak," Journal of the Optical Society of America A, vol. 30, no. 1, pp. 7-12, January 2013. [14] J. C. Nedelec, "Mixed Finite Elements in R3," Numerische Mathematik, vol. 35, no. 3, pp. 315-341, September 1980. [15] J. C. Nedelec, "A new family of mixed finite elements in R3," Numerische Mathematik, vol. 50, no. 1, pp. 57-81, January 1986. [16] E. M. Klopf, N. J. Šekeljić, M. M. Ilić, and B. M. Notaroš, "Optimal Modeling Parameters for Higher Order MoM-SIE and FEM-MoM Electromagnetic Simulations," IEEE Transactions on Antennas and Propagation, vol. 60, no. 6, pp. 2790-2801, June 2012. [17] M. M. Ilić, M. Djordjević, A. Ţ. Ilić, and B. M. Notaroš, "Higher Order Hybrid FEM-MoM Technique for Analysis of Antennas and Scatterers," IEEE Transactions on Antennas and Propagation, vol. 57, no. 5, pp. 1452-1460, May 2009. [18] G. Strang, Linear Algebra and Its Applications, 4 ed.: Brooks Cole, 2005. [19] G. Strang, Introduction to Linear Algebra, 4 ed. Wellesley, MA: Wellesley - Cambridge Press, 2009. [20] C. H. Wilcox, "An Expansion Theorem for Electromagnetic Fields," Communications on Pure and Applied Mathematics, vol. 9, no. 2, pp. 115-134, May 1956. [21] "WIPL-D Pro," 11.0 WIPL-D d.o.o., 2013 Available: http://www.wipld.com. [22] A. C. Woo, H. T. G. Wang, M. J. Schuh, and M. L. Sanders, "Benchmark Radar Targets for the Validation of Computational Electromagnetics Programs," IEEE Antennas and Propagation Magazine, vol. 35, no. 1, pp. 84-89, February 1993. [23] "FEKO," Altair Development S.A. (Pty) Ltd,, 2011 Available: http://feko.info/applications/RCS. http://www.wipld.com/ http://feko.info/applications/RCS