Convolution quadrature methods for 3D EM wave scattering analysis

Engineering Analysis with Boundary Elements 43 (2014) 50–55

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Convolution quadrature methods for 3D EM wave scattering analysis Amir Geranmayeh Continental Automotive GmbH (I ID RD EE EL ED), VDO-Str. 1, 64832 Babenhausen, Germany

art ic l e i nf o

a b s t r a c t

Article history: Received 8 September 2013 Received in revised form 15 March 2014 Accepted 16 March 2014 Available online 12 April 2014

The time-domain boundary integral equations describing the electromagnetic wave scattering from arbitrary three-dimensional metallic structures are solved by applying spectral domain finite-difference approximations while mapping from the Laplace domain to the z-transform domain. The validity of the results are verified through comparison with high-resolution finite difference time domain method results and the convergence rate of the introduced time-marching schemes is compared with the time basis functions expansion methods. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Boundary element methods Convolution quadrature methods Finite difference delay modeling Time-domain analysis

1. Introduction The usage of boundary element methods to solve the timedomain integral equations (TDIE) [1] for the transient analysis of the wave radiation and scattering problems [2] has been of continuous interest in computational electromagnetic and numerical modeling communities for almost half a century [3]. The temporal discretization of the (TDIE) is commonly accomplished by either the implicit marching-on-in-time (MOT) schemes using subdomain Lagrange polynomial interpolation [4] or the alwaysstable marching-on-in-order/degrees (MOD) of Laguerre entiredomain bases [5]. An alternative approach for discretizing the time convolution integrals in the TDIE, competitive to the time basis functions expansions in the MOT or MOD recipes, is the Lubich's convolution quadrature methods (CQM) [6], using the (first or) second order backward finite difference (BFD) approximations in the Laplace domain. As a great advantage, in modeling dispersive dielectric material where the relative permittivity ϵr ðsÞ and permeability μr ðsÞ and the Green's functions are functions of s (e.g., in Debye equation for the complex permittivity), the usage of CQM for time-discretization in the Laplace domain let the frequencydependent characteristics be directly incorporated into the timedomain solver [7]. The underlying physics describing the wave scattering process is time invariant [8], as the material properties do not change over time. The CQM, hence, can be utilized to transform continuoustime representation of the time-invariant integral kernel (system transfer function) to discrete-time domain while approximating the TDIE derivatives in the spectral domain. The CQM have been

E-mail address: [email protected] http://dx.doi.org/10.1016/j.enganabound.2014.03.008 0955-7997/& 2014 Elsevier Ltd. All rights reserved.

successfully applied to transient heat conduction [9], acoustic wave propagation [10], elastic and viscoelastic dynamic problems [11] as well as poroelastic formulations [12]. The CQM are called finite difference delay modeling (FDDM) when the scattering analysis of arbitrarily shaped three-dimensional (3D) structures is carried out in a marching style [13]. The FDDM is a provably stable method when the Lubich's CQM for the time discretization is used in conjunction with the Galerkin moment method for the spatial discretization [14]. In the mathematical literature, the method is called convolution quadrature when it applies to the single layer potential for the Helmholtz operator [15]. It is also called Tustin's method in digital signal processing and control theory to transform continuous-time representation (transfer function) of a linear time-invariant system to discrete-time domain. In the CQM or FDDM methods, a conformal mapping from the Laplace domain to the z-transform domain (bilinear transform) based on a finite difference formula accomplishes the discretization in the z-domain, and the result is inverse transformed to create a time-domain method. The bilinear transformation preserves the stability by exact mapping every point of the jω-axis in the s-plane onto the unit circle jzj ¼ 1 in the z-plane. Following this approach, each frequency response of the continuous-time system can be processed using a discrete-time filtering technique. The numerical solution of retarded functional equations can benefit from this frequency wrapping once the system's unit delays are replaced by first order all-pass filters z  1 in the discrete domain, as schematically depicted in Fig. 1. The system eigenvalue analysis of the FDDM on small-scale case studies reported in [16] revealed that the CQM render a symplectic energy-conserving time integrator for the numerical solution of the TDIE. The temporal discretization is carried out by

A. Geranmayeh / Engineering Analysis with Boundary Elements 43 (2014) 50–55

Fig. 1. Block diagram of the time discretization/integration procedures by the CQM or FDDM methods.

either first or second order BFD approximation when the transformed TDIE is mapped from the Laplace domain to the z-transform domain. Implicit Runge–Kutta schemes can also be applied for the temporal discretization when mapping from the Laplace domain to z-domain. Thereafter, s parameter is replaced with a matrix function of z and the inverse z-transform is computed numerically using the discrete Fourier transform (DFT) [17]. The BFD approximations of greater than second order are, however, never absolutely stable. The absolutely stable (A-stable) Radau IIA methods with two- and three-stage has third- and fifth-order convergence, respectively. Wang et al. [13] determine the time derivative of the combined field integral equation (CFIE) in conjunction with the higher-order spatial bases. The time derivative of the integral equations, however, seeks for derivative of the transient excitation that may not be available for impulsive pulses. The CQM can also be applied to solve the original form of the electric field integral equation (EFIE) containing a time integral [16]. In this paper, the FDDM is also adopted to the primary EFIE without additional time derivative as well as separately to the magnetic field integral equation (MFIE) with a single time derivative, using the linearly varying divergenceconforming space basis functions. Additionally, the calculations of recursive convolutions are accelerated here using the time-FFT algorithms on Toeplitz block aggregates of the retarded interaction matrices [14].

2. TDIE and CQM methods Let S denotes the surface of a perfect electric conducting (PEC) body that is excited by a transient electromagnetic field Ei ðr; tÞ. The total tangential electric field on S remains zero for all times. As a result, the induced surface current vector Jðr; tÞ satisfies the following time-domain EFIE: Z Z Z τ μ ∂ Jðr0 ; τÞ 0 ∇r ∇r 0 :Jðr0 ; t 0 Þ 0 0 dS  dt dS R 4π ∂t S R 4πϵ S  1 ¼ E i ðr; tÞ

ð1Þ

where E ðr; tÞ ¼ n^  ðn^  E ðr; tÞÞ, R ¼ jr  r j and the observation point r and the source point r0 indicate arbitrarily located points on the surface S and n^ denotes an outward-directed unit vector normal to S at field point r. The variable τ ¼ t R=c is the retarded pffiffiffiffiffiffi time, in which the speed of wave propagation c ¼ 1= μϵ is determined through the permeability μ and permittivity ϵ of the surrounding environment. The derivative counterpart of the EFIE (DEFIE) is also of interest to preferably avoid the laborious computation of charge accumulation in solving the original EFIE (1) by the MOT methods. The DEFIE is obtained by taking a time derivative from both sides of (1) Z Z μ ∂2 Jðr0 ; τÞ 0 ∇r ∇r 0 :Jðr0 ; tÞ 0 dS  dS 2 R 4π ∂t S R 4πϵ S i

¼

∂E i ðr; tÞ : ∂t

i

0

ð2Þ

Note that the Hertz vector potential can also be adapted to solve the EFIE [5] instead of direct solving for the unknown surface

51

current. Evidently, the Hertz approaches demand extra postprocessing stages for computation of some desired electromagnetic quantities. Excluding the temporal derivation on the excitation term, the Hertz approach results in formulations identical to the DEFIE (2). Considering S as a closed surface, one may also consider the time-domain MFIE  Z  Jðr; tÞ 1 1 ∂Jðr0 ; τÞ R R  n^   2 þ Jðr0 ; τ Þ  3 dS0 ¼ n^  Hi ðr; tÞ 2 4π S0 c ∂t R R ð3Þ where S0 denotes the surface S from which the contribution of the singularity at R¼0 has been removed [14]. To numerically solve the time-domain EFIE (1), DEFIE (2), and MFIE (3) or any linear combinations of them, they are first transferred to the Laplace domain  Z Z  sR=c  μs ~ 0 e  sR=c 0 ∇r ~ 0 ; sÞe Jðr ; sÞ dS  dS0  ∇r 0 :Jðr 4π S R 4πϵs S R tan i ¼ E~ ðr; sÞ

μs2 4π

Z S

~ 0 ; sÞe Jðr

 sR=c

R

ð4Þ dS0 

∇r 4πϵ

Z

 sR=c ~ 0 ; sÞe ∇r 0 :Jðr R S

  dS0 

tan

¼ sE~ ðr; sÞ i

ð5Þ

 Z  ~ sÞ 1 s~ 0 R ~ 0 R Jðr; Jðr ; sÞ  2 þ Jðr ; sÞ  3 e  sR=c dS0  n^  4π 2 R R S0 c ~ i ðr; sÞ: ¼ n^  H

ð6Þ

~ sÞ is approxiThe unknown induced surface current density Jðr; mately expanded using Bernoulli's separation of variables in the space and Laplace domains by M

~ sÞ ¼ ∑ I~k ðsÞf k ðrÞ: Jðr; k¼1

ð7Þ

where I~k ðsÞ are unknown weighting coefficients of the spatial vector basis functions f k ðrÞ. Substituting (7) in the EFIE (4), the DEFIE (5), and the MFIE (6) respectively gives Z μs M ~ f k ðr0 Þ  sR=c 0 e ∑ I k ðsÞ dS 4π k ¼ 1 R S Z ∇r M ~ ∇r 0  f k ðr0 Þ  sR=c 0 ~ i e ∑ I k ðsÞ  dS ¼ E ðr; sÞ ð8Þ R 4πϵs k ¼ 1 S Z

μs2 M ~ f k ðr0 Þ  sR=c 0 e ∑ I k ðsÞ dS R 4π k ¼ 1 S 

∇r M ~ ∑ I ðsÞ 4πϵ k ¼ 1 k

Z

i ∇r 0  f k ðr0 Þ  sR=c 0 e dS ¼ sE~ ðr; sÞ R S

" Z M s 1 M ~ 1 f k ðr0 Þ  R  sR=c 0 0 ∑ I k ðsÞf k ðr Þ  n^  ∑ I~k ðsÞ e dS 2k¼1 4π c R2 S k¼1 # Z M f k ðr0 Þ  R  sR=c 0 ~ i ðr; sÞ: e dS ¼ n^  H þ ∑ I~k ðsÞ R3 S k¼1

ð9Þ

ð10Þ

Performing the Galerkin's testing procedure in space using the same set of vector basis functions f m ðrÞ; m ¼ 1; 2; …; M, respectively, then gives Z Z μs M ~ f k ðr0 Þ  sR=c 0 e ∑ I k ðsÞ f m ðrÞ  dS dS 4π k ¼ 1 R S S Z Z M 1 ∇r 0  f k ðr0 Þ  sR=c 0 e ∑ I~ ðsÞ ∇r  f m ðrÞ þ dS dS 4πϵs k ¼ 1 k R S S Z i ¼ f m ðrÞ  E~ ðr; sÞ dS ð11Þ S

52

A. Geranmayeh / Engineering Analysis with Boundary Elements 43 (2014) 50–55

Z

μs2 M ~ ∑ I ðsÞ f m ðrÞ  4π k ¼ 1 k S

Z

f k ðr0 Þ  sR=c 0 e dS dS R S Z Z 1 M ~ ∇r 0  f k ðr0 Þ  sR=c 0 e þ ∑ I k ðsÞ ∇r  f m ðrÞ dS dS 4πϵ k ¼ 1 R S S Z i ¼ s f m ðrÞ  E~ ðr; sÞ dS S

Z 1 M ~ ∑ I k ðsÞ f m ðrÞ  f k ðr0 Þ dS 2k¼1 S " Z Z 1 M s~ f k ðr0 Þ  R  sR=c 0 ∑ I k ðsÞ f m ðrÞ  n^  e dS dS  4π k ¼ 1 c R2 S S # Z Z M f k ðr0 Þ  R  sR=c 0 ~ e dS dS þ ∑ I k ðsÞ f m ðrÞ  n^  R3 S S k¼1 Z ~ i ðr; sÞ dS: ¼ f m ðrÞ  n^  H S

ð12Þ

S

Z

μ M i ∑ ∑ I ðt Þ f m ðrÞ  4π k ¼ 1 j ¼ 0 k j S

m ¼ 1; 2; …; M:

ð14Þ

In the following, the temporal discretization mapping from the Laplace domain to the z-transform domain is explained. In the absence of inhomogeneous initial conditions, multiplication by s in the Laplace domain corresponds to the temporal differentiation. Thus, one may replace s with a finite difference approximation. Considering the unit delay property of the z-transform, the firstorder backward Euler and the second-order backward difference approximations to s are obtained respectively by s¼

1

1z Δt

ð15Þ

3  4z  1 þ z  2 : s¼ Δt

ð16Þ

~ The impedance matrices ZðsÞ are function of sp e  sR=c , p ¼  1; ~ 0; 1; 2. Substituting (15) and (16) to the functions ZðsÞ, 1

sp e  sR=c js ¼ 1  z  1 ¼ ∑ ωpn z  n Δt

ð17Þ

n¼0

Δt

ω0n ¼

n¼0

ξ e

ð19Þ

n!  n=2 qffiffiffiffiffiffi 1 ξ e  3=2ξ H n ð 2ξÞ n! 2

ð20Þ

ωpn ¼

1

Δt

1 ðωpn  1  ωpn   1Þ

p40

1 1 p1 ð3ωpn  1  4ωpn   1 þ ωn  2 Þ 2 Δt

ð21Þ p 40:

ð22Þ

Therefore, using the second-order backward difference formula to approximate the Laplace variable (differentiation operator) s, the inverse z-transform of the coefficient matrix appears as a function of ordered Hermite polynomials. For p ¼  1 [18], n ξι ξn  ι ¼ Δt e  ξ ∑ : ðn  ι Þ! ι¼0 ι ¼ 0 ι! n

ωn 1 ¼ Δt e  ξ ∑

ð23Þ

ð25Þ

ð26Þ

where t i ¼ iΔt, or in general when (14) is transferred back to the time-domain, the sequential calling of the coefficient matrices forms a block-triangular Toeplitz system matrix 1

i

j¼0

j¼0

∑ Z i  jIj  ∑ Z i  jIj ¼ Vi

0 r i o1:

ð27Þ

Having assumed that all the coefficients up to i 1 are known, they are sent to the right-hand sides of (27) and the retained coefficients associated with the present weight j¼ i are found for i ¼ 0; 1; 2; …; N sequentially,

j¼0

0 ri o N o 1:

ð28Þ

Using a proper cutoff strategy [19] before preforming the matrixvector products, the high storage cost is reduced for marching on smooth tails of the late transient response, i.e., (28) becomes Z 0Ii ¼ Vi 

n ξ

and for the higher order derivation p 4 0 terms in (17) and (18) respectively [13,15]

ωpn ¼

S

ð18Þ

and assuming ωpn ¼ 0 for n o 0, the Laurent series expansion gives the time sequence ωpn versus ξ ¼ R=cΔt respectively for (15) and (16) substitutions

ω0n ¼

Z 1 M ∑ I k ðt i Þ f m ðrÞ  f k ðr0 Þ dS 2k¼1 S " Z Z 1 M i f ðr0 Þ  R 0 ∑ ∑ I k ðt j Þ f m ðrÞ  n^  ω1i  j k dS dS  4π k ¼ 1 j ¼ 0 cR2 S S # Z Z 0 M i 0 0 f k ðr Þ  R ^ þ ∑ ∑ I k ðt j Þ f m ðrÞ  n  ωi  j dS dS R3 S S k¼1j¼0 Z ¼ f m ðrÞ  n^  Hi ðr; t i Þ dS:

i1

sp e  sR=c js ¼ 3  4z  1 þ z  2 ¼ ∑ ωpn z  n ;

ω2i  j

S

Z 0Ii ¼ Vi  ∑ Z i  jIj

1

Z

f k ðr0 Þ 0 dS dS R S Z Z 1 M i ∇r 0  f k ðr0 Þ 0 dS dS ∑ ∑ I k ðt j Þ ∇r  f m ðrÞ ω0i  j þ 4πϵ k ¼ 1 j ¼ 0 R S S Z ¼ s f m ðrÞ  E i ðr; t i Þ dS

ð13Þ

Each of the above three equations forms such a matrix equation in Laplace domain ~ ~ ~ ½ZðsÞ mk ½IðsÞk ¼ ½VðsÞm

Now, taking the inverse z-transform of (11), (12) and (13), the multiplications in the frequency-domain become convolutions in the time-domain [19], Z Z μ M i f ðr0 Þ 0 dS dS ∑ ∑ I k ðt j Þ f m ðrÞ  ω1i  j k R 4π k ¼ 1 j ¼ 0 S S Z Z 1 M i ∇r 0  f k ðr0 Þ 0 dS dS ∑ ∑ I k ðt j Þ ∇r  f m ðrÞ ωi1j þ 4πϵ k ¼ 1 j ¼ 0 R S S Z ð24Þ ¼ f m ðrÞ  E i ðr; t i Þ dS

i

∑

j ¼ maxð0;i  N g Þ

Z i  jIj;

0 r io N

ð29Þ

where the almost-zero late-appearing matrices Z n C 0, n 4 Ng can be eliminated from the construction of the matrix equation, after n marching on the solutions up to ξmax 5 n! [16]. Thus, the knowledge of the largest electrical dimension of the structure is sufficient to determine after how long marching on time samples, the newly appearing interaction matrices can be neglected without disturbing the late-time accuracy and stability. Additionally, by taking the advantages of the Toeplitz properties on the contribution of the retarded matrices Z i  j and using the inherently parallelizable time-FFT for the discrete convolutions, the computational complexity is reduced further. Note that the singularity may exist for several terms of ωpn 1=R, n ¼ 0; 1; 2; … in (19)–(22) when R-0. For instance, for the treatment of singularity associated with ω2i  j 1=R in (22) five terms j ¼ i; i 1; …; i  4 have to be handled separately (Figs. 2 and 3). Alternatively, taking unilateral inverse z-transform from the one-sided Laplace transform of the ideal sampled quantities in

A. Geranmayeh / Engineering Analysis with Boundary Elements 43 (2014) 50–55

53

3. Numerical results

Fig. 2. Coefficients ω0k for the first-order backward-difference approximation as a function of the electric distance.

Fig. 3. Coefficients ω0k for the second-order backward-difference approximation as a function the electric distance.

(8)–(10), one obtains convolutions with discrete time-domain samples Z μ M i f ðr0 Þ 0 dS ∑ ∑ I k ðt j Þ ω1i  j k R 4π k ¼ 1 j ¼ 0 S þ

∇r M i ∑ ∑ I ðt j Þ 4πϵ k ¼ 1 j ¼ 0 k

Z S

∇r 0  f k ðr0 Þ 0 dS ¼ E i ðr; t i Þ R

ð30Þ

∇r 0  f k ðr0 Þ 0 ∂ i dS ¼ E ðr; t i Þ R ∂t

ð31Þ

ωi1j

The validity of the schemes (24), (25), and (26) was checked on several 3D PEC objects under the plane-wave incidence. For the sake of brevity, three case studies are presented here, namely a sculpture, propeller and cube-spheres shapes. The exact geometrical dimension of the selected objects can be found in the opensource tool NETGEN, which is utilized here for 3D automatic triangular mesh generation, optimization, and hierarchical refinement. The presented discretization approaches’ solution converges when sufficiently large number of spatio-temporal unknowns are solved. This has been first checked by simultaneous space-mesh and time-resolution refinement on an electrically small cube as the benchmark problem without dictating the thresholding process that gradually takes energy out of the system [14]. To better distinguish which numerical scheme converges faster to the reference solution, the results are shown here when the CQM are about to converge after marching on a limited number of spatiotemporal degrees of freedom with the specified thresholding criteria. ^ The representative structure is illuminated by the x-polarized Gaussian plane-wave propagating along the z^ direction with the full-width half max of s, where cs ¼ 0:7071 light-meter. The governing integral equations are discretized in the spectral domain for the unknown induced surface current using the CQM procedure described in Section 2. The spatial variation of the current is represented by divergence-conforming vector bases f k ðrÞ with constant tangential and linear normal distribution, also known as triangular roof-top or Rao–Wilton–Glisson (RWG) functions, defined on each common edge joining two flat triangular panels [20]. The temporal sampling rates are determined according to the smallest electric distance between the mesh barycenters, Δt ¼ αRmin =c. The surface integrals are evaluated by adaptive refining quadruple quadratures [14]. For the numerical calculation of the vector and scalar potentials, respectively, four and sevenpoint quadrature methods are used on (partitioned) spatial subdomains. The impedance matrices Z n are not further calculated after n 4 N g . Fig. 4 considers solving the scattering problem of a PEC sculpture with an outer radius of 0.8 m meshed by 1720 triangular patches with M¼ 2580 common edges at Nt ¼260 time samples using the DEFIE (2) and discretized according to the CQM (25) with ^ α ¼ 5:65. The x-directed induced current density on the equator of the PEC sculpture, shown in Fig. 4, reveals that the higher-order BFD in CQM provides better agreement with the finite difference

Z

μ M i f ðr0 Þ 0 dS ∑ ∑ I k ðt n Þ ω2n k R 4π k ¼ 1 j ¼ 0 S þ

∇r M i ∑ ∑ I ðt n Þ 4πϵ k ¼ 1 j ¼ 0 k

Z

S

ω0n

1 M ∑ I ðt Þf ðrÞ 2k¼1 k i k " Z M i 1 f ðr0 Þ  R 0 ^ dS  n  ∑ ∑ I k ðt j Þ ω1i  j k 4π cR2 S k¼1j¼0 # Z M i f ðr0 Þ  R 0 dS ¼ n^  Hi ðr; t i Þ: þ ∑ ∑ I k ðt j Þ ω0i  j k 3 R S k¼1j¼0

-k

-z

ð32Þ

Now, the Galerkin's testing procedure in space on (30)–(32) gives the final matrix equations (24)–(26), respectively.

Fig. 4. The x-component of the transient induced surface current density on the equator of the PEC sculpture model for different CQM associated matrix equations.

54

A. Geranmayeh / Engineering Analysis with Boundary Elements 43 (2014) 50–55

H

H k y z

E

z

E

-z x

x

Fig. 5. The magnitude of the transient induced surface current density on the equator of the PEC sculpture model for different CQM system matrix equations in logarithmic scale.

k

y

-k

H E

Fig. 7. The x-component of the transient induced current density on the middle of the top surface of the two cylinders (propeller) with 0.5-m diameter and 2-m length obtained by the solution of different CQM systems.

of the response of the original form of EFIE in comparison with the results of DEFIE where different backward integrators are applied. ^ Fig. 7 illustrates the x-component of the induced current density on a 2-m long propeller surface meshed by 1704 triangle patches. The MFIE (3) is solved via (32) for M¼2556 unknown edges. The CQM parameters are set α ¼ 2:5, Ng ¼46, Nt ¼285 for the first-order BFD and α ¼ 3, Ng ¼60, Nt ¼238 for the secondorder BFD. Fig. 7 expresses agreement between the higher-order FDDM and the MOD results for the PEC propeller illuminated by s ¼ 0:7071-width incident pulse, even when only the early Ng computed interaction matrices are employed, i.e. Z i C 0 is assumed for i 4 N g .

y x

Fig. 6. The magnitude of the transient induced current density on the cube-sphere surface at point r ¼ ð0; 0; 0:02Þ in logarithmic scale. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.)

time-domain (FDTD) reference solution. It is worth mentioning that the solution procedure associated with first-order BFD terminates impedance matrices calculations after Ng ¼37 and the second-order BFD after Ng ¼50 samples. Fig. 5 shows the absolute values of the responses plotted in logarithmic scale. The resultant flat level for the late-time current in Fig. 5 demonstrates that using the first-order spatial basis functions in CQM yields stable results. Fig. 6 exhibits a PEC cube-sphere object with 1-m side lengths meshed with 592 triangles resulting M ¼888 edges. The surface current induced on the cube-sphere, obtained through solving the original form of the EFIE (1) by (30) using α ¼ 8:89, is plotted in logarithmic scale for Nt ¼178 time samples. Fig. 6 also demonstrates the late-time stability of solving the DEFIE (2) using (31) at two different time-frequency sampling steps, namely α ¼ 10 in red and α ¼ 14:125 in blue colors. Fig. 6 indicates that the properly chosen N g Z 64 for eliminating the almost-zero late-appearing matrices from the construction of the right side of the matrix equation do not harm at all the stability and accuracy of the transient analysis. Upon simultaneous start of excitation signal, as the time integral is approximated by summations on time sequences, up to a sampling rate time-shift is observed in the peak

4. Conclusions Alternative forms of CQM for the solution of boundary integral equations were investigated based-on spectral domain finitedifference approximation. They all apply a mapping from the Laplace domain to the z-transform domain. The CQM provide stable frameworks for transient wave scattering analysis. The computational cost for long time simulations can be considerably reduced by avoiding the calculation of impedance matrices associated with late-time interactions. The validity of the results was verified through comparison with FDTD and MOD results. References [1] Miller EK. An overview of time-domain integral equations. J Electromagn Waves Appl 1987;1(April (3)):269–93. [2] Martin R, Salinas A, Bretones A. Time-domain integral equation methods for transient analysis. IEEE Antennas Propag Mag 1992;34(June (3)):15–23. [3] Sayre EP, Harrington RF. Time domain radiation and scattering by thin wires. Appl Sci Res 1972;26(September (1)):413–44. [4] Pantoja M, Bretones A, Martin R. Time-domain analysis of thin-wire loaded antennas using integral equations. IEE Microw Antennas Propag 2000;147 (June (3)):203–6. [5] Ji Z, Sarkar TK, Jung BH, Yuan M, Salazar-Palma M. Solving time domain electric field integral equation without the time variable. IEEE Trans Antennas Propag 2006;54(January (1)):258–62. [6] Lubich C, Schneider R. Time discretization of parabolic boundary integral equations. Num Math 1992;63(December (1)):455–81. [7] Wang X, Weile DS. Electromagnetic scattering from dispersive dielectric scatterers using the finite difference delay modeling method. IEEE Trans Antennas Propag 2010;58(May (5)):1720–30. [8] Tijhuis A, Peng Z. Marching-on-in-frequency method for solving integral equations in transient electromagnetic scattering. IEE Proc H (Microw Antennas Propag) 1991;138(4):347–55.

A. Geranmayeh / Engineering Analysis with Boundary Elements 43 (2014) 50–55

[9] Abreu A, Canelas A, Mansur W. A CQM-based BEM for transient heat conduction problems in homogeneous materials and FGMs. Appl Math Mod 2013;37(February (3)):776–92. [10] Abreu A, Carrer J, Mansur W. Scalar wave propagation in 2D: a BEM formulation based on the operational quadrature method. Eng Anal Bound Elements 2003;27(February (2)):101–5. [11] Schanz M, Antes H. A new visco- and elastodynamic time domain boundary element formulation. Comput Mech 1997;20(October (5)):452–9. [12] Schanz M. Wave propagation in viscoelastic and poroelastic continua: a boundary element approach. In: Lecture notes in applied computer mechanics, vol. 2. Berlin: Springer Verlag; 2001. [13] Wang X, Wildman RA, Weile DS, Monk P. A finite difference delay modeling approach to the discretization of the time domain integral of electromagnetics. IEEE Trans Antennas Propag 2008;56(August (8)):2442–52. [14] Geranmayeh A. Time domain boundary integral equations analysis. Südwestdeutsche Verlag für Hochschulschriften: SVH-Verlag; 2011. [15] Hackbusch W, Kress W, Sauter SA. Sparse convolution quadrature for time domain boundary integral formulations of the wave equation. IMA J Numer Anal 2009;29(March (1)):158–79.

55

[16] Geranmayeh A, Ackermann W, Weiland T. Finite difference delay modeling of potential time integrals. In: Proceedings of IEEE antennas and propagation society international symposium no. 11. Toronto, ON; 2010. p. 1–4. [17] Wang X, Weile DS. Implicit Runge–Kutta methods for the discretization of time domain integral equations. IEEE Trans Antennas Propag 2011;59(December (12)):4651–63. [18] Kielhorn L, Schanz M. Convolution quadrature method-based symmetric Galerkin boundary element method for 3-d elastodynamics. Int J Numer Meth Eng 2008;76(June (11)):1724–46. [19] Banjai L, Sauter SA. Rapid solution of the wave equation in unbounded domains. SIAM J Numer Anal 2008;47(October (1)):227–49. [20] Wilton DR, Rao SM, Glisson AW, Schaubert DH, Al-Bundak OM, Butler CM. Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains. IEEE Trans Antennas Propag 1984(March (3)):276–81.