A Stable Higher Order Space-Time Galerkin Scheme for Time Domain Integral Equations

1.1 MB·

by  A. J. Pray

1
A Stable Higher Order Space-Time Galerkin
Scheme for Time Domain Integral Equations
A. J. Pray, Student Member, IEEE, Y. Beghein, Member, IEEE, N. V. Nair, Member, IEEE,
K. Cools, Member, IEEE, H. Bag˘cı Member, IEEE, and B. Shanker, Fellow, IEEE
4
Abstract—Stability of time domain integral equation (TDIE) functions can be effective in stabilizing TDIE solutions. The
1
0 solvers has remained an elusive goal for many years. Advance- smoothness of these functions allows for accurate evaluation
ment of this research has largely progressed on four fronts: (1)
2 of retarded potential integrals using conventional quadrature
Exact integration, (2) Lubich quadrature, (3) smooth temporal
rules. However, these functions are symmetric about their
n basis functions, and (4) Space-time separation of convolutions
a with the retarded potential. The latter method was explored interpolation point and, therefore, produce marching systems
J in [1]. This method’s efficacy in stabilizing solutions to the that are noncausal, requiring extrapolation. This necessitates

0 time domain electric field integral equation (TD-EFIE) was the use of very small time steps to maintain stability. Quasi-
1 demonstrated on first order surface descriptions (flat elements)
exact integration schemes transform surface integrals into
in tandem with 0th order functions as the temporal basis. In
volume integrals of the correlation function and the Green’s
] this work, we develop the methodology necessary to extend to
h higher order surface descriptions as well as to enable its use function,whicharethenevaluatedquasi-analytically.Thecrux
p with higher order temporal basis functions. These higher order of this procedure is in determining the correlation function
p- temporal basis functions are used in a Galerkin framework. A and its domain of support, which is difficult for higher order
number of results that demonstrate convergence, stability, and
m geometries. This motivates the twin goals of this work, which
applicability are presented.
are to develop a method that maintains temporal locality,
o
c Index Terms—Time-domain integral equations, higher order while achieving higher order temporal accuracy and being
. temporal basis, stability, time-domain analysis, marching-on-in- extensible to higher order geometric descriptions. As alluded
s time, space-time Galerkin method.
c to earlier, [1] was a purely numerical approach to stabilizing
i
s TDIEsolutions.Itis,therefore,easilyextendedtohigherorder
y I. Introduction descriptions of the geometry. However, achieving high order
h
accuracy in time is more challenging.
p Since the initial development of time domain integral equa-
Historically, collocation in time has been the preferred
[ tion (TDIE) solvers in the 1960s [2], their use in electromag-
method for time-marching schemes, although recent research
netic simulations has only recently been on the uptick. This
1
has shown that Galerkin methods yield more accurate and
v renewedinterestisthankstothesolutionofthecomputational
stable methods [15,16]. To achieve higher order accuracy in
5 complexity bottleneck [3,4] and increasingly sophisticated
3 approaches to addressing instability. While a plethora of time, in this work, we build on [16] and use a higher order
4 temporal basis, albeit with a small variation-the support of
stabilization schemes have been developed over the last two
2 these functions is restricted to one time interval. This permits
decades [5–10], the most promising schemes appear to be
1. the use of Galerkin testing to develop a causal marching
Lubich quadrature [11,12], smooth and bandlimited temporal
0 scheme. The stability of the resulting time marching scheme
basis functions [8], quasi-exact integration [13,14], and the
4 is further improved by using the purely numerical approach
1 separable expansion method [1]. Lubich quadrature relies on
developedin[1].Additionally,separationofspace-timeallows
: applying a Laplace transform to the appropriate TDIE and
v foreasyimplementationofhigherorderdiscretizationinspace
using a finite differencing scheme to map to the Z transform
i within this framework.
X domain. This discretized equation is then solved by marching
The paper is organized as follows: Section II formulates
r in time. These methods have yielded excellent results, but are
a the scattering problem, details its discretization, and presents
characterized by interactions with infinite temporal tails, lead-
modifications to impedance matrix elements using a separa-
ing to higher scaling in both computational cost and memory.
ble representation; and Section III shows the interpolation
Methods to overcome this cost scaling are current topics of
accuracy of the new temporal basis and presents a number
research.Alternatively,smoothandbandlimitedtemporalbasis
of scattering results, including convergence studies of farfield
scattering.
A.J.Pray,N.V.Nair,and B.ShankerarewiththeDepartmentofElectrical
and Computer Engineering, Michigan State University, East Lansing, MI,
USA,48824e-mail:[email protected]. II. Formulation
Y. Beghein is with the Department of Information Technology (INTEC),
GhentUniversity,Belgium
A. Time domain EFIE, MFIE, and CFIE
K.CoolsiswiththeElectricalSystemsandOpticsDivision,Universityof
Nottingham,NottinghamNG72RD,UK Consider a perfect electrically conducting (PEC) scatterer
H. Bag˘cı is with the Division of Computer, Electrical, and Mathematical residing in free space. Let ∂Ω represent the surface of the
SciencesandEngineering,KingAbdullahUniversityofScienceandTechnol-
ogy(KAUST),4700KAUST,Thuwal23955-6900,KingdomofSaudiArabia scatterer; the regions interior and exterior to the scatterer
2
are denoted as Ωi and Ωe = R3 \ Ωi, respectively; the functions Tl(t) are Lagrange polynomials, given as
j
unit vector normal to ∂Ω is denoted as nˆ. The scatterer is
illuminated by an incident plane wave {Einc(r,t),Hinc(r,t)}, Tlj(t)=Tl(t− j∆t)

efiwffehleidccthievxiescliyatesbssauncmuderlidrmentiottsebdJet(orv,astn)oimsohneinfg∂rleΩyq,usemwnhacliylchffomiranxt.itmuTrehnetg≤inecn0iedraeanntdet Tl(t)=(cid:96)0l(t) tot∈he[−rw∆its,e0]
(8)
scattered fields {Es(r,t),Hs(r,t)}, (r,t) ∈ Ωe × [0,∞) and (cid:89)p t−t
the total fields in Ωe are expressed as {Et(r,t),Ht(r,t)} = (cid:96)l(t)= t −ti
{Einc(r,t)+Es(r,t),Hinc(r,t)+Hs(r,t)}. The boundary condi- ii=(cid:44)0l l i
tions are given as t =(l/p−1)∆t, l=0,1,...,p .
l
(cid:12)
nˆ×nˆ×Et(r,t)(cid:12)(cid:12) =0 where ∆t is the time step. Here we note that, as opposed
(cid:12) r∈∂Ω (1)
nˆ×Ht(r,t)(cid:12)(cid:12) =J(r,t) . to the standard definition of basis functions in, e.g. [13],
r∈∂Ω these functions have support over one time step, ∆t. This is
Expressing {Es(r,t),Hs(r,t)} in terms of J(r,t) yields
essentialtomaintaincausalityintimemarchingwhenGalerkin
nˆ×nˆ×Einc(r,t)=L{J(r,t)} (2) testing is used in time. On the contrary, basis functions that
are typically used for MOT analyses [8,13] cannot be used
nˆ×Hinc(r,t)=J(r,t)+K{J(r,t)} (3) within a Galerkin framework, while maintaining causality, as
shown in the appendix. We note here that the temporal basis
where
(cid:40) (cid:90) functions in (8) are discontinous between time segments, i.e.
L{J(r,t)}=nˆ×nˆ× µ0∂t dr(cid:48)δ(t−R/c) ∗ J(r(cid:48),t) at t = j∆t. While this does not present any difficulties in
∇ (cid:90) δ(t4−πR/c∂)Ω (cid:90) t R t (cid:41) (4) dcaisncnroettizbiengutsheedTinD-dEiFscIrEetaiznidngTDth-eMtFimIEe, tdhiiffsetreesnttiinagte/bdasfiosrmsest
− dr(cid:48) ∗ dt(cid:48)∇(cid:48)·J(r(cid:48),t(cid:48))
4πε0 ∂Ω R t 0 of TD-EFIE and TD-MFIE because the correlation function
(cid:90) resulting from Galerkin testing is continuous, but not smooth.
1 δ(t−R/c)
K{J(r,t)}=−nˆ×∇× dr(cid:48) ∗ J(r(cid:48),t) (5) To discretize the time differentiated EFIE and/or MFIE, the
4π ∂Ω R t functions outlined in [16] should be used. This claim is
where ∗t denotes the convolution operation in t, R = |r−r(cid:48)|, justified in the appendix.
and the prime on ∇(cid:48) denotes that the differentiation is with
respect to r(cid:48). The integral in (5) is computed in the principle
value sense. Both (2) and (3) suffer from the well known Substituting (7) into (6), and employing Galerkin testing in
interior resonance problem for closed surfaces [17]. This can both space and time yields
be overcome by combining them as
−α/η nˆ×nˆ×Einc(r,t)+(1−α)nˆ×Hinc(r,t)=
0
−α/η0L{J(r,t)}+(1−α)(J(r,t)+K{J(r,t)})(cid:17)C{J(r,t)} (6) (cid:88)i Z I =V (9)
i−j j i
(cid:112)
where η = µ /ε is the intrinsic impedance of free space j=0
0 0 0
andαisarealnumberbetween0and1.Theresultspresented where
in Section III will use a discrete version of either (2) or (6). I =(cid:104)Ij,0 Ij,1... Ij,p(cid:105)TV =(cid:104)vi,0 vi,1... vi,p(cid:105)T (10a)
j 1 1 Ns i 1 1 Ns
vi,k =−α/η (cid:104)j (r)Tk(t),nˆ×nˆ×Einc(r,t)(cid:105)
m 0 m i (10b)
+(1−α)(cid:104)j (r)Tk(t),nˆ×Hinc(r,t)(cid:105)
m i
Zmn =(cid:104)j (r)Tk(t),C{j (r)Tl (t)}(cid:105) (10c)
i−j,kl m i n j−i
Bcu.rTrSoepnsatoculevs-eiTniegmiteahJeGs(rera(t,l2teo))r,f=k(si3np(cid:88))Na,Dstoiiajrsnlc((ar6ren))t,d(cid:88)iNzwatteet(cid:88)miobppneogIrinjna,llTbbljy(atsr)eisprfeusnecnttiionngst(ah7se) Zi−j = ZZZZiii111N−−−111......sjjj1,,,01p000 Zi1−1j,01 ...... ZZii11−−11...jj,,0ppp ...... ZZi1N−N...sjN,s0sp  (10d)
n=1 j=0 l=0 i−j,p0 i−j,pp
where Ij,l aretheunknowncoefficientstobedetermined.j (r) The dimensions of Zi−j are Ns(p+1)×Ns(p+1) as opposed
n n
are RWG vector basis functions [18]. The temporal basis to N ×N as in the case of a traditional temporal basis set.
s s
3
C. Evaluation of Inner Products
The integrand in (10c) is piecewise continuous and ap- Observation βc∆t
plication of standard quadrature rules to (4) and (5) does
not lead to accurate evaluation of the integrals. One needs R−ζ
ζ
to explicitly account for the discontinuities in the integrand.
r(cid:48)
This can be achieved via quasi-exact integration [13,14], r c∆t
i.e. by identifying regions where the integrand is smooth
Source
and analytically integrating over each subregion. However, as
discussedinSectionI,thisprocedureisrestrictedtofirstorder
spatial discretizations. This restriction leads us to adopt the
Fig. 1: Separable expansion parameters
approach of [1], which is summarized here for the sake of
completeness. To begin, we define the vector potential as
(cid:90) J(r(cid:48),τ)
A{J(r,t)}= dr(cid:48) (11)
∂Ω 4πR Next, we present the modifications to the matrix elements
where τ = t − R/c. Employing space-time Galerkin testing in (10c) using the approximation (13). Continuing from (12),
to the contribution of the vector potential due to a source the contribution from the time differentiated vector potential
jn(r)Tlj(t) yields is given by
(cid:104)jm(r)Tik(t),∂tA{jn(r)Tlj(t)}(cid:105) (cid:68)j (r)Tk(t),∂A(cid:110)j (r)Tl(t)(cid:111)(cid:69)≈ µ0 (cid:88)Nh a αq T˜q (15a)
(cid:90) (cid:90) (cid:90) i∆t ∂Tl(τ) m i t n j 4π q mn ij,kl
= drj (r)· dr(cid:48)j (r(cid:48)) dtTk(t) t j q=0
(cid:90)∂Ωm m (cid:90)∂Ωn n (i−1)∆t i 4πR (12) where
whe(cid:17)re ∂∂ΩΩmdirsjmth(er)s·upp∂Ωonrtdor(cid:48)fjnj(r(r(cid:48)))ψ.iNj,kelx(rt,,rw(cid:48))e make the approx- αqmn =(cid:90)∂Ωmdrjm(r)·(cid:90)∂Ωndr(cid:48)jn(r(cid:48))RPq(ξ) (15b)
n n
imation Next, we define the scalar potential as
(cid:90) i∆t (cid:18) ζ(cid:19) (cid:90) (cid:90) t ∇(cid:48)·J(r(cid:48),t(cid:48)−R/c)
ψ (r,r(cid:48))= dtTk(t)δ t− Φ{J(r,t)}=− dr(cid:48) dt(cid:48) . (16)
ij,kl (i−1)∆t i c ∂Ω −∞ 4πR
(cid:16) (cid:17)
δ t− R−ζ Using the expansion (13) the contribution of this quantity to
∗t 4πRc ∗t∂tTlj(t) (10c) is
= 4π1R(cid:88)q∞=0aqPq(ξ)T˜iqj,kl (13) (cid:68)jm(r)Tik(t),∇Φ(cid:110)jn(r)Tlj(t)(cid:111)(cid:69)≈ 4π1ε0 (cid:88)qN=h0aqφqmnTˆiqj,kl (17a)
≈ 1 (cid:88)Nh a P (ξ)T˜q where
4πR q=0 q q ij,kl φq =(cid:90) dr∇·j (r)(cid:90) dr(cid:48)∇(cid:48)·jn(r(cid:48))Pq(ξ) (17b)
where ξ =k (R−ζ)/c+k , a =k 2q+1, mn Ωm m Ωn R
1 2 q 1
(cid:90) i∆t (cid:18) ζ(cid:19)2 Tˆq =(cid:90) i∆t dtTk(t)δ(cid:18)t− ζ(cid:19)∗
T˜iqj,kl = (i−1)∆tdtTik(t)δ t− c ∗t ij,kl (i−1)∆t i c t (cid:90) t (17c)
Pq(k1t+k2)P0,β(t)∗t∂tTlj(t) . Pq(k1t+k2)P0,β(t)∗t dt(cid:48)Tlj(t(cid:48)) .
−∞
Here, Pq(k1t+k2) is a Legendre polynomial of order q with Lastly, the tested operator K{J(r,t)} becomes
support t∈[0,β∆t], k and k are real numbers providing the
1 2
mapping[0,β∆t]→[−1,1],ζ isthemaximummultipleofc∆t (cid:68) (cid:110) (cid:111)(cid:69) 1 (cid:88)Nh
j (r)Tk(t),K j (r)Tl(t) ≈ a κq T¯q (18a)
between r and the source triangle (see figure 1), and P0,β(t) m i n j 4π q=0 q mn ij,kl
is a window function defined as
 where
P0,β(t)=10 tot∈he[0rw,βi∆set]. (14) κmqn =(cid:90) drjm(r)·nˆ×(cid:90) dr(cid:48)jn(r(cid:48))
Ω Ω
The expansion in (13) enables accurate evaluation of the spa- m (cid:32)k1P(cid:48)q(ξ) − nPq(ξ)(cid:33)×Rˆ (18b)
tial integrals in (10c). It can be shown that the expansion (13) cR R2
isuniformlyconvergent.Foragivennumberofharmonics,Nh, (cid:90) i∆t (cid:18) ζ(cid:19)
anexpressionfortheerrorboundisgivenintheappendix(the T¯q = dtTk(t)δ t− ∗
same expansion can be used in computing the tested gradient ij,kl (i−1)∆t i c t (18c)
of the scalar potential). P (k t+k )P (t)∗ Tl(t) .
q 1 2 0,β t j
4
The prime on P(cid:48)q(·) denotes the 1st derivative of Pq(·) with p ksamp 5 10 20 40
respect to its argument. Returning to (17a), it can be shown 1 1/1/1 1/1/2 2/1/2 3/1/3
that, for (i− j)>ζˆ+β+1, 2 2/1/1 4/2/2 4/3/3 4/3/3
1 3 4/3/4 4/4/4 6/4/5 7/5/5
Tˆq = A δ
ij,kl (cid:32)a(cid:90)0 i∆klt q0 (cid:33)(cid:32)(cid:90) j∆t (cid:33) (19) TtiAonBaLcEcuIr:aLcioewse(runbdoiuffnedreonftiNathedn/e1esdteddertiovarteipvreo/idnutcegerianl)terpola-
A = dtTk(t) dtTl(t)
kl i j
(i−1)∆t (j−1)∆t
whereζˆ =ζ/(c∆t)andδmn istheKroneckerdelta.Substituting error is defined as Error = (cid:107)f − f (cid:107)/(cid:107)f (cid:107)
approximate analytic analytic
(19) into (10c) yields the modified marching system
where,
(cid:88)i Z I +(cid:88)i Zφ I =V (20a) fanalytic =[f(∆t−∆) f(2∆t−∆) ... f(Nt∆t−∆)]
j=0 i−j j j=0 i−j j i fapproximate =[f˜(∆t−∆) f˜(2∆t−∆) ... f˜(Nt∆t−∆)]
(23)
(cid:34)Ziφ−,mjn(cid:35) =04π1ε0φ0mnAkl ii−− jj>≤ββ++11++ζζˆˆ . (20b) f˜(t)=(cid:88)jN=t1(cid:88)l=p0 IljTlj(t) .
k,l
Therefore, The coefficients Il are computed by approximating an un-
j
shifted Gaussian and applying Galerkin testing. ∆/∆t is de-
(cid:88)i (cid:34)Zφ,mn(cid:35) Ij,l =i−(cid:88)β−2(cid:32)φ0mnAkl(cid:33)Ij,l finedassomeshiftbetween0and1,whichinthisexperiment
i−j n 4πε n is set to 0.906. The center frequency is set to f =200 MHz
j=0 k,l j=0 0 (20c) 0
and f = 300 MHz. As expected, the error can be seen in
φ0 A max
= mn klCi−β−2,l
4πε n
0
where
(cid:88)i
Ci,l = Ij,l . (20d)
n n
j=0
(20d) is used while marching in time to recursively compute
the charge terms Ci,l while retaining the O(NN 2) scaling of
n t s
the solver.
The integrals in (15b), (17b), and (18b) are of the form
(cid:90) j(r(cid:48))P (ξ)
dr(cid:48) q (21)
Ω Rα
n
Fig. 2: Scaling of interpolation error for various p and sam-
which can be evaluated to arbitrary precision using an ap-
pling frequencies
propriate order of integration rule. Obviously, if q is too high
thentheorderofintegrationwill,ofnecessity,beimpractically
high. The order of integration is determined by Nh, i.e. the figure 2 to scale as O(∆tp+1) for the different values of ksamp.
highest order harmonic required for an interaction pair. While In order to estimate the lower bound of N we repeat this
h
theoretical bounds on N can be obtained for a given error interpolation experiment for various basis function orders and
h
(see Section V-C), this is not a tight estimate. In practice N time step sizes. For each setup, we compute the number of
h
is significantly lower. This is elucidated further in the next harmonics needed to reproduce the interpolation accuracy of
Section. the same basis function order/time step size achieved without
the approximation (13). The value of β was chosen to mimic
III. Results a typical simulation on a tessellation consisting of elements
approximately λ /10 in length, i.e.
min
A. Interpolation (cid:38) (cid:39)
λ /10
Totesttheaccuracyofthesebasisfunctions,theyareusedto β= mcin∆t (24)
interpolate and shift an approximately bandlimited modulated (cid:108) (cid:109)
= k /5
gaussian pulse of the form samp
f(t)=cos(ω0t)e−(2t−σt2p)2 (22) wNherreequλimreind=toc/pfrmoadxu.cTeatbhlee eIrtraobrulleavteeslsthperomviidneidmuinmfivgaulruee2o.f
h
where σ = 3/(2πf ), t = 6σ, and f denotes the
max p max
frequency at which the power is approximately 160 dB below
the peak value at f = ω /2π. Figure 2 shows the order of B. Scattering Results
0 0
interpolation accuracy for various p. The time step is chosen The remainder of this Section presents scattering results
(cid:16) (cid:17)
as∆t=1/ 2k f ,wherek >0isarealnumber.The from a variety of scatterers. For each simulation, a PEC
samp max samp
5
scatterer is illuminated by a plane wave of the form is computed for θ = −180◦ and φ = 0◦. Figure 4 shows the
convergence in the farfield for the various simulation setups.
Einc(r,t)=uˆcos(2πf t)e−(t−r·kˆ/c−tp)2/2σ2 (25)
0 Similarly to the previous result, the rates of convergence for
where uˆ is the electric polarization and kˆ denotes the di-
rection of propagation. σ, t , and f are defined in Section
p 0
II-B. For all simulations, except where noted to the contrary,
∆t=1/(20f ) (k =10) and p=2. RCS comparisons are
max samp
made at f , f +∆f, and f −∆f, where f +∆f < f . The
0 0 0 0 max
error against a frequency domain solver (or analytical result
where available),
(cid:107)ζ −ζ (cid:107)
Error= tdie fdie (26)
(cid:107)ζ (cid:107)
fdie
is given at these three frequencies, while the plot is shown
for only f and f + ∆f in order to maintain a reasonable
0 0
range. The (cid:96)2 norm is used. Here RCS = 10 log (ζ )
tdie 10 tdie Fig. 4: RCS convergence for sphere
and RCS = 10 log (ζ ) denote the radar cross sections
fdie 10 fdie
obtained at a discrete set of angles using a TDIE and FDIE
variousvaluesofθinthe x−zplanewerefoundtoliebetween
solver, respectively. O(∆tp+1) and O(∆tp+2).
For our first result, we examine convergence in temporal
Our next scatter is a thin PEC box, discretized using 1,146
basisfunctionorderandsamplingfrequency.Aplateofdimen-
unknowns, with dimensions 100×50×10 m. The scatterer
sions 1m × 1m, discretized using 200 triangular elements, is
is excited by a plane wave with parameters f = 1.4 MHz,
illuminated by a plane wave of parameters f = 150 MHz, 0
0 f =2.7MHz,uˆ =yˆ,and xˆ=kˆ.Thisisadifficultscattering
f = 225 MHz, kˆ = zˆ, and uˆ = xˆ. The time domain max
max problem to stabilize due to the relatively small dimensions in
simulationwasperformedforthreevaluesofpandfourvalues
thez-direction,particularlysousingtheTD-EFIE.Thecurrent
of k . For this result, we do not compare against frequency
samp
domain results, but rather look at convergence in farfield
scattering. Here, convergence is defined for given values of
θ and φ as
Convergence = |Eθk+1(θ,φ)−Eθk(θ,φ)| (27)
k+1 |Eθk+1(θ,φ)|
where Es(θ,φ)≈φˆEφ (θ,φ)+θˆEθ (θ,φ) is the farfield approx-
k k
imation to the scattered field for a given sampling frequency,
indexed by k. The result is shown in figure 3 for an angle of
θ=−130◦ and φ=0◦. The rates of convergence were seen to
be nonuniform across observation angles θ in the x−z plane,
andwerefoundtolieroughlybetweenO(∆tp+1)andO(∆tp+2).
Fig. 5: Current coefficient on box
isobservedfor10,000timesteps(57transitsacrossgeometry)
and can be seen in figure 5 to remain stable throughout the
duration of the simulation. To validate the accuracy of the
Fig. 3: RCS convergence for plate
For our second result, we perform the same convergence
study, but for scattering from a sphere of radius 1 m, dis-
cretized with 576 unknowns. The sphere is illuminated by a Fig. 6: RCS of box
plane wave with f = 60 MHz, f = 90 MHz, uˆ = xˆ, and
0 max
kˆ = zˆ. The TD-CFIE is used with α = 0.5 and the farfield solution,theRCSiscomputedinthe x−yplaneandcompared
6
with a validated frequency domain EFIE solver. The radar
cross section of the box is obtained from the two solvers at
3 different frequencies. Figure 6 shows the RCS values at the
twohigherfrequenciesandcloseagreementisseen.Theerror
between the two solutions was found using (26) to be 0.40%,
0.36%, and 0.19% at the 0.2, 1.4, and 2.6 MHz, respectively.
Our next result is scattering from a cone-sphere of 7,965
unknowns discretized with 2nd order elements with its axis
along the z-direction. The spherical portion of the scatterer
has radius 1 m, while the height of the cone is 4 m. The
excitation has f = 80 MHz, f = 150 MHz, uˆ = yˆ, and
0 max
kˆ = x. The current for 4,000 time steps (80 transits across
geometry) is computed using the TD-CFIE with α = 0.5. Fig. 9: Current coefficient on almond shaped scatterer
Fig. 10: RCS of almond shaped scatterer
Fig. 7: Current coefficient on cone-sphere
is challenging given the sharp tip and the fact that it is very
thin in the z-direction. Validation is presented in figure 10 via
comparisonofRCSinthe x−zplanewithanFD-CFIEsolver.
The error is found to be 0.038%, 2.1%, and 0.59% at 11, 40,
and 69 MHz, respectively.
Our final result is solution to the EFIE for scattering from
a VFY218. This structure fits in a box of size 15×9×4 m
and is discretized using 6,498 flat elements. The excitation is
x-polarized with kˆ =−yˆ, f =60 MHz, and f =100 MHz.
0 max
In this result, the current is discretized using temporal basis
functions of order p = 0. This is an extremely challenging
Fig. 8: RCS of cone-sphere
This is a challenging geometry due to the sharp tip of the
scatterer, yet as is seen in figure 7, the current remains stable
throughout the simulation. The solution is again validated in
figure 8 via RCS comparison with an FD-CFIE solver. The
RCSiscomputedinthe x−yplaneand,again,closeagreement
betewenthetimeandfrequencydomainsolutionsisseen,with
theerrorsgivenas0.29%,1.7%,and0.49%at11,80,and149
MHz, respectively.
Our next result is a thin almond shaped scatterer of 1,668
unknowns,discretizedwith2ndorderelements.Thegeometry Fig. 11: Current coefficient on VFY218
fits within a box of 5×4.33×0.865 m. The incident wave has
parameters f = 20 MHz, f = 35 MHz, uˆ = xˆ, and kˆ = zˆ. problem due to the sharp corners and edges on the geometry,
0 max
The TD-CFIE with α = 0.5 is used. Figure 9 shows that yetthecurrentinfigure11isstablethroughoutthesimulation
the current remains bounded throughout the simulation (80 (7,000 time steps or 100 transits across the geometry). The
transits across geometry or 6,000 time steps). This geometry validity of this result is demonstrated in figure 12. The same
7
B. Smoothness of temporal basis functions
To justify this claim, we examine the time dependent terms
in the tested vector potential
(cid:104)U(t),∂2δ(t−R/c)∗T(t)(cid:105) (30)
t
where U(t) is a compactly supported, discontinuous, locally
integrable function and T(t) is compactly supported. (30) can
be written as
(cid:90) ∞ (cid:90) ∞
dtδ(t−R/c) dt(cid:48)U(t(cid:48))∂2T(t(cid:48)−t) (31)
t(cid:48)
−∞ −∞
Fig. 12: RCS of VFY218 Due to the compactness of U(t) and T(t)
(cid:90) ∞ (cid:90) ∞
dt(cid:48)U(t(cid:48))∂ 2T(t(cid:48)−t)=− dt(cid:48)∂ U(t(cid:48))∂T(t(cid:48)−t) (32)
t(cid:48) t(cid:48) t
problem was solved using a frequency domain EFIE solver −∞ −∞
and the two results are plotted against each other. Excellent which is finite if T(t) is at least weakly differentiable to first
agreement is seen. order, i.e. T(t) is continuous. For the undifferentiated TD-
EFIE, the quantity of interest is
IV. Conclusion
(cid:90) ∞ (cid:90) ∞
We have presented a higher order (in geometric discretiza- dt(cid:48)U(t(cid:48))∂t(cid:48)T(t(cid:48)−t)=− dt(cid:48)∂t(cid:48)U(t(cid:48))T(t(cid:48)−t) (33)
tion and temporal basis function order) space-time Galerkin −∞ −∞
schemeforTDIEsbasedonaseparableexpansioninspaceand which will be finite if T(t) is integrable.
timeoftheretardedpotentialGreen’sfunction.Theexpansion
yields smooth spatial integrands, which can be evaluated to C. Convergence of (13)
arbitrary precision through purely numerical means. This en- The error incurred through truncation (13) can be deter-
ablesthemethodtobeusedonhigherordertessellationswhile mined via passage to the Fourier domain in both t and r,
retaining high accuracy in the matrix elements. To extend the yielding quantities depending on variables ω and λ, respec-
accuracyofthesolver,ahigherordertemporalbasiswasused tively. Given maximum temporal and spatial frequencies of
which expands the current using multiple functions within a interest, ω and λ , the bound on this error can be shown
max max
single time step. This allows the solver to be used within a to be
svpaalicdea-tteimteheGamleertkhiondfrwameehwavoerkpwreitsheonutetdviocolantvinegrgceanucsealrietys.uTltos (ztzs)NhK1(12−ztzzsz ) −K2(2Nh+1)(cid:32)1−ztzzsz (cid:33)2 (34)
for the scattered farfield of a PEC sphere as compared to t s t s
the Mie series. To elucidate the stability properties of the where
solver we have presented stable results of scattering from a
eω e|λ |c
variety of geometries, the solution of each being difficult or z = max , z = max (35)
t k (N +3/2) s k (N +3/2)
impossibletostabilizeviaexistingmethods.RCScomparisons 1 h 1 h
have been presented in each example to verify the accuracy and K1 and K2 are some positive constants. This bound
of the solution. converges rapidly for values of ztzs <1.
V. Appendix VI. Acknowledgement
A. Causality of marching system The authors would like to acknowledge support under NSF
Let the surface current be approximated by (7) and let the grantCCF1018516andtheDoDSMARTProgramundergrant
support of Tl(t) lie on the interval t ∈ [a,b], for all j > 0, N00244-09-1-0081.
j
0 ≤ j ≤ p, where b > a. To examine the causality of the
References
marching system we look at the value of
[1] A.Pray,N.Nair,andB.Shanker,“Stabilitypropertiesofthetimedo-
(cid:104)Tik(t),Tlj(t)∗tδ(t−∆)(cid:105)=(cid:104)δ(t−∆),Tl(t)∗tTk(−t−(j−i)∆t)(cid:105) (28) mainelectricfieldintegralequationusingaseparableapproximationfor
theconvolutionwiththeretardedpotential,”AntennasandPropagation,
Causality is maintained if (28) is zero for j−i > 0, ∆ ≥ 0, IEEETransactionson,vol.60,pp.3772–3781,aug.2012.
and t≥0. It can be shown that [2] M. B. Friedman and R. Shaw, “Diffraction of pulses by cylindrical
obstacles of arbitrary cross section,” J. Appl. Mech., vol. 29, pp. 40–
Tl(t)∗Tk(−t−(j−i)∆t)(cid:44)0 46,1962.
t
(29) [3] A.Yilmaz,J.-M.Jin,andE.Michielssen,“Timedomainadaptiveinte-
for (j−i)∆t∈((a−b)∆t−t,(b−a)∆t−t)
gralmethodforsurfaceintegralequations,”AntennasandPropagation,
Considering the extreme case in which t = 0, the support of IEEETransactionson,vol.52,pp.2692–2708,oct.2004.
[4] A.Ergin,B.Shanker,andE.Michielssen,“Theplane-wavetime-domain
this function is seen to include values of (j−i) > 0 only for algorithmforthefastanalysisoftransientwavephenomena,”Antennas
casewhenb−a>1.ThisshowsthatthesupportofTl(t)must andPropagationMagazine,IEEE,vol.41,pp.39–52,Aug.1999.
[5] B.P.RynneandP.D.Smith,“Stabilityoftimemarchingalgorithmsfor
be only one time step. We note that this also holds true for
the electric field integral equation,” Journal of Electromagnetic Waves
the important cases of the integral and derivative of Tl(t). andApplications,vol.4,no.12,pp.1181–1205,1990.
8
[6] A. Sadigh and E. Arvas, “Treating the instabilities in marching-on-in-
time method from a different perspective,” Antennas and Propagation,
IEEETransactionson,vol.41,no.12,pp.1695–1702,1993.
[7] J.-L. Hu, C. Chan, and Y. Xu, “A new temporal basis function for
the time-domain integral equation method,” Microwave and Wireless
ComponentsLetters,IEEE,vol.11,pp.465–466,nov2001.
[8] D. Weile, G. Pisharody, N.-W. Chen, B. Shanker, and E. Michielssen,
“Anovelschemeforthesolutionofthetime-domainintegralequations
ofelectromagnetics,”AntennasandPropagation,IEEETransactionson,
vol.52,pp.283–295,Jan.2004.
[9] S.J.Dodson,S.P.Walker,andM.J.Bluck,“Implicitnessandstability
of time domain integral equation scattering analysis,” Applied Compu-
tationalElectromagneticsSocietyJournal,vol.13,no.3,pp.291–301,
1998.
[10] T.K.Sarkar,W.Lee,andS.M.Rao,“Analysisoftransientscattering
from composite arbitrarily shaped complex structures,” Antennas and
Propagation,IEEETransactionson,vol.48,pp.1625–1634,2000.
[11] C.Lubich,“Convolutionquadratureanddiscretizedoperationalcalculus.
i,”NumerischeMathematik,vol.52,no.2,pp.129–145,1988.
[12] X.Wang,D.S.Weile,andS.Member,“ImplicitRunge-KuttaMethods
for the Discretization of Time Domain Integral Equations,” Antennas
andPropagation,IEEETransactionson,no.c,pp.1–13,2011.
[13] B.Shanker,M.Lu,J.Yuan,andE.Michielssen,“Timedomainintegral
equationanalysisofscatteringfromcompositebodiesviaexactevalua-
tionofradiationfields,”AntennasandPropagation,IEEETransactions
on,vol.57,pp.1506–1520,May2009.
[14] Y. Shi, M.-Y. Xia, R.-S. Chen, E. Michielssen, and M. Lu, “Stable
electric field tdie solvers via quasi-exact evaluation of mot matrix
elements,” Antennas and Propagation, IEEE Transactions on, vol. 59,
pp.574–585,feb.2011.
[15] T.Ha-Duong,B.Ludwig,andI.Terrasse,“Agalerkinbemfortransient
acoustic scattering by an absorbing obstacle,” Int. J. Numer. Meth.
Engng,vol.57,pp.1845–1882,2003.
[16] Y.Beghein,K.Cools,H.Bag˘cı,andD.DeZutter,“Aspace-timemixed
galerkin marching-on-in-time scheme for the time domain combined
fieldintegralequation,”AntennasandPropagation,IEEETransactions
on,vol.PP,no.99,p.1,2012.
[17] B. Shanker, A. A. Ergin, K. Agyu¨n, and E. Michielssen, “Analysis
of transient electromagnetic scattering from closed surfaces using a
combined field integral equation,” Antennas and Propagation, IEEE
Transactionson,vol.48,pp.1064–1074,2000.
[18] S.M.Rao,D.R.Wilton,andA.W.Glisson,“Electromagneticscatter-
ing by surfaces of arbitrary shape,” Antennas and Propagation, IEEE
Transactionson,vol.30,pp.408–418,1982.