ANDRIULLI, F

Cools, Kristof

MICHIELSSEN, E

English

Ghent University Academic Bibliography

Time-Domain Caldero´n Identities and Preconditioning of theTime-Domain EFIEKristof Cools*1, Francesco P. Andriulli2, and Eric Michielssen21Department of Information Technology (INTEC), Ghent University2Department of Electrical Engineering and Computer Science (EECS),The University of MichiganIntroductionFrequency domain EFIEs discretized via the method of moments have long been usedto simulate time harmonic scattering from perfect electrically conducting (PEC) sur-faces. Unfortunately, at low frequencies, the resulting method of moments matricesbecome ill conditioned; this phenomenon is known as low frequency breakdown.Recently, an analytical preconditioner for the frequency domain EFIE [1, 2] thatalleviates this problem was developed by leveraging the Caldero´n identities [3]. Notsurprisingly, time domain EFIEs for analyzing transient scattering from PEC sur-faces that are discretized via marching on in time (MOT) recipes also suffer from lowfrequency breakdown; this phenomenon manifests itself in terms of ill-conditionedMOT equations for large time steps. Here, an analytical preconditioner for the time-domain EFIE is developed. First, a time-domain operator identity that enables theconstruction of an analytically preconditioned time domain EFIE is constructed.Second, a procedure for efficiently discretizing this new time domain EFIE is eluci-dated. Numerical results that demonstrate the performance of the new analyticalpreconditioner are presented.Equations and DiscretizationConsider a closed PEC surface S with external normal nˆ, which is illuminated by atransient electric field Ei(r, t). The “derivative form” of the time-domain EFIE forthe current on S, J(r, t), readsT˙J = M˙i. (1)Here M i(r, t) = −nˆ ×Ei(r, t) represents the excitation and the operator T mapsJ(r, t) to the scattered field nˆ×Es; T can be decomposed asTJ = ThJ + TsJ (2)whereTsJ = −1cnˆ×∫Sds′δ˙(t−R/c)4piR∗ J(r′, t) (3a)ThJ = cnˆ×∫Sds′∇H(t−R/c)4piR∗ ∇′ · J(r′, t). (3b)Here H is the Heaviside function, δ is the Dirac impulse, and “∗” denotes tem-poral convolution. The capacitive part Th of this operator is hypersingular with1-4244-0123-2/06/$20.00 ©2006 IEEE 2975a nontrivial nullspace. Discretizing equation (1) using the standard MOT recipeyieldsj−1∑k=0Zk · Ij−k + Vj = 0, j = 1...NT (4)where (Zk)m,n =〈b′ml′j |T˙ |bnlj−k〉and (Vj)m = −〈b′ml′j |M i〉; b′m and bn are spatialbasis and testing functions and l′j and lj−k are temporal ones. Successively solvingthese equations for j = 1, ..., NT is equivalent to solving the lower triangular blockmatrix equation Z0Z1 Z0Z2 Z1 Z0.......... . .I0I1I2...+V0V1V2... = 0. (5)When the speed of light times the timestep is much larger the spatial element size,the block matrix on the left side of (5) becomes ill-conditioned.In the frequency domain, the Fourier transforms T , Th, and Ts of the operatorsT , Th, and Ts are used in constructing the frequency domain EFIE. The followingCaldero´n identity holds [3]:T 2 = K2 − 14I. (6)Here, I is the identity operator on the space of time-harmonic current distributionsand K is the compact operator defined in [3]. It follows from (6) that T 2 is asecond kind operator. Letting this operator act on an arbitrary surface distributionyields an equality between functions. By computing the inverse Fourier transformof both sides, a time-domain counterpart of the Caldero´n identity is obtained. Toavoid Heaviside functions as convolution kernels, the second time derivative of thetime-domain identity is usedT˙ 2J = K˙2J − 14Iδ¨(t) ∗ J . (7)Leveraging techniques leading to the construction of analytically preconditionedfrequency domain EFIEs, this equation inspires the definition of the following newtime-domain EFIE:T˙ 2J = T˙M˙i. (8)This equation is of “the second kind”. Therefore, it is expected to be better condi-tioned than the original time-domain EFIE. Upon discretizing, (8) yields a systemsimilar to (5). The MOT recipe can be used to expand J(r, t) and to test theequation. Because discretizing the complete operator T˙ 2 is too cumbersome, thetwo factors will be discretized separately. This amounts to expanding the identityoperator “in between” in a suitable set of basis functions. Moreover, this has tobe done in such a manner that the good properties of the operator T˙ 2 are con-served. By inspecting the definitions (3), it is recognized that the spatial behaviorof the operators Ts and Th resembles that of their frequency domain counterparts.In particular, the actions of these operators on the curl- and divergence-free sub-spaces of surface distributions is the same as in the frequency domain. From this,2976it follows that in the discretization of the spatial part of the “identity in between”,techniques similar to those adopted in the frequency domain [1], can be used in thetime-domain. This means that the spatial part of the domain and range space ofT˙ 2 is well approximated by Rao-Wilton-Glisson basis functions fm. The interme-diate space is best described by rotated RWG functions nˆ × fm. This leaves uswith chosing a suitable set of temporal basis functions. From [4], it is clear thatLagrange interpolants Li(t) = L(t − i∆t) are good candidates to describe a causalbandlimited function of time. On the other hand, if we want to arrive at the strictlylower triangular form of (5), we need to test the two operators T˙ with a Dirac im-puls δj(t) = δ(t − j∆t). As a consequence, the temporal part of the intermediatespace is described with two different sets of basis functions. The range of the rightT˙ in (8) is described with Dirac impulse functions and the domain of the left T˙ isdecomposed using Lagrange interpolants Li. Since for small timesteps, the quotientLi/∆t limits to a Dirac impulse, this dual approach causes no problems and the twosets of basis functions can describe the same intermediate space. Finally, a genericelement of the discretized T˙ 2 can be cast as< fmδj |T˙ 2|fnLi >=< fmδj |T˙ |nˆ× fpLr > (G−1S )p,q(G−1T )r,s < nˆ× f qδs|T˙ |fnLi >(9)where summation over recurring indices is understood and (GS)m,n =< fm|fn >and (GT )i,j =< Li|δj > are spatial and temporal Gramm matrices. As can bereadily verified, the temporal Gramm matrix equals the identity. As can be seenfrom (9), the new system block matrix is the product of two different discretizationsof the T˙ operator W0W1 W0....... . . = Z0lZ1l Z0l....... . . ·G−1S · Z0rZ1r Z0r....... . . . (10)Here, the Zkl and Zkr blocks stem from the discretization of the left and right T˙operator and Wk = Z0l ·G−1S ·Zkr + ...+Zkl ·G−1S ·Z0r . The new system matrix is theproduct of two complete MOT system matrices. However, if our primary goal is tocure the condition number of the system (5), we only need to alter the block Z0.Indeed, in the MOT scheme, this is the only block that has to be inverted. Inspiredby this observation, the system matrix W0W1 W0....... . . = Z0lZ0l. . . ·G−1S · Z0rZ1r Z0r....... . . (11)is used. Now Wk = Z0l · G−1S · Zkr . The advantage of this approach (which willhenceforth be referred to as time-localized preconditioning) is that it retains thebenevolent effect on the condition number of the system block matrix in generaland Z0 in particular, without introducing the computational burden of computingcomplete expressions forWk. The eigenvalues of a triangular block matrix are com-pletely described by the eigenvalues of the blocks that are located on the diagonal.The new system (11) and (10) equally affect the condition number of the systemblock matrix. If an iterative solver is used, the decrease in the number of iterations2977Figure 1: Condition number and backscattered far field.will be the same. However, if we want to address the stability issue for the MOT, itis expected that in order to obtain benefits, the entire block matrix (10) has to beused, as all blocks Zk are involved.Numerical ResultsA simulation using the system (5) is compared to one using (11). The geometryis a sphere of radius 0.25m. This sphere is discretized using 170 triangular surfacepatches with 255 unknowns. The incident wave is a base-band GaussianEi(r, t) =4T√pixˆe−γ2(12)with γ = 4T (ct− ct0 − zˆ · r), T = 200 meter, and t0 = 300 lightmeter. In figure 1(a)the condition number is plotted versus the sampling frequency 1/dt. The conditionnumber of the preconditioned system is quasi independent of ∆t. In figure 1(b),the backscattered far field obtained from a EFIE simulation and a modified EFIEsimulation are plotted versus time. The timestep used is 5 lightmeter.References[1] Adams R. and Champagne N. A numerical implementation of a modified formof the electric field integral equation. IEEE Transactions on Antennas andPropagation, 52(9):2262–2266, september 2004.[2] Contopanagos H., Dembart B., Epton M., Ottusch J., Rokhlin V., VisherJ., and Wandzura S. Well-conditioned boundary integral equations for three-dimensional electromagnetic scattering. IEEE Transactions on Antennas andPropagation, 50(12):1824–1830, december 2002.[3] Hsaio G. and Kleinman R. Mathematical foundations for error estimation innumerical solutions of integral equations in electromagnetics. IEEE Transactionson Antennas and Propagation, 45(3):316–328, march 1997.[4] Aygu¨n K., Schanker B., Ergin A., and Michielssen E. A two-level plane wavetime-domain algorithm for fast analysis of emc/emi problems. IEEE Transac-tions on Electromagnetic Compatability, 44(1):152–164, february 2002.2978

Time-domain Caldéron identities and preconditioning of the time-domain EFIE

https://biblio.ugent.be/publication/356492/file/572030

Time-domain Caldéron identities and preconditioning of the time-domain EFIE

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography