Time-Domain Methods for the Maxwell Equations

160 Pages·2001·1.73 MB·English

by  Andersson U.

Chapter 1
troInduction
1.1 Computational electromagnetics
1.1.1 The ellMaxw equations
TheellMaxw equationserew (cid:12)rstulatedformyb James Clerk ell.Maxw They are:
r(cid:1) D = (cid:26) (Gauss 0 w)la
r(cid:1) B = 0 (Gauss 0 w)la
(1.1)
@ B = (cid:0)r(cid:2) E yarada(F 0sw)la
@t
@ D = r(cid:2) H (cid:0) J (Ampere 0sw)la
@t e
where E is the electric (cid:12)eld, D is the electric (cid:13)ux,ydensit H is the magnetic (cid:12)eld,
B isthemagnetic(cid:13)ux,ydensit J istheelectrictcurrenydensitand (cid:26) isthehargec
e
ydensit [ Che89 ].
The ellMaxw equations describe electromagnetic phenomena. This includes
micro-, radioandradar es.vaw TheellMaxwequationsarediscussedinmoredetail
in Chapter 3 . y’saradaF and (cid:18)Amp ere’s wsla constitute a (cid:12)rst-order ypherbolic
system of equations. The owt Gauss’ wsla can be edderiv from y’saradaF and
(cid:18)Amp ere’s wslavidedpro that the initial conditions ful(cid:12)ll the Gauss’ ws.la
These equations are linear and it yma hence appear to be rather easy to esolv
them.analytically er,evwHoboundary and terfacein conditions emak the ellMaxw
equations hard toesolv.analytically They can beedsolvanalytically only for a few
eryv simple shapes hsuc as a sphere or an in(cid:12)nite circular cylinder. Hence one has
torelyonamixofexptserimenandeximativapproand/orumericalnmethods. Nu-
mericalmethodsfortheellMaxwequationsareusuallyreferredtoascomputational
electromagnetics (CEM).
1
2 Chapter 1. ductionoIntr
Exptserimenand CEMtcomplemenheacother whenelopingdeva product. An
tageanadvof umericaln methods is that they emak it possible to test a large um-n
ber of tdi(cid:11)eren constructions without actually building them. CEM is also useful
when exptserimenare di(cid:14)cult and/or dangerous to perform. hSucan example is a
tninglighestrik on an aircraft in t.(cid:13)igh
Thefastariationsvintheelectromagnetic(cid:12)eldsemakitahallengectoconstruct
umericaln methods for the ellMaxw equations.
1.1.2 Applications
There is a wide range of applications for CEM. Some of the more imptortan ones
are:
(cid:15) electromagnetic ycompatibilit (EMC),
(cid:15) tennaan analysis and thesis,syn
(cid:15) radar cross section CS)(R calculations,
(cid:15) cellular umanphone{h body teraction,in
(cid:15) micr evawovo ens,
(cid:15) target recognition and
(cid:15) ybrid/monolithich micr evawoin tegrated circuits.
yMan of these applications are impossible to model in eryev detail. orFinstance,
the teriorin of a modern aircraft is (cid:12)lled with umerousn wires and other small
details that are impossible to eresolv in a computation. us,Th when modeling a
radarpulsestrikinganaircraftitisimpossibletoumericallyncomputetheinduced
tcurren in eryev cable. er,evwHo it tmigh still be possible to accurately predict
the radar cross section of the aircraft. Whether this is possible depends strongly
on the \electrical size" of the aircraft. By \electrical size" ew mean the relation
eenbwetsome appropriate length scale, for instance the length of the aircraft, and
theelengthvawof the radar e.vaw
A major task in industrial CEM is the creation of the computational grids.
The objects of a calculation are usually described with CAD models. Quite often,
the CAD (cid:12)les are enbrok (in the sense that the CAD surfaces are not connected
together properly), hwhic means that the geometry ustm be repaired. This thesis
will not address this imptortan aspect of industrial CEM: instead, ew will assume
that the computational grids exist.
eWwill also be rather brief on the imptortanaspect of postprocessing the com-
putational results. The eransw of a computation is seldom a simple \YES" or
\NO". Inymancases, large e(cid:11)orts need to be sptenon analyzing the results. This
isoftendoneybvisualizingthedata. Thevisualizationof3Delectromagnetic(cid:12)elds
is a non-trivial business, and is further discussed in Sections 4.8 and 6.6.3 .
1.1. Computationalomagneticsctrele 3
1.1.3 Numerical methods for the ellMaxw equations
The ellMaxw equations can be edsolv either in the time domain or the frequency
domain. urthermore,Fthe umericaln method can be applied either on the partial
tialdi(cid:11)eren equation (PDE) ulationform of the ellMaxw equations in ( 1.1 ) or on a
boundary tegralin ulation.form ableT 1.1 ysdispla examples of methods with this
classi(cid:12)cation.
Time Domain requencyFdomain
PDE ulationform FD-TD FEM
tegralInulationform MOT MoM
ableT1.1. Classi(cid:12)cationofumericalnmethodsfortheellMaxwequations.
The abbreviations inableT 1.1 :
(cid:15) FD-TD = Finite-Di(cid:11)erence Time-Domain
(cid:15) MOT = hing-On-in-TimeMarc
(cid:15) FEM = Finite tElemen Method
(cid:15) MoM = Method of tsMomen
ableT 1.1 lists only the most commonly used method in heac.category There are
of course umerousn other methods.
Time-domain methods canesolva problem foreralsevfrequencies in one single
calculationandtheycanalsowfollothepulseolutionevintime. Becausethisthesis
concernstime-domainmethodsewwilltrateconcenourdiscussiononthesemethods
and only brie(cid:13)y tcommen on frequency-domain methods. The latter methods are
of course best suited for applications where only a few frequencies are t.presen
requency-domainFmethods, tegralin ulationform
requency-domainFulationtegral-forminmethods,hsucasMoM[ an91W ],reducethe
olumetricv equations to surface equations and usth reduce the bumern of spatial
dimensions of the problem yb one. Anothertageanadvis that after solving a par-
ticular problem for one angle of incidence, it is elyrelativ easy to (cid:12)nd the response
for another angle of incidence. er,evwHo it is bcumersome to handle cases with
aryingvmaterial properties.
MoM results in a dense linear system of equations. Solving this system directly
with Gaussian elimination has the ycomplexit O ( N 3 ) if the size of the matrix is
N (cid:2) N . Assuming that eweepk thebumern of tselemen per elengthvawt,constan
N increases proportionally to f 2 , where f is the.frequency The orkw to esolv the
MoM system directly is therefore O ( f 6 ). Oneyawto diminish this orkloadw is to
esolv the system with eiterativ methods. eIterativ methods are usually based on
ectormatrix-vultiplication,mhwhichas aycomplexitof O ( N 2 ). Theorkwtoesolv
4 Chapter 1. ductionoIntr
the MoM system elyiterativ is O ( f 4 ) if the eiterativ method ergesvcon.nicely An
enev better ycomplexit can be edhievac yb Multipole methods. In this case the
linear system of equations can beedsolvwith O ( Nlog ( N )) arithmetic operations if
aelultilevmmethod is used [ W93CR , SC95 ].
Anotheryaw to reduce the ycomplexit of MoM is to use the so-called ysicalph
optics (PO) method. Here the wnsunkno on the surface are computed directly
from the tinciden (cid:12)eld. teractionIn eenbwettdi(cid:11)eren parts of the surface is hence
neglected. This is a high frequency ximation,appro PO and MoM egiv ticaliden
results as the frequency tends to.yin(cid:12)nit Using POewcan compute thewnsunkno
on the surface in O ( N ) arithmetic operations
requency-domainFmethods, PDE ulationform
One PDE ulationform of the ellMaxw equations in the frequency domain is the
ectorv Helmholtz equation. orFthe electric (cid:12)eld it is
1
r(cid:2) ( r(cid:2) E )= ! 2 (cid:15) E ; (1.2)
(cid:22)
where ! istheangularfrequencyand (cid:15) and (cid:22) arespace-deptendenmaterialproper-
ties. The most commonyawtoesolvthis is to use (cid:12)nitetselemen(FEM) [ CK98V ]
because of their geometric .y(cid:13)exibilit er,evwHo (cid:12)nite di(cid:11)erences are also used
[Lar00 ]. The widespread commercial code HFSS [ HFS ] uses FEM.
Time-domain methods, tegralin ulationform
Time-domain methods for the tegralinulationform of the ellMaxw equationsevha
not been widely used. er,evwHo in the last few earsy there has been an increase
in e(cid:11)orts on this subject. Most methods are so-calledhing-on-in-timemarc(MOT)
methods. The ycomplexit of original MOT methods is O ( N N 2 ), where N is the
t s t
bumern of timesteps and N is the bumern of surface hes.patc This ycomplexit
s
can be edvimpro yb using so-called plane- evaw time-domain (PWTD) [ ESM98 ]
methods.
PWTD methodsevhabeen created yb adapting ideas from ultipmole methods
described e.abvo The ycomplexit of the elo-levwt PWTD is O ( N N 4 = 3 log ( N )),
t s s
and the ycomplexit of theelultilevmPWTD is O ( N N log ( N )).
t s s
Some tagesanadv of tegralin equation methods as compared to PDE methods
in the time domain are:
(cid:15) They do not su(cid:11)er from dispersion errors.
(cid:15) They only discretize a surface.
(cid:15) No absorbing boundary condition is needed.
1.1. Computationalomagneticsctrele 5
Akwbacdraof MOT methods is that they are prone toyinstabilit[ RS90 ]. The
issue has been studied in detailyber’salkWgroup. They state that MOThemessc
for solving magnetic (cid:12)eld tegralin equations \can be stabilized for all practical
purposes" using implicit timestepping methods [ WB97D ].
Time-domain methods, PDE ulationform
In the time domain there are eralsev possible hniquestec for termediatein frequen-
cies,including(cid:12)nitedi(cid:11)erences(FD-TD)[ af00T ],(cid:12)niteolumesv(FV-TD)[ SHM89 ],
and(cid:12)nitetselemen(FE-TD)[ SF90 ]. ThetagesanadvandtagesandisadvofFD-TD,
FE-TD and FV-TD are thoroughly discussed in this thesis, in particular in Chap-
ter 7 . eWwill describe them brie(cid:13)y here.
In CEM, the ymacron FD-TD refers to a (cid:12)nite di(cid:11)erence ximationappro of
y’saradaF and (cid:18)Amp ere’s wsla using second-order accurate tralcen di(cid:11)erences in
timeandspaceonagridthatisstaggeredinspaceandtime. Thegridisillustrated
in Figure 1.1 yb wingsho one cell of the grid. This method asw troinduced in
1966 ybeeY[ ee66Y ] and asw further elopdeved ybeva(cid:13)oTin the 1970s. It is the
most commonly used time-domain method. It is conceptually easy to grasp and
eryv te(cid:14)cien for homogeneous domains. The major kwbacdra is its yinabilit to
handle edcurv boundaries.accurately The FD-TD method is described in [ af00T ]
H x H y
z E
z
H
z
E
y y
z
y
E
x
x
x
Figure1.1. ositionsPoftheelectricandmagnetic(cid:12)eldectorvcomptsoneninaunit
eeYcell.
and discussed in Chapter 4 .
ItispossibletoconstructFD-TDhemessconunstructuredgrids(seeforinstance
Chapter 4 of [ af98T ]). In this case it is eryv kytric to ehievac a stable method.
There are owt other hesapproac ailableva on unstructured grids, namely FV-TD
and FE-TD.
6 Chapter 1. ductionoIntr
Finite olumesV erew troinduced to CEM yb arShank [ SHM89 ] yb exporting
methodsfromcomputational(cid:13)uiddynamics(CFD).Hisearlyorkwusedstructured
grids, but lately he has turned to unstructured grids. His main reason for doing so
is the ydi(cid:14)cult of creating a global body-conforming grid for realistic geometries,
hsuc as a complete aircraft. This orkw is also described in Chapter 4 of [ af98T ].
RileytroinducedanotherypteofFiniteolumeVmethod[ T97R ]. Hishemescasw
basedonstaggeringtheelectricandmagnetics(cid:12)elds. Thisorkwhasbeenuedtincon
ybEdelvik[ Ede00 ],whoseorkwwithexplicit(cid:12)niteolumeverssolvisatalfundamen
part of the time-domainybridhcodes in the GEMS project [ GEM ] (see Section 1.2
for a description of the GEMS project). The (cid:12)nite olumev grid is illustrated in
Figure 1.2 .
E
4
n p
H i
r t d
j,k H
q
E
3 t p E
i,m 1
n d
j
E
2
Figure 1.2. Acellintheprimarygridandadualface.
Anothermethodthatisellwadaptedforunstructuredgridsisthe(cid:12)nitetelemen
time-domain method (FE-TD) [ LLC97 ]. The (cid:12)nite telemen method is based on a
ariationalvulationformofthePDEinsomesuitableHilbertspace. ximationsAppro
to the solution are then tsough in a (cid:12)nite-dimensional subspace.
AcommonhapproacfortheellMaxwequationsistodiscretizespacewithtetra-
hedra(trianglesin2D)anduseso-called\edgets"elemen[ Ned80 ]asbasisfunctions
for the (cid:12)nite-dimensional subspace. The ectorv basis function for an edge e in 2D
is plotted in Figure 1.3 . en(Evthough only one triangle iswnshoin Figure 1.3 , the
basis function actually has support on both the triangles that has e as an edge.)
j
edge e
Figure 1.3. Theectorvbasisfunction ’ e foredge e .
1.1. Computationalomagneticsctrele 7
This ectorv basis function is designed to (cid:12)t eryv ellw with the ysicsph of the
ellMaxwequations. It enforcestialtangenyuittinconbutwsallonormalu-tindiscon
yit of the (cid:12)elds. urthermore,Fit ful(cid:12)lls r(cid:1) ’ =0, where ’ is the basis function.
e e
This is in tagreemen with theowtGauss’ ws.la
Akwbacdrato unstructured grid methods as compared to FD-TD is that they
need more memory and (cid:13)oating ptoin operations per wn.unkno urthermore,Fthe
computer code for unstructured grid methods is usually erwslo than the code for
FD-TD, measured in (cid:13)oating ptoin operations per second. This is due to the
indirect addressing needed for unstructured grids.
eW think that the best happroac is to binecom FD-TD with unstructured
grid methods toin so-called ybridh methods. Unstructured grids are used near
edcurv objects and around small geometrical details, while structured grids are
used in the homogeneous parts of the computational domain. This binescom
the e(cid:14)ciency of structured grids with the geometric y(cid:13)exibilit of unstructured
grids. uW and Itoh [ WI95 ] erew (cid:12)rst to tpresen a binationcom of the FD-TD
methodandanimplicitFE-TDmethod. Theyevhabeenedwfolloyberalsevothers
[MM98 , KLI97 , SDPP98 , eu99Y , Ryl00 , RB00 , Ril01 ]. Abinationcomofanexplicit
FV-TDersolvand FD-TDaswproposedybRiley andurnerT[ T97R ] and has been
further estigatedvin and edvimpro in [ EL00 , Ede00 ]. The ybridh concept is illus-
trated in Figure 1.4 , hwhicysdispla a ybridh grid for a dielectric circular cylinder.
(This (cid:12)gure is a reproduction of Figure 1 in [ WI95 ], and is used with the tconsen
of the authors.)
Figure 1.4. Aybridhgridinowtdimensons.
One yaw to decrease the umericaln errors is to use higher-order methods. orF
homogeneous domains this is ratherard.tforwstraigh A fourth-order accurate FD-
8 Chapter 1. ductionoIntr
TD method is easily realized. er,evwHo staircasing of boundaries and terfacesin
ysdestro.accuracy Computingthescattered(cid:12)eldfromacircularperfectlyconduct-
ingcylinderouldwresultinlessthansecond-order,accuracybothforasecond-order
accurate discretization and a fourth-order accurate discretization. This is wnsho
for the second-order accurate discretization in Chapter 8 .
oT oidva staircasing errors, ew could again try to use an unstructured grid
close toedcurvboundaries andterfaces.in er,evwHoto get a fourth-order accurate
heme,scewouldwneedatleastathird-orderaccuratetationimplemenofthebound-
ary condition. This is not easy toe.hievac urthermore,Fthe terfaceineenbwetthe
structured and unstructured grids ustm be designed to support fourth-order accu-
racy and not cause umericaln.yinstabilit This is also di(cid:14)cult toe.hievac
It is our opinion that a fully fourth-order accurate method for industrial appli-
cations is not feasible in the near future. er,evwHohigher-order discretizations can
still be useful. orFinstance, ew could use a fourth-order accurate FD-TD hemesc
yawa fromthetransitionregionandsmoothlyertrevtotheeeYhemescclosetothe
transition region. This ouldw egiv a second-order heme,sc but with smaller error
than our tpresenybridh methods.
High-frequency methods
Inbothtime-domainandfrequency-domain methodsfortheumericalnxima-appro
tion of the ellMaxw equations, one needs at least ten mesh ptsoin per elengthvaw
for practical engineering.accuracy It wsfollo that for moderately high frequencies
a largebumern of mesh ptsoin is required to be able to eresolv the problem.
orF eryv high frequencies it becomes impossible to eresolv the problem using
time-domain methods. Here one has to use high-frequency methods, hsuc as the
geometrical theory of di(cid:11)raction (GTD) [ BK94 ] and uniform theory of di(cid:11)raction
(UTD)[ KP74 ]. High-frequencymethodsarebasedonanalyticalximationsapproof
the ellMaxw equations.
1.2 GEMS
The arallelP and ti(cid:12)cScien Computing Institute (PSCI) [ PSC ] is a tercen of excel-
lencefundedybanindustrialconsortium,theedishSwNationalBoardforIndustrial
andhnicalecTtelopmenDev (NUTEK), KTH and Uppsala.yersitUniv PSCI asw
created in 1995. One of the programs within PSCI is Computational Electromag-
netics (CEM). romF 1995 to 1998, the project \Large Scale FD-TD" [ Lar ] asw
conductedwithintheCEMprogramwiththeauthorasprojectleader. Duringthis
project, ewelopdeved a 3D FD-TD code, hwhicew called pscyee .
In 1998, \Large Scale FD-TD" asw succeeded yb another PSCI project, the
hucm more eextensiv General ElectroMagnetic ersSolv (GEMS) project [ GEM ].
ThisaswaedishSwearthree-ycodetelopmendevproject thataswsupportedyban
1.3. Outline and main esultsr 9
eextensivhresearcprogram. Atialsubstanpart of the fundingaswsuppliedybthe
National Aeronautical hResearc Program (NFFP).
The main obejectiv of the GEMS project asw to elopdev a arewsoft suite for
solving the ellMaxw equations. Thisarewsoftsuite aims to be state-of-the-art and
toformaplatformforfuturetelopmendevybedishSwindustryandacademia. The
code will be used in an industrialt.vironmenen
The core of the arewsoft suite is owt ybridh codes, one for the time domain
and one for the frequency domain. The time-domain code is aybridheenbwetFD-
TD, explicit FV-TD and implicit FE-TD. The frequency-domain code is a ybridh
eenbwetMoM, PO and GTD/UTD.
TheindustrialpartnersinGEMSareEricssonMicr evawo Systems(EMW),Saab
Ericsson Space (SES) and Ericsson Saab Avionics (Avionics). Code elopdevers
are PSCI, the edishSw Institute of Applied Mathematics (ITM) and the edishSw
Defence hResearc tEstablishmenA).O(F The industrial partners also etak part in
the codet.elopmendev
1.3 Outline and main results
kgroundBac
The next hapterc tainscon a ernacularv description of ym h.researc The owt fol-
winglohapterscegivkgroundbacinformation on thehresearctedpresenwithin this
thesis. Chapter 3 ersvcotheellMaxwequationsandChapter 4 addressestheFinite-
Di(cid:11)erence Time-Domain (FD-TD) method [ af00T ].
Chapter 5 is a brief description of the GEMS time-domain codes. Chapters 6
through 9 taincontheresultsofymh.researc Theorderofthesehapterscishrono-c
logical, though there has been considerableerlap.vo
arallelizationP
Chapter 6 ersvcoparallelizationoftheleap-frogupdateintheFD-TDmethod. Do-
main decomposition is used to distribute the computations on the nodes of a par-
allelhine,macandunicationcommisperformedusingthemessagepassingterfacein
(MPI)standard. vingHa p nodesofaparallelcomputer,ewsplitthecomputational
domain in p domains of almost equal size. The Cartesian topology yfacilit of MPI
is used to distribute these domains on the p nodes.
eW wsho that perfect scale-up can be edhievac on a parallel computer with
distributed.memory Ontheotherhand,perfectspeed-upisusuallynotpossibleto
obtain. The time to complete a time step onheacnode is proportional to n n n ,
x y z
while the time needed to unicatecomm is proportional to n n + n n + n n ,
x y x z y z
where n (cid:2) n (cid:2) n is the problem size on heac node. As this size decreases, the
x y z
unicationcomm time will no longer be negligible compared with the computation
time.
10 Chapter 1. ductionoIntr
The resultstionedmenin the previous paragrapherewedhievacon an IBM.SP
Similar results can be obtained of for instance a yCra T3E. An exception, "super-
linear speed-up", occurs on computers where hecac e(cid:11)ects are t.dominan orFex-
ample, this happens on a cluster of Dec Alpha computers.
Ontheparallelshared-memoryectorvcomputeryCraJ90,ewdemonstratethat
autotaskingesgivximatelyapprothesameperformanceastheMPItation.implemen
eWalsowshothatonaujitsuFectorvcomputer, itispossibletoehievacmorethat
50% of the peak performance. The performance on theujitsuFectorv computer is
moredeptendenontheproblemsizethanothercomputers. vingHaalargealuevon
the bumern of cell in the x-direction ( N ) will egiv the best performance because
x
it leads to long ectorv lengths.
eWwshothatourparalleltationimplemencanbeusedfortuangargancomputa-
tionsybperformingaone-billion-cellcomputationonanaircraft. Thiscomputation
aswedhievacwith 125 nodes with 160 MHz RS/6000 processors of an IBM.SP
Hybrid time-domain methods
Chapters 7 and 8 ervco a new hniquetec for ybridizationh of the (cid:12)nite-di(cid:11)erence
time-domain (FD-TD) method with methods for unstructured grids. On the un-
structuredgrids,eweitheruseanimplicit(cid:12)nitetelemen(FE)methodoranexplicit
(cid:12)niteolumev(FV) method. Theybridizationhis performedybvinghaa transition
erylaeenbwetthestructuredandunstructuredgrids,wherestructuredandunstruc-
tured cells erlap.vo In 2D, this region is half a cell k,thic and in 3D it is one cell
k.thic Chapter 7 tainscon 2D, and Chapter 8 tainscon 3D.
In Chapter 7 ew wsho that both the FD-FE and FD-FV ybridh for the trans-
ersevmagneticellMaxwequationsaresecond-orderaccurate. Thisiswnshofore(cid:12)v
tdi(cid:11)eren cases: a perfect electric conducting circular cylinder, a perfect magnetic
conducting circular cylinder, a dielectric circular cylinder ( (cid:15) = 4), a diamagnetic
r
circular cylinder ( (cid:22) =4) andacuum.v yStabilit is thoroughly studied ybumeri-n
r
cal tests. The FD-FVybridhis stable for all test cases,videdprothat theystabilit
condition is not violated. The FD-FE ybridh can be unstable when the Crank-
holsonNic method is used for timestepping. eWevha found examples where this
happens. In all these cases,ystabilitcould be restoredybmaking the timestepping
methodtlyslighmore implicit orybhingswitcto theowtstageardkwbacdi(cid:11)erence
ulaform (BDF-2) method. eWalso wsho that our ybridh methods perform ellw on
a test case with a ptoin source and a perfectly conducting allw with 45 degrees
inclination. This case is eryv similar to one of the test cases in the classical paper
yb Cangellaris andtrighW[ CW91 ], hwhic is the most tfrequen reference when the
problems of staircasing are discussed.
InChapter 8 ewwshothatourybridhmethodin3Dcanbeusedtoehievacgood
resultsonagenericaircraftmodelgeometryandtheNASAalmondmodelproblem
yb computing the radar cross section of these objects. The methods wsho super-
linear ergencevcon for a acuumv test case. er,evwHo they are not second-order
accurate. This iswnshoto be causedybtheterpinolation of diagonal compts,onen