An amplitude finite element formulation for electromagnetic radiation and scattering

Computers and Mathematics with Applications 71 (2016) 1364–1391

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

An amplitude finite element formulation for electromagnetic radiation and scattering Arup Nandy, C.S. Jog ∗ Department of Mechanical Engg., Indian Institute of Science, Bangalore 560012, India

article

info

Article history: Received 2 April 2015 Received in revised form 31 December 2015 Accepted 10 February 2016 Available online 8 March 2016 Keywords: Electromagnetic radiation and scattering Amplitude formulation Nodal finite elements High frequency

abstract Electromagnetic radiation and scattering in an exterior domain within the context of a finite element method has traditionally involved imposing a suitable absorbing boundary condition on the truncation boundary of the numerical domain to inhibit reflection from it. In this work, based on the Wilcox asymptotic expansion of the electric far-field, we propose an amplitude formulation within the framework of the nodal finite element method, whereby the highly oscillatory radial part of the field is separated out a-priori so that the standard Lagrange interpolation functions that are used have to capture a relatively gently varying function. Since these elements can be used in the immediate vicinity of the radiator or scatterer (with few exceptions which we enumerate), it is more effective compared to methods that impose absorbing boundary conditions at the truncation boundary, especially for high-frequency problems. The proposed method is based on the standard Galerkin finite element formulation, and uses standard Lagrange interpolation functions, standard Gaussian quadrature and the same degrees of freedom as a conventional formulation. We show the effectiveness of the proposed formulation on a wide variety of radiation and scattering problems involving both conducting and dielectric bodies, and involving both convex and non-convex domains with sharp corners. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction The traditional approach to solving electromagnetic radiation and scattering problems on unbounded domains within the context of a finite element method has been to truncate the domain at some finite distance, and imposing some suitable absorbing boundary condition (ABC) on the truncation boundary so as to reduce the reflection of the outgoing waves from it. Results using the most simple (first-order) ABC, known as Sommerfeld radiation condition, are presented in [1]. Some second and higher-order ABCs have been presented in [1–4]. A modified form of the second-order ABC (‘Joly–Mercier ABC’) has been presented in [5], and has been applied to the three-dimensional Maxwell equations in [6]. Another approach to model the Sommerfeld radiation condition accurately at the truncation boundary has been to combine the Boundary Integral method (BI) with the standard finite element method [7–15]. In these works, the interior domain is modeled with edge-based finite elements and the truncation surface with BI, with field continuity imposed between the two. Although BI terminations are more accurate compared to the conventional ABCs, a major drawback of this method is the fully populated matrix for the BI block, resulting in increased memory requirements and computational time. A relatively recent approach introduced by Berenger is the concept of a ‘perfectly matched layer’ [16–23], where a ‘lossy material’ which surrounds the domain is used to damp out the propagating wave. Another recent approach to

∗

Corresponding author. E-mail addresses: [email protected] (A. Nandy), [email protected] (C.S. Jog).

http://dx.doi.org/10.1016/j.camwa.2016.02.013 0898-1221/© 2016 Elsevier Ltd. All rights reserved.

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

1365

model the non-reflecting boundary conditions is the use of infinite elements [24–26], which the authors have shown to be more efficient compared to a conventional formulation based on an ABC. Conventional elements are required between the radiator/scatterer and the infinite elements in most of the infinite element formulations, while the proposed amplitude formulation can be used in the full domain provided the origin is not part of the domain. Special shape functions are required for the infinite elements, while in the proposed formulation standard Lagrange interpolation functions are used. A comparison of the performance is carried out in Section 5.4. In most of the existing literature, electromagnetic radiation and scattering problems have been solved using edge elements, for example, see [1–28]. These elements impose only tangential continuity of the primary field variable across interelement boundaries, whereas the normal component is allowed to be discontinuous. There are relatively very few works that use nodal finite elements either for interior or exterior electromagnetic problems, see for example, [29–37]. While the use of nodal finite elements is especially desirable while solving multiphysics problems where structural or fluid flow variables need to be coupled with the electromagnetic field variables, the main difficulty has been in modeling sharp corners and edges where a regularized formulation (which is generally used to eliminate the spurious modes) is unable to model the singular fields correctly [1,29,32]. Otin and coworkers [33–35] have suggested a remedy to this problem, namely, that of using a zero penalty in regions near to the singularity. We use a slightly modified version of this strategy in our work in conjunction with the new proposed strategy that we outline below for solving scattering problems involving nonconvex domains with sharp corners. In the current work, we propose an amplitude formulation that is based on the Wilcox asymptotic expansion of the electric field at large distances from the radiator/scatterer [38,39]. Although a similar idea was used in [40–42] in the context of acoustical problems, the presence of singularities, material discontinuities and the existence of spurious modes in the current problem (to name just a few of the additional complications) makes an extension of that strategy to the current problem nontrivial. A similar idea in the electromagnetics context was also proposed in [43]. However, they restrict their formulation to two-dimensional problems, do not consider scatterers with nonconvex boundaries/sharp corners, and require the solution of an equation for the phase. So in some sense, the current work may be considered as a generalization of their work to solve more realistic problems. If r denotes the radial coordinate of a spherical coordinate system, and k the wave number, the main idea in the amplitude formulation is to a-priori separate out the highly oscillatory e−ikr part from the interpolation function so that the finite element interpolation function has to capture only a relatively gently varying function which it does quite effectively, as we will demonstrate with several examples where the solutions obtained using the proposed strategy are compared with analytical solutions and existing numerical strategies. The superior performance of the proposed amplitude formulation over a conventional strategy that uses an absorbing boundary condition is more evident at higher frequencies. Furthermore, unlike the PML where new elements are introduced in a ‘lossy’ layer surrounding conventional elements, the AMP elements (elements based on the proposed amplitude formulation) can be used in the immediate vicinity of the radiator/scatterer (except when the origin is part of the domain, in which case, we use a combination of conventional elements and AMP elements with the conventional elements used near to the origin); following this principle, all the numerical examples in Section 5 which deal with radiation and scattering from conductors use only AMP elements, while those dealing with dielectrics use a combination of conventional and AMP elements. One aspect that we would like to emphasize at this stage is that although the amplitude formulation have been developed within the context of a nodal finite element formulation in this work (since our goal is to ultimately couple it with structural or fluid mechanics variables where nodal finite elements, which have a very different data structure compared to edge elements, are generally used), the same ideas can be used to develop amplitude formulation in edge elements. Thus, there is no intrinsic dependence of the amplitude formulation on the nodal finite element framework. The outline of the remainder of the paper is as follows. In Sections 2–4, we develop the variational formulation and the finite element formulations for the conventional and the AMP elements; the formulations for both radiation and scattering from conductors and dielectrics are presented. Section 5 presents a series of examples demonstrating the superior performance of the amplitude formulation over the conventional strategy that uses an absorbing boundary condition by comparing the solutions obtained against either analytical solutions or other numerical results such as edge element results obtained using HFSS [44]. Besides considering scattering from ‘smooth’ bodies such as spheres or ellipsoids, we also consider examples where the scattering body has sharp corners such as cubes or circular cylinders, where, as is well-known, the (unmodified) regularized nodal finite element method fails. 2. Mathematical formulation Under harmonic excitation, the governing differential equation is given as

 ∇×

 k2 ∇ × E − E = −iωj , µ µ 1

(1)

√

√

where E is the electric field, µ is the magnetic permeability, i = −1, j is the current density, k = k0 µr ϵr is the wave √ number of the medium, k0 = ω2 /c 2 is the wave number of vacuum, ω is the frequency of excitation, c = 1/ ϵ0 µ0 is the speed of light in vacuum, ϵr = ϵ/ϵ0 and µr = µ/µ0 are the relative permittivity and relative permeability, and ϵ0 and µ0

1366

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

are the permittivity and permeability for vacuum. In the absence of charges, the electric field is subject to the constraint

∇ · (ϵ E ) = 0.

(2)

For exterior electromagnetic problems, we use the following first-order absorbing boundary condition on the truncation boundary Γ∞ [1]: n × (∇ × E ) = −ikn × (n × E ) ,

(3)

where n is the unit normal to the boundary Γ∞ . 2.1. Conventional potential formulation Similar to the strategy followed in [31,37], we replace E by A + ∇ψ , where A and ψ are vector and scalar potentials, i.e., E = A + ∇ψ, H =

i

µω

(4a)

∇ × A.

(4b)

The above potential-based formulation ensures tangential continuity of the electric field at the interface in an inhomogeneous medium, and also ensures that the condition ∇ · (µH ) = 0 is automatically satisfied. The governing differential equation (1) and the constraint equation (2) in terms of the potentials is given by

 ∇×

 k2 ∇ × A − (A + ∇ψ) = −iωj , µ µ 1

(5a)

∇ · (ϵ A) + ∇ · (ϵ∇ψ) = 0,

(5b)

while Eq. (3) gets modified to n × (∇ × A) = −ikn × (n × A) ,

(6)

where we have neglected the ∇ψ term since it is of the order 1/R lower compared to the term that is retained (this was also verified numerically for all the examples), where R is the radius of the truncation surface Γ∞ . For the scalar field ψ , we impose the following far-field condition:

∇ψ · n = −

ψ R

.

(7)

2.1.1. Variational formulation The variational statements corresponding to Eqs. (5a) and (5b) which incorporate a regularization term corresponding to the Coulomb gauge condition ∇ · A = constant, and the first-order absorbing boundary conditions at Γ∞ are

 2  k 1 (∇ × Aδ ) · (∇ × A) dΩ − Aδ · (A + ∇ψ) dΩ + (∇ · Aδ )(∇ · A) dΩ Ω µ Ω µ Ω µ    k ¯ dΓ − iω +i Aδ · H Aδ · j dΩ , (Aδ × n) · (A × n) dΓ = −iω µ Γ∞ Γh Ω       ϵ ϵ∇ψδ · ∇ψ dΩ + ϵ∇ψδ · A dΩ + ψδ ψ − ϵψδ A · n dΓ − ϵψδ (∇ψ · n + A · n) dΓ = 0. 

1

Ω

Ω

Γ∞

R

(8a) (8b)

Γh

Γh denotes the part of the boundary Γ where H × n is specified as H¯ . On a conducting surface, we specify the essential boundary condition E × n = 0. The variables A and Aδ belong to the function space H (curl, Ω ), while ψ and ψδ belong to H 1 (Ω ). 2.1.2. Finite element formulation We discretize A, Aδ , ψ , ψδ as

ˆ, A = NA

ˆ δ, Aδ = N A

∇ × A = BAˆ , ∇ × Aδ = BAˆ δ ,

∇ · A = Bp Aˆ , ∇ · Aδ = Bp Aˆ δ ,

ˆ ψ = Nψ ψ, ψδ = Nψ ψˆ δ ,

ˆ ∇ψ = Bψ ψ, ∇ψδ = Bψ ψˆ δ ,

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

1367

where N1 0 0 N2 0 0 · · · N = 0 N1 0 0 N2 0 · · · ,  0 0 N1 0  0 N2 · · · Nψ = N1 N2 N3 · · · ,   ∂ N2 ∂ N2 ∂ N1 ∂ N1 0 − ··· 0 − ∂z ∂y ∂z ∂y     ∂ N1 ∂ N2 ∂ N2   ∂ N1 B =  0 − 0 − · · · ,   ∂z ∂x ∂z ∂x   ∂N ∂N ∂ N2 ∂ N2 1 1 − 0 − 0 ··· ∂y ∂x ∂y ∂x   ∂ N1 ∂ N1 ∂ N1 ∂ N2 ∂ N2 ∂ N2 ··· , Bp = ∂x ∂y ∂z ∂x ∂y ∂z  ∂N ∂N  1 2



Bψ



(9)

···  ∂x ∂x   ∂N ∂N  2  1  · · · . =   ∂y ∂y    ∂ N1 ∂ N2 ··· ∂z ∂z

We express the unit normal n to a surface parametrized by natural coordinates (ξ , η) as ∂x × ∂ξ =   ∂x  ∂ξ ×

n1 n2 n3

  n=

∂x ∂η

.

∂x  ∂η 

It follows that

ˆ, A · n = nT N A

ˆ δ, Aδ · n = nT N A

ˆ, A × n = nTmat N A

ˆ δ, Aδ × n = nTmat N A

where nmat is the skew-symmetric matrix given by

−n 3

0 n3 −n 2

 nmat =

0 n1

n2 −n 1 . 0



We can now write Eqs. (8) as



KAA Kψ A

KAψ Kψψ

    ˆ F A = A , Fψ ψˆ

(10)

where

 KAA =

 1  T B B + BTp Bp − k2 N T N dΩ + i

µ

Ω

k2

 KAψ = −

 Kψ A =

Ω

Ω

µ

k Γ∞

µ

N T nmat nTmat N dΓ ,

N T B ψ dΩ ,

ϵ BTψ N dΩ −

 Γ∞

ϵ NψT nT N dΓ −



ϵ T ϵ BTψ Bψ dΩ + Nψ Nψ dΓ − R Ω Γ∞    ¯ dΓ + FA = −iω NT H N T j dΩ ,

Kψψ =

Fψ = 0.



Γh

Γh Γh

ϵ NψT nT N dΓ , ϵ NψT nT Bψ dΓ ,

Ω

ˆ and ψˆ , we obtain E and H using Eq. (4) as After solving for A ˆ + Bψ ψ, ˆ E = NA i ˆ. BA H = µω

1368

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

2.2. Amplitude formulation According to the expansion theorem of Wilcox [38,39], we can represent E as ∞ e−ik|x|  an (θ , φ)

E=

|x|

n =0

|x|n

,

(11)

where |x| is the distance of a point from the origin, or, alternatively, the radial coordinate in the spherical coordinate system (r , θ, φ). As in [40–43], we use this expansion to separate out the rapidly-varying part e−ik|x| /|x|, so that the finite element interpolation functions only have to capture a gently-varying part as discussed in the following section. 2.2.1. Variational formulation In the amplitude formulation that we now propose, we assume (A, ψ) and their variations (Aδ , ψδ ) as A = Aδ =

ψ = ψδ =

1

|x| 1

|x| 1

|x| 1

|x|

¯ (x)e−ik|x| , A ¯ δ (x)eik|x| , A (12)

¯ x)e−ik|x| , ψ( ψ¯ δ (x)eik|x| .

¯ and ψ¯ only have to capture a Since the rapidly-varying part has been separated out, the interpolation functions for A relatively gently-varying function. Note that as in [42], the main idea is to replace the spherical radial coordinate by |x|, which allows the use of a Cartesian framework (i.e., no angular coordinates, which cause problems of nonuniqueness, need to be discretized), and hence, allows the use of standard Lagrange interpolation functions, standard Gaussian quadrature and so on. In fact, the same data structure used for the conventional formulation is used in our implementation of the amplitude ¯ and ψ¯ as formulation. The physical variables of interest, namely, E and H , can be expressed in terms of A    ∇ ψ¯ ik 1 ¯ e−ik|x| , − + ( x ψ) |x| |x| |x|2 |x|3     ¯ i ∇×A ik 1 ¯ ) e−ik|x| . H = − + ( x × A µω |x| |x|2 |x|3 

E =

1

¯+ A

(13)

Using Eqs. (12), we have

   ∇ · A¯ ik 1 ¯ − + 3 (x · A) e−ik|x| , |x| |x|2 |x|     ¯ ik 1 ∇ · Aδ ¯ δ ) eik|x| , + − ( x · A |x| |x|2 |x|3     ∇ × A¯ ik 1 ¯ ) e−ik|x| , − + ( x × A |x| |x|2 |x|3     ∇ × A¯ δ ik 1 ¯ δ ) eik|x| , + − ( x × A |x| |x|2 |x|3     ∇ ψ¯ ik 1 ¯ − + 3 (xψ) e−ik|x| , |x| |x|2 |x|     ¯ ∇ ψδ ik 1 ¯ δ ) eik|x| . + − ( x ψ |x| |x|2 |x|3 

∇·A = ∇ · Aδ = ∇×A = ∇ × Aδ = ∇ψ = ∇ψδ =

Substituting Eqs. (12) and (14) into Eqs. (8a) and (8b), we have



  ¯ δ ) · (∇ × A¯ ) + ik − 1 (x × A¯ δ ) · (∇ × A¯ ) (∇ × A µ |x|2 |x|3 |x|4    2   ik 1 ¯ δ ) · (x × A¯ ) + k + 1 (x × A¯ δ ) · (x × A¯ ) dΩ − + (∇ × A |x|3 |x|4 |x|4 |x|6      2 k 1 1 ¯ δ · A¯ dΩ + ¯ δ )(∇ · A¯ ) + ik − 1 (x · A¯ δ )(∇ · A¯ ) A (∇ · A − 2 |x|2 |x|3 |x|4 Ω µ|x| Ω µ 1

Ω



1

(14)

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391





k2



ik

|x|

3

+

4

(∇ · A¯ δ )(x · A¯ ) +

1369



1 + 6 (x · A¯ δ )(x · A¯ ) dΩ |x| |x|4 |x|     2 ¯  ¯ k Aδ · ∇ ψ 1 ik k ¯ δ · x)ψ¯ dΩ + i − + ( A (A¯ δ × n) · (A¯ × n) dΓ − 2 |x|2 |x|4 |x|3 Ω µ Γ∞ µ|x|   ik|x|  eik|x|  e ¯ δ × n · H¯ dΓ − iω ¯ δ · j dΩ , = iω A A Γh |x| Ω |x|       1 ik 1 ik 1 ¯ δ · ∇ ψ¯ − ¯ δ · x)ψ¯ + ϵ ∇ ψ + (∇ ψ − ψ¯ δ (x · ∇ψ) |x|2 |x|3 |x|4 |x|3 |x|4 Ω  2       ¯ ¯  k 1 ∇ ψδ · A ik 1 ¯ δ ψ¯ dΩ + ¯ δ x · A¯ dΩ + + ψ ϵ + − ψ |x|2 |x|4 |x|2 |x|3 |x|4 Ω         ψ¯ δ (∇ ψ¯ · n) ϵ ψ¯ δ ψ¯ ik 1 ¯ · n dΓ − ¯ δ ψ¯ dΓ ϵ + − A − + ( x · n ) ψ 2 R |x|2 |x|3 |x|4 Γh Γ∞ |x|  ϵ ψ¯ δ A¯ · n − dΓ = 0. |x|2 Γh

−

1



(15)

(16)

It may appear that the load terms in Eq. (15) are now oscillatory. However, the j term is nonzero only in scattering problems involving dielectrics on the domain occupied by the dielectric. Since we use conventional elements in that region, no problem arises. Similarly, the eik|x| term on Γh also does not cause a problem since this term varies gently on the surface of the radiator, e.g., if the surface is spherical, this term is a constant. To summarize, the highly oscillatory terms have been eliminated so ¯ and ψ¯ have to capture only relatively gently-varying functions. that the interpolation functions for A If the origin is part of the domain, then we obviously cannot use Eqs. (12). In such a case, we divide the domain into two parts using a spherical surface Γ1 with center at the origin and having radius r1 . Conventional elements are used inside this spherical surface and AMP elements outside it (up to the truncation surface Γ∞ ). Now instead of Eq. (12), we use the following interpolations: A= Aδ =

ψ= ψδ =

r1

|x| r1

|x| r1

|x| r1

|x|

¯ (x)eik(r1 −|x|) , A ¯ δ (x)e−ik(r1 −|x|) , A ¯ x)eik(r1 −|x|) , ψ( ψ¯ δ (x)e−ik(r1 −|x|) .

Note that the above fields are defined in a manner so as to maintain continuity at the interface between the conventional and the AMP elements. Corresponding to the above modified fields, Eqs. (15) and (16) get modified as follows. All terms ¯ ψ¯ and their derivatives in Eqs. (15) and (16) are multiplied by r1 eikr1 . Similarly, all the terms having A¯ δ , ψ¯ δ involving A,      ¯ , A¯ δ , A¯ , ψ¯ δ , ψ, ¯ A¯ δ , ψ, ¯ ψ¯ δ or and their derivatives, are multiplied with r1 e−ikr1 . Thus, the terms having combinations A corresponding derivatives, are multiplied by r12 . For problems where the origin is not part of the domain such as radiation problems or scattering problems involving conducting surfaces, we use AMP elements throughout the domain. 2.2.2. Finite element formulation ¯ A¯ δ , ψ¯ , ψ¯ δ and the fields associated with them as follows: We discretize A,

¯ = N Aˆ¯ , A

¯ δ = N Aˆ¯ δ , A

∇ × A¯ = BAˆ¯ , ∇ · A¯ = B Aˆ¯ ,

∇ × A¯ δ = BA¯ˆ δ , ∇ · A¯ = B Aˆ¯ ,

p

¯ · n = nT N Aˆ¯ , A ¯ · x = xˆ T N Aˆ¯ , A ¯ × n = nTmat N Aˆ¯ , A ¯ × x = xTmat N Aˆ¯ , A ˆ¯ ψ¯ = Nψ ψ,

δ

p

δ

¯ δ · n = nT N A¯ˆ δ , A ¯ · x = xˆ T N Aˆ¯ , A δ

δ

¯ δ × n = nTmat N A¯ˆ δ , A ¯ × x = xT N A¯ˆ , A δ

mat

ψ¯ δ = Nψ ψˆ¯ δ ,

δ

1370

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

ˆ¯ ∇ ψ¯ = Bψ ψ,

∇ ψ¯ δ = Bψ ψˆ¯ δ , ∇ ψ¯ δ · x = xˆ T Bψ ψ¯ˆ δ , ∇ ψ¯ · n = nT B ψˆ¯ ,

ˆ¯ ∇ ψ¯ · x = xˆ T Bψ ψ, ˆ¯ ∇ ψ¯ · n = nT B ψ, ψ

δ

ψ

δ

where x y , z

  xˆ =

−z

0 z −y

 xmat =

y −x . 0



0 x

Then we can write Eqs. (15) and (16) in finite element form as



KAA Kψ A

KAψ Kψψ

 ˆ    ¯ F A = A , ˆ F ψ ¯ ψ

(17)

where

 KAA =

1



µ 

1

k2

|x|4 

−

1 4

|x| 

+i KAψ = −



1

+



ik

+

3

|x|

k

µ|x|2  2

1

|x|4 T



ik

−

|x|3 

xmat xTmat N

BTp xˆ T N +



k2 4

|x|

N T xmat B +



dΩ +

+

1

 6

|x|

1

|x|4 

µ

Ω

1



1

|x|

+

ik

|x|3

BT B 2 p p

N T xˆ xˆ T N





BT xTmat N

 +

dΩ −

ik

|x|3

−

1

|x|4

k2

 Ω

µ|x|2

ϵ

1

µ 

|x|2 ik

N T xˆ Bp

N T N dΩ



T

N Bψ −

1

|x|4

+

ik

|x|3



T



N xˆ Nψ

dΩ ,

ϵ T T N n N dΓ − 3 4 2 ψ |x| |x| |x| Ω Γ∞ |x|       ik ik 1 1 T 1 T ˆ B + B x N + − NψT xˆ T Bψ Kψψ = ϵ B − ψ ψ 2 ψ ψ 4 3 3 4 | x | | x | | x | | x | | x | Ω  2    k 1 ϵ T T + + N N Nψ dΓ N d Ω + ψ ψ |x|2 |x|4 | x |3 ψ Γ∞      1 1 ik T T T − ϵ N + (ˆ x · n ) N n B − N dΓ , ψ ψ ψ |x|2 ψ |x|3 |x|4 Γh    ik|x| e eik|x| T T ¯ N H dΓ + N j dΩ , FA = −iω Ω |x| Γh |x| 

Kψ A =



N T nmat nTmat N dΓ ,

k Ω

N

|x|6

Γ∞





|x|

Ω

+

BT B + 2

−

1



NψT xˆ T N +

1

BT N 2 ψ



dΩ −



 Γh

ϵ T T N n N dΓ , |x|2 ψ

Fψ = 0.

ˆ¯ and ψˆ¯ , we recover E and H using Eqs. (13) as After solving for A    ˆ¯ Bψ ψ ik 1 ˆ ˆ ¯ ¯ NA + − + 3 (xNψ ψ) e−ik|x| , E= |x| |x| |x|2 |x|    ˆ¯  ik i BA 1 ˆ ¯ H = + + 3 (xmat N A) e−ik|x| . µω |x| |x|2 |x| 

1

3. Axisymmetric formulation An axisymmetric formulation, if possible, is obviously more cost-effective and accurate than a full three-dimensional analysis. Let (r , θ , z ) denote a cylindrical coordinate system. An axisymmetric formulation is possible if the excitation and

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

1371

geometry are such that Er = Ez = 0 and Eθ = Eθ (r , z ). Corresponding to this, we have Hθ = 0, and Hr and Hz are function of r and z. Under the above conditions, ∇ × E has following simple expression in the r–θ –z system:

∂ Eθ ∇×E =− eˆ r + ∂z



∂ Eθ + r ∂r

Eθ



eˆ z .

(18)

Since the only nonzero component of E is tangential to the mesh in the r–z plane, there is no discontinuity in the normal component of E at the surface of an inhomogeneity, and ψ can be set to zero, i.e., E = A in this case. Furthermore, since ∇ · E = 0 we do not need to include the regularization term. A general solution to Eq. (1) (assuming j = 0) obtained using separation of variables is given by Eθ =

∞ 

 



h(n2) k r 2 + z 2 Pn(1)

 √

n=1

z



r 2 + z2

.

(19)

Note that

∂ Eθ Eθ = lim . r → 0 ∂r r

lim

r →0

(20)

3.1. Conventional formulation The various terms are discretized as follows: Eθ = N Eˆ θ ,

ˆ, H = N¯ H ∇ × E = BEˆ θ , ˆ Eˆ θ , E × n = nN where N = N1 N2 · · · ,  N1 0 N2 0 · · · ¯ N = , 0 N1 0 N2 · · ·  ∂ N1 ∂ N2





−

−

···



 ∂z ∂z , ∂ N1 N2 ∂ N2 + + ··· ∂r  r  ∂r r nz ˆ = n . −n r

B = 



(21)

N1

Hence, we can write the finite element equations as Kp Eˆ θ = f , where

 Kp =

 1  T rB B − k2 rN T N dΩ + i



rk

µ Γ∞ µ    T ¯ T f = −iω rN H dΓ + rN jθ dΩ . Ω

Γh

N T N dΓ ,

Ω

After obtaining Eˆ θ , we can calculate H as H =

i

µω

BEˆ θ .

From Eq. (20), it follows that for points on the axis r = 0, we should use the following B matrix instead of that given by Eq. (21):

 ∂N 1 −  ∂z B= ∂N 2

1

∂r

∂ N2 ∂z ∂ N2 2 ∂r

−

··· ···

  .

(22)

1372

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

3.2. Amplitude formulation In place of Eq. (12), we now have Eθ =

(Eθ )δ =

1

Gθ e−ik|x| ,

|x| 1

(Gθ )δ eik|x| .

|x|

With G = Gθ eθ , the discretizations are given by Gθ = N Gˆ ,

∇ × G = BGˆ , ˆ Gˆ , G × n = nN G × x = xˆ N Gˆ , where

 xˆ =



z

−r

.

The finite element equations are given by K Gˆ = f , where

   ik ik 1 T T ˆ − N + BT xˆ N x B + 4 3 4 3 µ | x | | x | | x | | x | | x | Ω   2    1 irk k rk2 T T T ˆ ˆ + N N N d Ω + N T N dΓ , + x x N d Ω − 2 2 |x|4 |x|6 µ| x | µ| x | Γ∞ Ω    reik|x| T reik|x| T ¯ dΓ + f = −iω N H N jθ dΩ . Γh |x| Ω |x| 

r

K =



1

BT B + 2



1

ˆ we can recover Eθ and H using the relations After obtaining G, Eθ = H =

1

|x|

ˆ −ik|x| , N Ge

i

µω



BGˆ

|x|

 +

ik

|x|2

+

1

|x|3



 (ˆxN Gˆ ) e−ik|x| .

4. Modification of the formulation for scattering from conducting surfaces/dielectrics When an incident wave Einc impinges on a conducting body then the scattered field Escat in addition to satisfying all the Maxwell equations also satisfies the following equation on the conducting surface

(Escat + Einc ) × n = 0.

(23)

We solve for Escat in the domain outside the conducting body. The total field E is obtained by adding Escat to Einc . We implement the above constraint using a Lagrange multiplier technique in the elements adjacent to the scattering boundary. The Lagrange multiplier degrees of freedom are suppressed for elements not adjacent to the scattering boundary. 4.1. Conventional formulation We express Escat as A + ∇ψ , so that Eq. (23) becomes

(A + Einc ) × n = 0

(24)

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

1373

with all the ψ degree of freedoms prescribed to be zero over the scattering boundary. Eq. (10) gets modified to [37] 1



     

KAA

KAψ

bKψ A

Kψψ

−1 1

1

N T t T Nλ dξ dη

−1

0

1



Nλ tN dξ dη T

−1

−1



0

0



FA Fψ

    Aˆ    ψˆ  =      λ −

1

−1





1

−1

NλT t Eˆ inc dξ dη

 , 

(25)

where Eˆ inc are nodal values of Einc on the conducting surface, t is a 2 × 3 matrix containing along its rows two linearly ¯ (ξ , η) are natural coordinates that parametrize independent vectors t1 and t2 that are both perpendicular to the normal n, the surface of the element, and

¯= n

∂x ∂x × , ∂ξ ∂η  N1 0

Nλ =

0 N1

N2 0

0 N2

 

∂x Note that the surface Jacobian  ∂ξ

 ··· . ···  ∂x  × ∂η  gets canceled and hence does not appear in Eq. (25). Since the Lagrange multiplier

technique cannot handle abrupt changes in slope of the continuum (e.g., edges of a cube), the entire A vector along edges with a discontinuous normal n is set to zero, and the corresponding Lagrange multipliers are also suppressed. 4.2. Amplitude formulation Using the form of A given in Eq. (12), Eq. (17) gets modified to





     

KAA

KAψ

Kψ A

Kψψ

1

−1 1

1



−1

e−ik|x|

−1

|x|

NλT tN dξ dη

0



1

−1

eik|x|

|x|



 N t Nλ dξ dη   ˆ ¯  A    0  ψˆ¯  =    T T

0



λ

FA Fψ 1





1

− −1

−1

NλT t Eˆ inc dξ dη

  . 

(26)

4.3. Scattering from dielectric bodies When an incident wave Einc impinges on a dielectric body, then Eq. (23) is modified to

(Escat + Einc ) × n = Et × n,

(27)

where Escat is the scattered field in the outer domain (outside the dielectric) and Et is the transmitted field in the inner domain (inside the dielectric). Each of these fields individually satisfies the Maxwell equations. Following the scattered field formulation [1], similar to the outer domain, we split the field inside the dielectric Et into two parts Escat and Einc . The advantage of this splitting is that one can impose tangential continuity between the inner and outer scattered fields at the surface of the dielectric. Assuming j = 0, Eq. (1) can be written on the inside domain as

     k2 1 k2 ∇ × Escat − Escat = − ∇ × ∇ × Einc − Einc . µ µ µ µ Writing Escat as A + ∇ψ , we get       1 k2 1 k2 ∇× ∇ × A − (A + ∇ψ) = − ∇ × ∇ × Einc − Einc . µ µ µ µ 

∇×

1

(28)

(29)

Comparing Eqs. (5a) and (29), and using the fact that Einc satisfies Eq. (1) with k = k0 , µ = µ0 and j = 0, we see that scattering from dielectric bodies can be treated as an excitation problem with the following equivalent current density: jeqv = iωϵ0

=0



ϵr −

1

µr



Einc ,

(inner domain) (outer domain).

After obtaining Escat for the inner domain, Et can be obtained by adding Einc to Escat . As already mentioned, we use only conventional elements within the dielectric region, and hence there is no amplitude formulation for the inner domain. We use amplitude formulation only in the outer domain where the equivalent current density is zero, and hence, they can be used using the same formulation as discussed earlier.

1374

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

4.4. Scattering from domains with sharp corners and edges It is well-known that the nodal finite element based on the A–ψ potential formulation fails to capture the singular fields in case the domain is non-convex and has sharp corners and edges. Hence, some modification of the formulation is required in order to capture these singular fields accurately. Costabel and Dauge [33] use a strategy where the penalty factor associated with the regularization term varies with the distance from the corner or edge, while Otin [34,35] uses a zero penalty in the layers adjacent to corners and edges. We use a slightly modified version of the latter strategy where we surround the entire geometry by a thin layer of AMP elements where zero penalty is used, and use the usual AMP elements with nonzero penalty in the remaining domain. Our numerical experiments on a variety of problems shows that the thickness of the zero penalty layer is not very critical as long as it is small compared to the dimensions of the scatterer. We have found that even if a smooth scatterer such as a sphere is surrounded by a thin zero penalty layer (although this is not required), the results are almost the same as with the regular amplitude formulation. This is significant since if we have a scatterer which has a combination of smooth surfaces and sharp corners, then we can simply surround the entire body with a thin layer of zero penalty elements without having to identify particular zones where such elements need to be used. Note that this strategy is applied only to problems involving domains that have sharp edges/corners and are non-convex. For scattering from a dielectric cube for example, where the domain has sharp corners but is convex, we do not need to use this thin layer of zero penalty elements. 5. Numerical examples In this section, we present a wide variety of examples involving both radiation and scattering, involving both convex and non-convex geometries with sharp corners, and involving both conductors and dielectrics, to demonstrate the good performance of the amplitude formulation in comparison with the conventional strategy based on an ABC condition as described in Section 2.1. Although there is a significant improvement for both scattering and radiation problems, the improvement is much more significant for the latter, since in the former there is a significant variation of the fields along the polar and azimuthal directions, while the dominant term in the Wilcox asymptotic expansion (on which the AMP element is based) is in terms of the radial coordinate alone. For some scattering problems, we have also compared our solutions with the HFSS solution [44], which uses a combination of edge elements and a PML. Both yield accurate answers when compared against an analytical solution (where available), or results that are in agreement on more complicated problems that involve non-convex geometries. A precise comparison in terms of computational cost with the PML is difficult since they are based on different kinds of finite elements (nodal versus edge), but our experience is that both strategies use roughly the same total number of degrees of freedom (which is the primary factor that decides the computational cost). Thus, the current proposed method may be thought of an alternative to the PML method, albeit formulated under a different framework (nodal finite elements, although as already mentioned, it can be implemented with edge elements also), and based on a different logic (the Wilcox expansion). Our experience with acoustics problems shows that a conventional strategy based on an ABC can yield an accurate magnitude of the pressure although the real and imaginary parts themselves have significant errors, especially at high frequencies. Hence, in the current work, we have mainly plotted the real and imaginary parts of the field variables instead of say the radar cross section (RCS), although in one example, we present the RCS plot also; if the real and imaginary parts of the field variables are accurate, then it follows that the RCS which depends on them is also accurate. In the examples, by ‘far field’ we mean the field at a large finite distance from the scatterer, and not the field at infinity. 5.1. Radiation problem based on an axisymmetric analytical solution An axisymmetric analytical solution to the Maxwell equations is given by Er = Eθ = 0, Eφ = h(n2) (kr ) Pn(1) (cos θ ) , Hr = −

(2)  inhn (kr )  (1) (1) Pn (cos θ ) cos θ − Pn+1 (cos θ ) , µωr sin θ

(30)

(1)  iPn (cos θ )  (2) Hθ = − c (n + 1)h(n2) (kr ) − ωrhn+1 (kr ) , µωrc Hφ = 0, (2)

(m)

where hn denote the spherical Hankel functions of the second kind, and Pn In the limit as θ tends to 0 or π , we have Eφ = Hθ = 0 and Hr |θ→0 =

−in(n + 1) (2) hn (kr ) , µωr

denote the associated Legendre polynomials.

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

1375

Table 1 Comparison of near- and far-field response for different k0 a and k0 R∞ between conventional and proposed strategies (axi-symmetric case). k0 a = 1, k0 R∞ = 5, Mesh = 4 × 4 (r × θ), θ = π/4 Field variables

Eφ (real) Eφ (imag)

k0 r = 1.5

k0 r = 4

Analytical

Conventional

Current

Analytical

Conventional

Current

−0.2801 −0.4925

−0.2687 −0.4706

−0.2809 −0.4745

−0.0821

−0.0826

−0.0826

0.1627

0.1560

0.1605

k0 a = 10, k0 R∞ = 20, Mesh = 1 × 4 (r × θ), θ = π/4 Field variables

k0 r = 10 Analytical

Eφ (real) Eφ (imag)

−0.0555 0.0444

k0 r = 15 Conventional 0.0255 0.0533

Current

Analytical

Conventional

Current

−0.0549

−0.0379 −0.0283

−0.00001 −0.0091

−0.0379 −0.0283

Conventional

Current

0.0448

k0 a = 100, k0 R∞ = 300, Mesh = 1 × 6 (r × θ), θ = π/4 Field variables

k0 r = 100 Analytical

Eφ (real) Eφ (imag)

0.0061 0.0035

k0 r = 200 Conventional 0.0001

−0.0002

Current

Analytical

0.0061 0.0035

0.0017 0.0031

−0.00001 0.00001

0.0017 0.0031

k0 a = 100, k0 R∞ = 105, Mesh = 1 × 6 (r × θ), θ = π/4 Field variables

k0 r = 100 Analytical

Eφ (real) Eφ (imag)

0.0061 0.0035

k0 r = 105 Conventional 0.0106 0.0105

(a) Eφ (real).

Current 0.0061 0.0035

Analytical

Conventional

Current

−0.0016

−0.0084

−0.0016

0.0066

0.0047

0.0065

(b) Eφ (imaginary).

Fig. 1. Variation of Eφ as a function of θ for k0 r = 100, 200, φ = π/6, k0 a = 100, k0 R∞ = 300 and 1 × 6 (r × θ ) mesh for the axi-symmetric problem.

Hr |θ→π =

−(−1)n in(n + 1) (2) hn (kr ) . µωr

Consider a sphere of radius a in a spherical domain of radius R∞ . We prescribe H as per Eq. (30) with n = 1 on the surface r = a of the sphere. The results for k0 a = 1, 10, 100 are presented in Table 1. As expected, for the low-frequency case k0 a = 1, the AMP element is only marginally better than the conventional element (since both are close to the analytical solution), but at the higher frequencies k0 a = 10 and 100, the AMP element performs significantly better. Fig. 1 shows the variation of Eφ with θ , where again we see that the amplitude formulation performs significantly better.

1376

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

5.2. Transverse Magnetic (TM) solution of Maxwell equations The transverse magnetic (along r) solution of the Maxwell equations is given by iωn(n + 1) (2) hn (kr ) sin(mφ)Pn(m) (cos θ ) , kr Eθ = t1 t2 t3 , Er = −

(m)

Eφ = −

im cos(mφ)Pn

(cos θ) t3

r sin θ

Hr = 0,

,

(31a) (31b) (31c) (31d)

ωm (2) Hθ = hn (kr ) cos(mφ)Pn(m) (cos θ ) , c µ sin θ ω sin(mφ) (2) Hφ = hn (kr ) t2 , c µ sin θ

(31e) (31f)

where t1 =

i sin(mφ) r sin θ

, (m)

t2 = (n + 1) cos θ Pn(m) (cos θ) + (m − n − 1)Pn+1 (cos θ ) , (2)

t3 = c (n + 1)h(n2) (kr ) − ωrhn+1 (kr ) . We have the following limiting values for m = n = 1:



(2)

2 (kr ) − h(12) (kr ) , kr



(2)

2 (2) h (kr ) , kr 1

Eθ |θ→0 = −iω sin φ h2

Eφ |θ →0 = −iω cos φ h2 (kr ) −





−iω(1 + ikr ) cos φ [cos(kr ) − i sin(kr )] , µck2 r 2 iω(1 + ikr ) sin φ [cos(kr ) − i sin(kr )] , = µck2 r 2 = −Eθ |θ →0 , = Eφ |θ→0 , = Hθ |θ→0 , = −Hφ |θ→0 .

Hθ |θ →0 = Hφ |θ →0 Eθ |θ →π Eφ |θ →π Hθ |θ →π Hφ |θ →π

We take the domain to be a hollow sphere with inner radius a and outer radius R∞ . We implement the Sommerfeld condition at the outer radius, while at the inner radius we prescribe the magnetic field H as per Eq. (31) (for the case m = n = 1). We choose k0 a = 1, 5, 10 and 100 for our simulations. Table 2 presents the near- and far-field E values for three different combinations of k0 a and k0 R∞ where θ = π /4 and φ = π /4. For all cases considered, the AMP element exhibits better performance than the conventional element. With the increase of k0 a, the superiority of the amplitude formulation, both for the near- and the far-field responses, becomes more evident. In fact, as the frequency is increased, the amplitude formulation continues to yield an extremely accurate prediction even at very high frequencies such as k0 a = 100 (see Fig. 2), where the conventional strategy using the same mesh breaks down completely (as seen from the ‘flat’ curves). 5.3. Transverse Electric (TE) solution of the Maxwell equations The transverse electric (along r) solution is given by ETE = −µHTM , HTE = ϵ ETM ,

(32)

where ETM and HTM are given by Eq. (31). In this case, we again consider the same spherical domain as in the TM case, and prescribe H on the inner surface of the domain as per Eq. (32) (with m = n = 1). To show that our amplitude formulation works even if the medium is not vacuum, we choose ϵr = 1.5 and µr = 2.0. We carry out simulations for k0 a equal to 1, 10 and 100. Table 3 again shows the

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) Ex at φ = π/4 (real).

(b) Ex at φ = π/4 (imaginary).

(c) Ey at θ = π/4 (real).

(d) Ey at θ = π/4 (imaginary).

(e) Ez at φ = π/4 (real).

(f) Ez at φ = π/4 (imaginary).

Fig. 2. Variation of electric field components for k0 a = 100, k0 R∞ = 500 and 8 × 8 × 16 (r × θ × φ) mesh for the TM case.

1377

1378

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

Table 2 Comparison of near- and far-field responses for different k0 a and k0 R∞ between conventional and proposed strategies (TM case). k0 a = 5, k0 R∞ = 45, Mesh = 6 × 4 × 8 (r × θ × φ), θ = π/4, φ = π/4 Field variables

Ex Ex Ey Ey Ez Ez

(real) (imag) (real) (imag) (real) (imag)

k0 r = 15

k0 r = 35

Analytical

Conventional

Current

Analytical

Conventional

Current

−2.199e7 −1.224e7

1.819e7 −1.608e7 −4.970e7 3.024e7 2.319e7 −1.284e7

−2.226e7 −1.212e7

−0.927e7

−0.736e7

−0.941e7

0.541e7 2.891e7 −1.404e7 −1.310e7 0.765e7

0.573e7 2.129e7 −1.525e7 −1.060e7 0.831e7

0.547e7 2.862e7 −1.384e7 −1.297e7 0.777e7

Conventional

Current

5.798e7 4.743e7 −3.110e7 −1.731e7

5.705e7 4.726e7 −3.142e7 −1.688e7

k0 a = 10, k0 R∞ = 50, Mesh = 8 × 8 × 16 (r × θ × φ), θ = π/4, φ = π/4 Field variables

Ex Ex Ey Ey Ez Ez

(real) (imag) (real) (imag) (real) (imag)

k0 r = 15

k0 r = 45

Analytical

Conventional

Current

−4.398e7 −2.448e7

−4.649e7 −6.015e7

−4.398e7 −2.443e7

11.595e7 9.487e7 −6.220e7 −3.463e7

9.611e7 17.917e7 −6.572e7 −8.508e7

11.579e7 9.489e7 −6.222e7 −3.455e7

Analytical 0.780e7

−1.474e7 −2.595e7 4.273e7 1.103e7 −2.085e7

−2.049e7 −1.206e7 6.027e7 3.967e7 −2.897e7 −1.702e7

0.779e7

−1.474e7 −2.594e7 4.266e7 1.102e7 −2.087e7

k0 a = 100, k0 R∞ = 500, Mesh = 8 × 8 × 16 (r × θ × φ), θ = π/4, φ = π/4 Field variables

Ex Ex Ey Ey Ez Ez

(real) (imag) (real) (imag) (real) (imag)

k0 r = 125

k0 r = 475

Analytical

Conventional

Current

Analytical

Conventional

Current

4.814e7 3.582e7 −14.208e7 −11.051e7 6.808e7 5.066e7

−0.159e7

4.876e7 3.548e7 −14.378e7 −10.932e7 6.902e7 5.013e7

−1.280e7

0.024e7 0.019e7 −0.064e7 −0.033e7 0.026e7 0.016e7

−1.285e7

0.276e7 0.196e7 −0.342e7 −0.191e7 0.333e7

0.925e7 3.855e7 −2.753e7 −1.810e7 1.308e7

0.931e7 3.867e7 −2.770e7 −1.818e7 1.319e7

Table 3 Comparison of near- and far-field responses for different k0 a and k0 R∞ between conventional and proposed strategies (TE case). k0 a = 1, k0 R∞ = 9, Mesh = 4 × 4 × 8 (r × θ × φ), θ = π/4, φ = π/4 Field variables

Ex Ex Ez Ez

(real) (imag) (real) (imag)

k0 r = 2

k0 r = 8

Analytical

Conventional

Current

Analytical

Conventional

Current

3.0298e−1 −2.0886e−1 −2.1424e−1 1.4768e−1

2.2836e−1 −1.3764e−1 −1.6177e−1 0.8780e−1

2.8976e−1 −2.0750e−1 −2.0501e−1 1.4668e−1

−1.8362e−2 8.6695e−2 1.2984e−2 −6.1303e−2

−6.8496e−2 5.2047e−2 5.0289e−2 −3.5200e−2

−1.6570e−2 8.3696e−2 1.1717e−2 −5.9278e−2

k0 a = 10, k0 R∞ = 50, Mesh = 8 × 8 × 16 (r × θ × φ), θ = π/4, φ = π/4 Field variables

Ex Ex Ez Ez

(real) (imag) (real) (imag)

k0 r = 15

k0 r = 45

Analytical

Conventional

Current

Analytical

Conventional

Current

−2.9821e−1 3.6554e−1 2.1087e−1 −2.5848e−1

−3.7755e−1 0.3788e−1 2.6635e−1 −0.2660e−1

−2.9879e−1 3.6561e−1 2.1124e−1 −2.5850e−1

1.3104e−1 8.6738e−2 −9.2661e−2 −6.1333e−2

−0.0920e−1 −0.0179e−2 0.0773e−2 −0.0194e−2

1.3115e−1 8.6720e−2 −9.2729e−2 −6.1290e−2

k0 a = 100, k0 R∞ = 500, Mesh = 8 × 8 × 16 (r × θ × φ), θ = π/4, φ = π/4 Field variables

Ex Ex Ez Ez

(real) (imag) (real) (imag)

k0 r = 125

k0 r = 475

Analytical

Conventional

Current

Analytical

Conventional

Current

5.4683e−1 1.4484e−1 −3.8667e−1 −1.0242e−1

0.0243e−1 0.0540e−1 −0.0172e−1 −0.0382e−1

5.5044e−1 1.3401e−1 −3.8989e−1 −0.9481e−1

−138.69e−3 −54.10e−3 98.07e−3 38.25e−3

0.0004e−3 0.0002e−3 −0.0003e−3 0.0003e−3

−140.54e−3 −52.78e−3 99.63e−3 37.41e−3

superior performance of the amplitude formulation, especially at high frequencies, both for the near- and far-field response compared to the conventional strategy for three different combinations of k0 a and k0 R∞ . The variations of Ex and Ez (which are the dominant components) against θ for k0 a = 100 are plotted in Fig. 3; we see that the amplitude formulation performs well even if the medium is not vacuum.

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) Ex (real).

(b) Ex (imaginary).

(c) Ez (real).

(d) Ez (imaginary).

1379

Fig. 3. Variation of electric field components against θ for k0 r = 200, 500, φ = π/4, k0 a = 100, k0 R∞ = 500, and a 8 × 8 × 16 (r × θ × φ) mesh for the TE case.

5.4. Scattering from a conducting sphere An electromagnetic plane wave Einc = E0 e−ikz ex is incident on a conducting sphere of radius a. We are interested in the scattered field Escat . On the spherical surface, the scattered field is subject to the boundary condition (Einc + Escat ) × n = 0. Denoting the relative permeability by µ, the scattered fields are given by [45], Escat = E0

∞   2n + 1  r o an mn + ibrn nen , (−i)n n ( n + 1 ) n =1

Hscat = −

∞   2n + 1  r e bn mn − iarn non , (−i)n ωµ0 µ n=1 n(n + 1)

kE0

where jn (ka) arn = − (2) , hn (ka)

(33)

1380

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391 Table 4 Comparison of near- and far-field responses for different k0 a and k0 R∞ for the scattering from a conducting sphere problem. k0 a = 1, k0 R∞ = 5, θ = π/4, φ = π/6 k0 r = 1.25

Field variables

Ex Ex Ey Ey Ez Ez

(real) (imag) (real) (imag) (real) (imag)

Conventional

HFSS

Current

−0.0369 −0.3944

−0.0347 −0.3652

−0.6936

0.2782 −0.3804 0.4906 −0.8418

0.2787 −0.3695 0.4932 −0.8160

0.3519 0.5944 −0.8779 0.3115 −0.7362

−0.0450 −0.3842

Analytical

Conventional

HFSS

Current

0.0513 0.0452 −0.0447 −0.0158 −0.0759 0.0241

0.0518 0.0399 −0.0419 −0.0167 −0.0736 0.0221

0.0618 0.0603 −0.0511 −0.0106 −0.0801 0.0244

0.0542 0.0425 −0.0440 −0.0166 −0.0767 0.0233

Analytical

Conventional

HFSS

Current

−0.0383

0.0492 0.2578 −0.0770 0.2308 −0.8094 −0.4358

−0.0374

−0.0048

Analytical

Conventional

HFSS

Current

−0.0353

0.0077 0.0473 0.0019 0.0730 −0.0224 −0.0521

−0.0375

−0.0322

0.0161 0.0758 −0.0397 0.0282 0.0088

0.0437 0.1007 −0.0292 0.0027 −0.0117

0.2747

−0.3761 0.4829

−0.8316

k0 r = 4.75

Field variables

Ex Ex Ey Ey Ez Ez

Analytical

(real) (imag) (real) (imag) (real) (imag)

k0 a = 10, k0 R∞ = 40, θ = π/4, φ = 2π/9 k0 r = 11.6667

Field variables

Ex Ex Ey Ey Ez Ez

(real) (imag) (real) (imag) (real) (imag)

0.7193

−0.2612 −0.1307 −0.2974 −0.0411

(real) (imag) (real) (imag) (real) (imag)

0.0442 0.0941 −0.0461 −0.0021 −0.0186

(n + 1)jn (ka) − kajn+1 (ka)

brn = −



mon men

mon men





(n + 1)h(n2) (ka) − kah(n2+)1 (ka)

= 0,

 = θ

 o

mn men

1 sin θ

h(n2) (kr ) Pn(1) (cos θ )

(2)

φ

 o nn nen

= −hn (kr )

 =

(1) d Pn (cos θ) 

d Pn

=

sin φ , cos φ

dθ



kr sin θ (1)

=



Pn(1) (cos θ )

(1)

nPn+1 (cos θ ) − (n + 1)Pn (cos θ ) cos θ

In the limit as θ → 0 and θ → π , we have

cos φ

dθ

 (n + 1)h(n2) (kr ) − krh(n2+)1 (kr )

φ

 (cos θ )

dθ



kr



(1)



   (n + 1)h(n2) (kr ) − krh(n2+)1 (kr ) d Pn(1) (cos θ )  sin φ 

θ

 o 

cos φ , − sin φ



r

nn nen



n(n + 1) (2) sin φ = hn (kr ) Pn(1) (cos θ) , cos φ kr

 o

nn nen

,

r





0.6763

−0.2144 −0.1772 −0.2399 −0.1084

k0 r = 38.3333

Field variables

Ex Ex Ey Ey Ez Ez

0.6189

−0.3038 −0.1238 −0.3670 −0.0879

sin θ

.



cos φ , − sin φ



,

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) Real.

(b) Imaginary. Fig. 4. Variation of Ex along φ for k0 r = 1.25, θ = π/4, k0 a = 1, k0 R∞ = 5 for the scattering from a conducting sphere problem.

(a) Real.

(b) Imaginary.

Fig. 5. Variation of Ex along θ for k0 r = 38.3333, φ = 2π/9, k0 a = 10, k0 R∞ = 40 for the scattering from a conducting sphere problem.



mon men



 o θ→0 nn nen

mon men



 θ→0

 o θ→π nn nen

θ→π

     − cos φ ˆ sin φ ˆ = t1 θ+ φ , sin φ cos φ      sin φ ˆ cos φ = t2 θ+ φˆ , cos φ − sin φ      cos φ ˆ sin φ ˆ = (−1)n t1 θ+ φ , − sin φ cos φ      sin φ ˆ − cos φ ˆ = (−1)n t2 θ+ φ , cos φ sin φ

where n(n + 1) (2) hn (kr ) , 2   n(n + 1) (n + 1) (2) (2) t2 = hn+1 (kr ) − hn (kr ) . 2 kr

t1 =

1381

1382

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391 Table 5 Mesh specifications for the scattering from a conducting sphere problem. k0 a

1 10

Conventional and current

HFSS

No. of elements r × θ × φ

No. of unknowns

No. of unknowns

8 × 6 × 12 9 × 12 × 18

16 698 61 138

18 750 75 108

(a) k0 r = 12, θ = 2π/9 (real).

(b) k0 r = 12, θ = 2π/9 (imaginary).

(c) k0 r = 24, φ = π/4 (real).

(d) k0 r = 24, φ = π/4 (imaginary). Fig. 6. Variation of Ex for the scattering from a conducting ellipsoid problem.

We assume E0 = 1 and the medium to be vacuum so that k = k0 , and consider the cases k0 a = 1 and k0 a = 10. In scattering problems, there is a strong dependence of the fields on θ ; since the AMP element separates out only the highly oscillatory radial part a-priori, the performance of the AMP element compared to the conventional element is not as spectacularly superior as in radiation problems. But, nevertheless, is still significantly better at higher frequencies. The results found using HFSS [44], which uses a combination of edge elements and the PML, along with the results obtained using the conventional and AMP elements are presented in Table 4 and Figs. 4 and 5. We have used (almost) equivalent meshes to ensure a fair comparison (see Table 5 for the mesh specifications). It is evident from the presented results that at the low frequency, the conventional and amplitude formulations perform equally, while HFSS underperforms. For the high frequency, AMP performs slightly better than HFSS, while the conventional strategy fails. Thus, the amplitude formulation works quite well at both the low and high frequencies in contrast to the conventional strategy or HFSS. The same problem has been solved using infinite elements in [26] for k0 a = 6 (see their Fig. 3), and the number of equations (i.e. number

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) k0 r = 10, φ = π/4 (real).

(b) k0 r = 10, φ = π/4 (imaginary).

(c) k0 r = 30, φ = π/4 (real).

(d) k0 r = 30, φ = π/4 (imaginary).

1383

Fig. 7. Variation of Ex with θ in the scattering from a conducting cube problem.

of unsuppressed degree of freedoms) that they report is 3,11,100. For this case, using our amplitude strategy, we get an almost perfect match with the analytical solution using the fine mesh (61138 equations), and a solution that is very close to the analytical solution using the coarse mesh (16698 equations). This shows the cost-effectiveness of the proposed method. 5.5. Scattering from a conducting ellipsoid This example is chosen to show that the amplitude formulation works well with a non-spherical geometry. Let the ellipsoid have dimensions 2a, 2a and 2c along the X , Y and Z axes, respectively. We present our results for k0 a = 12 and c = a/4. The domain is truncated at a spherical surface given by k0 R∞ = 24. The incident wave is of the form e−ik0 z ex . The converged HFSS result obtained using a very fine mesh (1,49,530 unknowns i.e. unsuppressed degree of freedoms) is used as the benchmark solution. Fig. 6 presents the real and imaginary parts of the dominant solution component, namely Ex , as a function of φ (near field) and as a function of θ (far field). These plots are along arcs of constant radii whose values are stated in the captions. Although results for the Ey and Ez components have not been presented, we have ensured that they are also accurate. The conventional and AMP element results have been obtained using a mesh of nr × nθ × nφ = 12 × 20 × 16 elements. The mesh uses a combination of 18 node wedge and 27 node brick elements. The total number of unknowns is 1,23,018. We once again observe that the amplitude formulation results are significantly better compared to the conventional ones.

1384

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) φ = 0.

(b) φ = π/2. Fig. 8. RCS as a function of θ in the scattering from a conducting cube problem.

5.6. Scattering from a conducting cube In this example, we demonstrate the performance of the amplitude formulation when the scatterer (in this case, a cube) has sharp corners and edges. Let the cube have dimension a. We present results for k0 a = 10 and k0 R∞ = 30. As discussed in Section 4.4, we surround the cube surface by a thin layer of thickness k0 t = 1 within which the regularization term (∇ · Aδ )(∇ · A) is suppressed, while in the remaining domain it is not. Let a plane wave Einc = eik0 z ex be incident on the cube. Once again, we use the HFSS results obtained using a very fine mesh (813 758 equations) as the benchmark solution. Fig. 7 presents the near- (k0 r = 10) and far-field (k0 r = 30) results for Ex , solved using tetrahedral elements with 112 442 equations. In the figures, ‘CON’ and ‘Current’ denote the conventional and amplitude formulations, respectively. Once again, we see that in both the near- and far-fields, the AMP element is considerably better than the conventional formulation, and has an almost perfect match with HFSS results; in the far-field the conventional formulation breaks down completely. The bistatic RCS plots at φ = 0, π /2 as a function of θ are shown in Fig. 8, where again we see the almost perfect match with the benchmark solution of the proposed method, and the large errors in the conventional formulation. 5.7. Scattering from a conducting cylinder This problem is different from the previous one in that both curved surfaces and sharp edges are present in the scatterer. We denote the radius and the height of the cylinder by a and h, respectively, and take k0 a = 3, k0 R∞ = 30 and k0 h = 6. As in the previous example, the regularization term is suppressed in a thin layer of thickness t, where k0 t = 0.5. A plane wave Ei = eik0 z ex is incident on the cylinder. The converged HFSS result, obtained using 267 572 unknowns, is taken as the benchmark solution. We have presented the results obtained using a tetrahedral mesh with 68 738 unknowns in Fig. 9, which again shows the superior performance of the amplitude formulation over the conventional one both in the near(k0 r = 6) and far-fields (k0 r = 30). 5.8. Scattering from a conducting L-shaped domain An L-shaped domain is constructed from a rectangular parallelepiped [−a, a] × [−b, b] × [−a, a] by removing the rectangular parallelepiped [−c , a] × [−b, b] × [−c , a]. A plane wave Ei = eik0 z ex is incident on the domain. The results are shown in Fig. 10 for k0 a = 5, k0 b = 2 and k0 c = 3. We have suppressed the regularization term in a layer of thickness t, where k0 t = 1. The benchmark solution is taken as very fine mesh HFSS result with 524 310 unknowns. The results for conventional and amplitude formulation (along with the benchmark solution), obtained using a tetrahedral mesh with 124 536 unknowns is shown in Fig. 10. The superiority of the amplitude formulation results is again evident from the plots. 5.9. Scattering from a dielectric sphere Consider a plane wave Einc = E0 e−ikz ex incident on a dielectric sphere. Let the relative permittivity and relative permeability of the dielectric sphere and the surrounding medium be denoted by (ϵ1 , µ1 ) and (ϵ2 , µ2 ), respectively. The

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) k0 r = 6, φ = π/4 (real).

(b) k0 r = 6, φ = π/4 (imaginary).

(c) k0 r = 30, φ = π/4 (real).

(d) k0 r = 30, φ = π/4 (imaginary).

1385

Fig. 9. Variation of Ex with θ for the scattering from a conducting cylinder problem.

analytical solution has been presented by Stratton [45]. The fields Escat and Hscat have the same expressions as in Eq. (33) with µ replaced by µ2 and k replaced by k2 , while Et and Ht are given by E t = E0

∞  n=1

(−i)n

2n + 1 

n(n + 1)

atn mon1 + ibtn nen1 ,



∞  k1 E0  2n + 1  t e Ht = − (−i)n bn mn1 − iatn non1 . ωµ0 µ1 n=1 n(n + 1)

The terms arn and brn are given by

µ1 jn (k1 a) [(n + 1)jn (ka) − kajn+1 (ka)] − µ2 jn (ka) [(n + 1)jn (k1 a) − k1 ajn+1 (k1 a)]   , µ1 jn (k1 a) (n + 1)h(n2) (ka) − kah(n2+)1 (ka) − µ2 hn(2) (ka) [(n + 1)jn (k1 a) − k1 ajn+1 (k1 a)]  2 µ1 jn (ka) [(n + 1)jn (k1 a) − k1 ajn+1 (k1 a)] − µ2 kk1 jn (k1 a) [(n + 1)jn (ka) − kajn+1 (ka)] brn = −  2  , µ1 h(n2) (ka) [(n + 1)jn (k1 a) − k1 ajn+1 (k1 a)] − µ2 kk1 jn (k1 a) (n + 1)h(n2) (ka) − kah(n2+)1 (ka) arn = −

(34)

1386

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) k0 r = 10, φ = π/4 (real).

(b) k0 r = 10, φ = π/4 (imaginary).

(c) k0 r = 30, φ = π/4 (real).

(d) k0 r = 30, φ = π/4 (imaginary).

Fig. 10. Variation of Ex with θ for the scattering from a conducting L-shaped domain problem.

(2)

atn = btn =

jn (ka) + arn hn (ka) jn (k1 a)

,

µ1 jn (ka) + µ1 brn h(n2) (ka)   . µ2 kk1 jn (k1 a)

The terms mon , men , non , nen and their limit values have the same expressions as in the conducting sphere case. The expressions (2)

for mon1 , men1 , non1 , nen1 and their limit values are obtained by replacing hn (kr ) and k by jn (k1 r ) and k1 , respectively, in  mon , men , non , nen . For the purpose of the numerical simulation, we assume the dielectric sphere to be in vacuum, i.e., ϵ2 = µ2 = 1. For the dielectric sphere, we have assumed ϵ1 = 1.5 and µ1 = 1. We take ka = 8 and kR∞ = 24. The incident wave is of the form e−ik0 z ex . For both formulations, we use a mesh of 29 413 ten-node tetrahedral elements with a total of 117 652 unknowns. Since the origin is part of the domain, we use conventional elements to mesh the dielectric sphere, and AMP elements to mesh the vacuum part of the domain. HFSS results were obtained using a total of 130 370 unknown degrees of freedom. Fig. 11 presents the results for the scattered electric field Ex inside the dielectric sphere at k0 r = 6.4, at a point just outside the dielectric sphere (k0 r = 9.6), and a point in the far-field (k0 r = 16) respectively. Other components are less dominant, and hence are not plotted. For all the three values of k0 r, the AMP elements performs better than conventional elements. Both AMP and edge elements closely match the analytical results.







A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) k0 r = 6.4, φ = π/6 (real).

(b) k0 r = 6.4, φ = π/6 (imaginary).

(c) k0 r = 9.6, φ = π/4 (real).

(d) k0 r = 9.6, φ = π/4 (imaginary).

(e) k0 r = 16, φ = π/6 (real).

(f) k0 r = 16, φ = π/6 (imaginary). Fig. 11. Variation of Ex with θ for the dielectric sphere problem.

1387

1388

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) k0 r = 4, φ = π/4 (real).

(b) k0 r = 4, φ = π/4 (imaginary).

(c) k0 r = 10, φ = π/4 (real).

(d) k0 r = 10, φ = π/4 (imaginary).

(e) k0 r = 30, φ = π/4 (real).

(f) k0 r = 30, φ = π/4 (imaginary). Fig. 12. Variation of Ex with θ for the dielectric cube problem.

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

(a) k0 r = 4.8, φ = π/4 (real).

(b) k0 r = 4.8, φ = π/4 (imaginary).

(c) k0 r = 7.2, φ = π/6 (real).

(d) k0 r = 7.2, φ = π/6 (imaginary).

(e) k0 r = 18, φ = π/3 (real).

(f) k0 r = 18, φ = π/3 (imaginary). Fig. 13. Variation of Ex with θ for the two-layer dielectric sphere problem.

1389

1390

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

5.10. Scattering from a dielectric cube A plane wave Einc = e−ik0 z ex with k0 a = 10 is incident on a dielectric cube with material properties ϵr = 1.5, µr = 1, and the surrounding medium is vacuum. A converged HFSS result with 310 220 unknowns is taken as the benchmark solution. We have used 10-node tetrahedral elements with 110 744 unknowns to obtain our results. Although the dielectric cube has sharp edges and corners, the domain is convex, and so we do not use the strategy outlined in Section 4.4. Similar to the previous problem, since the origin is part of the domain, we use conventional elements in a region bounded by a sphere of radius r with k0 r = 10, and AMP elements beyond that. Fig. 12 presents the scattered electric field Ex inside the cube (k0 r = 4), at a point just outside the cube (k0 r = 10), and in the far-field (k0 r = 30). The superior results of the proposed formulation are again evident. 5.11. Scattering from a two layer dielectric sphere We consider scattering from a two-layer dielectric sphere, with ϵr = 1.75 − 0.3i up to k0 r = 5.4, ϵr = 1.25 − 1.25i from k0 r = 5.4 to k0 r = 6, and vacuum beyond k0 r = 6 (a similar problem has been considered in [7]). For all the three domains, µr is taken as 1. We truncate the domain at k0 R∞ = 24. The incident wave is of the form eik0 z ex . A converged HFSS solution with 176 898 unknowns is used as the benchmark solution. For our formulation, we use 28 000 ten-node tetrahedral elements with a total of 112 000 unknowns to model the domain. The plots for Ex at k0 r = 4.8 (inside the dielectric sphere), k0 r = 7.2 (just outside the dielectric sphere), and k0 r = 18 (far-field) are presented in Fig. 13. Once again, we see that the performance of AMP elements is significantly better than conventional elements. 6. Conclusion and future work In this work, we have developed an Amplitude formulation for exterior electromagnetic radiation and scattering problems within the framework of the nodal finite element method. This strategy is based on the Wilcox asymptotic expansion for the far-field electromagnetic field. Based on this expansion, the highly oscillatory part of the solution is separated out a-priori, so that the finite element interpolation functions only have to capture a relatively, much-more gently varying function. Because of the outward wave-favoring nature of the assumed solution, each element absorbs radiation, as opposed to the conventional strategy where the elements allow both incoming and outgoing waves, and the Sommerfeld condition is imposed only through an absorbing boundary condition. For radiation problems, the amplitude formulation far outperforms the conventional strategy based on a first-order absorbing boundary condition, especially at high frequencies. For scattering problems, there is sharp variation of the electric field even along the angular directions, and for this reason, the performance of the amplitude formulation is not as dramatically superior as in radiation problems, although it is still significantly better. Since the amplitude formulation can be used in the immediate vicinity of the radiator/scatterer, the strategy is efficient. We have shown that the proposed strategy works for non-spherical geometries, and, with a slight modification of the formulation, yields good results even for scatterers with sharp edges and corners provided a thin layer of zero penalty elements is used in the vicinity of the scatterer. In problems where the origin is part of the domain, a small region of conventional elements in the vicinity of the origin with a surrounding layer of AMP elements needs to be used. The computational cost of the proposed strategy is only marginally higher compared to conventional elements due to some additional terms in the element stiffness matrix, with the global number of degrees of freedom the same. Although we have implemented this strategy within the framework of nodal finite elements, it can also be implemented within the framework of edge elements. Acknowledgment We gratefully acknowledge many helpful discussions with Mr. Neeraj Kumar of the ECE Department, IISc. References [1] A. Chatterjee, J.M. Jin, J.L. Volakis, Edge-based finite elements and vector ABC’s applied to 3-D scattering, IEEE Trans. Antennas and Propagation 41 (2) (1993) 221–226. [2] Bruno Stupfel, Numerical implementation of second- and third-order conformal absorbing boundary conditions for the vector-wave equation, IEEE Trans. Antennas and Propagation 45 (3) (1997) 487–492. [3] V.N. Kanellopoulos, J.P. Webb, A Numerical study of Vector absorbing boundary conditions for the Finite element solution of Maxwell’s equations, IEEE Microw. Guid. Wave Lett. 1 (11) (1991) 325–327. [4] J.P. Webb, V.N. Kanellopoulos, Absorbing boundary conditions for the finite element solutions of the vector wave equation, Microw. Opt. Technol. Lett. 2 (10) (1989) 370–372. [5] P. Joly, B. Mercier, A new second order absorbing boundary condition for Maxwell’s equations in dimension 3, INRIA Res. Report 1989. [6] F. Assous, E. Sonnendrucker, Joly–Mercier boundary condition for the finite element solution of 3D Maxwell equations, Math. Comput. Modelling 51 (2010) 935–943. [7] K. Sertel, J.L. Volakis, Multilevel fast multipole method solution of volume integral equations using parametric geometry modeling, IEEE Trans. Antennas and Propagation 52 (7) (2004) 1686–1692.

A. Nandy, C.S. Jog / Computers and Mathematics with Applications 71 (2016) 1364–1391

1391

[8] S. Alfonzetti, G. Borzi, N. Salerno, FEM analysis of unbounded electromagnetic scattering by the Robin iteration procedure, Electron. Lett. 32 (19) (1996) 1768–1769. [9] X.Q. Sheng, J.M. Jin, J. Song, C.C. Lu, W.C. Chew, On the formulation of Hybrid finite-element and boundary-integral methods for 3-D scattering, IEEE Trans. Antennas and Propagation 46 (3) (1998) 302–311. [10] P. Soudais, H. Steve, F. Dubois, Scattering from several test-objects computed by 3-D Hybrid IE/PDE methods, IEEE Trans. Antennas and Propagation 47 (4) (1999) 646–653. [11] J. Liu, J.M. Jin, A novel hybridization of higher order finite element and boundary integral methods for electromagnetic scattering and radiation problems, IEEE Trans. Antennas and Propagation 49 (12) (2001) 1794–1806. [12] X.Q. Sheng, E.K.N. Yung, Implementation and experiments of a Hybrid algorithm of the MLFMA-enhanced FE-BI method for open-region inhomogeneous electromagnetic problems, IEEE Trans. Antennas and Propagation 50 (2) (2002) 163–167. [13] M.N. Vouvakis, S.C. Lee, K. Zhao, J.F. Lee, A symmetric FEM-IE formulation with a single-level IE-QR algorithm for solving electromagnetic radiation and scattering problems, IEEE Trans. Antennas and Propagation 52 (11) (2004) 3060–3070. [14] K.C. Donepudi, J.M. Jin, W.C. Chew, A higher order multilevel fast multipole algorithm for scattering from mixed conducting/dielectric bodies, IEEE Trans. Antennas and Propagation 51 (10) (2003) 2814–2821. [15] L.E. Garcia-Castillo, I. Gomez-Revuelto, F.S. de Adana, M. Salazar-Palma, A finite element method for the analysis of radiation and scattering of electromagnetic waves on complex environments, Comput. Methods Appl. Mech. Engrg. 194 (2005) 637–655. [16] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 114 (1994) 185–200. [17] J.P. Berenger, Three-dimensional perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys. 127 (1996) 363–379. [18] Y. Xiao, Y. Lu, The prolate and oblate spheroid perfectly matched layers, IEEE Trans. Magn. 38 (2) (2002) 669–672. [19] Y. Xiao, Y. Lu, Combination of PML and ABC for scattering problem, IEEE Trans. Magn. 37 (5) (2001) 3510–3513. [20] U. Pekel, R. Mittra, An application of the perfectly matched layer (PML) concept to the finite element method frequency domain analysis of scattering problems, IEEE Microw. Guid. Wave Lett. 5 (8) (1995) 258–260. [21] J.H. Bramble, J.E. Pasciak, Analysis of a finite element PML approximation for the three dimensional time-harmonic Maxwell problem, Math. Comp. 77 (261) (2008) 1–10. [22] S.A. Cummer, A simple, nearly perfectly matched layer for general electromagnetic media, IEEE Microw. Wirel. Compon. Lett. 13 (3) (2003) 128–130. [23] R. Mittra, U. Pekel, A new look at the perfectly matched layer (PML) concept for the reflectionless absorption of electromagnetic waves, IEEE Microw. Guid. Wave Lett. 5 (3) (1995) 84–86. [24] L. Demkowicz, M. Pal, An infinite element for Maxwell’s equations, Comput. Methods Appl. Mech. Engrg. 164 (1998) 77–94. [25] W. Cecot, W. Rachowicz, L. Demkowicz, An hp-adaptive finite element method for electromagnetics. Part 3: A three-dimensional infinite element for Maxwell’s equations, Internat. J. Numer. Methods Engrg. 57 (2003) 899–921. [26] W. Rachowicz, A. Zdunek, An hp-adaptive finite element method for scattering problems in computational electromagnetics, Internat. J. Numer. Methods Engrg. 62 (2005) 1226–1249. [27] M.M. Ilic, S.V. Savic, B.M. Notaros, First order absorbing boundary condition in large-domain finite element analysis of electromagnetic scatterers, in: Telecommunication in Modern Satellite Cable and Broadcasting Services, 10 th International Conference, 2011; pp. 424–427. [28] T.V. Yioultsis, T.D. Tsiboukis, Development and implementation of second and third order vector finite elements in various 3-D electromagnetic field problems, IEEE Trans. Magn. 33 (2) (1997) 1812–1815. [29] W.E. Boyse, A.A. Seidl, A Hybrid finite element method for 3-D scattering using Nodal and Edge elements, IEEE Trans. Antennas and Propagation 42 (10) (1994) 1436–1442. [30] J. Tang, K.D. Paulsen, S.A. Haider, Perfectly matched layer mesh terminations for Nodal-based finite-element methods in electromagnetic scattering, IEEE Trans. Antennas and Propagation 46 (4) (1998) 507–516. [31] I. Bardi, O. Biro, K. Preis, Finite element scheme for 3D cavities without spurious modes, IEEE Trans. Magn. 27 (5) (1991) 4036–4039. [32] I. Bardi, O. Biro, R. Dyczij-Edlinger, K. Preis, K.R. Richter, On the treatment of sharp corners in the FEM analysis of high frequency problems, IEEE Trans. Magn. 30 (5) (1994) 3108–3111. [33] M. Costabel, M. Dauge, Weighted regularization of Maxwell equations in polyhedral domains, Numer. Math. 93 (2002) 239–277. [34] R. Otin, Regularized Maxwell equations and nodal finite elements for electromagnetic field computations, Electromagnetics 30 (2010) 190–204. [35] R. Otin, ERMES: A Nodal-based finite element code for electromagnetic simulations in frequency domain, Comput. Phys. Comm. 184 (2013) 2588–2595. [36] R. Otin, L.E. Garcia-Castillo, I. Martinez-Fernandez, D. Garcia-Donoro, Computational performance of a weighted regularized Maxwell equation finite element formulation, Prog. Electromagn. Res. 136 (2013) 61–77. [37] C.S. Jog, Arup Nandy, Mixed finite elements for electromagnetic analysis, Comput. Math. Appl. 68 (2014) 887–902. [38] C.H. Wilcox, An expansion theorem for electromagnetic fields, Comm. Pure Appl. Math. 9 (2) (1956) 115–134. [39] R. Mittra, O. Ramahi, A. Khebir, R. Gordon, A. Kouki, A review of Absorbing boundary conditions for two and three dimensional electromagnetic scattering problems, IEEE Trans. Magn. 25 (4) (1989) 3034–3039. [40] R.J. Astley, Wave envelope and infinite elements for acoustical radiation, Internat. J. Numer. Methods Fluids 3 (1983) 507–526. [41] R.J. Astley, W. Eversman, Finite element formulations for acoustical radiation, J. Sound Vib. 88 (1983) 47–64. [42] C.S. Jog, An outward-wave-favouring finite element-based strategy for exterior acoustical problems, Int. J. Acoust. Vibr. 18 (1) (2013) 27–38. [43] C. Geuzaine, J. Bedrossian, X. Antoine, An amplitude formulation to reduce the pollution error in the finite element solution of time-harmonic scattering problems, IEEE Trans. Magn. 44 (6) (2008) 782–785. [44] ANSYS Academic Research HFSS 2014.0. ANSYS, Inc. USA. [45] Julius Adams Stratton, Electromagnetic Theory, John Wiley & Sons, Chicago, 2007.