See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/224857925
An Empirical Approach for Design of Wideband,
Prob-Fed, U-Slot Microstrip Patch Antennas on
Single-Layer, Infinite, Grounded...
Article in Applied Computational Electromagnetics Society Journal · November 2003
CITATIONS
READS
20
42
2 authors, including:
Deb Chatterjee
University of Missouri - Kansas City
85 PUBLICATIONS 194 CITATIONS
SEE PROFILE
Some of the authors of this publication are also working on these related projects:
Sommerfeld Integrals in Layered Media View project
Characteristic Mode Analysis of wideband microstrip patch antennas View project
All content following this page was uploaded by Deb Chatterjee on 05 March 2014.
The user has requested enhancement of the downloaded file.
APPLIED
COMPUTATIONAL
ELECTROMAGNETICS
SOCIETY
JOURNAL
Editor-in-Chief
Atef Z. Elsherbeni
November 2003
Vol. 18 No. 3
DISTRIBUTION STATEMENT A
Approved for Public Release
Distribution Unlimited
ISSN 1054-4887
The ACES Journal is abstracted in INSPEC, in Engineering Index, and in DTIC.
The first, fourth, and sixth illustrations on the firont cover have been obtained from the Department of
Electrical Engineering at the University of Mississippi.
The third and fifth illustrations on the front cover have been obtained fix)m Lawrence Livermore National
Laboratory.
The second illustration on the front cover has been obtained from FLUX2D software, CEDRAT S S
France, MAGSOFT Corporation, New York.
GENERAL PURPOSE AND SCOPE: The Applied Computational Electromagnetics Society (ACES)
Journal hereinafter known as the ACES Journal is devoted to the exchange of information in
computational electromagnetics, to the advancement of the state-of-the art, and the promotion of related
technical activities;' A primary objective of the information exchange is the elimination of the need to "reinvent the wheel" to solve a previously-solved computational problem in electrical engineering, physics, or
related fields of study. The technical activities promoted by this publication include code validation,
performance analysis, and input/output standardization; code or technique optimization and error
minimization; innovations in solution technique or in data input/output; identification of new applications
for electromagnetics modeling codes and techniques; integration of computational electromagnetics
techniques with new computer architectures; and correlation of computational parameters with physical
mechanisms.
SUBMISSIONS: The ACES Journal welcomes original, previously unpublished papers, relating to
applied computational electromagnetics. Typical papers will represent the computational electromagnetics
aspects of research in electrical engineering, physics, or related disciplines. However, papers which
represent research in applied computational electromagnetics itself are equally acceptable.
Manuscripts arc to be submitted through the upload system oi ACES web site http://aces.ee.olemiss.edu
Sec "Information for Authors" on inside of back cover and at ACES web site. For additional information
contact the Editor-in-Chief:
Dr. AtcfEIshcrbeni
Department of Electrical Engineering
The University of Mississippi
University, MS 386377 USA
Phone: 662-915-5382
Fax:662-915-7231
Email: atef(5).olemis.edu
SUBSCRIPTIONS: All members of the Applied Computational Electromagnetics Society who have paid
their subscription fees are entitled to receive the ACES Journal with a minimum of three issues per
calendar year and are entitled to download any published journal article available at
http://accs.ce.olcmiss.cdu.
Back issues, when available, are $15 each. Subscriptions to ACES is through the web site. Orders for
back issues of the ACES Journal and changes of addresses should be sent directly to ACES Executive
Omcer:
Dr. Richard W. Adicr
ECE Department, Code ECAB
Naval Postgraduate School
833 Dyer Road, Room 437
Monterey, CA 93943-5121 USA
Fax:831-649-0300
Email: rwa(3).altslohaI.net
Allow four week's advance notice for change of address. Claims for missing issues will not.be honored
because of insufficient notice or address change or loss in mail unless the Executive Officer is notified
within 60 days for USA and Canadian subscribers or 90 days for subscribers in other countries, from the
last day of the month of publication. For information regarding reprints of individual papers or other
materials, sec "Information for Authors".
LIABILITY. Neither ACES, nor the ACES Journal editors, are responsible for any consequence of
misinfonnation or claims, express or implied, in any published material in an ACES Journal issue. This
also applies to advertising, for which only camera-ready copies are accepted. Authors are responsible for
infonnation contained in their papers. If any material submitted for publication includes material which
has already been published elsewhere, it is the author's responsibility to obtain written permission to
reproduce such material.
APPLIED
COMPUTATIONAL
ELECTROMAGNETICS
SOCIETY
JOURNAL
Editor-in-Chief
Atef Z. Elsherbeni
November 2003
Vol. 18 NO. 3
ISSN 1054-4887
The ACES Journal is abstracted in INSPEC, in Engineering Index, and in DTIC.
The first, fourth, and sixth illustrations on the front cover have been obtained from the Department of
Electrical Engineering at the University of Mississippi.
The third and fifth illustrations on the front cover have been obtained from Lawrence Livermore National
Laboratory.
The second illustration on the front cover has been obtained from FLUX2D software, CEDRAT S.S.
France, MAGSOFT Corporation, New York.
THE APPLIED COMPUTATIONAL ELECTROMAGNETICS SOCIETY
http//:aces.ee.oIeiniss.edu
ACES JOURNAL EDITORS
EDITOR-IN-CHIEF/ACES/JOURNAL
AtefElshcrbcni
University of Mississippi, EE Dept.
University, MS 38677, USA
ASSOCIATE EDITOR-IN-CHIEF
Alexander Yakovlcv
University of Mississippi, EE Dept.
University, MS 38677, USA
MANAGING EDTTOR
Richard W. Afflcr
833 Dyer Rd, Rm 437 EC/AB
NPS, Monterey, CA 93943-5121, USA
EDITORIAL ASSISTANT
Matthew J. Inman
University of Mississippi, EE Dept.
University, MS 38677, USA
EDITOR-IN-CHIEF, EMERITUS
Ducnn C. Baker
EE Dept. U. of Pretoria
0002 Pretoria, South Africa
EDITOR-IN-CHIEF, EMERITUS
Robert M. Bcvcnsec
Box 812
Alamo, CA 94507-0516, USA
EDITOR-IN-CHIEF, EMERITUS
David E. Stein
USAF Scientific Advisory Board
Washington, DC 20330, USA
EDITOR-IN-CHIEF, EMERITUS
Allen Gllsson
University of Mississippi, EE Dept.
University, MS 38677, USA
EDITOR-IN-CHIEF, EMERITUS
Ahmed Kbhk
University of Mississippi, EE Dept.
University, MS 38677, USA
ACES JOURNAL ASSOCIATE EDITORS
Giandomcnico Amcndola
Univeisita' della Calabria
Rcnde, Italy
Chun-Wen Paul Huang
Anadigics, Inc.
Warren, NJ, USA
Edmund K. Miller
LASL
Santa Fe, NM, USA
John Bcggs
NASA Langley Research Center
Hampton, VA, USA
Todd H. Hubing
University of Missouri-Rolla
Rolla, MO, USA
Krishna Naishadham
Wright State University
Dayton, OH, USA
John Braucr
Ansoft Corporation
Milwaukee, WI, USA
Nathan Ida
The Univeisity of Akion
Akron, OH, USA
Giuseppe Pclosi
University of Florence
Florence, Italy
Magda El-Shcnawce
University of Arkansas
Fayetleville AR, USA
Yasusbi Kanai
Niigata Institute of Technology
Kashiwazaki, Japan
Vicente Rodriguez
ETS-Lindgren
Cedar Park, TX, USA
Pat Foster
Microwave & Antema Systems
Gt. Malvcrn, Wore. UK
Leo C. Kempel
Michigan State University
East Lansing MI, USA
Harold A. Sabbagh
Sabbagh Associates
Bloomington, IN, USA
Cynthia M. Furse
Utah State University
Logan UT, USA
Andrzcj Krawczyk
Institute of Electrical Engineering
Warszawa, Poland
John B. Schneider
Washington State University
Pullman, WA, USA
Christian Hafncr
Swiss Fedci<il Inst. of Technology
Zurich, Switzerland
Stanley Kubina
Concotdia Univeisity
Montreal, Quebec, Canada
Abdel Razck Sebak
University of Manitoba
Winnipeg, MB, Canada
Michael Hamid
University of South Alabama,
Mobile, AL, USA
Samir F. Mahmoud
Kuwait University
Safat, Kuwait
Amr M. Sharawce
American University
Cairo, Egypt
Andy Harrison
Radiance
Huntsville, AL
Ronald Marhcika
Ohio State Univeisity
Columbus, OH, USA
Norio Takahashi
Okayama Univeisity
Tsushima, Japan
THE APPLIED COMPUTATIONAL ELECTROMAGNETICS SOCIETY
JOURNAL
Vol. 18 No. 3
November 2003
TABLE OF CONTENTS
"Bistatic Scatterer and Antenna Imaging"
S.Inge
136
"Reducing Electromagnetic Coupling in Shielded Enclosures using a Genetic AlgorithmFinite-Difference Time-Domain Solver"
R. I. Macpherson and N. J. Ryan
144
"Impact of Some Discontinuities on the Convergence of Numerical Methods in
Electromagnetics"
M.M. Bibby
151
"Wireless Propagation in Non Line of Sight Urban Areas using Uniform Theory of Diffraction"
G. A. Ellis
162
"Automatic Calculation of Band Diagrams of Photonic Crystals Using the Multiple Multipole
Method"
J. Smajic, C. Hafiier, and D. Emi
172
"Analysis of Conventional and Novel Delay Lines: A Numerical Study"
O.Ramahi
181
"An Empirical Approach for Design of Wideband, Probe-Fed, U-Slot Microstrip Patch Antennas
on Single-layer, Infinite, Grounded Substrates"
V. Natarajan and D. Chatterjee
191
"Utilizing Particle Swarm Optimization in the Field Computation of Non-Linear Magnetic
Media"
.
-.
_
A.A.AdlyandS.K.Abd-El-Hafizr.....-f.
202
"An Approach to Multi-Resolution in Time Domain Based on the Discrete Wavelet Transform"
C. Represa, C. Pereira, A.C.L. Cabeceira, IBarba, and J. Represa
210
© 2003, The Applied Computational Electromagnetics Society
136
ACES JOURNAL, VOL. 18, NO. 3. NOVEMBER 2003
BISTATIC SCATTERER AND ANTENNA IMAGING
Steve Inge (in-j)
Broadcast Communications Ltd.
P.O. Box 98, Wellington, New Zealand
ABSTRACT: Two computing techniques to create
radiating centers on a conducting body. The currents
an image of the radiating centres on an antenna or
on the body are found using either a computer
scatterer using Fourier optics is presented.
program or measurements, for a particular excitation.
The excitation may be a plane wave being scattered
by the body, or sources on the body for an antenna.
Steve Inge ilieil on October 24, 2002 at the age of 61 after this
This paper describes two simple approaches to
article was submitted for publication. The reviewers accepted
imaging a radiator.
this article subject to revisions. Since that is not possible, this
article is published as submitted.
Brief Overview
The late Ian McEnnis, Chief Antenna Engineer at Broadcast
Communications Limited in New Zealand, and one of Steve's
colleagues, notified us of his death and had the following
remembrances of Steve, which we print below.
The radiation emanating from any antenna (or
scatterer) may be visualised by simple operations
upon the antenna currents. First choose the far field
Steve had been with our organization far 42 years. His immense
direction from which to look at the antenna. Create a
intellect enabled him to become a "guru" in any area he focused
plane surface at the antenna, normal to the look
on. Most of BCL's in house engineering software has been
developed by Steve. BCLIPPS (BCL's Interactive Planning and
direction. Project all of the antenna currents onto this
Propagation Software) was Steve's crowning. It has formed the
surface, adjusting the phase as necessary, to maintain
basis of our coverage work for the past 10 years and has enabled
the same phases and polarisation in the far field.
us to be a leader in this area. Steve was the engineer who solved
Filter the current image on the surface by removing
the most difficult of problems and he was the one his colleagues
turned 1o when they were having difficulty. His mathematical
all
spatial
frequencies above
ability was awesome and once he got his teeth into a problem he
wavelength.
would not let it go until he had it sorted.
There is nothing Steve
brightaess in a 2D graphic image, "which will show
Uked more than to share his knowledge with others. Having said
which parts of the antenna are "glowing" when
all that, in the end Steve was just a "good bloke" who will be
sorely missed by many.
Display
the
one cycle per
current
intensity
as
viewed from this direction. This is repeated for both
polarisations, which may be combined into a single
total radiation image if desired.
Introduction
Bistatic k-space (BSKS) imaging from antenna
Imaging ttieory
currents was introduced by John Shaeffer et al (Ref
1). It is a technique to graphically show the active
How does it work? By flattening the currents on to a
surface we create a two dimensional current plane. If
1054-4887 © 2003 ACES
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
137
we now take the Fourier transform (FT) of this
The image is a view of the far field antenna radiation
cuirent surface, we create a new description of the
toward the look direction. For scatterers, the look
currents. The FT creates a set of sinusoidal waves
direction may be the same
equivaknt to the original currents. Since the currents
illumination, which is the monostatic (single antenna)
occupied a finite area on the surface, the FT range is
radar case. Similarly, when the look direction is
infinite. Each point on the FT surface represents a
different to the illumination direction, we have the
particular uniform linear phase advance, and thus
bistatic (separate antennas) radar case.
as the incident
represents a plane wave propagating in a single
direction. This is known as a plane wave spectrum of
Imaging process
the radiation. The Fx, Fy positions on the spectrum
are the direction cosines of the plane wave direction.
1. Given the locations and currents of wire segments,
is
from a NEC (Ref 4) analysis for example, calculate
Fz = J\-(Fx^+Fy^) The centre of the plane (0,0) is
the contribution to far field (Efar) for each
The
z
constant
component
phase,
of
which
propagating toward the
the
creates
radiation
a
wavefront
look direction.
When
Fx^ + Fy^ = \, the wave direction is normal to the
segment, in a chosen (theta, phi) radiation
direction. The two polarisations are treated
separately. The contribution from a single segment
is
look direction. All points on this unit circle
correspond to a phase advance of one cycle per
wavelength. The spectrum outside this circle gives
direction cosines > 1, indicating imaginary angles.
AEp =y30A/exp(yAR.S){(R«AL)R-AL}.p (1)
where p is the chosen polarisation unit vector, / is
This energy does not radiate, but represents
evanescent waves which eventually return to the
the current, AL is the segment length vector, R is
the far field direction unit vector.
antenna as reactive energy. (Ref 3)
Sis the
coordinate vector of the segment centre and it is the
The antenna image is created by projecting the
wave number In IX. Note that the 1/r distance
currents onto a normal plane, and removing spatial
factor is omitted, making the units in volts. This is
frequencies greater than one cycle per wavelength.
the same operation as calculating the far field
The latter process may be achieved by setting the FT
radiation pattern, without summing the voltages."
spectrum to zero at all points outside the one cycle
The code can be extracted from a suitable antenna
per wavelength circle and inverse FT, or equivalently
program far field radiation pattern routine.
by
convolving
the
current
surface
with
the
J\{7ip)l{7rp) point spread function (PSF), where
J] 0 is the Bessel function of the first kind and order
2. Each segment is replaced with an isotropic source
in the same location, with the magnitude and phase
of the Efar contribution to the chosen polarisation.
one, and rho is the radius. This point spread fiinction
is known from optics as the Airy disk.
By placing this light source at each point, we create
an image, viewable firom any direction. At this
point all the detail can be seen, so that wires with
138
S. Inge: Bistatic Scatterer and Antenna Imaging
rapidly varying phase from a particular viewpoint
the currents from broadside, we can see the
are visible. This does not agree with radiation
alternating pattern. However if we apply the blur
filter, we find that at all internal points along the
reality.
wire, the cyclic currents integrate to zero. Only at the
3. Filter the image to remove spatial frequencies
ends is the last current pattern not cancelled by its
above one cycle per wavelength. This physical
neighbour, so we see a bright spot centred on each
limit is due solely to the wave nature of light. Note
wire end. This illustrates the maxim "radiation is due
that even the tiniest point source appears as a
to the acceleration of charge". See the first example
circular blob, surrounded by faint rings, known as
below.
the Airy disk. The Airy disk has the form of
y, (M) /(«), where Jj is the Bessel fiinction of
Taper windows
the first kind and order 1. The resolution limit is
the Raleigh criterion of 0.61 wavelength, which is
where the maximum of one Airy disk overlaps with
Tapering of the spatial window illumination reduces
both resolution and sidelobe levels. The taper
fiinction used is (1 — /? )^ where rho is the fraction
the first minimum of another.
of the window radius, andj? is the taper exponent. In
Starting with the Efar elements, we Fast Fourier
Transform (FFT) the (ID, 2D, or 3D) space,
2D,
the
resulting
J_+l («)/«''
blur
fiinction
form
is
giving both reduced resolution and
remove all spatial frequencies above one cycle per
wavelength, and inverse FFT. The wanted spatial
frequencies are found in the line, circle or sphere
sidelobe levels as p increases, p = 0 for uniform
illumination, p = 1 for parabolic taper, etc.
sectors in the comers of the transformed
Example radio wave photographs
matrix, the remainder are set to zero.
4. Now the image appears as it would in our
"microwave
microscope".
Note
that
optical
microscopes have the same limit to magnification,
The grey scale is in decibels. The number at top left
of each image is the maximum "brightness" across
the image. The look direction for all the images is
where a point is reached where more magnificati(m^
from the zenith, 0 = 0, (p =1". The small-tp offset is
just gives more blur. Hence the move to higher
to ensure that the polarisation is properly defined. All
frequencies and electron microscopes. These
these images used NEC-2 to calculate the current
images may be named radio wave photographs,
distributions. Unless specified, horizontal (Ephi)
since they show the luminosity of each part of the
polarisation is used. The faint grid lines are spaced
antenna or scatterer.
one wavelength apart. A multidimensional mixed
radix FFT is used (Ref 5), with each image as
Example Consider a long wire with a sinusoidal
200x200 pixels. In the largest image spaces the pixel
standing or traveling wave current pattern. If we view
size is Xy20.
139
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
Table 1. Window Taper Performance
Tapsf Fxpfi'ioni )• ."
0
!
2
i
4
5
6
7
Rd;iii\!.>heas« width
1
1 lO
I 21
lU
1 47
1 62
1.78
1.54
-17
-24
-.10
-J(^
-41
-46
-50
• 54
S)delr»heiC'.cliiB
Long v\iri:, end fed
Th\<, iningc is trikon \mn bTo:Mh\ik to a ?
ur.\c!\:ni>ih wire, icd near ihc K^iiom ond.
'['ho biiMdi-iLk" ludiiiiiiMi puiiL-rit is i'onr.ai in
llv jrucifercncc of ihc radtaiMn Trom tlic tw*^
end piMIV.":.
l.i.Mii; uiri.'. end fed. radmliui! in mam lobe.
Ths.N iiiMgo s.h<vAN ilic radiaiian in tiio
direction of the main lobe, froivi nil nlonc ijie
mrc.
The coiurilniuon dcctcuses from the
rccd |>oinl u.nMiid liic liir end.
I lnuii;i.' si/e II.','. .'-i.niiirci
S. Inge: Bistatic Scatterer and Antenna Imaging
LoiiL' wire, centre fed. broadsitlc racJiiUMHi.
"['hu: i.s ilic .same 5 wtivcIcngSh wire, biil
centre foci
"Hie current disamtiiuiity causes
ladinliflii i'mai (he feed poini as wdl as the
end |>ointe. t lma<ic size -= 10). square)
Loop, raidiulkm towards far riglU.
This is a ."> x 3 /. squnre liwp anlcnna, in the
Iwrizonlal plane.
feed |X)inl. and
'I'hc r.ndiafion from the
[TOUT
Ihe finif coriier.s h
visible. (Image size 6'/.)
Loop, radiation roward bottom ()f page
In ibis direction there is ntdiatioii orilj- front
ihe corKer.<;.
f Imaiie size 6Xi
140
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
141
S^iiiiirc p\:\\t:. cdjio rioniinl lo illiimiirjiion
[1)1.'' !.'> ;i Tv. stiiinrc Pl-X." piaSc in
lii>n/<iii!ai pkmc
tAdiiallv niodcilcJ a.s ii
ynd (M* wires ai 5.10 siKtciiii;;
fiuiu.v
tho
radar
IIK'
cro?.<.
The imaiic
scciion
uhon
illiinii!utL\l wiih .1 plane w:nc friMn iltc iiiili!.
The brii;hl line
.-IH^AS >,\1ICIC
the \sa\c K i'lcld
is shi>ncd out h> the condiiciiivj pblc.
AKo
ihorc IS MHTiv' ladiniion Iroiri the irnilmii
ftirners. Kntc llie iiumlvr ol' spurious iiniii^c
puinis. LHcn h;?\on;l ilie plmc. duo to ihe
isntiperA! I cycle per ^ filter dmnge <iyo < '/.)
Sqiiiirc pkile. iHuminaicd daiiionnllv.
T'liifi slunv: ihe rndnr cmss section v.lien ihe
same ?/. square pliitc is iliiimimUod h a pliiiie
^\a ve iirn\ ing dMut>nidl> ih>m llie riyht.
The
iviain cornrihution to ilic radar echo comes
iVoiii
' llic
(Iniaje si/e6/.l
outer
• corners.
142
S. Inge: Bistatic Scatterer and Antenna Imaging
Arr^himedian spiral, Htheta ixMnrisahon
Because; lliC itiiagc
liitcring pnvcss is
separated fmni t\\<: radiaiion iiUeyral. i! is
prtSKible to make near field iniases, a?^ well as
far ik'lii. 'Iliis inuigc is a near field innagc of
ati horixonlal Arcliinicdian spiiat aiiteiiriu
x'ie\wd from ateve. Tin's \% a near ficfd view
of the image formed from ihc far field
radiation. Hie spiral diameter overatl is Imlf a
\vas elcniilh.
Resolution 0.!/„<.yde.
latni: tiulci rings are
The
filiering iutcfacLs.
(fmaue mze 2k)
Airhinisdian spiral, Epiii |.x>larisition
riK same antcnmi, but EpIii. 'fogelhcr these
l'.u> imaucs shsHv ilie circular |x>larised
radiaiion ooining from the spiral turu.i wheri!
llsi"
turn
one wavelength.
t'Imauesfze2A'(
circumrerence
et|iKils
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
143
Radio Camera
short course, and also for the many conversations
since that helped further my understanding of the
An alternative imaging process is to emulate a
subject.
camera. In this process a plane is created in the far
field of the antenna/scatterer, and the E field is
Conclusion
calculated at many points on the plane (usually in a
rectangular grid). This array of complex field values
Two practical radiator imaging techniques have been
is a numerical hologram. A reference wave is not
described. These are thought to be new, but
needed, as the phase is recorded in the complex E
similarities to near field measuring techniques and
field value. This numerical hologram may be
holography suggests that they have simply been
reconstructed using a camera lens, or equivalently by
rediscovered. However, it is timely to bring imaging
taking the FT of the 2D complex field values. The
to the attention of CEMists as a useful tool in
image is formed at one focal length behind the
understanding radiation. The author has used the first
camera lens. The radiation integral has already
technique to examine currents induced in a support
removed the evanescent energy, so it is not necessary
tower, which caused an unwanted null in the
to filter the image any further. Note that the point
radiation pattern. Viewing the currents only confused
spread function may be larger than the theoretical
the situation, since most of the currents did not
Airy disk, since the finite lens aperture (or grid size)
radiate in the null direction.
can only worsen the resolution. In microscopy this
reduction is known as the "numerical aperture".
References
Using NEC for example, we can evaluate the far field
over a rectangular grid of points on a plane normal to
the look angle. The difference from the previous
1. "Bistalic it-Space Imaging for Electromagnetic Prediction Codes
for Scattering and Antennas", John F. Shaeffer, Kam W. Horn,
Craig Baucke, Brett A. Cooper, and Noel A. Talcoft, Jr. NASA
imaging process is that the plane is in the far field.
The antenna/scatterer should be rotated (with a GM
Technical
Paper
3569
(1996).
http://techreports.larc.nasa.gov/ltrs/PDF/NASA-96-tp3S69.pdf
card) so that the desired look angle is moved to the x
axis. The near field command (NE) is used to create a
2. High Resolution Radar Imaging, Dean L. Mensa, Aitech House
Inc., 1981. Describes a related technique to BSKS imaging.
planar grid of field values in the YZ plane. There are
still two-polarisations required. After gathering the.
output values into a suitable matrix form, the FT
3. "Introduction to Fourier Optics", by Jos^h W. Goodman,
McGraw Hill, 1968. Chapteis 2 & 3.
image may be formed using Matlab or Mathematica
4. "Numerical Electromagnetic Code (NEC) - Method of
etc.
Moments", Naval Electronic Systems Command Technical
Document 116, Part III Users Guide, GJ Burke and AJ Poggio,
Acknowledgement
1977.
5. Subroutine FFT, multivariate complex Fourier transform,
I wish to thank John Shaeffer and Edmund Miller for
computed in place using mixed-radix fast Fourier transform
introducing the concept of imaging at their ACES
algorithm. R. C. Singleton, Stanford Research Institute, Sept.
1968.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
144
Reducing Electromagnetic Coupling in Shielded Enclosures
using a Genetic Algorithm - Finite-Difference
Time-Domain Solver
Russell Iain Macpherson' and Nick J Ryan^
Department of Engineering
King's College
University of Aberdeen
Aberdeen AB24 3UE
United Kingdom
r.macpherson@ieee.org', n.j.ryan@eng.abdn.ac.uk^
Abstract: Comprehensive shielding in modern electronic
equipment may lead to resonant behaviour within the
equipment enclosure. This paper presents a method for
optimising the placement of sources of electromagnetic
(EM) energy and susceptors to this EM energy within an
enclosed resonant cavity. The source and susceptor sire
placed on a dielectric slab within a perfectly conducting
enclosure to reduce the level of EM coupling between the
two. Optimisation is facilitated using a genetic algorithm
coupled with a finite-difference time-domain solver. Results are presented for single objective optimisation and
multi-objective optimisation cases.
There is obviously an optimal size and shape of enclosure and component layout however, exhaustively placing
components inside enclosures and then computing the resultant fields is a task that that would require a massive
undertaking on behalf of the designer. A far better approach to component placement is to use some kind of optimisation method to place the components to achieve a
desired radiation level. This paper describes the novel use
of a genetic algorithm (GA) and a finite-difference timedomain (FDTD) solver to optimise source and susceptor
placement inside a perfectly conducting structure, building upon initial work completed in [2].
Genetic algorithms Eire briefly introduced in section 2 emd
then the Finite-Difference Time-Domain method is introduced in section 3. Section 4 describes the setup of the
The explosive increase in the use of electronic equipment computer simulation and section 5 discusses the results
in the information age has led to electromagnetic com- of these simulations, finally section 6 presents conclusions
patibility (EMC) and electromagnetic interference (EMI) firom the work.
issues becoming more important to designers of electronic
equipment. These issues must be considered at the time
of design and not once designs are at the prototype stage. 2
Genetic Algorithms
Higher clock speeds and faster switching transitions lead
to greater levels of electromagnetic emissions, while higher GAs are a stochastic search mechanism with their operacomponent integration and lower power demands lead to tion firmly rooted in natural selection and survival of the
greater sensitivity. So at once systems are becoming more fittest, [3] and [4]. GAs have been shown to provide roprone to generating and also more sensitive towards EMI. bust search and optimisation in complex spjices, [3] and
[5]. They use simple operations on a population of individAn area of interest to the authors is electronic enclosures
uals, which lead to an emergent evolution of an individual
and the placement of devices inside these enclosures. Modor individuals. Each individual in the population repreern electronic items have many different sources of EMI.
sents a potential solution to the given problem scenario
These sources are often placed inside a rectangular shaped
and as such is evaluated. After an individual has been
metal box to limit the EM emissions from the equipment.
evaluated a figure of merit (FoM) is attributed to the inCare is taken to ensure the shielding is as comprehensive as
dividual. This FoM is a measure of how well the individual
possible thus apertures are kept to a minimum and covered
solves the problem. GAs can lead to the optimal solution
in metallic blanking plates. Such a shielded, rectangular
to a problem, more often however, they are used to openclosure has a good chance of forming a resonant cavity
timise a problem towards an improved solution facilitatso that field strengths generated by the circuit may well
ing a trade-off between excessive computational time and
be enhanced once the circuit is mounted within the casing.
meaningful results. The GA process can be summarised
Normally the circuit may have a number of sources within
graphically as shown in figure 1.
the enclosure, sdl of which add to the EM fields within;
changing the positions of the various sources will change The GA methodology used in this work is the micro
the amount of resonance and constructive interference, [1]. genetic algorithm (/iGA). The /JGA technique has been
There are also likely to be a number of susceptors at mul- shown to reach the optimal area of multi-modal solution
tiple locations throughout the circuit.
spaces earlier them conventional GA methods, with min-
1
Introduction
1054-4887 ©2003 ACES
145
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
of local optima. Once the population has converged to a
determined optimal value then further GA operating is no
longer needed and the GA can be halted. Due to the small
OEHERATEIHITIAL POPULATION
population size used in the fiGA premature convergence is
often seen to be happening, to prevent this a restart mechEVALUATEFoM FUNCTION
anism is used. This restart mechanism involves checking
for convergence in the current generation of individuals
and then restarting the next generation with only the elite
REPRODnCrlON
individual and new randomly created individuals. The use
of a restart operator ensures random search takes place
and leads /iGA away firom local optima. Convergence is
checked for in this application by measuring the changes
in the chromosomes, the individuals, between generations.
Once there "are few changes between generations then the
population can be considered to be converged, the setting
EVALUATE FoM FUNCTION
of this limit is a parameter that the user has control over.
An important component of any GA code is the device by
which random numbers are generated, or to be more correct pseudo-random numbers. This is usually a portable
Figure 1: GA Process Flowchart.
pseudo-rsuidom number generator (PRNG) that produces
the same sequence of random numbers for a given seed
imal computing time, [6]. It has already been demon- value. The quality of the PRNG is an important factor.
strated that the /iGA can be successfully linked with The area of PRNG research is vast and not going to be, for
an FDTD solver, [7], for optimi.sation in EM problems. brevity, discussed here. Common references on the subject
Initially in a GA a fixed size population is created and are [9], [10], [11] and [12].
populated with randomly generated individuals of a fixed
length, this length is determined by the encoding of the A clear difference can be seen here between optimisation
problem parameters. The population size for the /JGA is based on calculus methods and GA based optimisation.
generally kept small e.g. five, this is different from the GAS make use of interim performance in the optimisation
classic GA where the population size is typically much problem, calculus based methods are only concerned with
larger, in the order of 100s. These individuals are all optimisation toward an overall optimal point and do not
po.ssible solutions to the same given problem. They axe "remember" any interim performances.
as.se.ssed and each one given a score, the FoM, which is
then used to determine the likelihood of reproduction into
3 The FDTD Method
the next generation. The FoM is analogous to fitness in
the natural world where fitter individuals survive longer The FDTD method, [13], [14], has become one of the most
and hence have a better chance of continuing their ge- popular methods for solving Maxwell's equations. It is a
netic lineage. The particularly successful, higher scoring volumetric domain decomposition technique that is second
individuals may be reproduced more than once into the order accurate in space and time, easily implemented in
next generation. In the fiGA used here tournament se- software and accurately models the physical world. It is a
lection is employed, this is perhaps the most effective for widely used method for EMC/EMI work, [15], as well as
many application types as it has been shown to provide rcidar, bioengineering and antenna analysis. The FDTD
better convergence towards a solution in the initial stages method is described widely in the literature, [16], and so
of optimisation, [8]. Tournament selection works by ran- only a very brief description is given here.
domly choosing N members of the population and in this
instance as in many others N = 2. These individuals are Meixwell's equations, in differential form, equations (1) then pitched against each other to determine which has the (2), are simply modified to central-difference equations,
better FoM, the winner of the tournament being selected discretised, and implemented in software. The equations
to be a member of the new population; this is repeated un- are solved in a leap-frog manner, [14], i.e. tKe electric field
til the new population is filled. In conjunction with tour- is solved at a given instant in tiiue, then the magnetic field
nament selection elitist reproduction is also employed, this is solved at the next instant in time, and the process is
guarantees that the best individual from a population is repeated over and over again. In equations (1) and (2) H
present in the next population, hence preserving the best and E have their usual meanings.
current solution to the given problem and maintaining a
good genetic stock. After reproduction the individuals undergo crossover. Here two randomly picked individuals are
mated together, swapping information between their chromosomes. In a classic GA mutation is also applied, howdE
1_ „ (T „
ever in the f(GA mutation is turned off. Mutation serves
-5-=-V X H --£;.
(2)
at
i
e
to alter the genetic makeup of these new offspring with a
small fixed probability. Once again mimicking nature, muUsing a three dimensional Cartesian co-ordinate system
tation can lead to either a detrimental or beneficial effect
we can now write out the vector components of the curl
on performance. Mutation allows, in a classic GA, random
operator. Given below is only one of these components,
search to take place and hopefully leads to the avoidance
namely the x component of the H field,
f = -(ivx^-.jH),
(1)
Macpherson and Ryan: Reducing EM Coupling in Shielded Enclosures
146
4.1
-400-
Physical Problem Description
The problem geometry is that of a simplified metal box,
which has perfectly electrically conducting (PEC) walls,
an internal structure and DS representation. The internal
structure is modelled as two PEC sheets in the top corner
of the geometry. The DS can be thought of as a representation of the substrate of a printed circuit board (PCB).
The problem geometry is shown in figure 2. The DS has a
relative permittivity, fn of 4 and possesses unity magnetic
permeability fi^.
4.2
y
^
PEC Sheets
nelcc-trlc StQb
Figure 2: Enclosure geometry showing the representations
of the internal structure and dielectric slab, all sizes are
in millimetres, mm.
dih _ 1 fdEy dE,
, \
at ~ n\dz ~ dy -P"')^
(3)
where n is magnetic permeability and />' is a magnetic loss
parameter. It is the complete set of coupled partial differential equations, six in total, that are the fundamental
basis of the numerical FDTD algorithm. The discretised
form of equation (3) is shown below:
m
=
MEDIAa,\i^^,,,
HAljl
=
D«('")«x|.,-|
+D\,{m) \
,p ,n
^
p ,„
^ I,
(4)
where n is the time step under consideration and i,i,k
are the three orthogonal spatial co-ordinates. The integer
array MEDIA defines material conditions for each field
vector component. This allows the medium at each point
to be mapped out. The arrays Da,Di, are magnetic field
update coefficients, there are corresponding electric field
update coefficients also. A full explanation of all of these
equations is presented in [13].
4
Simulation Setup
The initial results are obtained from using the GA to place
a source and susceptor point relative to each other on the
surface of a dielectric slab (DS) to achieve the lowest electromagnetic coupling between the two. The complete simulation is implemented in Fortran 77 compiled on a SUN
ENTERPRISE server with 512MB of RAM utilising 4 of
8 processors.
pGA and FDTD Setup
The ^GA used is based on the implementation by Carroll,
[17], modified to accommodate the given task of moving
the source point and susceptor point on the DS. The two
points are specified in a three dimensional Cartesian system with the y ordinate being held constant as this represents the surface of the DS. Binary encoding is used in
the ^GA and this leads to a chromosome string length
of 24 bits, 4 bits per i, y and z co-ordinate of the source
and susceptor. It should be noted that it would be possible to omit the y ordinate from the chromosome however the memory saving would be insignificant especially
when compared to the coding effort required to convert
the GA for this one specific case. Care must be taken to
avoid the y ordinate being altered during reproduction and
crossover, this is achieved using a uniform crossover operator. Uniform crossover, [18], has also been found to give
faster convergence than single point crossover in a /iGA,
[6], [19]. The population size is maintained at five. These
co-ordinate values are passed to the FDTD solver which
computes the resulting field distribution inside the enclosure. The peak electric field at the susceptor point is returned as the FoM to the ^GA, no fitness scaling is used as
tournament selection is the method of selection employed
in this GA implementation, [18]. The fiGA attempts to
minimise this FoM, i.e. minimise the peak electric field at
the susceptor point.
The FDTD solver is based on the equations given in [13].
The DS is one cell thick and has one cell of free space
between itself and the enclosure wall. The PEC sheets
are modelled as infinitely thin. A soft Et direited sinusoidal, continuous wave, ideal Hertzian dipole of 2.Q5GHz
frequency and amplitude of lOVm"* is uSfed as the source
in the FDTD simulation. The domain is discretised into
dimensions of 7.2727jnm in the x and z directions and
7.2mm in the y direction. Based on these dimensions and
the material within the computational domain the time
step size for the solver is set at the Courant limit, in this
instance approximately 14ps.
Initially two optimisation cases are run; a non-lossy and
then a lossy dielectric slab, in the lossy case the conductivity is set to 0.045m~*. After this a multi-objective simulation is run where two objectives are to be optimised,
these being the reduction in coupling between source and
susceptor points at two specific frequencies.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
147
5
Results
5.1
Single Objective Optimisation
The progress of the optimisation process can be seen in
figure 3 where the results of both the non-lossy and lossy
simulations are presented. This figure shows the evolution
of the best individual In the simulation, the elite individual. While this individual remains the best the FoM
remains at the same value, when a superior one is bred, a
step change occurs indicating the presence of a new elite
value.
"
'
1
,
K
rfo n_l mijf
^
u
the DS are indicated; for the non-lossy case these are at
(2,16) and (3,35) respectively, and for the lossy case the coordinates are (2,35) and (23,11). These co-ordinate values
given are of course on the x and z axes as the y ordinate
is constant. These false colour plots clearly show that the
^GA has indeed found good solutions but not the best
solutions, the global minimum point in figure 4 is at (2,2)
with a value that is 2.ldB down on the position found by
the fiGA. In figure 5 the global minimum is at (12,14) with
a value that is lAdB down the /tGA position. The final
values found by the ;tGA are 57.GdB and 64.6dfl down on
the global maximum for the non-lossy and lossy cases respectively. It should be noted that the patterns presented
in these figures are ones of over 1 million possible patterns
due to source and susceptor placement, finding the global
minimum on each of them would require direct evaluation
of each pattern, the /iGA vastly cuts down on the number
of evaluations required to arrive at an acceptable result.
11
\^
>\ ^
^*^in.m
0
20
«
60
N
GeKnritm Nuirber
Figure 3: Generation optimisation of both the non-lossy
and tossy cases.
The final value for the non-lossy cajse is 0.0541Vm"' from
an initial value of 0.1172Vm"' and the final value for the
lossy case is 0.02401^m"' from a value of 0.0685Vm-'.
Both of the curves are normalised to the value of
0.1172Vin~' for ease of comparison. As expected the values of field strength for the lossy case are lower than those
for the non-lossy case as in the lossy case energy is dissipated in the dielectric slab. Figure 3 also shows us that
the lossy case reaches its best FoM at generation 66, earlier than the non-lossy case which reaches its best FoM
at generation 91. This is a random dynamic of the GA
and attributable to the solution surface topology and the
choice of PRNG used with the /iGA. The number of generations was limited to 100 as experience, with this code
has shown that there is a minimal return on computing
time by exceeding this point. Testimony to this is given
by the fact that running the code to 200 generations gives
only a final value of 0.0507Vm~^ for the non-lossy case
and 0.0204P'm~' for the lossy case, a marginal increase
in performance for a doubling of computation time. This
embraces one philosophy of using a GA, naihely to produce an acceptable impro\'ement toward an optimal point
for a limited amount of effort.
X co-Drdinato nuntof
Figure 4: False colour plot of final field distribution on
surface of dielectric slab after non-lossy simulation. The
X and z axes are marked in FDTD cells. Source and susceptor positions are indicated on the figure.
5.2
Multi-objective Simulations
Multi-objective optimisation(MO) is often sought in practice as often compliance in one particithir objective upsets
the compliance of another objective. MO aims to achieve
a solution that can not be improved upon.-simultaneously,
in each of the objectives. This is called a Pareto optimal
solution and the set of all these solutions is called the
Pareto optimal set, [20]. For the MO simulations some
changes had to be made to the GA-FDTD code. As two
objectives were being optimised a method of measuring
their fitness at one and the same time is required, to accomplish this the method of objective weighting is used,
The final field distribution on the surface of the DS can [20]. This method is explained mathematically by equabe seen in figure 4 for the non-lossy case and in figure 5 tion (5),
for the lossy case. These plots are generated by applying
a peak hold to the Ez component at each location on the
surface of the DS throughout the simulation.
Z = y^Wifi{x) where x € X.
(5)
The final positions of the source and susceptor points on
Macpherson and Ryan: Reducing EM Coupling in Shielded Enclosures
148
before the optimisation process begins, the solid line in
the figure displays the final DFT after completion of the
optimisation process. A clear diiference in the two curves
at the respective objective frequencies can be seen. Also,
it can be clearly seen that the field values at the objective
frequencies, on the final DFT, are considerably lower than
at other frequencies, within the same response, excepting
of course at the tails of the response where there is minimal energy from the source. The values of the objectives
on the final response are 54dS down for the l.002GHz objective on the max value in the response and the value at
2.00iGHz is 44dB down. Relative to the initial DFT the
final response in down 27dB at the 1.002GHz objective
and down 35dB at the 2.004Gffz objective.
xn-fl(diiut«numlMr
Figure 5: False colour plot of final field distribution on
surface of dielectric slab after lossy simulation. The x and
z axes are marked in FDTD cells. Source and susceptor
positions are indicated on the figure.
In equation (5) Z \s a. scalar valued variable that is a
weighted sum of the individual objectives, where there are
N objectives in total. The function / is the function that
returns the FoM, in this case the FDTD solver, and x represents the parameters of the function, co-ordinate values
in this case. The feasible region of solutions is denoted by
X. The sum of the individual weights adds to 1 and each
weight lies in the range 0 to 1 i.e. 0 < u;; < 1. It is the
scalar value Z that is optimised by the GA. The method
of objective weighting is an easily implementable scheme
that produces a solution from the Pareto optimal set.
Figure 6: Time series of the pulsed source, upper graph,
and its DFT, lower graph. This is the source characteristic
The problem setup for the MO optimisation involves the
used in the multi-objective optimisation simulations.
seime geometry as previously with the only change being
the source. A Gaussian pulse modulated onto a sine wave
is now used as the source. This gives a symmetrical spectrum about the carrier frequency, no DC component, and g
ConclusionS
a bandwidth controlled by the length of the pulse. Its
mathematical form given in [13] is described by equation A method of using a /iGA in conjunction with a FDTD
(6),
solver to facilitate electromagnetic optimisation has been
shown. The conjoining of the two codes allows a difficult
design problem to be tackled, namely that of radiating
E: = Eoe ^"• '"'• "Jsin(27r/o(n-no)At).
(6) and susceptible component placement within an arbitrary
metallic-structure to minimise coupling. The /JGA techThe Gaussian pulse is centred ziround frequency /o at step nique is used as this gives a significant reduction in comno with a 1/e characteristic decay of nitcny steps, and At' putation time with little loss in performance of the final
is the step size. The time series and frequency response optimised result. Ensuring that the ^GA finds the global
of the source is shown in figure 6. The centre frequency optimum in difficult optimisation problems is an area that
for the source is chosen at l.SOOGifz and the objectives needs further attention but the ability of the technique
for optimisation were chosen as the minimisation of the to reach a near optimal solution has been demonstrated.
1.002Gifz and 2.004(?i?z components in the frequency Both single %ni multiple objectives for optimisation have
response which is computed on the fly as the FDTD sim- been presented with promising results from each. It is
ulation progresses as detailed in [13]. These frequencies worth noting that if an exhaustive search were to be comwere chosen as they tie in with equipment used in a related pleted on the problems presented then over one million
measurement study. Only non-lossy simulations were run simulations would be required taking over one year to complete; the GA based search, however, took approximately
in this setup.
four hours.
The resulting frequency response from the E: component
at the susceptor point the simulations can be seen in fig- The design optimisation cases presented are rather simure 7. The dashed line in figure 7 shows the initial DFT plistic but do prove the concept of the hybrid code. More
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
149
[7] J. Jiang and G. P. Nordin, "A rigorous unidirectional
method for designing finite aperture diffractive optical elements," Optics Express, vol. 7, no. 6, pp. 237242, 2000.
[8] G. J. E. Rawlins, ed.. Foundations of Genetic Algorithms. Morgan Kaufmann, 1991.
[9] D. E. Knuth, The Art of Computer Programming,
vol. 2, Semi-numerical Algorithms. Reading, Massachusetts: Addison-Wesley, 2 ed., 1981.
10] P. Hallekalek, "Good random number generators are
(not so) easy to find," Mathematics and Computers
in Simulation, vol. 46, pp. 485-505, 1998.
11] P. L'Ecuyer, "Uniform random number generation,"
Ann. Oper. Res., vol. 53, pp. 77-120, 1994.
Frequency (Hz)
Figure 7: DFTs of E. at tlie susceptor point of the
multi-ohjective simulation with objective of minimising the
l.002GHz and 2.004GHz frequency components of the response. Shovin are the initial DFT before optimisation
begins (dashed line) and the final DFT after optimisation
(solid line).
challenging problem geometries can easily be accommodated in the FDTD code with little burden to computing
time a-s once the domain has been discretised it is the
number of re-sulting cells that chiefly determines the computation time. It is envisaged that this tool can be used to
provide design rules for component placement within an
enclosed re.TOnant cavity and that it can can also be used
to assess specific designs.
References
[1] J. Mix, G. Haussmann, M. Piket-May, and
K. Tliomas, "EMC/EMI design and analysis using
FDTD," in IEEE Int. EMC Symp., vol. 1, (Denver,
CO), pp. 177-181, 1998.
12] P. L'Ecuyer, Random Number Generation. New York:
In Handbook of Simulation, Jerry Banks (ed.) Wiley,
1997.
13] A. Taflove, Computational Electrodynamics: The
Finite-Difference Time-Domain Method. Boston,
MA: Artech House, 1995.
14] K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in
isotropic media," IEEE Transactions on Antennas
and Propagation, vol. 14, pp. 302-307, Mar. 1966.
15] B. Archambeault, O. M. Ramahi, and C. Brench,
EMI/EMC Computational Modeling Handbook.
Kluwer Academic Publishers, 1998.
16] J. B. Schneider and K. Shlager, "Finite-Difference
Time-Domain literature database." WWW. Available at http://wwu.fdtd.org Last Accessed 20"'
April 2002.
17] D. L. Carroll, "Fortran Genetic Algorithm driver."
WWW. Available at http://cuaerospace.com/
carroll/ga.html Last Accessed 20"" April 2002.
18] M. Mitchell, An Introduction to Genetic Algorithms.
MIT Press, 5 ed., 1999.
19]
[2] R. I. Macpherson and N. J. Ryan, "The use of a genetic algorithm and the FDTD method for limiting
resonance in shielded enclosures," in Proceedings of
EMC Europe &000, Brugge, A"' European Symposium
on Electromagnetic Compatibility, vol. 2, pp. 203-208,
September 2000.
[20]
[3] D. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning. Reading, MA:
Addison-Wesley, 1989.
D. L. Carroll, Genetic Algorithms and Optimizing
Chemical Oxygen-Iodine Lasers, vol. XVIIl of Developments in Theoretical and Applied Mechanics,
pp. 411-424, H. Wilson, R. Batra, C. Bert, A. Davis,
R. Scliaper, D. Stewart, and F. Swinson eds. Guntersville, AL: The University of Alabama, 1996.
N. Srinivas and K. Deb, "Multiobjective optimization
using nondominated sorting in genetic algorithms,"
Evolutionary Computation, vol. 2, no. 3, pp. 221-248, .
1994.
Russell Macpherson graduated with a
MEng (Hons) degree in electrical and
electronic engineering degree from the
University of Aberdeen, Scotland, U.K.
[5] 3. Holland, Adaption in Natural and Artificial Sysin 1999. Currently he is pursuing his
tems. The University of Michigan Press, 1975.
PhD studies at the same institution.
His research interests include genetic
[6] K. Krishnakumar, "Micro-genetic algorithms for staalgorithms, artificial neural networks
tionary and non-stationary function optimization," and fuzzy logic and how these methods can be applied
SPIE: Intelligent Control and Adaptive Systems, to electromagnetic optimisation problems.
vol. 1190, pp. 289 296,1989.
[4] R. L. Haupt and S. E. Haiipt, Practical Genetic Algorithms. John Wiley & Sons, 1998.
Macpherson and Ryan: Reducing EM Coupling in Shielded Enclosures
Nick Ryan received the BEng degree
in electrical and electronic engineering
from the University of Bristol, U.K.
in 1990 and the PhD degree from the
department of electronic and electrical engineering, University of Sheffield,
U.K. in 1998. He joined the department of engineering in the University
of Aberdeen, U.K. in 1997 where he currently lectures
safety and reliability in electronic systems. His research
interests include computational electrodynamics, electromagnetic compatibility, and renewable energy systems.
150
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
151
Impact of Some Discontinuities on tlie Convergence
Of Numerical Methods in Electromagnetics.
Malcolm M. Bibby
Gullwings, Weston, MA 02493
e-mail: mbibbv@qullwinqs.com
Abstract. The effect of discontinuilles at edges
and at feed-points of antennas on numerical
convergence rates is investigated. In the case
of edges, higher order representations of the
edge-mode, expressed with the aid of Hermite
splines, are shown to provide improved
convergence in both global and local measures.
When using a magnetic frill to excite an antenna,
it is shown that when the current representation
allows for the "charge jump" across the frill, then
convergence is accelerated. The use of both
sub-domain and entire-domain functions is
explored.
Introduction.
The treatment of edge conditions in
electromagnetic problems is generally thought to
be well understood, particularly after the
publication by Meixner [1] of the relationship
between 'edge angle' and the power of the
distance from the discontinuity. An extensive
review by Van Bladel [2] brought together much
of the evidence in support of the phenomenon's
existence. Nevertheless, the model proposed
by Meixner, while widely accepted, has not
benefited from further development.
More
recently, Shen, et al., [3][4] reported on an
analytical solution for the current and charge
distributions at the ends of a tubular dipole.
Their work showed that the Meixner method
applied only at the very ends of such a dipole.
While extending our knowledge in this area, their
work
is,
nevertheless,
restricted
to
current/charge distributions on a tubular dipole.
A different form of discontinuity arises at the
feed-point of an antenna. This fact is less well
recognized.
The author knows of only onereport that discusses it quantitatively.
In this
report, Popovic, Dragovic and Djordjevic [5, p38]
pointed out that "as a consequence of the
annular magnetic-current frill excitation, the
current derivative has a discontinuity at the frill
location".
They then went on to quantify the
charge jump that occurs in this situation.
Unfortunately, it is not clear, from the balance of
their writing, how they consistently incorporated
this observation into their work in a general
manner.
Recently, the author published reports [6][7] that
examined the use of higher order basis functions
in modeling a tubular dipole antenna/scatterer.
The results presented in these reports did not
support the expectation that higher-order basis
functions would provide faster convergence
properties.
The models used in these
investigations employed splines up to fifth order.
They also incorporated the Meixner edge
condition. Subsequent analysis has shown that
the application of the simple form of the edge
condition actually inhibited the convergence
properties of the higher order models.
The
simple form of the Meixner edge condition has
only one degree of freedom - it's amplitude. In
order for the splines covering the rest of the
antenna/scatterer to meet their potential, the
number of degrees of freedom in the edge
model were increased (in this study) to match
those of the adjacent spline.
As part of the analysis referred to in the
preceding paragraph an appreciation for the
problem associated with the charge jump across
the magnetic frill developed.
The problem
exhibits itself through a slow convergence
and/or oscillation in the_ ioput admittance of the
-antenna as the model size is increased. ' The
slow convergence problem was found to occur
for both entire-domain and sub-domain basis
functions.
The Basis Functions.
The sub-domain basis functions used on the
main portion of the antenna/scatterer are all
splines and are shown in Table I. They are the
same as those used in [6][7]. They include the
linear spline, which acts as a reference for the
performance of the higher-order functions. As
the results when using the 3-term-sin/cos
1054-4887 ©2003 ACES
Bibby: Impact of Some Discontinuities on the Convergence of Numerical Methods
152
Expansion Function
3-term sin/cos
Mathematical Description
liy) = 4 + B,sm(Jc{y - >-,)) + C,(cos(*(y - y,) - 1)
where y, < y < y,^,
Linear Spline
Ky) = «, + My -Vi) y, ^ y ^ VM
I(y) = a, + bXy - >',) + c,(y - yf yi < y < y^
Quadratic Spline
Cubic Spline
liy) = a, + bXy - y,) + c^{y - y>f + d^y - y^y
where y, < y < y^
Quartic Spline
liy) = fl, + bXy - yi) + c,iy - yf + d^y - y,y + e^y - y>y
where y, < y ^ :v,^,
Table 1. Definitions o F thie splines considered in this study.
N
N-
i=
(=1
c,A; £ A, + i.X)A- + ^1
functions were almost identical to tfiose found
when using the quadratic spline, only the latter
are reported.
The reader should be aware
however that the 3-term-sin/cos functions were
studied.
When using splines such as those
found in Table I, one asserts as much continuity
at the junction of contiguous cells as the
functions allow.
For example, a quadratic
spline, which is of order 3 and degree p = 2,
provides C' continuity. Hence, in the case of a
quadratic spline both the amplitude and the first
derivative can be matched. By this means the
number of variables is reduced. In the case of
the quadratic spline, for N cells, one needs to
determine
only
the
free
variables
N
a,,ip^Cj,i;,^,,ajy^.,.
The amplitude and first
derivative at one extremity are related to those
at the other extremity through the relationships
given in equations (1a) and (1b) above.
The sub-domain edge/end basis functions used
in this study are based on the form given in [2,
p120]:
£, = (/"(Oo + a,rf + a^d^ + a,d^ +
) (2)
where E^ is the electric field nomnal to the
edge/end, d is the distance from the edge/end,
V is a number greater than -1.0 and the series
Oj.fli.aj.aj
are constants to be determined.
Expression (2) needs to be matched up, at Its
cell boundary, with each of the functions shown
in Table I. Thus derivatives of this function are
required. Straightfonward differentiation of the
(la)
6^„ = 22cA + h
(1b)
function can be perfonned no more than once'
while still maintaining a function that can be
integrated, and thus might be employed with the
quadratic spline.
However, it would fail for
higher order splines.
An alternative, and
numerically superior, method is to apply the
concepts found in the development of Hermite
splines.
Hemiite splines are different from most splines
found in computational electromagnetic studies
which are concerned only with amplitudes. The
splines in Table I are examples of this latter
situation. The constants at the open ends of
these splines are the source of much writing on.
ways to handle them. - see, for example, {8].
The derivatives of these splines are functions of
Uiese constants. In the development of Hermite
splines the issue of derivatives is tackled
directly.
Hermite splines are typically of odd'
degree, including p = 1, and are defined in
temns of their' amplitude and (p - l)/2
derivatives at each end.
The process for
developing Hermite cubic splines is covered in
[9] and the more general case, with applications
to CEM, is described in [10]. In the case of
equation (2) only one end Is open and the
development for a two term expression will be
described next.
Two terms will provide
amplitude,
f(1),
and
derivative,
f'(1),
inforanation that can be used directly with the
quadratic spline described eariier.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
153
Equation
(2)
restated
IS
as
(A - A) < z' < A, A is the width of the cell and
/» is identified with a cell boundary.
//(< ) is a
polynomial of order 2 selected so that //,(1) = 1
and //,'(«) = 0. Similarly, H^(u) is a similar
polynomial selected so that //j(l) = 0 and
H^ '(1) = 1. When these conditions are applied,
the solution is given by:
P{u) = M"{(1 + V) - vu)f(l)
(3)
+«-(! - «)A/'(1)
This provides us with what we are looking for - a
function that separates the amplitude and the
first derivative at the open end and one that can
be integrated in a straightforward manner. This
concept can be extended to functions of higher
orders, and is done in this study for use with the
cubic and quartic splines of Table I.
When investigating the discontinuity at a
magnetic frill, a structure that possesses no
other discontinuities is useful for investigative
purposes. A circular loop is such a structure
and is examined here.
The 'natural' basis
functions to use in this situation are Fourier
series. Due to symmetry, only cosine temis are
typically used. The expansion series is then:
I = ^O;C0s(;V) 0 < <fi < 2/r
(4a)
This expansion has been in use for decades and
the results of its use are well documented [11].
When one examines the reported results, the
lack of finality in the convergence of the
admittance at the feed-point is notable. Zhou
and Smith [12] attributed this poor convergence
to the use of a delta-function source.
Their
results using
a
magnetic frill source
demonstrated that indeed the frill source
performed better than the delta-function source.
However the admittance curves still did not fully
converge until U in (4a) was extremely large.
This is surprising given the 'natural' suitability of
the functions in (4a). The work reported here
shows that the model in (4a) is deficient for use
when dealing with the magnetic frill as it does
not account for the 'charge jump' across the frill.
The derivatives of (4a) are all continuous, and
what is needed is a set of functions that allows
for this 'charge jump'.
Instead of discarding
(4a), one can modify it as in (4b) below.
I = 0,+ f,(a,cos(i^) + Z.,sin((2/ - 1)^/2)) (4b)
The derivatives of this series are such that a
'charge jump' can be modeled properly.
As
discussed in [10], in the immediate vicinity of the
frill it is the odd derivatives of the current that
must be in anti-phase on opposite sides of the
frill.
Solution Tools
The equations used in the analysis that follows
are based on the Electric Field Integral
Equation, EFIE, formulation. In particular, In the
case of the linear dipole, the equation is Hallen's
derivation [13], shown in (5) below. In (5) / is
the desired current, G is the Green's function.
R^p-zy + (2asin(^))'.
the radius is a, C
is a constant to be determined and E, is the
incident excitation.
The dipole, which is an
open circular cylinder with an infinitesimally thin
wall thickness, is located with its center at the
origin of a cylindrical coordinate system with its
length aligned with the z axis.
When investigating the loop, the derivation due
to Mel [14], specialized for a loop is used - see
(6) below. The center of the loop is located at
the origin of g cylindrical coordinate system, with
s = Rt^. The loop itself lies in the p,j* plane
with its feed-point at j^ = 0. hi this case, the
expression
for
R
in
G
is
R = a^ + IRlQ - cos(*') + 20^4(1 - cos(!>')cos^'
, where ^' is a local variable around the circular
cross-section of the loop element.
Here the
1 V,
f/(z')G(z, z')(/z'+C cos(fa) = — f £Xz')sin(A(z - r • ))£&
in '
Jli
(5)
J
IIJ{S')
1 't.
G{s,s') - kjsmkis - ^)G(^.s')(.\ - I • s')d^ ds'+ Ccosks = — UEI • i)smk(s -^)d^
2". e-;«
where G(z,z')=— | — df
(6)
Bibby: Impact of Some Disconlinuities on the Convergence of Numerical Methods
Extended Boundary Condition is used and the
net electric field on the axis of the loop is zero.
The numerical solutions for / in equations (5)
arid (6), which are expressed in the usual
ZI = V system, are obtained with the aid of the
Boundary Residual Method, BRM, [15][16].
This is a Least Squares Method variant that has
two important elements. The first significant
element is that the entries in the Z and V
matrices are calculated as they would be for a
point-matching method. However the location
of the matching points is not arbitrary.
Specifically they are located at the nodes of an
appropriate integration rule such as GaussLegendre. Furthermore, the number of sample
points exceeds the number of free variables,
typically by a ratio of 2:1 or higher. This means
the Z matrix is rectangular and the matrix
equations are solved using QR or SVD methods
[17]. The second important element is that the
rows in the Z and V matrices are weighted by
the weights from the same integration rule used
to establish the nodes mentioned earlier. This
has the important result of minimizing the
residuals integrated over the entire surface
being studied. Formally, this statement is min
pfj = \ZI - vf^.
This is an important
measure. It reports on the performance of the
solution over the complete surface and hence it
is referred to here as a global measure. For
comparative purposes, it is best to nonnalize
this function. The normalized residual is:
H
(7)
The Dipole as a Scatterer.
In order to calculate the currents on a straight
dipole excited by a plane wave equation (5) is
evaluated in conjunction with the basis functions
shown in Table I.
The boundary condition,
Ii±h) = 0, is satisfied in one of three ways. In
the first method a, and a„^^ are set to zero. In
the second approach, the two end cells support
the
basic/conventional
end
function
/(z') = a^^ih - |Z'|)/A
and then enforce
flj = fl, = a^^^. The derivatives, i, and b^^^
and higher, are left unconstrained. In the third
approach, the special Hermite end splines
defined earlier are used and all applicable
constraints are enforced. In the case of the
linear spline, the second and third approaches
154
utilize the same function.
For the quadratic
case, a two term end spline is used, for the
cubic case a three term end spline Is used and
the quarlic spline employs a four term end
spline. The plane wave is polarized parallel to
the axis of the dipole and is normally incident.
Figure 1a shows the results for the current at the
center of a half-wave dipole, with a radius of
0.007;i. As the order of the spline is increased
the results without the inclusion of any special
end sections show only a small benefit from the
inclusion of higher order temns. The inclusion of
the conventional end spline provides an
improvement, particularly for the lower order
splines. Interestingly, the convergence curves
for all four splines look very similar. It Is only
when the Hermite type end splines are added
that the real benefit of using high order splines is
observed.
Figure lb reports the results for the Normalized
Residual. It is observed that the conventional
end function improves the performance of the
lower order basis functions.
However, the
performance of the Hermite end functions is
consistently better.
The Loop Antenna.
As mentioned above, the circular loop is used to
investigate the discontinuity associated with a
magnetic frill. The basis functions defined in
equations (4a) and (4b) are used in conjunction
with equation (6) to evaluate the currents at the
feed-point of the antenna.
When using
equation (4b) the value of N is only half that of
the con-esponding value used with equation
(4a).
The circumference of the loop is one
wavelength and
n = 2]ii(2ftRJ a) = lO.O,
where R,, Is the radius of the loop and a is the
Tadius otthe loop element.
The convergence curves for the central current
are shown in Figure 2a and the associated
normalized residuals are shown in Figure 2b.
The results show clearly the superiority of a
series that accommodates the 'charge jump' at
the feed-point, thereby confimning the presence
of the phenomenon.
Only 32 temns are
considered in these graphs. If the number of
terms were increased substantially, the results
for the cosine series only would eventually
converge to those observed with the use of the
supplemented basis functions.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
155
3 4 56710''
3 456710''
2 3 456
2
3 456
_L
O without end cell
«i with conventional end cell
n with Hermite end cell
3 4 567.]Q1
2
3 456
3 4567 .|Q1
2
3 456
# of cells on a half-wavelength dipoie antenna
Figure la. Plots of the current oh a dipoie excited by a plane wave modeled with four splines
3 456710''
2 3 456
3 456r10''
2 3 4 56
O witout end cell
A with single end cell
D With Hermite end cell
3 4567 .JQ1
2
3 456
3 4567 .|Q1
2
3 456
#of cells on dipoie scatterer
Figure 1b. Plots of the normalized residual on a half-wave dipoie excited by a plane wave.
Bibby: Impact of Some Discontinuities on the Convergence of Numerical Methods
156
10.0
9.0
8.0
to
|7.0H
E
^6.0
5.0
O Cosine series only
£t Cosine series + single sine term
n Cosine series + sine series
4.0
3.0
T
2
3
456789 .^Q1
2
3
# of temis in a Fourier series expansion on a one wavelength circular loop antenna
Figure 2a. Plots of the input admittance of a circular loop excited by a magnetic frill.
-1.0
O
-3.0 m
o ^
O - O -
D
°-°-^^--o.
Wen
g> -5.0
D .
TJ
A.
D
N
1-7.0
D.
oooc
'n. ^"--A AA^A,
^n.,
o
-11.0
O Cosine series only
A Cosine series + single sine term
n Cosine series + sine series
^flflxmnWj^
-13.0
6
7
8
9
10^
# of terms in a Fourier series expansion on a one wavelength circular loop antenna
Figure 2b. Plots of the normalized residual for a circular loop eXcited by a magnetic frill.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
157
/ = tjj + '^a^cos(izz'/h)^h - \z\
{8a)
I = a,+ J^{a,cos(i!rzph) + 6,sin((2( - l);r|2'|/2/0yA - |z'|
(8b)
The Dipole as an Antenna.
Discussion.
The dipole antenna incorporates botii of the
types of discontinuity already discussed. Here,
it has radius, a, of 0.007/1. The frill has
dimension of 6/a=2.3 where b is the outer radius
of the frill.
The combination of these two
discontinuities is explored initially using subdomain functions.
The sub-domain functions
are the same as those used earlier with a dipole
scatterer. However, at the center of the dipole,
the phase in the tenns representing the first and
third derivatives is shifted 180° when these
derivatives exist.
The impact on the
convergence of the current at the feed-point can
be seen in the results reported in Figure 3a.
The results for the associated normalized
residuals are shown in Figure 3b. It is clear,
from both sets of data, that the phase reversal of
the odd order derivatives at the frill contributes
significantly to an improved model for the
current.
End Effects.
When considering the results
reported in Figure 1a, one must bear in mind
that, for the linear spline, the conventional end
temi and the Hennite end term are one and the
same. Adding a single end term to any of the
splines produces an observable improvement in
the convergence of the current at the center of
the dipole compared with not using one at all.
This is particularly true for the lower order
splines. As mentioned earlier, the convergence
curves obtained when using a single
conventional end term all appear to have the
same convergence characteristics.
This
suggests that the end term is limiting the
convergence process. When the higher order
end tenns are included, to the maximum order
possible, the convergence process is speeded
up significantly.
When viewed from a global
perspective, as reported in Figure 1b, the
conventional end term is important in the case of
the linear spline. For the higher order splines
the benefits become less as the order of the
spline is increased. However, when the higher
order end ternis are included, the global results
are distinctly better for all splines.
Entire-domain functions have long been used to
model the current on a linear dipole antenna. A
study that compared ten such functions was
published recently [18]. That study concluded
that entire-domain functions performed poorly on
the linear dipole.
Historically, however, the
choice for an entire-domain function on a linear
dipole has been one of various forms of a
Fourier series.
In the absence of a better
alternative, such a series is used here. The
particular form of the Fourier series used is
shown in (8a) above. In order to accommodate
the "charge jump" at the magnetic frill, a
modified form of the latter was investigated this is shown in (8b) above.
The results, for the input admittance of a halfwavelength antenna, using these two series are
shown in Figure 4a.
The corresponding
findings for the normalized residuals are shown
in Figure 4b. When using equation (8b) the
value of ;\^ is only half that of the corresponding
value used with equation {8a).
Both figures
clearly demonstrate the utility of incorporating
terms that address the issue of discontinuity
across the frill.
Magnetic Frill. The importance of properly
incorporating "charge jump" properties into
models of current flowing on an antenna excited
by a magnetic frill is dramatically illustrated in
Figures 2-4. The significant oscillations seen
when employing a Fourier series containing only
cosine terms is an IndicatLon of its lack of.
suitahility - for both the l<56p and the linear
dipole! To be sure that the oscillations were the
result of the feed mechanism and nothing else,
the excitation was changed from the magnetic
frill to a plane wave.
The oscillations
disappeared in the case of the loop and were
highly attenuated on the linear dipole.
This
allowed the conclusion that the oscillations are
indeed a result of using the magnetic frill
excitation and not something else. The addition
of the sine terms in the current model, while
improving the convergence results significantly,
does raise the condition number of the
impedance matrix from the region of ~ (0* up to
Bibby: Impact of Some Discontinuities on the Convergence of Numerical Methods
158
567810''
9.8
quadratic
2
3
456
cubic
quartic
O without phase reversal of odd derivatives
£> with phase reversal of the odd derivatives
9.6 -
# of cells on a half-wavelength dipole antenna
Figure 3a. Plots of admittance on a dipole excited by a magnetic frill for three splines.
34567 10''
2
3456
3 4 567 .^Q1
2
3 4 5 6
#of cells on a half-wavelength dipole antenna
Figure 3b. Plots of the normalized residual on a dipole with magnetic frill for three different splines.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
159
O Cosine series only
O Cosine + sine series
D Cosine series + abs(z')
9.7 -
I 9.4
E
E
-O-O-O-D-
^9.1 -^
8.8 -
:
8.5
1
,
\
\
>
\
<
\
<
\
<
I
r-
4
8
12
16
20
24
28
# of terms in a Fourier series expansion on ah half-wave dipole antenna
32
Figure 4a. Plots of the admittance of a half-wave dipole excited by a magnetic frill
-2.0
O Cosine series only
O Cosine + sine series
D Cosine series + abs(2')
-4.0
(0
T3
S-6.0
a
N
2 -8.0 -
MO.O
o
0>(/X>^
-12.0
--—-^.^^
-14.0
T
0
4
8
12
16
20
24
28
# of terms in a Fourier series expansion on a half-wave dipole antenna
32
Figure 4b. Plots of the normalized residual on a half-wave dipole excited by a magnetic frill
Bibby: Impact of Some Discontinuities on the Convergence of Numerical Methods
- 10'° when more than 12 terms are used in the
series.
The linear dipole problem was also solved using
a Method-of-Moments solution, using the
standard Pocklington equation. This was done
as a check on the testing method and the
integral equation formulation.
The testing
functions were pulses of equal dimensions, the
basis functions were again the Fourier series of
(8a) and (8b) and the Pocklington equation was
regularized per the procedure described in [19].
The results of using (8b) converged significantly
faster than when using (8a) while the oscillations
noted in connection with Figure 4 were almost
completely attenuated. These findings provided
further confirmation of the need to accommodate
the discontinuity at the feed-point.
The presence of a high condition number in the
entire domain solutions that include a series
expansion of the absolute value of sine temis
makes them unattractive.
Consequently, the
effect of incorporating a single term only was
Investigated. In the case of the loop this term
was sin([(*'|/2). In the case of the dipole it was
simply |z'|.
The associated results are shown
in Figures 2 and 4 respectively.
Use of just
these single terms contributes significantly to
alleviating the convergence problem.
The
condition numbers were comparable to those
observed when using the cosine temis only.
This suggests that one judiciously selected term
could significantly speed up the convergence
rate in these situations.
However, the ones
tested here are not those functions - othenwise
they would exhibit normalized residual curves at
least as steep as those shown by the addition of
the sine series.
The basic EFIE for a loop is shown in equation (9a).
After some manipulation, this assumes
the fomn shown in (9b). When the current is
represented by (4a) the fonn In (9c) results.
The last term on the right-hand side of this latter
jafi
160
equation evaluates to zero. But this is the term
that represents the "charge jump". Therefore
the representation of (4a) is inadequate.
A
temi. or terms, that permit(s) the current
derivative at the origin to be discontinuous is
necessary.
The inadequacy of the
representations of (4a) and (8a), is believed to
be the main source of the oscillations seen in
the relevant convergence curves.
In 1980, Richmond [20] reported that the
incorporation of a term that represented an
outgoing wave on an infinite dipole improved
convergence considerably on a finite dipole.
This was confirmed in [18], where it was also
noted that the condition numbers were high.
This outgoing term is essentially a log term at
the feed-point and consequently presents a
discontinuity there.
It is possible that this
discontinuity had more to do with the
improvement realized by Richmond than the
fonm of the temi itself.
Conclusions.
This study was motivated by the apparent failure
of high-order basis functions to improve
conviergence rates when compared with basis
functions of lower order [6][7]. As a result of
this work, it is believed that the reasons for this
failure originate with inadequate current
modeling in the vicinity of discontinuities.
Specifically, it has been demonstrated that to
achieve superior convergence properties one
must bear in mind three things.
1.
Current models in the vicinity of
discontinuities must support the same level of
continuity as that supported by the models
employed away from the discontinuities.
2.
When an Isolated excitation source,
such as a magnetic frill, is part of the system,
the current model(s) must have "sufficient
flexibility to accommodate issues such as
"charge jump" in the excitation region.
3.
The above comments apply to both subdomain and entire-domain models.
= Jcos(<# - ^')I(^')Gi^J')Rdr + -J-5- ]^^Rd^'
d'l \Gd^'+ ^G
dl
J ^ j ]][cos((!> - (^')k^R^I(^') + -^
1*1 cos{(p')k'R'a„ + ^(cosmeR' - f)ajCosU</>') Gdf+ J^-JajSmU^yj
I
1—
-EM =
(9a)
K
(9b)
N
M
{9c)
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
161
Acknowledgements.
The author thanks Professor Andrew Peterson
of the Georgia Institute of Technology for his
encouragement and useful suggestions in
connection with the preparation of this paper.
References.
[I]
L.IVleixner,
"The
Behavior
of
Electromagnetic Fields at Edges", IEEE
Trans, on Ant. & Prop., Vol. 20, No. 4,
pp.442-446,1972.
[2]
J. Van Bladel, Singular Electrdmagnetic
Fields and Sources, Chapter 4, Clarendon
Press, 1991.
[3]
H. Shen and T. T. Wu, "The Universal
Current Distribution Near the End of a
Tubular Antenna", J. Math. Phys., 30(11),
Nov., 1989.
[4]
H. Shen, R. W. P. King and T. T. Wu, "The
Combination of the Universal End-Current
and the Three-Term Current on a Tubular
Dipole", J. of Electromagnetic Waves and
Applications, Vol. 4, No. 5, 189-200,1990.
[5]
B. D. Popovic, M. B. Dragovic and A. R.
Djordjevic, Analysis and Synthesis of Wire
Antennas, Research Studies Press, 1982.
[6]
M. M. Bibby, "A Comparative Study of
Expansion Functions Using the Boundary
Residual Method on a Linear Dipole - Part
II: Sub-Domain Functions", Journal of the
Applied Computational Electromagnetic
Society, Vol. 17, No. 1, pp. 54-62, March
2002.
[7]
M. M. Bibby, "An adaptive Approach to
Faster Convergence in the Solution of
Electromagnetic Problems by Using nonUniform Cell Sizes", Proc. 1^. Annual
Review
of
Progress
in
Applied
Computational Electromagnetics, pp. 457463, Monterey, CA., March, 2002.
[8]
C. F. Gerald, Applied Numerical Analysis,
2"". Ed., pp. 477-482, Addison-Wesley
Publishing Co., 1978.
[9]
J. F. Epperson, An Introduction to
Numerical Methods and Analysis, pp. 171181,Jo hn Wiley & Sons, 2002.
[10] M. M. Bibby, "Hermite Splines and Their
Application in Electromagnetic Problems",
in preparation.
[II] R. W. P. King, & G. S. Smith, Antennas in
Matter, pp. 529-537, The MIT Press, 1981.
[12] G. Zhou & G. S. Smith, "An Accurate
Theoretical Model for the Thin-Wire
Circular Half-Loop Antenna", IEEE Trans,
on Ant. & Prop., Vol. 39, No. 8, pp11671177,1991.
[13] E. Hallen, "Theoretical Investigations into
the Transmitting and Receiving Qualities
of Antennae", Nova Acta Regiae Soc. Sci.
Upsaliensis, Ser. IV, No. 4, pp. 1-44,
1938.
[14] K. K. Mei, "On the Integral Equations of
Thin Wire Antennas", IEEE Trans, on Ant.
& Prop., Vol. 13. pp. 374-378, May 1965.
[15] K. J. Bunch & R. W. Grow,, Numeric
Aspects of the Boundary Residual
Method", Int. J. of Numerical Modeling:
Electronic Networks, Vol. 3, pp.57-71,
1990.
[16] K. J. Bunch & R. W. Grow, "On the
Convergence of the Method of Moments,
the Boundary Residual Method and PointMatching Method with a Rigorously
Convergent Formulation of the PointMatching Method", Journal of the Applied
Computational Electromagnetic Society,
Vol. 8, No. 2, pp. 188-202,1993.
[17] C. L. Lawson & R.-J. Lawson, Solving
Least-Squares Problems, Prentice-Hall,
Inc.. 1974.
[18] M. M. Bibby, "A Comparative Study of
Expansion Functions Using the Boundary
Residual Method on a Linear Dipole - Part
I: Entire-Domain Functions", Journal of the
Applied Computational Electromagnetic
Society, Vol. 17, No. 1, pp. 42-53, March
2002.
[19] D. B. Miron, "The Singular Integral
Problem in Surfaces". IEEE Trans, on Ant.
& Prop., Vol. Ap-32, No. 3. May 1983.
[20] J. H .Richmond, "On the Edge Mode in the
Theory of Thick Cylindrical Monopole
Antennas", IEEE Trans, on Ant. & Prop.,
Vol. 28, No 6, pp. 916-920, Nov. 1980.
Malcolm M. Bibby received the B.Eng. and
Ph.D. degrees in Electrical Engineerinjg
from the University of Liverpool, England in
1962 and 1965 respectively. He has been
interested in the
numerical aspects
associated with antenna design for the last
twenty years.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
162
Wireless Propagation in Non Line of Sight Urban Areas using Uniform Tlieory of
Diffraction
Grant A. Ellis
Department of Electrical and Computer Engineering
National University of Singapore
Abstract: This paper describes a three-dimensional
electromagnetic propagation model for signal
power prediction in a non line of sight urban area
located in Singapore.
The model, which is
implemented using the Ohio State NEC-BSC V42
basic scattering code and is based on ray theory
and uniform theory of diffraction (UTD) takes into
account first-order, second order, and third order
effects including triple reflections and diffractions.
The simplified propagation model of Raffles Place,
the central business district of Singapore, uses
more than 360 geometrical structures and is
compared with data measured at 937.6 MHz along
a drive route at the site. The results from the
propagation model produced reasonable agreement
with the measurements showing that a simplified
model using UTD can be used to simulate the gross
features of electromagnetic scatter in an urban area.
It is shown that use of penetrable dielectric
building walls are also necessary for accurate
prediction of radio wave propagation in urban areas
under non line of sight conditions. It is also found
that third order scattering effects can be dominant
in non line of sight situations and may be necessary
for accurate predictions.
1
B<TRODUCnON
As cellular mobile communications systems
become more widely used, there is an increasing
need to develop reliable planning tools for base
station deployment. Part of this growth is due to
frequency reuse whereby the same frequencies_are
utilized in geographically distinct areas or c3ts
such that signal path loss over distance can be
-exploited.
Frequency reuse provides for
significantly greater capacity.
One of the
drawbacks with frequency reuse however, is the
potential for interference between nearby cells. As
cells become smaller in size, the use of reliable
planning tools becomes even more necessary. In
these environments, accurate signal strength
prediction plays an important role in controlling
intercell interference and providing coverage for
low-power handheld units. The use of reliable
planning tools means that the use of expensive and
time-consuming measurements in the field can be
minimized. These tools can also provide a cost
effective means of anticipating problems at the
early stages of base station deployment and of
optimizing parameters such as antenna location and
orientation.
Several techniques have been developed [1-6]
for analysis of multipath and other propagation
effects in macrocellular and microcellular
environments. In several of these approaches, ray
tracing or launching techniques are combined wifli
GTDAJTD (Geometrical Theory of Diffraction
AJniform Theory of Diffraction) to analytically
determine the sum of all incident, reflected, and
diffi-acted fields at the receiver location. The use
of GTD/UTD is appropriate, since the frequencies
used in mobile communications (e.g. > 900 MHz)
are usually high enough such that the assumption
that the features of most building structures are
much larger than the wavelength used becomes
valid. In [2], a computer based propagation model
is used to obtain a time delay representation of the
received signal. In [3, 4], a multislit waveguide
with randomly distributed gaps between the sides
of buildings is considered as a model of straight
streets under line of sight conditions. In [5], a ray
launching technique is used to detect the
intersection of a ray with an object and the
diffi'acted and reflected rays are then computed.
This process is repeated until a maximum number
of reflections is exceeded. In [6], 2D FDTD is
corrected to account for spherical wave spreading
and applied to a simple outdoor environment. Each
building in the model consists of a perfectly
. conducting core with a dielectric coating.
In this paper, a particular urban area was
selected, a three-dimensional model was created
and the propagation characteristics were simulated
using NEC-BSC V4.2 (Numerical Electromagnetic
Code-Basic Scattering Code) [7].
The model
presented includes contributions to the received
signal from all possible propagation paths,
including ground and wall reflections from
diffracted and specularly reflected signals both in
the LOS (Line of Sight) and NLOS (Non Line of
Sight) regions using GTD/UTD. This code is
structured such that all of one type of scattered
1054-4887 © 2003 ACES
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
163
field or interaction is computed at one time to
increase efficiency. Buildings having an arbitraiy
layout can be modeled even if the streets are not
rectilinear and the buildings walls on each side of
the streets are not coplanar.
The purpose of this research is to investigate
how accurately UTD could predict electromagnetic
scattering in a non line of sight urban environment
using a simplified model consisting of plates and
curved surfaces.
In order to evaluate the
limitations of the model, predictions from the
model considering f, t^, and 3"* order scattering
mechanisms were compared to CW test data
measurements from a private cellular provider.
This paper is organized as follows: Section 2
describes the urban propagation model, Section 3
compares the measured vs. predicted results, and
Section 4 summarizes the results and provides
some directions for future research.
2
URBAN PROPAGATION MODEL
The goal of this investigation was to determine
how well UTD could be used to predict
electromagnetic scattering in a heavily urbanized
environment using a simplified model of the site.
Raffles Place, the central business district of
Singapore was chosen as an appropriate site for
this study due to its high building density. Other
reasons for choosing Raffles Place were:
1.
2.
3.
The building layout does not lie completely in
a rectangular grid. This allows for
investigations to be carried out in an area
without a rectangular grid layout and having a
somewhat random building orientation.
Drive routes can be investigated which arc
mainly non-line-of-sight (NLOS).
The building heights in this area range from 5
to 282 meters. The effects of diffraction from
the edges of rooftops and sides of buildings
can be investigated.
~
Figure 1 shows a three-dimensional view of
the geographical layout of Raffles Place used in the
UTD model with the drive route shown. The drive
route begins along Cecil Street and follows onto
Collyer Quay up to Fullerton Square then turns left
onto Battery Road where it ends on Chulia Street.
The area covered is about 500 m by 500 m and the
length of the drive route is 1.12 km. The data
predicted along this route is compared with the
measured drive data collected by MobileOnc, a
private cellular network provider in Singapore. In
all field measurements, GPS was used to record the
location of the measiired signal power along the
route. The following input model parameters are
required in the NEC-BSC UTD code [7]:
1.
2.
3.
4.
5.
6.
Coordinates of the plates, curved surfaces, and
other structures necessary to simulate building
exterior and roofs.
Estimate of the electrical permittivity,
conductivity, and thickness of the building
exteriors.
Transmitter and receiver locations.
Transmitter and receiver antenna patterns.
Operating frequency used.
Types and order of UTD interactions to be
included in the simulation.
The coordinates of the building exteriors were
extracted from an aerial map by using an overlay
having a fixed origin selected on the map. The
comers of the buildings were measured relative to
this origin. Most building shapes were modeled
using only plates, which were assumed to be flat
surfaces having average electrical permittivity and
conductivity. Structures having low curvature
sections were modeled using multiple plates to
form a piecewise linear approximation of that
section. When the radius of curvature of the
structure was substantial, cylinders or truncated
cones were used. Elliptic cylinders were also used
to model some building structures.
The ground plane was simulated using a single
homogeneous half space of a relative dielectric
pcmiittivity £r = 10 and conductivity CT = O.OOI
S/m, which can be considered reasonable for a
typical ground. Preliminary studies suggest that
the predicted signal power near the ground is
somewhat insensitive to the ground dielectric
constant and conductivity. The seawater surface
was also included in the model and was modeled as
a dielectric coated plate with a relative permittivity
of Er= 80 and-electrical conductivity o = 4 S/m
The model database used in the UTD
propagation model consists of 360 plates and 5
cylinders to represent the 64 buildings and other
structures located at Raffles place. The number of
plates in the model could easily exceed 500 if the
entire area is to be accurately modeled. However,
following the assumption that effects from
structures like trees, lamp posts, telephone boxes.
Mass Rapid Transit station buildings, etc can be
ignored, these structure were not included in the
model [8], Every attempt was made to define the
model in more detail in areas where it was
illuminated strongest by the antenna. Because of
Ellis: Wireless Propagation in Non Line of Sight Urban Areas
the wide variation in building heights, the detailed
description in the model is not limited to the area
local to the transmitter. Since the model used
omnidirectional antennas for the transmitter and
receiver, the entire area including but not limited to
the region between the transmit antenna and the
receiver field point was modeled. As a result of the
above, the number of plates required for modeling
was quite large to account for all the possible ray
paths.
The building exteriors were modeled as
dielectric plates to approximate materials, such as
glass windows, absorbing panels, concrete, etc.
Research is still underway to produce useful
approximations for the thickness, permittivities and
conductivities of these delectric layers [8]. The
choice of a good approximation for the electrical
properties is of utmost importance, since the
Fresnel reflection coefficients depend on these
values. These coefficients dictate the reflection
coefficient of the materials influencing the
modeled results.
The buildings located at Raffles Place were
constructed using a variety of materials, such as
concrete, glass, and steel. It is difficult to identify
the electromagnetic properties of each individual
building exterior and it is assumed that all the
building walls have electrical property values
averaged over those of glass and concrete. The
building walls were modeled using plates and
cylinders coated with dielectric on one side to
simulate the building exterior and a perfect electric
conductor on the other. The material properties
were chosen as Ei- = 15, a = 7 S/m, and thickness =
0.5 meter [9, 10]. The accuracy of the model is
limited by the lack of accurate building and street
databases and the approximation of building
exteriors to be "smooth" flat surfaces having an
average relative permittivity and conductivity.
- Both the transmitter and receiver antennas
were simulated using vertically polarized halfwave dipoles. The transmitter input power was
32.4 dBm. The transmitter dipole antenna was
placed 3 m above the machine room on the rooftop
of Ocean Building, a building of about 120 m in
height. The receiver was placed 2 m above the
ground to simulate mounting of a vertical
monopole on top of a car. The measured data was
taken at intervals of between 28 - 59 m usmg a
GPS based measurement system. The simulated
data is sampled in 20 m intervals with and without
spatial averaging. AH UTD scattering mechanisms
were selected during the simulation. The NEC-
164
BSC UTD code however, allows for the individual
scattering terms to be selected to investigate their
behavior along the route.
3
COMPARISON OF THE MEASURED
VS. PREDICTED RESULTS
Three sets of simulation cases were
considered. The first simulation case considers
only the l" and 2"'' order interaction terms listed in
Table 1. The second simulation case considers up
to and including the P order interactions. The
additional S^ order interaction terms are listed in
Table 2. For these first two cases, the building
walls only reflect and diffract the radio waves
using the material dielectric properties given. The
third simulation case considers f, T^, and 3^
order interactions listed in Tables 1 and 2 using
penetrable dielectric plates as building walls. For
this case, radio waves are both reflected and
transmitted through the plates to simulate
propagation through building structures.
Figure 2 shows a comparison between the
predicted signal powers resulting from the first two
simulation cases (1*' and 2°'' order interactions
only, and l", ?'' and 3"* order interactions) and the
measured fields. The f' and T^ order interaction
terms listed in Table 1 result in a total of 5015 ray
paths along the drive route. Table 2 lists the Sf"
order interaction terms resulting in an additional
42850 ray paths. The curve showing all 1st, 2!"*
and i order interactions appears to give better
agreement with the measured field data.
The greatest discrepancy between the
measured and predicted data for cases 1 and 2
occurs over the 0.2 km stretch located immediately
before point D near Fullerton Square (Figure 1).
The predicted results for the case of up to 3"* order
interactions underestimate the measured data by
more than 20 dB. "This is difficuU to explain using
the model because there are no nearby structures
located in this area on which the rays can get
scattered. Along the drive route to the right of
Collyer Quay is an open area immediately adjacent
to Marina Bay. One possible explanation for the
discrepancy between the neasured and predicted
data is the lack of radio wave propagation through
buildings since modem office-building interiors are
composed of mostly air. Studies have shown that
radio waves may penetrate into buildings with only
between 9-11 dB attenuation at 900 MHz [11].
Our model vvas modified to allow the building
walls to be penetrable by radio waves such that the
total attenuation through the building structures is
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
165
between 12-15 dB. The properties of the building
walls were changed to Er = 15, o= 0.0023 S/m, and
thickness = 1 m and were simulated as penetrable
dielectric plates.
Shown in Figure 2 is the
measured vs. predicted signal power along the
drive route in Figure 1 using penetrable walls. The
correspondence between measured and predicted
power is much better than for the cases where the
building walls only provide reflection and
diffraction.
It can be observed that use of
penetrable building walls results in 11 dB less
variation in the predicted field along the drive route
than for the case of 3"* order scatter only. Table 3
lists the RMS error and standard deviation for each
simulation case.
Additional observations can be made from the
measured vs. predicted signal power results shown
in Figure 2. At location A on Cecil Street, the peak
in the measured received field may be attributed to
3'^ order interactions since this peak cannot be seen
in the 2"'' order results. For instance, in the
prediction, a 3"" order reflection occurs between the
Singapore Land Tower building, the OUB Center,
and a second reflection from the Singapore Land
Tower (Figure 3). Multiple interactions between
the Ocean Building and the Ocean Towers also
occur before the ray arrives at the receiver location.
At location B, the measured signal power is at
the lowest level along the drive route and the
number of 2"' and 3'^ order ray paths appears to be
minimum. At location B the only contribution to
the signal power at the receiver location appears to
be from the diffraction from the rooftop of the
Ocean building. This location is also the point
closest to the transmitter along the drive route.
At location C, the predicted results that include
up to 2"'' order interactions only indicate a drop of
nearly 50 dB in the received signal strength due to
obscuration by an overpass located along the route.
For the predicted results that include 3^ order
interactions, the effect of the overpass is minimal.
The measured signal power also indicates
essentially no effect due to the overpass. The
effects of the 3"* order interactions appear to be
dominant directly under and immediately beyond
the overpass. Multiple rays that originate at the
transmitter location may be diffracted from the
rooftop edge on top of the Ocean Building and then
diffract and reflect fom the north side of the John
Hancock building. These rays are then reflected
from the sides of both the Union Overseas
Shopping Center and a restaurant and then finally
arc reflected from the side of the Hitachi Center
(see Figure 4) to the receiver location immediately
beyond the overpass location.
Probably the most interesting observation in
the measured and predicted data is shown at
location D. Here, the measured and predicted
signal powers are at their maximum value over the
entire route.
Location D is in front of the
Singapore Land Tower on Battery Road. The
received signal power is about 12 dB above the
average signal power along the entire route. At
this location, the transmit antenna on the Ocean
Building is obstructed by several buildings. The
received signal appears to come from several ray
paths that first begin as reflections from the OUB
Center and UOB Plaza and then are reflected from
the Singapore Land Tower and Standard Chartered
buildings to the receiver location. This area of
Battery Road appears to behave like a parallel-plate
waveguide so that the rays from several locations
are guided into a relatively narrow street bounded
by the Singapore Land Tower on one side and
Standard Chartered Bank and May Bank on the
other side (see Figure 5).
At the end of the route and furthest from the
transmitter on Battery street (location E) is the
location of the largest discrepancy between the 2"
and 3"* order predicted data. The measured signal
power in the area around location E however, is
nearly constant. The predicted data however,
shows almost a 50 dB variation for the l" and 2"''
order interactions and nearly 27 dB of variation for
the case of up to 3"* order interactions. The
predicted data shows a 2"'' order interaction that is
the Esult of a specular reflection from the City
House building and then diffraction between the
edges of the Exchange Building and Golden Shoe
Car Park (Figure 6) to the receiver location. The
measured power increase is only about 1 dB at this
location. Rays reflected from the Republic Plaza
and the UOB Plaza channel into this area resulting
in a peak in the received signal power. Like the
results shown at location D, this area may behave
like a parallel-plate waveguide where multiple
reflections occur between UOB Plaza on one side
and the Sinsov building on the other side across
Battery Street.
4
CONCLUSION
In this paper, a simplified model using UTD
with ray tracing and consisting of 360 plates and 5
cylinders was used to predict high-frequency
electromagnetic wave propagation in an urban
environment. Taking into account single, double,
Ellis: Wireless Propagation in Non Line of Sight Urban Areas
and third order UTD terms interacting between
building walls and comers, the gross features of the
electromagnetic
response
were
accurately
predicted.- The three-dimensional model was
developed using arbitrary building layouts modeled
using simple structures constructed using plates to
represent flat surfaces and cylinders for curved
surfaces.
Several sets of UTD simulations were
performed and the resulting predicted signal
powers were compared with measurements along
the drive route located in Raffles Place, the central
business district of Singapore.
Comparisons
between the predictions and measurements suggest
that:
1.
When there is no line of sight (NLOS),
diffraction cannot be neglected. The effects of
diffraction from the edges of the rooftop
appeared to be significant at the transmitter
location, f and 2!"* order interactions alone
may not be sufficient for accurate predictions.
Inclusion of higher order scattering is
important for the accurate prediction of the
scattering in NLOS situations. Inclusion of 3"*
order interactions also result in 9 dB less
variation than 2'^ order interactions in the
predicted signal power along the drive route.
For NLOS cases, the fields at the receiver
location often result from scattering on
structures local to the receiver location.
These fields may also result from reflections
from large structures or surfaces far away
from the receiver location. This can result in
propagation paths that may not be intuitively
obvious. These higher order scattering terms
(3"* order and above) may require a larger
model to account for all possible ray paths.
Lack of sufficient structures in the
propagation model and inaccuracies in
antenna positions, route locations, may
introduce errors into the predicted signal
powers. Another source of potential error is
scattering from cars, trucks, or other vehicles
located near the receiver and not included in
the model. It was also assumed that the
dielectric constant and conductivity of the
building exteriors throughout the model was
constant. These simplifications introduce
extra error.
4.
UTD can also produce artifacts in the
predicted fields.
At some locations
166
considering up to f* order interactions only
result in artifacts not observed in the measured
data. At those locations, the 3"* order scatter is
dominant over the f and 2"* order scatter
effects.
5.
At two locations along the drive route, it was
found that the peak in the Bceived signal
power might be the result of a parallel-plate
waveguide effect caused by multiple
reflections between adjacent buildings
separated by a street or road.
6.
Effects of radio wave propagation through
building walls and structures are necessary in
resolving differences between measured and
predicted data under NLOS conditions. Use of
penetrable buildmg walls in the model reduces
the RMS error in the prediction by 13 dB.
Future work may include incorporation of
higher order (e.g. 4"' order) interaction terms in the
propagation study. This may be usefiil for the
prediction of the signal power near Fullerton
Square. The dielectric properties of the building
exteriors are also very important and may have a
large effect on the magnitudes of reflections from
buildings.
Future work includes additional
modeling of propagation through buildings.
Analysis of these effects will be helptiil in
resolving differences between measured and
predicted signal power data. Also of interest are
the cross-polarized fields scattered in non line of
sight locations, for the study of polarization
diversity antenna schemes.
ACKNOWLEDGMENTS
The author wishes to thank the following for their
contributions:
1. MobileOne for providing the CW_drive test
~' measurement data.
2. BuHding and Construction Authority of
Singapore for providing a detailed map of
Raffles Place.
3. To my student, Cheong Si Neng for the
Raffles Place model.
4. Ronald J. Marhefka at The Ohio State
University Electroscience Laboratory for
strong support and help in using the NECBSC
code.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
167
REFERENCES
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
F. Ikegami, T. Takeuchi, S. Yoshida,
"Theoretical prediction of mean field
strength for urban mobile radio", IEEE
Trans, on Ant. and Prop. Vol. 39, No. 3,
pp. 299-302, March 1991.
J. P. Rossi, J.C. Bic, A. J. Levy, Y.
Gabillet, M. Rosen, " A ray launching
method for radio-mobile propagation in
urban area", IEEE Am. and Prop. Symp.,
Vol. 3, pp. 1540-1543, June 1991.
N. Blaunstein, M. Levin, "Propagation
Loss in the Urban Environment with
Rectangular Grid-Plan Streets", lO"'
International Conference on Antennas
and Propagation, pp. 2.178-2.181, 14-17
April 1997.
N. Blaunstein, R. Giladi, and M. Levin,
"Characteristics' Prediction in Urban and
Suburban Environments", IEEE Trans.
Veil. Tech., \b\. 47, No. 1., pp. 225-234,
Feb. 1998.
S. Y. Tan, H. S. Tan, " Propagation
model for microcellular communications
applied to path loss measurements in
Ottawa city streets", IEEE Trans, on Veil.
Tech.. Vol. 44, No.2, pp. 313-317, May
1995.
J. Schu.ster and R. Luebbers, "Comparison
of GTD and FDTD Predictions for UHF
Radio Wave Propagation in a Simple
Outdoor Urban Environment," IEEE APS
Int. Symp. Dig., Montreal, Canada, Vol. 3.
pp. 2022-2025, July 1997.
R. J. Marhefka, J.W. Silvestro, "Nearzone basic scattering code. User's manual
with space station applications" Technical
Report 716199-13, March 1989.
M. F. Catedra, J. Perez, F. Sacz de Adana,
and O. Gutierrez, "Efficient ray-tracing
techniques for three-dimensional analyses
" of propagation in mobile ommunications:
application to picocell and microcell
scenarios", IEEE Antennas and Prop.
Magazine, Vol. 40, No. 2, pp. 15-28,
April 1998.
[9]
[10]
[11]
S-C Kim, B. J. Guarino Jr., T. M. Willis
111, V. Erceg, S. J. Fortune, R. A.
Valenzuela, L. W. Thomas, J. Ling, and J.
D.
Moore,
"Radio
propagation
measurements and prediction using threedimensional ray tracing in urban
environments at 908MHz and 1.9GHz",
IEEE Trans, on Veh. Tech.. Vol. 48, No.
3, pp. 931-946, May 1999.
O. Landron, M. J. Feuerstein, and T. S.
Rappaport, "In Situ Microwave Reflection
Coefficient Measurements for Smooth and
Rough Exterior Wall Surfaces", in Proc.
43"" IEEE Veh. Technol Conf., pp. 77- 80,
1993.
A. M. D. Turmani, A. F. de Toledo,
"Modeling of radio transmissions into and
within multistory buildings at 900, 1800,
and 2300 MHz", lEE Proc. I, 140, (6), pp.
462-470,1993.
Grant A. Ellis was bom in
Portland, OR in 1960. He
received the B.S. and
M.S. degrees in electrical
engineering
from
Washington
State
University, Pullman, WA,
in
1583
and
1986,
respectively, and the
Ph.D. degree in electrical engineering from the
University of Washington, Seattle, in 1995. From
1983 to 1985 he was a staff engineer with the
Boeing Aerospace Company, and from 1987 to
1998 he was a principal engineer with the Boeing
Defense and Space Group. He was an assistant
professor of electrical engineering with the
National University of Singapore from 1998 to
2001. He is presently a member of the research
staff of HRL Laboratories in Malibu. His current
research interests are in microwave circuit design
including monolithic MJC design, active antenna
design, and propagation analysis.. Dr. Ellis is a member of Tau Beta Pi, Eta Kappa
Nu, ACES, and the IEEE.
168
Ellis: Wireless Propagation in Non Line of Sight Urban Areas
Republic
Plaza
Chulia Street
Location E
Cecil Street
(Location A)
OUB Plaza
Battery Road
(Location D)
Location B
Transmitter location
(Ocean Building)
Overpass
(Location C)
CoUyerQuay
Fullerton Square
Drive Route
Figure 1.
UTD model of Raffles Place and drive route.
Table 1.1"' and 2"'' order interactions considered in case 1.
RP
RPRP
DP
VP
RPDP
RPVP
DPRP
VPRP
DPDP
DN
RN
DC
RC
Single Reflection from plates
Double reflection from plates
Single edge diffraction
Single comer diffraction
Plate Reflection- Edge Diffraction
Plate Reflection - Comer DifiRraction
Edge Diffraction- Plate Reflection
Comer Diffraction- Plate Reflection
Edge Diffraction - Edge Diffraction
Diffiaction from End Cap Rim
Reflection from End Cap
Diffraction from Curved Surface
Reflection from Curved Surface
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
169
Table 2.3"* order interactions considered in case 2.
Triple reflection from plates
Double Plate Reflection- Edge Diffraction
Double Plate Reflection - Comer Diffraction
Plate Reflection- Edge Diffraction- Plate
Reflection
Plate Reflection- Comer Diffraction- Plate
Reflection
Edfie Diffraction- Double Plate Reflection
RPRPRP
RPRPDP
RPRPVP
RPDPRP
RPVPRP
DPRPRP
Measured vs. Predicted Signal Power
-2D
+
—
—
o
-40
measured
2nd order (case 1)
3rd order (case 2)
3rd order with penetrable walls (case 3)
D
-BO
S
-BO
9
-100
-120
-140
-1B0
0.2
Figure 2.
0.4
0.8
O.B
Distance (kmj
Measured and predicted signal-power along drive route.
Table 3. RMS error and standard deviation between measured and predicted power along drive route.
Order of interactions
included in model
RMS Error
Standard Deviation
1
js.jnd
33.0 dB
13.2dB
2
,s. 2nd 3rd
22.5 dB
13.1 dB
9.2 dB
8.1 dB
Case
3
]", 2"', 3'" with penetrable
walls
Ellis: Wireless Propagation in Non Line of Sight Urban Areas
170
OUB
Center
Singapore
Land Tower
Exchange
Building
Ocean Towers
Republic
Plaza
Transmitter
(Ocean Building)
Cecil Street
Location B
Location A
Figure 3.
Major Ray Contributions at locations A and B along Cecil Street.
Collyer Quay
Hitachi
Center
Overpass
(Location C)
Restaurant.
Transmitter
(Ocean Building)
Union
Overseas
Shopping
Center
John Hancock
Building
Figure 4.
Major Ray contributions at Oveipass (Location C).
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
171
Standard Chartered
Bank
UOB Plaza
Battery Road
Singapore Land
Tower
OUB
Center
Caltex Tower
Transmitter
Republic
Plaza
Figures.
Location E
Major Ray Contributions at location D.
UOB Plaza
Chulia Street
Sinsov
Building
Transmitter
Location
City House
Figure 6.
Major Ray contribution at Location E.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
172
AUTOMATIC CALCULATION OF BAND DIAGRAMS
OF PHOTONIC CRYSTALS USING THE
MULTIPLE MULTIPOLE METHOD
Jasmin Smajic, Christian Haftier, Daniel Emi
Laboratory for Electromagnetic Fields and Microwave Electronics
Swiss Federal Institute of Technology
ETH-Zentrum, Gloriastrasse 35, CH-8092 Ziirich, Switzerland
E-mail: ismaiic@ifh.ee.ethz.ch
ABSTRACT - In the framework of photonic crystal's band
structure calculations, we present a novel way - based on
several advanced techniques for searching and tracing
eigenvalues with the multiple multipole program - to compute
these diagrams automatically, efficiently, and with a high
accuracy. Finally, we validate the results for some well known
test cases.
Keywords - band diagram; multiple multipole program; photonic
crystal; eigenvalue analysis.
I. INTRODUCTION
PHOTONIC Crystals (PhCs) were proposed in 1987 by E.
Yablonovitch [1] at the University of California, as an
optical counterpart to semiconductors, i.e., PhCs should
provide a photonic bandgap in the same way that a
semiconductor possesses an electronic bandgap. In fact, PhCs
are rarely found in nature. Exceptions are opals and butterfly
wings. However, thanks to nano-technology it has become
possible to fabricate artificial PhCs in the last decade. These
PhCs essentially consist of a periodic assembly of dielectric
scatterers, i.e., there is a strong link to the well-known
structures of grating theory. One of the important differences
between PhCs and semiconductors is the size of the unit cell.
For a semiconductor, one has a 3D grid consisting of identical
atoms, i.e., the lattice constant in all three directions of the
crystal is in the order of the diameter of an atom, whereas the
cell size of a PhC is in the order of half an optical
wavelength, i.e., much larger. From this fact, one expects that
the Photonic Integrated Circuits (PICs) based on PhC concept
must be much larger than the traditional semiconductor ICs.
But - because of the macroscopic size of the PhC's unit cell one has much more freedom in introducing and fabricating
defects in a PhC than in a semiconductor. Note that a
semiconductor becomes interesting from the technological
point of view only when a few impurities or defects are
introduced and the same holds for PhCs. Doping traditional
semiconductors is a rather statistical way of introducing
defect atoms in a semiconductor and therefore, the building
blocks of semiconductor are relatively large blocks of material
with a specific doping. These blocks obviously consist of
many atoms. When designing a "doped" PhC, one can
precisely position and fabricate all defects with a high degree
of freedom - at least it is expected that this can be done in the
near future [2].
Although the variety of PhC structures that might be
fabricated one day seems to be almost infinite and although
many interesting structures were already proposed (various
types of waveguides, sharp bends in waveguides without any
reflection, couplers, resonators, etc.) or even fabricated on a
prototype level, one currently cannot say what kind of PhC
structures will be favored. At the moment, one can neither
know the materials that are best suited for PhCs - it is well
known that a large dielectric contrast is required for obtaining
a bandgap, which somehow limits the materials that may be
used, but there is no unique choice at all - nor what kind of
geometry (2D crystals or 3D crystals [3-5], symmetry, shape
of the scatterers) is most appropriate. Thus, tfiere is a strong
need for theoretical investigations and simulations of
potential structures. The first step of such investigations
consists in the computation of the band diagrams of perfect
PhCs without any defects. The goal is to fmd structures that
may easily be fibricated and exhibit a broad band gap, i.e., a
frequency range where no electromagnetic waves are allowed
to propagate within the crystal. In order to find the band gap,
one must compute the band diagram of the lowest order
modes of the PhC. This is essentially an eigenvalue problem
that exhibits several special cases that may cause difScult
mmierical problems, especially when one is designing a
procedure for the autmnatic, efficient, accurate and reliable
computation of the~"complete band diagrams for arbitrary
structures.
Currently, the most frequently used approach is the Plane
Wave Method (PWM) that mainly approximates the
electromagnetic field by a superposition of plane waves [610]. It is well known, that this method has a problematic
convergence [11-13, 10]. Other methods that were used for
PhCs are the Finite Difference Time Domain (FDTD)
[14,15], the transfer matrix method [16], the Finite Element
Method (FEM) [17], and the Boundary Element Method
(BEM) [18]. In the following, we apply the latest version of
the Multiple Multipole Program (MMP) [19] impleniented in
the MaX-1 software [20].
1054-4887 © 2003 ACES
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
173
In order to obtain efficient, reliable, and accurate results, we
carefully analyze the numerical problems that may occur and
introduce several new techniques. For reasons of simplicity
we focus on the 2D case.
The remainder of the paper is organized as follows: A commonly
used representation of PhCs in terms of their band diagram is
elucidated in Section II. In Section III we briefly explain the
core of our photonic crystal calculations, when MMP is
considered. The proper framework of the eigenvalue search is
reported in Section IV, whereas in Section V a successful
automation of such search procedure is proposed. A
validation of our band structure calculation by means of
various test examples is given in Section VI. And finally, we
conclude this contribution with a short summary in Section
VII.
II. DEFINITION OFTHE BAND DIAGRAM
As an introductory example let us consider the simple case
of a 2D PhC consisting of dielectric rods arranged in a square
lattice and embedded in, e.g., air. For periodic structures it is
possible to apply some fiindamcntal theorems from solid state
physics. The original lattice for this crystal is given on the leH
hand side of Fig. 1. For the dielectric constant we can write
e{r) = e(r + R) where R is one of the original lattice vectors.
According to Bloch's theorem [6, 7] for the modal field inside
the crystal we write:
(1)
If we want to construct the reciprocal lattice, we can use [7]:
/,=2;r
e,(ejXe,)
/,=2;r
5,(ejX5,)
/3 = 2;r
e,(ejXej)
These equations are derived firom the definition of the
reciprocal lattice vector space. For 2D crystals (cylindrical
geometry), the vector/, is omitted and the vector e, is the
unit vector e. along the cylinder axis.
From the equations above, we can conclude that the discrete
translational symmetry of a photonic crystal leads to the fact
that modes with the wave vector k and modes with the wave
vector k + G are identical, i.e., we have periodicity also in the
reciprocal space. A special representation of the primitive cell
for this periodicity is called the first Brillouin zone (l" BZ). It
can be defined as a zone around any lattice point in the
reciprocal space with points that arc closer to this lattice point
than to any other lattice point.
o
o
o
-^
f.
-*o
/l
o
o
Fig. I: The original (left) and reciprocal (right) lattice for a 2D photonic crystal
(square lattice). Construction details for reciprocal lattice are given in the text.
(2)
Note that (1) hold.s not only for the electric but also for the
magnetic field. Bloch's tlieorcm may be proven in classical
electro-dynamics [6]. Important consequences of this theorem are
[6,7]:
1. e'" = \, i.e., k-R = N-2}[, where N is an integer - the
wave vector space (reciprtxal space) is discrete,
2- E- {r) = E- - (;"). '-C- the reciprocal space is
periodic. G is one of the reciprocal lattice vectors.
This allows us to define the so-called reciprocal lattio:"Space,
spanned by the reciprocal lattice vectors. We first define the
original lattice vectors as follows:
R=;7,e,+7;jej+7jC,
(3)
where e,, Cj. e, are three independent lattice vectors and j^^,
Fig. 2; Construction of the 1" Brillouin zone (solid square), its iircducible part
(triangle P-X-M) and characteristic points for band stniclure computation (F,
X.andM).
i7j, ;;, arc integer numbers. Note that e, is missing in 2D
crystals. Similarly, we write for the primitive reciprocal
lattice vectors:
G = lcJ,+lcJ, ^K-J,
(4)
The Brillouin zone construction (using Bragg's planes dashed lines) for the square lattice is shown in Fig. 2.
Because of the high degree of symmetry, wc need to analyze
only a small part of the 1" BZ. This part is called the
irreducible BZ (IBZ), [6, 7]. In the case of periodic strucUires,
Smajic and Hafher: Automatic Calculation of Band Diagrams of Photonic Crystals
2E-6-
174
2E-6-
- ,'' \ Ipcriodic boiindaryl
Jr'/'
1E-6-
'\ Iperiodicboundaiyl
f+ei'
fictitious
excitation
fr+T
2E-6
Bessel, and multipole
expansions
II I I 11 I I
1E-6
2E-6
Bessel, and multipole
expansions
Fig. 3: The unit cell of the photonic crystal with dielectric rods arranged on a
square lattice.
Fig. 4: The unit cell of the photonic crystal with dielectric rods airanged on a
hexagonal lattice.
it is sufficient to perform the modal field analysis in the area
of the IBZ. As illustrated in Fig. 2 the IBZ for a square lattice
is a triangle with the corners T, X, and M. Since the maxima
and minima of the eigenvalues (resonance frequencies) are
supposed to be on the borders of the IBZ, it is sufficient to
trace the eigenvalues along the sides of the IBZ in order to get
the photonic bandgaps. Therefore, the standard band diagram
consists of three sections: T-X, X-M, and M-F (see Fig. 5).
For other lattices, the procedure is essentially the same [21,
22]. Assume that an arbitrary point in the reciprocal space is
considered. This point essentially defines a wave vector. For
the periodicity of (3) we then obtain for the field in the
original space:
subdivided into cells by an appropriate grid of fictitious
boundaries (dashed lines in Fig. 3 and Fig. 4). Assume that
the field in one of the infinitely many cells is known, then,
the field in all other cells is easily obtained fi-om the
periodicity conditions (5), i.e., the Floquet theorem [7].
The geometric shape of the original cell depends on the
crystallographic structure (i.e. the crystal symmetry), but it is
not unique for a given crystallographic structure at all,
because the fictitious boundaries we have introduced, are
quite ambiguous. For example, in Fig. 3 we used straight
lines between the circular rods. We could replace these lines
by curved, periodic lines and we could move these lines to
any other position in space. Since we will impose so-called
periodic boundary conditions along the fictitious boundaries
of the original cell, we have to minimize the numerical
problems when we select the fictitious boundaries in such a
way that the electromagnetic field along them is as well
behaved as possible. Therefore, straight lines in the middle
between neighbor rods are most reasonable when the rods are
circular or rectangular. When the geometric shape of the rods
is more complicated, it may be advantageous to use curved
lines.
Once, the original cell is isolated by introducing fictitious
boundaries, we can derive boundary conditions for the field
along them. In 2D PhCs, the original cell is boundeJby two
pairs of parallel lines. For example, wheit F is a point on the
left border of the -original cell in Fig. 3, r-^^e^ is the
corresponding point on the right border, where e,
#(F + ^) = F(F) • e' * * = F(r) e' (C,tl,|«,| + C,il!ftl)
(5)
where F stands for the electric as well as for the magnetic
field. In the MMP implementation of the periodic boundary
problem C\. Q are parameters that characterize the point in
the reciprocal lattice space. As a consequence, it is sufBcient
to know the field in a unit cell (as an equivalent
representation of the primitive cell) of the original space. Let
us call .this the original cell. Note that neither the shape nor
the location of the original cell is unique, but for both the
square and the hexagonal lattice we may simply use
quadrangular original cells as shown in Fig. 3 and 4.
For the square lattice, the relation between the periodic
constants (C/, Q) and the position in the IBZ is very straight
forward, i.e., these are the Cartesian components (C„ Cy) of
the wave vector k . For the hexagonal lattice, the situation is
a bit more complicated [23-25].
III. THE MMP SOLUTION OF PERIODIC PROBLEMS WITH
FICTITIOUS BOUNDARIES
Any software for computing band diagrams must handle both
eigenvalue problems and periodic structures. The MMP
implementation of MaX-1 contains a simple concept for
handling arbitrary periodic structures: First, the structure is
corresponds to one of the primitive lattice vectors. Because of
the periodicity, we obtain from (5):
F(? + e,) = F(r)-e'"^'l''l'
(5-)
This condition holds for both the electric and the magnetic
field in every point along the right boundary of the original
cell. We call this the periodic boundary condition that is
imposed on the right border of the original cell. Similarly, we
can introduce a periodic boundary condition for the upper
border.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
175
Fig. 5: The band diagram of the photonic ciystal with dielectric rods on a square lattice (for H-polarizalion).The algorithms used within the
eigenvalue .search procedure arc labeled correspondingly.
Having defined the original cell and its periodic boundary
conditions, one has to set up the MMP model of the scattering
body in the lattice point: We approximate the field in each
domain by a superposition of multipole expansions and sometimes by additional, analytic solutions of Maxwell's equations
(in the frequency domain). The amplitudes or parameters of
the resulting scries expansions are then computed with the
generalized point matching technique, i.e., b>' minimi2dng a
weighted error fimction defined along all natural and fictitious
boundaries. For example, for the simple geometry in Fig. 4
we use the following expansions:
E:=Y, ^.^-'«("') cosCn?)) + B,^J,(w)sin(n^)
(6)
i>=0
B"J„(Kr)sm(nip)
(7)
£, = 5] C^Hl (w) cos(,np) + D^H',(,io) sm{n p)
(8)
H,
^J^A"J.(K>)cos{n<p) +
it=0
-
11=0
—
H, =Y.C" Hl{Kr)cos(np)+ D'J Hl(Kr)sm{n</>)
(9)
where J„ is Bessel function of order n, HJ is Hanlcel function
of first kind and order n, vis transverse propagation constant
and (r,^) are polar coordinates with respect to the origin at
the position of the corresponding expansion. Expansion (6)
(Bessel expansion) is used in the case of E-polarization and
expansion (7) is used in the case of H-polarization. These
Bessel expansions are used for the domain inside of the
dielectric rod because these' functions have no singularity at
origin. Furthermore, these expansions are sufficient because
the domain is simply connected. The background domain is
not simply connected, because it contains a hole. Therefore,
we need at least two different expansions, namely a multipole
expansion (8) or (9) and Bessel expansion (6) or (7). Note
that the moltipole expansion essentially accounts for the field
scattered at the inner boundary, whereas the Bessel expansion
accounts for the outer, ficitious boundaries. This means that
the Bessel expansion simulates the field that comes from all
rods outside the original cell. According to Vekua [24], our
set of expansions is complete in the sense that the error of the
field is below an arbitrarily small value e provided that the
highest orders are big enough and provided that the
amplitudes {A. B, C. D in (6)-(9)) are computed correctly.
IV. THE MMP-MAS EIGENVALUE SOLVER
For obtaining the band diagram of a PhC, it is necessary to
solve an eigenvalue problem, because there is no excitation
like in scattering problems. This means that we only obtain
non-trivial solutions (i.e. frequencies) for an arbitrary point of
the IBZ_{i.e.^ for a given set of complex values C/, Ci). Thus,
we essentially have a periodic resonator problem to solve. Thesearch of resonance frequencies in the NfMP code MaX-1 is
somehow different fi-om many other numerical methods
because MMP uses a fiill rectangular system matrix obtained
from the generalized point matching technique. For such type
of matrix it is very demanding to obtain accurate results while
avoiding problems with the condition number [25]. Note that
condition number problems are especially crucial when one is
solving eigenvalue problems. If this is not properly done, one
can obtain a "noisy" behavior near the eigenvalues and this
can heavily disturb the numerical eigenvalue search
procedures. However, the standard MMP eigenvalue search
procedure first defines a real valued, positive eigenvalue
search function
Smajic and Hafner: Automatic Calculation of Band Diagrams of Photonic Crystals
*'-§S
(10)
where e is the eigenvalue (i.e. the resonance frequency), E is
the weighted residual, and A is an amplitude that may be
retrieved from any field component in a specific test point (or
an integral over some field profile). For the band gap
computation, it is most reasonable to define J^ as the total
electromagnetic energy within the original cell. According to
(10) the desired eigenvalues are characterized by the minima
of the search function 7. Analyzing the shape of j] near the
minima provides additional information on the accuracy of
the solution.
Although more reliable results are obtained when the
1E13-r
i|iiiijiiii|iiiijiiii|riiiiiiji|iiii
1.2E14 1.6E14
frequent^ (Hz)
2E14
2.4E14
Fig. 6: The behavior of the eigenvalue search function (in the F point) of one
single model with four different "last" expansions (the order of the last
expansion is labelled in the figure)
176
amplitude is defined by an appropriate integral, the definition
in one or a few test points is sufficient for most cases. Since
the numerical intergration may be time-consuming, one
usually prefers the simpler test point method. However, it is
important to note that the definition of the search fimction is
not unique. By defining different search functions, one can
gain even more intrinsic information providing a good error
estimation and for validation purposes. As depicted in Fig. 6,
even for a single model (fixed amplitude definition and fixed
multipole expansion), one can address the different minima of
the same eigenvalue search fiinctions simply by rearranging
the columns of the MMP matrix. In feet, in the Givens update
algorithm [25], which was used for solving the MMP matrix
equation, the last expansion somehow plays the role of an
excitation. When it happens that the spatial symmetry of such
excitation is not contained in the symmetry of the searched
eigenmode, this mode will not be "excited", hence, the
corresponding minimum of the eigen-value search fimction is
suppressed. Although, it may be desirable to suppress some
modes in applications where not all modes must be
considered, this is usually inconvenient for the automatic
computation of the complete band structure. We therefore
look for an alternative technique.
Remember that we have introduced fictitious boundaries for
handling the periodic problem. Similarly, we now can
introduce a fictitious excitation that is defined in such a way
that all modes are excited (Fig. 7). This concept mimics the
measurement of resonance frequencies, where one always
needs an excitation of the resonator and a test point (or port)
where the signal is measured. By sweeping the fitequency of
the exitation, the peaks of the amplitude A in the test point
can be readily assigned to the resonance frequencies of the
different modes. This procedure was first introduced by the
Method of Auxiliary Sources (MAS) [26] and a similar
method was used by Sakoda [27]. Finally, the method was
adapted to MMP by Moreno [28]. MAS uses eigenvalue
\-S
ranclDin I
position I
1E-11
TTT
4E'13
ii|iiii|iiiijiiii|iiiijnii|iiiiiiiii|Mii
8E13
1.2E14 1.6E14
frequeniy (Hz)
2E14
2.4E14
Fig. 7: The behavior of the eigenvalue search function (for a k-vector in the F
point) using the fictitious excitation in a random and symmetric position
respectively.
^^ '' '*"''' " " " " I'' " I'' T M " " I' M 11 I I I I I I
Z.136E14
2.139E14 2.142E14
frequoni^ (Hz)
2,145E14
Fig. 8: The "twin minima" phenomenon, behavior of the eigenvalue search
fiincuon within the eigenvalue search procedure using a randomly located
fictitious excitation.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
177
search fiinctions |i such as the energy density ^ at the test
point are used. The eigenvalues are then obtained from the
maxima of ji. The analysis of H near the maxima has yielded
a strange "double peak" phenomenon that disturbs the
numerical search procedure. The standard MMP-MAS
eigenvalue solver searches for minima of the eigenvalue
search function r\ = Vil = \lA^, ie., one obtains "twin
minima" instead of double peaks, as shown in Fig. 8. The
"double peak" phenomenon and the "twin minima" are caused
by a very sharp peak of the residual £ at the correct
eigenvalue position. Note that this peak is not obtained in the
standard MMP approach without fictitous excitations. Of
course, the residual peak may also be used for defining the
eigenvalues. Since these peaks are extremely sharp, it is very
likely that one of the eigenvalues is missed by the rough
search routine that searches for all eigenvalues. In order to
overcome these problems, one can define more complicated
eigenvalue search functions ;; as proposed in Fig. 8. This allows
one to suppress the double peak phenomenon. Unfortunately,
one may encounter numerical underflow problems in some
applications. Therefore, the current MaX-1 eigenvalue solvers
uses three different "competing" eigenvalue search functions;
1) A complicated one with user-definable exponent n, 2) the
inverse of the amplitude, and 3) the proper residual. Using all
of these three functions, the code is capable to detect the
correct locations of the eigenvalues. An alternative to
overcome the twin minima problems is the introduction of
"fictitious losses" that smoothen the resulting search function
nSince one oilen considers a broad frequency range, it is not
reasonable to find the eigenvalues by plotting the eigenvalue
search function over the entire range with a very high
resolution. It is much more efficient to subdivide the search
process into two steps: 1) Rough detection of all eigenvalues
and 2) fine search, i.e., accurate computation of the
eigenvalues. The first step seems to be trivial as soon as the
problems mentioned above have been solved. The second step
requires a fast minimum search procedure for real functions.
The algorithm used in MaX-1 is mainly based on a parabolic
interpolation because the search function near the minima is
usually almost parabolic - provided that the double peak
phenomenon has been removed.
ilefhition of the set 3f starting pcints
:Si|S, = (C,,C,)efliz,M,2,3)
X
Sc =
SI
, current sterling p3inl
*
Irequenci' renge definlior tor th3 rougi eigenv. search
(f 1, f d, Nt - lumber o' steps
M, - ni-mter 3f noces (eigenvaljes)
5
Ec= 1 ,'ine eijenvBlie search inte*
Ef = r or K or M, endrg point defiiition
a
'.S
•a
• SH
.1U
av
rougi e genvdue searcn
V
5
5
line search lor Ec
pare meter eslimati en techniqje;PET)
Cx.Cy tedefintion
Ec++
Sc++
Enc
Fig. 9: The algorilhm for Ihc band strucluri; computalion using MMP.
! from gam ma to x, lett part
read win square_gjt_m_g_middle.win
drawwmdowl In
reodturctionm.fun
drawvtundiool 2 1
readlundlonxjun
draiw fund ion 121
rirBwt«t173 4Q01"x"
draw text 348 4001 "m"
drawtext 0 4001 "gamma"
drawlKt 5B0 4001 "gamma"
witefundion • tieflder 01
vwite function / test "rC
wrtefunction /terf "triltouin zone"
write function /text "fat"
wrtefun !
read mmp bandOOl .mmp
se* eigenvalue fine 1
end pre prccessing drectives
!
-»• set period constant 1S7079B .3000
set VBr41
wrtefunction + header 0 2
«rHefunctian/tod"rr
wrtefunction /text "triltouin zone"
wrtefunction /ted "fat'
end outer loop pre dredfi^es
- !
mmp solve
wrtefun / var
wrtefun/Ire rea
increase period constant -39269.9075 0 0 0
increase ver -1
end inner loop direcfwes
wrtefun I
read function/
process fundion c)ivCi/2.3e14)
draw fund kin 1 31
increase eigenvElue 1
end outer loop post diredives
-Vi
read dir sijuarej j(_m^jiiiddle_2.dir
Fig. 10: The algorithm for band diagram compulation written in the MaX-1
saipt language.
178
Smajic and Hafner Automatic Calculation of Band Diagrams of Photonic Crystals
TABLE I
CONVERGENCE CHARACTERJSTICS, COMPUTATION FOR 1" AND 6" ElGENFREQUENCY AT X POINT OF IBZ
Number of
unknowns
Elgenftequency 1
Frequency (Hz)
I.0223585el4
1.009S9SSeI4
1.0078338el4
1.0074678el4
1.0075153el4
20
36
S2
94
164
Enor(%)
1.473
0.206
0.032
0.005
0.000
Eigenfrequency 6
Field mismatch. (%)
9.620748e-0
4.204639e-0
0.288115e-0
4.465747e-2
4.72972 Ie-7
Having a closer look to typical band diagrams (Fig. 3), we see
different situations which can cause problems for both the
rough search and the line search. Mainly at the T and the M
point we usually observe degenerate modes. Furthermore, we
have areas with almost degenerate modes and points where
different lines seem to cross each other, where the modes are
(accidentally) degenerated. When the rough search is
performed to degenerate points, it usually cannot detect all
modes involved. Even if the search procedure is started in a
close vicinity to such degeneracies, it will be too timeconsuming to iterate into all eigenmodes. In order to
overcome these problems, it is reasonable to start a rough
search in a domain where all eigenvalues could be easily
tracked down (e.g. the interval between r and X in the band
diagram of Fig. 5). Once this has been done, one can trace
each eigenvalue by moving a small step either to the left or
right side within the band diagram, and repeating this
U.O .
0.6-i
Frequency (Hz)
2.3465308el4
2.2567438el4
2.2S12072el4
2.2519421el4
2.2S20785el4
Enor(%)
4.194
0.207
0.039
0.006
0.000
Field mismatch. (%)
2.306028e+I
4.449937e+0
0.82436Se4
O.I28581e4)
3.551224e^
procedure until the border of the diagram is reached. For each
such step, only a fine search must be performed. Depending
on 1) the desired accuracy, 2) the step size, and 3) special
properties of the model, several iterations are required. The
number of iterations could be drastically reduced when using
the Eigenvalue Estimation Technique ^ET) implemented in
MaX-1 [19]. This technique uses the information of previous
eigenvalue solutions for the extrapolation of the current
eigenvalue's search interval. Typically, 4-12 iterations per
step are sufficient for obtaining an eigenvalue with a high
precision. For example, for tracing the first mode in Fig. 3,
280 search steps were performed and 5 iterations per step
were required in the average.
V. AUTOMATIC EIGENVALUE SEARCH
Referring to e.g. Fig. 5 a standard band diagram consists of
three different intervals corresponding to the three sides of
^
\^
^• ""^--S
^==-
eoa
^
2xc
/
n-
X
"\
r
M
Fig. 11: The band diagram of the photonic crystal with dielectric rods and square
lattice, H-poIarization, the first 6 modes.
X
M
r
Fig. 13: The band diagram of the photonic crystal with dieleciric rods and
hexagonal lattice, H-polarization, the fust 6 modes.
U.D .
- '
'~~"**-*^
0.6-:
ma
2jn:
—
A
0.4-
=> -^
0.2^
rFig. 12: The band diagram of the photonic crystal with dielectric rods and
square lattice, E-polarization, the first 6 modes.
3c
Ii
r
Fig. 14: The band diagram of the photonic crystal with dielectric rods and
hexagonal lattice, E-polarization, the first 6 modes.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
179
the IBZ. When the rough search is started somewhere in the
middle of such an interval (e.g. in the area between Fand X
in the band diagram), it must be repeated three times. After
each rough search the fine search must be repeated for each
obtained eigenvalue and, finally, the fine search routine must
nm for each eigenvalue once towards the left and once
towards the right side of the band diagram, as depicted in Fig.
5. MaX-1 contains a script language that allows one to define
complicated procedures such as the search procedure mentioned
above. The set of MaX-1 directives for the automatic
generation of a band diagram from the point in the middle
between T and X to the P point, is given in Fig. 10, and the
complete algorithm for this procedure is given in Fig. 9. It is
obvious that the algorithm is not simple and the overall
procedure relies on fast computer resources.
automatic algorithm was developed and evaluated along
several examples. The eigenvalue searching procedure in the
frequency domain has been performed using a fictitious
excitation. Optimal eigenvalue search fijnctions have been
found while evaluating the total eigenvalue spectrum for kvalues at three preferable points on the IBZ. The three
resulting sets of eigenvalues are evolved into a full band
diagram using a highly efficient Eigenvalue Estimation
Technique (EET). The overall algorithm performs photonic
band diagram calculations at a very high level of accuracy
and at reasonable computational costs. This algorithm is
easily extendable for applications involving localized defect
mode analysis [30], various PhC defect waveguide types
(supercell approach [31]) and photonic waveguide
discontinuities [31], as well.
VI. NUMERICAL VERIFICATION
ACKNOWLEDGMENT
We have applied MMP to various PhC lattices. Internal tests
show excellent convergence. Therefore high accuracy may
easily be obtained. Table 1 shows the MMP estimate of the
mismatching errors along the boundary for the model
outlined in Fig. 3 with different maximum orders of the
multipolcs and BesscI expansions, i.e., with different numbers
of unknowns. Note that the computation time typically is
proportional to the cube of the number of unknowns because
we use a brute-force flill matrix solver (Givens update
scheme). Despite of this, the computation time remains
reasonably short because the matrices obviously are much
smaller than the matrices used in other methods. For example
Fig. 11 was obtained with 3 rough-search routines, 100
frequency step.'; each. The total number of 1656 plotted points
required were then computed with 8280 MMP evaluations of
;;, i.e. approximately 5 iterations per point in the diagram
were performed. The total calculation time was 40 minutes on
a Pentium 4, 2GHz. Because of the excellent convergence, we
also can estimate the accuracy of the eigenvalues by
comparing them with a very accurate MMP model. As one
can sec from Table I, one only obtains one more digit when
doubling the number of unknowns.
In order to validate this algorithm, several calculations were
performed and results were compared with the results of MPB
package developed at the MIT [29]. For the PhC with square
lattice and dielectric rods (Fig. 3), a band diagram calculation
was performed for different field polarizations and the results
are given in Fig. 11 (H-polarization) and Fig. 12 (Epolarization). The results for the hexagonal lattice case (Fig.
4), are depicted in Fig. 13 (H-polarization) and Fig. 14 (Epolarization). These two types of PhC rely on the same lattice
data: A dielectric rod with radius r = 0.3a and a dielectric
constant of £ = 11.56, the lattice is embedded in air and the
lattice constant is a = I0"'(m)• From Figs. 11-14 we deduce a
perfect agreement with the MPB results documented in [29].
This work was supported by the Swiss National Science
Foundation in the fi-amework of the project NFP-2000065102.01/1 and the research initiative NCCR Quantvim
Photonics.
REFERENCES
[I]
(2)
[3]
[4]
[5]
(6)
[7]
(8)
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
VII. CONCLUSION
We have presented an efficient method for band structure
calculation for 2D dielectric PhCs. In this framework a fully
[17]
E. Yablanovich. "Inhibited spontaneous emission in solid-.state physics and
elcctronicii", Phys. Kcv. Leu., 58, pp. 2059-2062, 1987.
J. Mills, "Photonic ciystals head toward the marketplace", Nov. 2002,
httn://optics.onr/anides/ole/7/l I/I/I
K. Busch, S. John, "Liquid-ciystal photonic-band-gap materials: the
tunable electromagnetic vacuum", Pbys- Rev. Letters. 83, pp. 967-970,
1990.
A. Figotin, Y. A. Godin, "tw'o-dimensienal tunable photonic crystals",
Pliys.Rev.B. 57,pp.2841-2848,1998.
R. L. Sutherland, V. P. Tondiglia, L. V. Natarajan, S. Chandra,
"Swiichable orthorombic F photonic crystals formed by holographic
polymerization-induced phase sepatution of liquid ci>'star. Optics Express,
10, pp. 1074-1082,2002.
K. Sakoda, "Optical Propeities of Photonic CiyslaIs",Springcr, Berlin, 2001.
J. D. Joannopoulos, R. D. Mcade, J. N. Winn, "Molding the Flow of
Light", Princeton University Press, 1995.
K. M. Ho, C. T. Chan, and, C. M. Soukolis, "Existance of a Photonic Gap
in Periodic Dielectric Structures", Phys. Rev. Leu., 65, pp. 3152-3155,
1990.
S. G. Johnson, and J. D. Joannopoulos, "Block-iterative frequency-domain
methods for Maxwell's equations in a plancwave baids". Opt. Ejipress, 8,
pp. 173-190.2001.
R. D. Meadc, A. M. Rappe, K.D.Brommer.J. D. Joannopoulos, and O. L.
Alerhand. "Accurate theoretical analysis of photonic band-gap materials",
PAM. Rev. B, 48, pp. 8434-8437,1993.
H. S. Souzucr, J. W. Haus, and R. Ingiiva, "Photonic bands: Convergence
problems with the planc-witve method", Pbys^ Rev-. B, 45, pp. 1396213972,1992.
D. Hermann, M. Frank, K. Busch and P. Wollle, "Photonic band stnicture
compulation.s". Opt. Express. 8, pp. 167-172,2001.
K. Ohtaka, T. Ueta, and K. Amemiya, "Calculation of photonic bands
using vector cylindrical waves and reflectivity of light for an airay of
dielectric rods". Phys. Rev B. 57, pp. 2550-2568,1998.
C. T. Chan, Q. L. Yu, and K. M. Ho, "Order-N spectral method for
electromagnetic waves", Phys. RevB, 51, pp. 16 635-16 642,1995.
M. Qui and S. He, "A non-orthogonal finite-difference time-domain
method for computing the band structure of a two-dimensional photonic
crystal with dielectric and metallic inclusions", /. Appl. Phys., 87, pp.
8268-8275,2000.
J. B. Pcndry and A. MacKinnon, "Calculation of Photon Dispersion
Relations", Phys. Rev. Lett.. «9, pp. mi-mi, 1992.
W. Axmann and P. Kuchment, "An efficient finite element method for
computing spectra of photonic and acoustic band-gap materials", J.
Compui. Phys., 150, pp. 468-481,1999.
Smajic and Hafner: Automatic Calculation of Band Diagrams ofPhotonic Crystals
[18] P. A. Knipp, and T. L. Reinecke, "Boundary-element calculations of
electromagnetic band-structure of photonic crystals", Physica £, 2, pp.
920-924,1998.
[19] Ch. Hafiier, "Post-modem Electromagnetics Using Intelligent MaXwell
Solvers", John Wiley & Sons, 1999.
[20] Ch. Hather, "MaX-l: A visual electromagnetics platfonn", John Wiley &
Sons, 1998.
[21] P. R. Villeneuve and M. Pichc, "Photonic band gaps in two-dimensional
square and hexagonal lattices", Phys. Rev. B, 46, pp. 4969-4972,1992.
[22] M. Plihal and A. A. Maradudin, "Photonic band stracture of twodimensional systems: the triangular lattice", Phys. Rev. B, 44, pp. 85658571,1991.
[23] Ch. Hafner, J. Smajic, The Computational Optics Group Web Page (IFH,
ETH Zurich), httn://alDhard.ethz.ch/.
[24] I. N. Vekua, "New methods for solving elliptic equations", North-Holland,
Amsterdam, 1967.
[25] G. H. Golub and C. P. Van Loan, "Matrix Computations", John Hopkins
University Press, Baltimore, 1996.
[26] F. G. Bogdanov, D. D. Karkashadze, and R. S. Zaridze, in "Generalized
Multipole Techniques for Electromagnetic and Light Scattering", edited by
T. Wriedt, pp. 143-172, Elsevier, Amsterdam, 1999.
[27] K. Sakoda, N. Kawai, T. Ito. A. Chutinan, S. Noda, T. Mitsuyu, and K.
Hirao, "Photonic bands of metallic systems. I. Principle of calculation and
accuracy", Phys. Rev. B, 64, pp. 045116,2001.
[28] E. Moreno, D. Emi and Ch. Hafiier, "Band structure computations of
metallic photonic crystals with the multiple multipole method", Phys. Rev.
A 65, pp. 155120:1-10,2002.
[29] S. G. Johnson and J. D. Joannopoulos, The MIT Photonic-Bands Package
borne page, http://ab-initiQ.niitedu/mph/.
[30] R. D. Meade, K.I D. Brommer, A. M. Rappe and J. D. Joannopoulos,
"Photonic band states in periodic dielectric materials", Phys. Rev. B, 44.
pp. 13 772-13 774,1991.
[31] E. Moreno, D. Emi and Ch. Hafner, "Modeling of discontinuities in
photonic crystal waveguides with the multiple multipole method", Phys.
Jiev.fi, 66, pp. 036618,2002.
180
Jasmin Smajic was bom in Tuzla, Bosnia and
Herzegovina, in 1971. He received a Dipl. El.Ing. degree fix)m Faculty of Electrical
Engineering, Tuzla in 1996. Since 1996 he has
been working at die Faculty of Electrical Engineering
in Tuzla on numerical calculation of electromagnetic
field, numerical mathematics and optinuzaticn. He
got a M.Sc. degree in 1998 from the Faculty of
Electrical Engineering and Computing in Zagreb,
Croatia, for the analysis of magnetic field in nonlinear material. He received a PhD. in 2001 from
Faculty of Electrical Engineering and Computing
in Zagreb, Croatia, for numerical calculation of
time-varying field in non-linear material and electrical machine design
optimization. Since 2002 he is a postdoctoral research fellow in the Computational
Optics group of the Laboratory for Electrcmagnaic Fields and Mioowave Electionics
at the ETH Zuridi. His current research interest includes numerical field cakubtion
and desi^ optimization of photonic crystal devices.
Christian Hafner was bom in Zurich,
Switzerland, in 1952. He received a Dipl. El.-l'ng.
Degree, Doctoral Degree, and Venia Legendi
from Swiss Federal Institute of Technology
(ETH), Ziiridi in 1975, 1980, and 1987
respectively. In 1999 he was given the title of
Professor.
Since 1976 he has been working at the ETH on Ibe
devebpnmt of methods for computational
electromagnetics and for optimization problems. He
has developed the Multiple Multipole Program
(MMP), the MaX-1 package, die Generalized
Genetic Programnung (GGP) code, and various
optimization codes. He worked on various applications (electrostatics, EM scattering,
antenna, waveguides and waveguide discontinuities, gratings, diiral media, etc). His
cun-ent focus is on photonic crystals, rrucrostructuFed optical fibers, and Scanning
Nearjeld Optical Microscoes (SNOM). In 1990 he obtained the second prize of the
Seymour Cray avard for scientific computmg and in 2001 he has been awarded tite
2000 Outstaiiding Journal Paper Award by the Applied Computational
EkctiDmagnetics Sodely. He is member of the Electromagnetics Academy.
Daniel Erni was bom in Lugano, Switzerland, in
1961. He received an El.-Ing. HTL degree from
Inlerkantonales Tcchnikum Rapperswil HTL in
1986, and a Dipl. El.-Ing. degree from Swiss
Federal Institute of Technology (ETH), Ziirich in
1990, bodi in electrical engineering. Since 1990
he has been woridng at die Laboratory for
Electromagnetic Fields and Microwave Ekaronics
(ETH) on nonlinear wave propagation, laser dk}de
modeling (multi-secUon DFB and DBR lasers,
VCSELsX computational electromagnetics and on
the design of non-periodic optical vraveguide
gratings e.g. by means of evolutionary algorithms.
He got a Ph.D. degree in 1996 for die investigation of non-periodic waveguide
gratings and non-periodic coupled cavity laser concepts. His current research
interests includes highly multimode optical -signal transmission in optical
interconnects (i.e. in optical backplanes with extremely large waveguide cross
sections) as well as alternative waveguiding concepts for dense integrated optical
devices like e.g. photonic crystal devices, couplers and WDM filter structures. In
2001 he has been awarded die 2000 Outstandng Journal Paper Award by the ApplM
Computational Bedromagnetics Society for a contribution on die application of
evolutionary cptimizalion algorithms in computational opti'cs. Dr. Emi is the head of
die Communication Photonics Group at ETH Zurich ftvww.Dhotonics.ee.ethz.chV He
is also a member of the Swiss Physical Society (SPS), of die German Physical
Society (DPG), of the Optical Society of America (OSA), and of the IEEE.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
181
Analysis of Conventional and Novel Delay Lines: A Numerical Study
Omar M. Ramahl
Mechanical Engineering Department, Electrical and Computer Engineering Department,
and CALCE Electronic Products and Systems Center
2181 Glenn L Martin Hall
A. James Clark School ofEngineering
University of Maryland
College Park. MD 20742. U.S.A.
oramahi@calce.umd.edu
www, enme. iimd. edu/EMCPU
Abstract
Delay lines are convenient circuit elements used
to introduce delay between circuit board components
to achieve required timing. Serpentine or meander
lines are the most common of delay lines. These lines
introduce delay but also introduce spurious
dispersion that makes the signal appear as if it is
arriving earlier than expected. The cause of such
spurious speed-up or skew is analyzed qualitatively.
Previous work found that owing to the periodicity
inherent in the serpentine line structure, the crosstalk
noise accumulates synchronou.ily, thus creating a
higher potential for triggering false logic. Numerical
simulations arc performed using the Finite-Difference
Time-Domain (FDTD) method to corroborate the
qualitative prediction with physical behavior. Based
on the understanding of the coupling mechanism in
periodic serpentine lines, a qualitative prediction can
be made of the behavior of novel delay lines such as
the spiral line. A new delay line, the concentric Cs
delay lines, is introduced. The design of the new line
is based on forcing the crosstalk noise to spread over
time, or to accumulate asynchronously, thus
enhancing the integrity of the received signal.
large. As a consequence, a delay line used to
introduce precise delay (based on the length of the
line and board material) between circuit board
elements can no longer be considered to have TEM
or quasi-TEM propagation behavior, and to be
electromagnetically isolated from neighboring lines
that fall within its proximity. With the timing budget
shrinking to the picosecond regime, any nonpredictable behavior of delay lines can potentially
cause a timing imbalance at the receiver end.
Kev Words: Delay Lines, Serpentine Lines.'Meander
Lines, Spiral Lines, FDTD, Numerical Simulation,
Crosstalk Noise.
When serpentine lines arrused in high-speed
digital circuit applications, the time delay through a
single serpentine line can be much longer than the
rise time of the signal pulse. Under such conditions,
serpentine lines have been found to introduce a
dispersion that makes the signal appear as if it is
arriving earlier than would be expected based on the
exact electrical length of the line [1]. This spurious
speed-up in the signal can be characterized as a type
of dispersion. This dispersion is primarily caused by
the topology of the line and is directly related to the
crosstalk between the adjacent transmission line
sections. More specifically, the dispersion is related
to two fundamental geometrical parameters: the first
I. Introduction
It can be argued that timing problems are
amongst the most serious plaguing high-speed digital
circuit-board design. As the clock speed increases,
the wavelength shrinks. For instance, for a clock
speed of 1.5 GHz, the pulse harmonics containing
sufficient energy reach beyond 10 GHz, making a
typical motherboard or daughter cards electrically
Two mechanisms are employed to achieve
required signal delay between circuit components.
The first is achieved through internal electronic
circuitry. The second mechanism, which is the most
common and least expensive, is achieved through
meandering a transmission line as shown in Fig. 1.
The meandered line, commonly referred to as the
serpentine line, consists of a number of closely
packed transmission line segments.
The only
objective from meandering the line is to achieve high
density (of transmission line) per square inch of
circuit board while obtaining signal delay that is
directly proportional to the length of the line.
1054-4887 ©2003 ACES
Ramahi: Analysis of Conventional and Novel Delay Lines
182
coupling between non-orthogonal segments. We use
the FDTD method to analyze the spiral line and show
that more inclusive performance can be predicted.
Based on the qualitative insight gained from the
analysis of the serpentine and spiral line, we
introduce a new delay line that we refer to as the
concentric Cs delay line. Qualitative analysis of the
new lines is given followed by fiill-wave threedimensional simulation.
Fig. 1. A matched seven-section seqjentine delay
line. The Characteristic impedance of each identical
section is R^.
is the length of each serpentine section, and the
second is the spacing between adjacent sections,
which influences the capacitive and inductive
coupling [1], [2].
Based on these findings, Wu and Chao
introduced the flat spiral delay line to force the
crosstalk noise to accumulate asynchronously [3].
The flat spiral line is distinguished from the classical
spiral delay line, which requires three-dimensional
topology and therefore cannot be printed on a single
board layer. In this work, the flat spiral line will be
referred to as the spiral line.
The spiral delay line is considered a significant
improvement over the serpentine line. In [3], the
motivation behind the spiral line was introduced and
analysis was given based on the wave-tracing
technique which did not take into account higherorder mode coupling and comer effects. Also in [3],
quantitative analysis was also performed to include
the effects of multiple and feedback couplings
between lines. To incorporate high-order effects into
the modeling of serpentine line, especially the
coupling between the non-orthogonal lines, the threedimensional Finite-Difference Time-Domain (FDTD)
method was used to analyze serpentine and spiral
lines [4]. In a subsequent work, the FDTD method
was used to analyze the spiral line for Gaussian pulse
excitation [5]. More recently, the Method of
Moments (MoM) was used to present a fiill-wave
model for the serpentine line including physical
losses [6].
In this work, we use the same qualitative wavetracing methodology that was adopted in [2] to
analyze the spiral line. We show that the wavetracing methodology can only give qualitative
analysis of the line and that a full-wave model is
needed to accommodate the effect of comers and
The organization of this paper is as follows. In
section 2, the mechanism of coupling between two
transmission lines is introduced with application to
serpentine lines. Sections 3 and 4 present a
qualitative analysis followed by full-wave numerical
simulation of practical real-world structures for the
serpentine and spiral lines, respectively. Section 5
introduces the Concentric Cs delay lines. We note
that the FDTD method will be used purely as a
numerical simulation tool without giving details of
the simulation parameters involved such as the cell
size and time step, ...etc. The paper concludes with
critical observations that can be used in the design of
predictable-delay lines.
IT. Weak Coupling Crosstalk
a. Analysis
A typical serpentine delay line is composed of
transmission line sections closely packed as shown in
Fig. 1. Let us isolate two adjacent sections as shown
in Fig. 2. After examination, we notice that the two
isolated sections resemble the simplest case of two
parallel transmission lines [2]. After ignoring higherorder effects, the two-segment serpentine line shown
in Fig.2 is equivalent to the two transmission lines
matched at both ends and shown in Fig. 3. Once this
observation is made, using conventional transmission
line analysis, one can predict the crosstalk induced on
the second line due to the launched signal on the first
line.
Fig. 2. Segment of a serpentine line.
Considering Fig. 3, the near-end crosstalk,
shown as voltage V2, is proportional to the mutual
capacitance and mutual inductance characterizing the
two lines. The coupling coefficient k^/^ is given by
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
183
^ur ~
I \Cn
4
C,.
-12
+•
^12
(1)
-^22 J
where the crosstalk at the near end is given by
y,E(')=KE[uo-vA'-2^.)]
To demonstrate these two scenarios, we use the
FDTD method to simulate the behavior of two
parallel strip lines (similar topology to Fig. 3) with
cross section shown in Fig. 4.
(2)
where Cn and Ln are the self capacitance and
self inductance, respectively, of the lines, and Cu
and L|2 arc the mutual capacitance and mutual
inductance, respectively, of the coupled lines [1]. Let
us assume that the rise time is smaller than the line
delay, tj. If a step function is launched on line 1, then
the crosstalk at the near end of line 2 will be a pulse
with a duration twice the line delay (i.e., equivalent
to the round-trip time t<j). The crosstalk at the far end
has a much smaller duration and it equals zero when
the capacitive and inductive coupling coefficients are
equal (as in the case when the medium is
homogeneous as in strip lines) [l]-[2]. Since the
pulse on the active line propagates to the right, the
voltage at the near end of the passive line is in effect
due to a leftward propagating wave. Similarly, the
voltage at the far end of the passive line is the
accumulation of a rightward propagating wave
(which is zero in the case of homogeneous medium.
'd
Fig. 3. Parallel transmission lines with matched
terminations. A pulse transmitted on the upper line
generates crosstalk on the lower line that propagates
in a manner equivalent to the case of the serpentine
line segment shown in Fig. 2.
If line 1, in Fig. 3, is excited by a pulse starting
at t= to and of duration T assumed to be shorter than '
the length of the line U, then the near-end crosstalk
will consist of two opposite polarity pulses separated
by 2trt-T,and each of duration T. This can be shown
using the principle of superposition.
By
decomposing the finite duration input pulse into two
opposite polarity step functions separated by T, the
near-end crosstalk will be the sum of two pulses. The
first pulse will have positive polarity and is excited at
time to. The second pulse will have negative polarity
and will start at time to+T. The sum of the two pulses
yields a waveform consisting of two pulses, each of
duration T and separated by a time interval of width
2td-T.
W=0.]mm, L-.4mm, H=0.4mm, €^ = 5
Fig. 4. Strip line cross section.
Figure 5 shows the near-end voltages on the
active and passive lines when the excitation is a step
function. The far-end voltage at the active line is also
shown in Fig. 5. One can observe that the duration of
the near-end crosstalk is T.
Figure 6(a) shows the near-end voltage on the
active and passive lines when the excitation is a pulse
of duration 250 picoseconds. We observe that the
near-end passive line experiences two opposite
polarity pulses, each of the same duration as the
original pulse. This performance is in full agreement
with our earlier prediction. For completion, we show
the far-end voltage on the active and passive lines in
Fig. 6(b). Notice that the bump appearing in Fig. 5 is
due primarily to the numerical dispersion of the
FDTD method and the physical dispersion caused by
the propagation in the strip line configuration.
When forming a serpentine line by closely
connecting identical transmission line sections, a
similar crosstalk mechanism to that discussed above
• takes effect, except now one should also account for
coupling between non-orthogonal lines and the
coupling between orthogonal line segments. Let us
consider the serpentine line composed of seven
transmission line segments (shown in Fig. 1).
Suppose" Vj„ is a step voltage with a magnitude of one
volt..
Once this voltage waveform is launched at the
transmitter end, a voltage of magnitude 1/2 appears
on line 1. The voltage at line 1 will induce leftward
traveling waves at the near ends of lines 2, 4 and 6.
These signals continue to travel on lines 3,5 and 7
respectively. Notice that the launched signal does not
excite any forward propagating waves on any of the
lines as per the discussion given above. (Throughout
this work, forward propagating cross talk refers to
crosstalk propagating towards the receiver. Similarly,
Ramahi: Analysis of Conventional and Novel Delay Lines
184
"backward propagating crosstalk refers to crosstalk
propagating towards the source.)
NE-activa J ;
j/
i : FE-actn«
signals arrive ahead of the main signal by 2td, 4l<i and
6td, respectively. It is important to realize that as the
main signal hops from one line to the next, the
crosstalk induced is additive. In other words, the
crosstalk noise accumulates synchronously. This
process of coupling continues until the main signal
arrives at the receiver. In fact, by adding the effect of
all induced crosstalk, the waveform obtained at the
receiver end becomes a laddering waveform, where
each of the ladder levels is of duration It^.
/
NE-pasBive
I
/.
«
time (nanosecond)
n
Fig. 5. Wave forms on the active and passive lines
due to a step function excitation. VI is the input pulse
at the active line, V2 is the near-end voltage at the
passive line, and V3 is the far-end voltage at the
active line.
Assuming that multiple coupling is negligible
(multiple coupling refers to the feedback crosstalk
between lines), equation (1) can be used to find the
magnitudes of the crosstalk induced on any line by
the main signal as it propagates on any other line.
The coupling between any two lines m and n, where
m 9t n, can be written as
I
time (nanosecond)
(a)
•
i
a
oi
t
_ M I ^mn I , ^mn
'"-'^"41 C
~
(3)
I
—
where Cmm and Lmm are the self capacitance and
self inductance, respectively, of line m, and Com and
Lnn are the mutual capacitance and mutual
inductance, respectively, of the coupled lines m and
n. Let us assume that the rise time is much smaller
than the pulse width. Let us consider the leftward
traveling wave excited at line 2, which has a
magnitude of ki/2, and duration iXt,. Furthermore, this
wave arrives at the receiver end ahead of the main
signal by IX^. Similarly, backward traveling waves
will appear on lines 4 and 6 having duration of 21^
and magnitudes k3/2 and ks/l respectively. These
signals will arrive ahead of the main signal by 4t<i and
6td, respectively.
When the naain signal starts to propagate along
line 2, it induces forward propagating crosstalk on
lines 3, 5 and 7. The magnitudes of these crosstalk
signals are ki/2, Y-Jl and Y^l, respectively. These
-_
J
^p«««-«~
time (nanosecond)
(b)
Fig. 6. Wave fomis on the active and passwe sections
of a two-line ¥ystem obtained using the FDTD
method, (a) Near-end voltage (NE-a = Near end at
active line, NE-p = Near end at passive line), (b) Farend voltage (FE-a = Far end at active line, FE-p = Far
end at passive line).
Therefore, when the seven-section serpentine
line (Fig. 1) is excited with a step function starting at
to, we can expect the receiving end waveform to have
a laddering waveform. When the excitation is a finite
duration pulse, then the shape of the received signal
will depend on whether T is shorter or longer than the
round trip 2t<i. When T is longer than 2td, the
qualitative behavior of the received signal is depicted
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
185
in Fig. 7(a), showing degradation of the pulse
magnitude once the main signal is received. On the
other hand, when T is shorter than 2td, the received
waveform consists of a train of pulses, each of
duration T and separated by 2U - T, as shown in Fig.
7(b).
Tlic wave tracing analysis adopted here can be
further generalized to serpentine lines composed of
many sections. It can be shown that for a general
serpentine delay line consisting of 2N+1 sections, a
laddering waveform is generated and it arrives before
the main signal. The laddering wave includes N
ladders; each ladder is of width 2t<i and level /^iCN-iHi
fori=ltoN.
In light of the above analysis, several
observations can be made:
1. The laddering wave creates the impression
that the pulse is arriving earlier than
intended. In reality, however, the waveform
arriving at the receiver is composed of two
signals: the crosstalk laddering waveform
(noise) and the unadulterated signal.
2. The crosstalk accumulates synchronously at
the receiving end. The synchronous
accumulation is due to the periodicity of the
serpentine line structure.
3. The highest level of the laddering waveform
can reach several times the level of crosstalk
induced between two adjacent lines. In fact,
the cumulative magnitude can be higher than
the level of the main signal.
4. The magniftide of each ladder is
proportional to the capacitive and inductive
coupling between the lines.
5. The spreading of the laddering wave, i.e.,
the width of each ladder, is directly
proportional to the length of the serpentine
line segments.
6. When the pulse duration is greater than 2td,
where t^, is the delay through -one section,
the cumulative effect is less pronotmccd,
producing ladders of short duration. The net
effect is that the pulse arrives at the
receiving end as a dispersed wave with an
apparent rise time longer than that of the
original signal. In such case, the serpentine
line can be interpreted as a low-pass filter.
7. When the pulse duration is smaller than 2t;i,
where tj, is the delay through one section,
the received waveform consists of a train of
pulses, each of duration T and separated by
2tci-T.
b. Numerical Experiments
Several assumptions were made in the wave
tracing analysis presented above. These were:
1. Negligible non-adjacent line coupling.
2. The backward propagating crosstalk was
assumed to be zero. In the general case of
non-homogeneous media, there is some
backward induced crosstalk.
3. The small transmission line segments that
connect the longer sections (orthogonal
sections) have been assumed to have zero
delay (zero physical length).
4. Negligible multi-modal propagation.
Considering the above simplifications, the
qualitative model discussed earlier serves only to
give an understanding of the primary crosstalk
contributors. For a more accurate prediction of the
performance of the serpentine line, the threedimensional FDTD method will be used since it fiilly
integrates the constraints listed above.
—
The first design tested is that of the sevensection serpentine line discussed earlier (see Fig. 1).
A cross section of two adjacent lines of this geometry
is shown in Fig. 4. The total electrical length of the
line is 114.6 mm (4.51 in). The length of each section
is 16.6 mm (0.65 in). Figure 8 shows the signal at the
receiver end. Comparison is made to the control line.
(In all the experiments discussed in this paper, the
control line is a straight line with electrical length
equivalent to that of the delay line under
examination. The length of a delay line is considered
to be the length of a line that runs through the center
of the line, including bends.)
Observing Fig. 8, we notice that the wave
tracing analysis predicted the ladders that precede the
main signal only. While such prediction is quite
important, it did_not account for the unexpected high
signal level received after the arrival of the main
signal (overshoot). Although the high level occurs
beyond the logic-switching instant, nevertheless, it
can have an adverse effect when-the signal switches .
to zero as false logic could be triggered.
Ramahi: Analysis of Conventional and Novel Delay Lines
1S6
longer than Case A. However, this difference is of
minor importance since the low-to-high switching
occurs at approximately the same time.
Received
waveform
tune
i
Ojtputfiom
neiative polarity
step ftmctiDO
IS
I
(a)
time (nanosecond)
Fig. 8. Received wave form for the seven-section
serpentine line.
tune
10 mm
Ouptcfrom
neptive polarity •
step functnn
0>)
Fig. 7. Wave form at the receiving end of a
serpentine line due to a pulse with a duration of T. (a)
T > 2td. (b) T < 2td.
For the second test, we consider a more realistic
serpentine line, where the total electrical length of the
line is. 405.4 mm (15.96 in). Two variations of this
line are considered, as shown in Figs. 9 and 10, and
will be referred to as Case A and Case B,
respectively. The difference between the two
tepologies considered is that one has shorter sections
than the other. Notice that in both cases, the lines
have equivalent length, have identical line separation
between adjacent segments and both topologies
occupy equal board area. The excitation waveform is
a pulse having a finite duration of 1.4 nanoseconds
and rise time of 100 picoseconds.
Figure 11 shows the response of these two lines
in comparison to the reference line. The serpentine
line designated as Case B has longer sections and,
consequently, its receiver signal had ladders that are
Fig.-9. 19-section serpentine line.
187
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
rAA/V-i
the transmission line density (or etch density) to
remain unchanged, the separation between the lines
must not be changed. An alternate design, which
forces the crosstalk to accumulate asynchronously,
was proposed in [3] and is shown in Fig. 12. This
alternate design, the flat spiral delay line, minimizes
the periodicity in the transmission line routing as
much as practicable.
tOmm
B
Ro
10 mm
*
us
. LS
L26
:.
L2I
U
L9
L22
W
LU
Ll6
LIB
Fig. 10. 17-section serpentine line.
L2
LI?
U3
U2
22.2
L6
LU
L20
Ljsue
Lll
.
. '
1
U9
U
-.
refcicncc
•
Case A
- - - Case B
US
.< i'
u
.«»- 4-;-_'
.
•• •
•
'• : . '•
\
TJ 03 4
8. "
'
1
\
3 Dl
. .'• '
,
k
f 1 .«.«
\*. ^
.
Ikne (norosecond)
Fig. 11. Response of the 15.96 in long serpentine
lines in comparison to the reference line.
HI. The Spiral Delay Line
In the serpentine line, the crosstalk was found to
accumulate synchronously. This accumulation can be
significant enough to trigger false logic. Ideally, this
crosstalk needs to be eliminated; however, since the
crosstalk is inversely proportional to the separation
between the lines, the only way to eliminate or
reduce the crosstalk would be to increase the
separation between the lines and/or to increase the
coupling between the line and the reference plane.
Increasing the separation of the segments
unfortunately requires larger circuit board area,
which can be either expensive, or impossible in light
of the density requirements. Therefore, given that for
Fig. 12. Spiral delay line.
The most prominent feature of the spiral delay
line is the spreading, over time, of the crosstalk noise.
To demonstrate how the coupling mechanism works
in the spiral line, we will initially assume that
coupling between non-adjacent lines is negligible.
Let us consider the line shown in Fig. 12, and let us
assume that the signal arrives at point A at time t=0.
Once the signal arrives at A, forward propagating
crosstalk- is induced at point B.
This crosstalk
arrives before the main signal by the time delay"
needed to propagate through the entire length of the
line minus the segment L26. In other words, the
induced forward propagating crosstalk at B arrives
only after the time needed to travel the segment L26.
(Henceforth, the time delay through segment LN will
be denoted by ttN)- So the first crosstalk induced at B
by the arrival of the signal at A arrives at t=tu6When the main signal arrives at point C, forward
propagating crosstalk is induced at D, arriving at the
receiver at approximately t=t[,i+tL25+tL26Similarly, when the main signal arrives at E, it
induces at point F a forward crosstalk which arrives
Ramahi: Analysis of Conventional and Novel Delay Lines
188
at the receiver at t= ti+ tL2+ tL26+ ti^ tL24- Next, the
signal arrives at G, inducing forward propagating
crosstalk at H, which arrives in t= ti,i+ tL2+ tL3+ tL26+
tL25+ tL24 + tL23. Howcver, the signal arriving at H
will also induce forward propagating crosstalk at I.
The crosstalk induced at I arrives at the receiver at t=
tLi+ tL2+ tL3. As the signal propagates further towards
the receiver, it continues to induce crosstalk in the
adjacent lines.
The above wave-tracing analysis ignores
coupling between non-adjacent lines. However, even
if this coupling is included, it is easy to see that the
crosstalk noise will be distributive and will not
accumulate S5mchronously. When the voltage
excitation is a pulse of finite duration, the opposite
polarity crosstalk induced by the trailing edge of the
excitation pulse will help in reducing the final
cumulative crosstalk.
To demonstrate performance of the flat spiral
delay line, we construct a line with electrical length,
line separation, and total board area, all equal to the
serpentine lines shown in Figs. 9 and 10. The crosssection parameters are given in Fig. 4, and the rise
time and pulse duration are as before (pulse width of
1.4 nanoseconds and rise time of 100 picoseconds).
The dimensions of this spiral line are shown in Fig. 12.
Figure 13 shows the signal at the receiving end
of the spiral line in comparison with the control line.
The outstanding performance of the spiral line is
clearly visible. We observe that the crosstalk noise
was spread over time ahead of the main signal
resulting in a high fidelity signal. In fact, for the case
considered, we notice that the maximum crosstalk
stays at or below 10% of the signal level, thus
considerably reducing the potential for triggering
false logic. We note here that slight shift between the
signal due to the control line and the unadulterated
signal due to the spiral Hrie is due to -a slighrdifference in line lengths.
The singular feature of the spiral line is its
ability to spread.the noise over a larger duration. A
key observation is that the signal arriving at the
receiving end is composed of the unadulterated pulse
in addition to the asynchronous noise, This implies
that the spiral line can lead to an almost perfect
isolation of the accumulative noise from the
unadulterated pulse.
Mm* (nanosteonif)
Fig. 13. Response of the 15.96 in long spiral line in
comparison to the reference line.
IV. The Concentric Cs Delay Line
From the analysis and simulation of the spiral
line, the concept of asynchronous coupling is
emerging as a powerful mechanism to reduce the
apparent dispersion of a delay line. This section
introduces a novel design that exploits this
mechanism of asynchronous coupling. The new delay
line, shown in Fig. 14, will be referred to as the
concentric Cs delay line. (The name concentric Cs
was chosen since this new transmission line
resembles different lines shaped as the letter C and
centered at one location.) The design of the new line
was motivated by the theme of reducing the
periodicity in the topology of the structure in order to
minimize the accumulative or synchronous coupling.
Using wave-tracing analysis, one can
qualitatively predict the performance of the Cs line.
Similar to the spiral line, the Cs line distributes the
coupling over time. Consider a pulse entering the line
at the beginning of section LI (see Fig. J4). This
- ^ulse will induce a crosstalk that will precede the
primary (or unadulterated) signal by a time period
corresponding to the propagation over segments LI,
L2, L3, L4, L5 and L6. Notice that in contrast to the
spiral line, this (first-appearing) crosstalk will
precede the unadulterated signal by a relatively short
period, however, the topology of the Cs line prevents
synchronous accumulation. One can qualitatively
describe the periodicity in the Cs line structure as
lying in between that of the serpentine line and the
spiral line. The relative advantage of the mid-level
periodicity in the structure becomes more apparent
when observing the full-wave behavior as shown
below.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
189
To demonstrate the behavior of the Cs hne, we
design a line, equivalent in electrical length, line
spacing and total board area, to the spiral and
serpentine lines studied earlier. The dimensions of
this design are shown in Fig. 14. Figure 15 presents
the result of the FDTD simulation, under previous
excitation conditions, showing the signal at the
receiver end.
6
"^^^,,^
20 mm
In comparison to the spiral line, the Cs line
keeps the noise closer to the main signal, whereas the
spiral line distributes the noise over longer time
duration. From this perspective, the spiral line is
superior to the Cs line. However, if one were to
consider an excitation composed of a train of pulses,
as would be the case in practical scenarios, it might
be more advantageous to force the crosstalk noise to
accumulate closer to the pulse that generates it and
not interfere with other pulses.
V. Summary
u
L2 22.2 mm
LS
L4
U
Fig. 14. The concentric Cs delay line.
Comparison is made to the spiral, serpentine
and reference lines. Clearly visible from Fig. 15 is the
concentration, or accumulation of the noise in the
close proximity of the unadulterated signal. We
further observe that, for the particular topology and
line dimensions considered, that the receiver signal
remains at or below 35% of the peak amplitude of the
transmitted pulse.
time (nanosecond)
Fig. 15. Response of the concentric Cs line in
comparison lo the serpentine and spiral lines of equal
length.
This paper presented a qualitative analysis of
serpentine and flat spiral delay lines based on the
simply, yet powerful, ray tracing technique. The
three-dimensional FDTD method was used to predict
the full-wave performanc.e of these delay lines, thus
accounting for higher-order modes, multiple and
feedback coupling and the effect of orthogonal
segments. Despite the strength and completeness of
the FDTD method, the ray tracing method facilitated
an understanding of the coupling mechanism in
complex-shaped delay lines and lead to the
introduction of novel delay line such as the flat spiral
line. A new line, the concentric Cs delay line was
introduced based on the concept of reducing the
periodicity in the structural topology, thus forcing the
crosstalk noise to accumulate asynchronously.
Numerical simulation using the three-dimensional
full-wave FDTD method showed that the new
designs result in receiver waveform that is less
susceptible to triggering false logic than in the case of
the serpentine line.
References
[1] H. B. Bakoglu, Circuits, Interconnections, and
Packaging for VLSI, Addison Wesley, Reading,
MA. 1990.
[2] R. B. Wu and F.L. Chao, "Laddering wave in
serpentine delay-line," IEEE Trans. Components
Packaging Manufacturing
Tech.
Part-B
Advanced Packaging, Vol. 18, No. 4, pp. 644650, Nov. 1995.
[3] R. B. Wu and F.L. Chao, "Flat spiral delay line
design with minimum crosstalk penalty," IEEE
Trans. Components Packaging Manufacturing
Tech. Part-B Advanced Packaging, Vol. 19, No.
2,pp. 397-402, May 1996.
[4] O. M. Ramahi, "FDTD analysis of conventional
and novel delay lines," Presented at the 16"'
Annual Review of Progress in Applied
Computational
Electromagnetics
meeting,
Monterey, CA, March 20-25,2000.
[5] N. Orianovic, R. Raghuram and N. Matsui,
"Characterization of microstrip meanders in PCB
Ramahi: Analysis of Conventional and Novel Delay Lines
interconnects," Electronic Components and
Technology Conference, Conference Proc, May
21-24,2000.
[6] B. J. Rubin and B. Singh, "Study of meander line
delay in circuit boards," IEEE Trans. Microwave
Theory Tech., Vol. 48, No. 9, pp. 1452-1460,
Sept. 2000.
Omar
M.
Ramahi
received the BS degrees in
Mathematics, and Electrical
and Computer Engineering
(highest honors) from
Oregon State University,
Corvallis, OR in 1984. He
received his M.S. and
Ph.D. in Electrical and
Computer Engineering in
1986
and
1990,
respectively from the University of Illinois at
Urbana-Champaign. From 1990-1993, Dr. Ramahi
held a visiting fellowship position at the University
of Illinois at Urbana-Champaign. From 1993 to 2000,
he worked at Digital Equipment Corporation and then
Compaq (presently, Hewlett Packard), where he was
member of the Alpha Servers Product Development
Group. In August of 2000, he joined the faculty of
the A. James Clark School of Engineering at the
University of Maryland at College Park, as an
Assistant Professor in Mechanical Engineering,
Assistant Affiliate Professor of Electrical and
Computer Engineering, and a Faculty Member of
CALCE Electronic Products and Systems Center.
Dr. Ramahi served as a consultant to several
companies. He was instrumental in developing
computational techniques to solve a wide range of
electromagnetic radiation problems in the fields of
antennas, high-speed devices and circuits and
EMI/EMC. His interests include high-speed devices,
packaging, and interconnects, experimental and
computational EMI/EMC studies, biomedical
applications of electromagnetics, novel optimization
techniques, RF MEMS, interdisciplinary studies
involving interaction of novel materials and
electromagnetic energy, and high-impedance
electromagnetic bandgap surfaces. He has authored
and co-authored over 120 journal and conference
papers and presentations, and has given numerous
short courses to the industry. He is a co-author of the
book EM/EMC Computational Modeling Handbook
(Kluwer Academic, 2"^ Ed., 2001). Dr. Ramahi is a
member of Eta Kappa Nu, Tau Beta Pi, Senior
Member of IEEE and a member of the
Electromagnetics Academy.
190
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
191
An Empirical Approach for Design of
Wideband, Probe-Fed, U-Slot Microstrip
Patch Antennas on Single-layer, Infinite,
Grounded Substrates
V. Natarajan and D. Chatterjee
Abstract—Wideband microstrip antennas willi relatively
simple topologies continue to attract attention for design of
compact, higli-pcrformancc communication systems. The
coaxially-fcd, rectangular patch U-slot has recently been
investigated numerically and experimentally, and shown to
yield 10 dB return-loss bandwidlhs in excess of 20%. However there are no analytical models, nor any systematic design procedures currently available that can aid realizing
these conBgurations. To that end, based on extensive CAD
simulation results for a wide range'of commercially available microwave substrates (e, = 2.94 to 10.2), an empirical
design methodology is derived and illustrated by examples.
It is shown that the present empirical design technique,
with its attendant limitations, generate wideband U-Slot
designs that are optimized using CAD tools such as IE3D
within a few iterations, resulting in substantially reduced
overall process cycles.
I. INTRODUCTION
The theory and design of a wide variety of probefed, microstrip patch antennas for various applications
has been well documented [1]. For portable phone systems there is a need for low-profile (embedded) antennas
with a 10 dB return-loss bandwidth > 10% in addition
to other desirable electrical characteristics [2, p. 312].
(The return loss, for this paper, is taken as -201ogio |r|
in dB, where F is the reflection coefficient.) As found
in [I, ch. 9], studies on wideband microstrip antenna
designs emphasize techniques such as multi-layer substrates, parasitic elements and apcrture-ccnpled excitations. However such approaches obviate the realization
of low-profile, compact antenna topologies, or may complicate the fabrication process due to the need for sophisticated feed clement design(s) [3].
Interestingly the U-sIot antenna, first reported in [4],
was a new form in ultrawideband microstrip antenna design since it could generate « 40% bandwidth by maintaining very simple feed and patch designs on singlelayer foam substrates. For wideband applications it thus
appears that the U-slot design is a pioneering concept as
it is indeed a very formidable alternative to the existing
Manuscript received February 18,2003; revised June 4,2003
Division of Computer Science and Electrical Engineering, 570-F
Flariihcim Hall, University of Missouri Kansas City (UMK.C), 5100
Rockbill Road, KC, MO 64110, USA; e-mail:chald(Siunikc.edu
wideband patch topologies [1],[3]. The subject of this
paper is to further explore some advancements in design
of wideband, probe-fed, U-slot patches on single-layer
substrates.
Most results for this novel design are available for air
(fr = 1) or foam (tr »» 1) substrates [4]-[7]. One finds
appropriate results fore, = 2.33 in [8], inclusive of finite ground-plane and substrate truncation effects. Currently there are no analytical models, nor empirical design relations available to initiate an U-slot design from
some nominal specifications. In [6] an attempt has been
made to use a [Y]3x3 matrix representation for the Uslot, but the analytical expressions for the diagonal elements Yii,22,33 are unavailable. (It appears from [6]
that their determination was done using iterative experimental techniques). Explicit formulas for the two resonant frequencies of the U-slot, with validation results
for foam substrate are available in [9],[10] that are good
only as important checks in a simulation.
Reduced U-slot patch size topologies have been reported in [11]-[14] with an average bandwidth of ~
30%. Furthermore, 10 dB return loss bandwidths of
44% and 50% have been reported in [15] and [16], respectively. Dual-band designs have been reported in [17]
for substrates with permittivity Cr ~ 4.4, and in [18]
wideband patch designs on multi-layer substrates have
been reported. Results for U-Slot performance on jnicrowave substrates have been summarized in [19], but
without any design information. Recently, a different
design procedure for U-Slot on single-layer microwave
substrates has been reported in [20]. However, there are
significant differences between [20], and the proposed
approach in this paper. These differences will be identified later in this paper.
The information gleaned from [21] suggests avoiding
use of moment-method [22],[23] based CAD tools like
IE3D [24] for initiating single-element patch designs,
due to prohibitively high computational cost. Thus efficient design processes, despite their heuristic/empirical
nature, can still help the overall simulation cycle be costeffective. Normally such empirical procedures reduce
substantial savings in computational resources by requiring fewer iterations in the final CAD optimization of the
antenna topology [21]. Since no formal, systematic pro-
1054-4887©2003ACES
Natarajan and Chatterjee: Design of Wideband U-Slot Microstrip Patch Antennas
cedures are currently available for the design of U-Slot,
the purpose of this investigation is to develop guidelines
and present empirical formulas to aid in realizing such
goals.
The empirical formulas developed in this paper apply
to probe-fed designs on single-layer grounded substrates
that are (ideally) infinite in extent. In addition, empirical design formulas, obtained via moment-method based
parametric simulations, apply to U-slots patch designs
with definite geometrical symmetry as elaborated later in
this paper. The contents of the paper are outlined next.
To that end, following [25]-[27], validation results for
IE3D against appropriate measured data for microwave
substrates (Cr = 2.33) from [8], are included. The
IE3D code validation results are followed by a careful analysis of the various U-Slot designs studied earlier
[8],[28], resulting in various dimensional invariance relationships that are crucial for U-Slot design on microwave
(and foam/air) substrates. Selected results demonstrating the validity of the empirical formulas are included
from [29],[30]. Finally the major observations are summarized, and a list of relevant references is included.
II. PROBLEM DESCRIPTION
In Fig. 1 the geometry ofa probe-fed, U-Slot patch on
a single-layer substrate is shown with all the dimensions
indicated therein. This topology is a simple modification
to a probe-fed rectangular patch antenna, the latter being
generally a narrowband radiating element [1].
w
f
1
1
T
!
-• It.**.
B t
il
Metallic
/
1 m1
1
t>l ^
w.
F
X
L
DielecMc ReQJon
Raifalina
Patch
Cr .tan5
/
PEC Ground
V"
'CosdalFeed
Fig. 1. Physical topology of a coaxially-fed, single-layer, rectangular
patch U-Slot microstrip antenna
In absence of any analytical, i.e., cavity or transmission line models [1, chs. 4,5], one can still investigate the effects on a performance characteristic (such
as impedance, gain, etc.) due to variations in substrate/patch geometry via careful measurements or rigorous, fiill-vrave CAD simulations [21]. Development of
rapidly iterative design procedures could involve heuristic/empirical approaches, subject to further refinements
via CAD optimizations [21].
192
For the investigation reported here, the main aim is
to examine how parameters such as substrate thickness,
overall patch dimensions, slot width, probe location and
radius, as shown in Fig. 1, affect the wideband performance. The generic nature of the impedance chaiacter-
Fig. 2. Typical impedance loci characteristics. The performance of
any wideband design (or modifications) is desired such that the loop
of the impedance loci 1,2 & 3 encircle die center (VSWR = 1) of the
Smith Chart as in locus 4.
istics is shown via a Smith Chart in Fig. 2. The desired
characteristic is depicted in # 4. The wideband behavior
of the antenna will be superior if the size of the loop for
the impedance locus # 4 shrinks to the VSWR = 1 point
on tbe Smith Chart. (In addition, most of the frequencies
of interest has to lie on that loop.) For practical applications, the size and location of the loop in an impedance
loci should be such that the VSWR < 2, correq)onding to
a return loss of 10 dB.
Generally, the impedance loci will be far removed
from the desired behavior shown in # 4, i.e., it will be
more like the loci # 1, 2 or 3. The wideband problem
then reduces to the study of how the changes in various
dimensions in Fig. 1 could transform loci # 1, 2 & 3
such that a loop could be-obtained in the impedance loci
meeting the criterion VSWR < 2.
It is po'Ssible to analyze the impedance behavior using
analytical (cavity model) expressions, to a very good degree of accuracy [1, chs. 4 and 5]. Such would facilitate
rapid parametric simulations, prior to any computationally intensive, fiill-vrave analysis. Since there exists no
such analytical model for U-Slots, having recourse to an
alternate route for rapid parametric studies appears critical before any CAD-based optimization. To that end, it
is necessary to examine the capability of the IE3D code
for numerical modeling of U-Slot geometries. The results for the IE3D code validation are shown next.
III. VALIDATION RESULTS FOR THE IE3D CODE [24]
As mentioned in [25]-[27], a CAD tool needs to be
validated for an ensemble of appropriate test cases that
are closest to the topology being studied. Furthermore,
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
193
since code validation is an open-ended process, a judicious selection of the test cases forms an important part
of such investigation. To that end, it was decided to examine the capabilities of the IE3D code against the measured radiation pattern data (in 4 = 0°, 90° planes) for
U-Slots fabricated on substrate e^ = 2.33 as given in [8,
Fig. 5] corresponding to a frequency of 3.56 GHz. To
the best of the knowledge of the present investigators,
reference [8] is the only source for which measured and
computed data are available for U-SIots on finite ground
planes for microwave substrates. (Most of the data, as
mentioned earlier, is for foam (or air) substrates for USlot topologies.)
The dimensions of the antenna B as in [8, Table I]
arc: W = 3.6cms,L = 2.6cms,W, = 1.4cms,L, =
1.8cms,b = 0.4cms,t = 0.2cms,Xp = 0.0 cms, Yp =
0.0 cms, and F = 1.3 cms, referring to Fig. 1 herein. The
radius of the probe could not be found in [8]; after several trials, it was found that dpr„i„ = 0.127 cms in Fig.
1 provided the best agreement with the data in [8]. The
Fig. 3. Measured e - e - data [8. Fig. 5); computed -D - O(IE3D8.0)and-< - . - (1E3D 9.0), for total fields in 0 = 0° plane
IE3D results in Figs. 3 and 4 arc for a 12 x 10 cms rectangular ground plane. The agreement between measured
[8] and simulated results are reasonably acceptable at all
angular regions except near S -* 90°, 270°. The reason(s) for the discrepancies are explained below.
The actual topology analyzed in [8, Fig. 1] was a
U-Slot located on truncated, rectarigular substrate on a
finite rectangular ground plane. In the IE3D simulations, the radiating patch was located on an infinite substrate backed by a finite rectangular ground plane of the
same dimensions. This model cannot account for the
surface wave diffraction by the trancated dielectric substrate. Since surface waves arc dominant near the airsubstrate interface, (i.e. 9 -+ 90°, 270°), the radiation
behavior is not accurately predicted near this region as
seen in Figs. 3 and 4. This present limitation in the
Fig. 4. Measured - • e- e - data [8, Fig. 5]; computed -D - D(IE3D 8.0) and - < — * - {IE3D 9.0), total fields in ^ = 90° plane.
IE3D code [24] is due to lack of implementation of the
appropriate microstrip Green's fimction that can account
for diffractions due to dielectric and ground-plane truncations. In [8], since the FDTD technique was used to
calculate the radiation pattern, such edge diffractions are
considered in situ. Thus, computed patterns in [8] do not
show any ficticious discontinuities near d -> 90°, 270°,
unlike the IE3D results presented here.
The foregoing results demonstrate the limitations of
the IE3D code when applied to modeling of antennas
on grounded, truncated substrates. However the 2:1
VSWR bandwidth of antenna B, computed via IE3D,
was around 24% - in good agreement with measurements
[8, Table II]. Also, the IE3D code had been used to
replicate the results for microstrip antennas on infinite,
grounded substrates with the test cases chosen from various topologies in [1). In all these cases the agreements
were very good with published data. Since the scope of
this present investigation is limited lo infinite, grounded
substrales-the type of discrepanciesTn Figs. 3 and 4 are
not likely to affect the results.
IV. DIMENSIONAL INVARIANCE IN U-SLOT DESIGN
The key to the development of the empirical design
procedure is the establishment of the dimensional invariance of the U-Slot studied in [8] and [28]. These results
are summarized below in table I from [29],[30]. In table
I one finds that the only parameter which changes with
substrate Cr is ^, and all other dimensional ratios remain almost invariant. Consequently, to design a U-Slot
on an infinite, grounded microwave substrate the determination of ^ for a specific substrate (er and h) and
resonant/design frequency, fr, is the key step. One can
then use the information in table I to derive the topology
of the patch as shown in Fig. 1. From columns 3,4 and
7 in table I one can easily deduce that ^ ss 1.38. Inter-
Natarajan and Chatterjee: Design of Wideband U-Slot Microstrip Patch Antennas
194
TABLE I
DIMENSIONAL INVARIANCE IN U-SLOT DESIGNS.
Cr
I
^
^
¥
w:
W
1.0 8.168^ 1.515 0.835 4.237 0.13 3.203
2.33 4.49
1.445 0.777
4.5
0.144 2.573
4.0
3.87
1.443 0.776 4.51
0.144 2.573
9.8
2.87
1.442 0.777 4.48 0.144 2.574
2.33 5.624 1.444 0.777
4.5
0.143 2.571
The dat > have been obtained from 28] and 8], for the first
four and last row, respectively. For all the cases cited here,
the minimum and maximum bandwidths were 15% and
42%, respectively. The data for e, = 2.33 in the second
and fifth rows refer to U-Slot topologies from [28] and [8],
corresponding to 900 MHz and 3.26 GHz, respectively.
estingly, this fact appears to have been confirmed for the
U-Slot data presented in [19, table lb].
At this stage it is important to distinguish between the
approach in [20] and this paper. It is noted that [20], like
this paper, doesn't contain any full-wave mathematical
analysis for U-Slot. One of the main differences, in context of table I, is the determination of U-SIot dimensions
[20, sec. III]. The underlying assumption in [20, sec.
Ill] is the existence of four different resonant frequencies of the impedance loop (as shown in locus # 4 in Fig.
2) for the U-Slot. For high e, substrates such multiple
resonances may not occur, but an impedance loop could
still form, as shown in Fig.2, away from the zero reactance (jX = 0) line on the Smith Chart. Apparently, this
restricts the technique in [20, sec. Ill] primarily to low
Er substrates. In contrast, the dimensional invariances
(table I) apply to low, medium and high permittivity substrates.
Since it is important how the various dimensions in
Fig. 1 could affect the bandwidth, a detailed study was
undertaken to examine such effects, for which salient
features are shown in the following section.
Fig. 5. Effect of probe location on the impedance behavior of U-Slol:
«, = 4.5,tan* = 0.002, W= 4.89, L= 3.538, h= 1.27, L, =
2.45, W, = 1.9,1= 0.274, a= 0.777, b= 0.311, (i.e.,* = 2.5),
dprot. = 0.127 and Xp = 0.0 - all in cms;
(Y, = 0.4
cms), + - -I- - + (Yp = 0.2 cms), and o - o - o (Yp = -0.2 cms),
refciring to the dimensions shown in Fig. 1.
V. PARAMETRIC MODELING STUDIES VIA IE3D CODE [24]
The primaiy objective of the parametric simulations
is to examine the nature of the input impedance variations as shown in Fig! 2 in section IL Assuming that
the initial topology of the U-Slot has been designed using the infonnation in section IV, it is still possible that
the desired bandwidth may not have been achieved. This
implies that the initial design needs further optimization,
which in view of Fig. 2 implies that the impedance loop
should be shrunk to encircle the vicinity of the center of
the Smith Chart. The parameters that exercise significant
control on the impedance loop size and location are critical to the optimization process. Results for low and high
permittivity substrates are available in [29], and only selected results for e, = 4.5 are shown here from Figs.
j . n Ti. j .
• 1 j j • it • j. .j . o
5 to 9. The data are mcluded m the mdtvidual figures
captions, and hence are not repeated here to avoid tedium. In Fig. 5, Yp = -0.2 and 0.4 cms refer to probe
locations below and above the origin of the coordinate
Fig. 6. Effect of probe radius on the impedance behavior of U-Slot:
(r = 4.5,tan* = 0.002, W= 4.89, L= 3.538, h= 1.27, L, =
2.45, W. = 1.9, t= 0.274, a=: 0.777. b= 0,311, (i.«.,| = 2.5),
Xp =Yp = 0.0 - all in cms; o - o - o (d„„ij = 0.08636 cms),
(dp,„i, = o.i27 cms), and H- - + -+(dp„i, = 0.2 cms).
referring to the dimensions shown in Fig. I.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
195
system, respectively, with Xp = 0, as shown in Fig. 1.
The result indicates the trend that as the probe is moved
away from the edge of the slot, the impedance loop becomes more inductive and its size decreases. Similar
trends were observed for other substrate cases in [29].
Fig. 6 shows the effects of the probe radius (=
Idproir) on the input impedance. In contrast to the result in Fig. 5, variations in probe radius doesn't shrink or
expand the size of the loop. The comparison fiarther indicates that the dominant effect of the probe on the U-Slot
input impedance is determined by its location, and not
radius. Control of the probe radius could thus be viewed
as resulting in a 'fine-tuning' mechanism in order to obtain the desired wideband behavior.
Fig. 8. Effects of slot width, t, on the impedance behavior of USlot: tr = 4.5, tan i! = 0.002, W= 4.89, L= 3.538, Lj = 2.45,
W, = 1.9, a= 0.777, b= 0.311, (i.e.,| = 2.5), dp^te = 0.127,
Xp =Yp = 0.0-all in cms;
(t= 0.2 cms), o-o-o(i= 0.274
cms), H
1
1-0= 0.35 cms) referring to the dimensions shown in
Fig. 1.
-1.0
Fig. 7. EfTccts of substrate thickness on the iinpedance behavior of
U-Slol: c, = 4.5, tan i = 0.002, W= 4.89, L= 3.538, L, =
2.45, W, = l.D, 1= 0.274, a= 0.777, b= 0.311, (i.f.,J = 2.5),
dprotc = 0.127, X, =Y, = 0.0 - all in cms; o - o - o (h= 0.0
eras),
(h= 1.27 cms), and H---f--Kh= 1.5 cms), referring
10 the dimensions shown in Fig. 1.
In Fig. 7 the trends in impedance behavior with increase in substrate thickness, h, are shown. As the thickness increases from 0.6 cms to 1.5 cms, the impedance
loop decreases in size and becomes more capacitive in
character. Forh = 1.27 cms, a loop in the impedance behavior is formed closest to the center of the Smith Chart
- indicative of wideband behavior.
Variation in the U-slot width, t, results in impedance
changes (Fig. 8) similar to the one observed for probelocation variations in Fig. 5. As the U-slot width increases from 0.2 cms to 0.35 cms, the impedance loop
changes from being inductive to capacitive, and for a
slot-width t = 0.274 cms the loop is located close to
the center of the Smith Chart.
Increase in the j ratio from 0.5 to 4.5 doesn't cause
any significant change in the location of the impedance
loop, but results in shrinking of its size. This can be
Fig. 9. Effects of f ratio on the impedance behavior of U-Slot: «, =
4.5, tan J = 0.002, W= 4.89, L= 3.538,L, = 2.45, Wj = 1.9,
t= 0.274, dp^ie = 0.127, Xp =Yp = 0.0 - all in cms;
(a = 4.5), o - o - o (a = 1.0), -)- - -I- - H-Cf = 0.5), referring to
the dimensions shown in Fig. 1.
Natarajan and Chatterjee: Design of Wideband U-Slot Microstrip Patch Antennas
inferred from Fig. 9.
The information gleaned from the most important
parametric simulation results, as shown in Figs. 5 to 9,
suggests the following optimization guidelines for wideband U-Slot design:
(a) change the slot width, t, probe location Yp, and
substrate thickness h, such that the impedance
loop encircles the region in the close vicinity of
the center of the Smith Chart;
(b) following step (a), if the size of the loop is undesirably large or small, increase or decrease the f
ratio to reduce the loop size without affecting the
location of the loop to achieve larger bandwidth
(c) one may, optionally, change the probe radius to
move the impedance loop such that it encircles, or
is close to the Smith center (VSWR = 1)
The preceding simulation results were obtained for configurations where the U-slot is symmetric about the x
axis and the probe is located such that X, = 0, as
shown in Fig. 1. Again, it is important to distinguish
between [20] in contextof the parametric simulations in
the present investigation. This will be followed by the
last part, i.e., development of the empirical design equations, in section VI.
The parametric simulations in [20, sec. 11] focussed
mainly on the variation in resonant frequencies. (The USlot geometry studied in [4] was for air €r = 1 which
was scaled to Cr = 2.2 in [20]). The technique in [20,
sec. ni] is based on the availability of limited data. Furthermore, as stated in [20, sec. V], the effects of substrate and (probe) feed were not investigated in detail.
The information gleaned from the overall comparisons
between [20] and the present paper suggests that the results included here have broader scope of applicability
compared to [20].
196
for each case were noted. (If there were several resonant frequencies in the 2:1 VSWR range, an average estimate of /r was taken [29].) Each discrete pair of ^
and fr values, corresponding to an individual design,
were plotted. The MATLAB software (version 6.1, release 12.1) was used to obtain a quadratic relation ('best
fit') for these ^ vs. /, plots. For a any specific e,,
and various h, several such equations were obtained, as
shown in table 11. In this process, it was observed that
a > 20% fractional bandwidth for the 2:1 VSWR range
was obtained for those designs obeying 3.5 <x^ 5-5The same phenomenon also corresponded to the range
0.13 < i^ < 0.18 for the U-Slot topologies designed
and simulated on Cr = 2.94,4.5, and 10.2. (Here A corresponded to the resonant frequency /,, determined from
the Smith Chart.) The details can be found in [29],[30]
and are omitted here for brevity.
TABLE 11
EMPIRICAL (QUADRATIC) EQUATIONS FOR DESIGN OF U-SLOT
h
(cms)
0.635
£r = 2.94
Er = 4.5
«r = 10.2
f =0.088/? 1 = 0.063/?
-1.5/r -1- 8.5 -1.7/.-1-9.1
1.0
^ = 0.19/? f = -0.78/?
-2.1/^-1-7.8 -0.3/.-1-8.8
1.216
f = 0.085/? ^ = -0.21/?
-1.7/r-f6.8 -2.2/.-1-8.1
1.80
^ = 0.89/? f = -0.89/?
f = 1.8/'?
-8.6 fr +13 -4.8/,.-f 8.5 -|-0.11/r-^4.5
Here f, is the design resonant frequency in GHz, and h and e,
are the substrate thickness and permittivites, respectively.
-1.8/r -1- 10
^ = 0.2/if
-2.4/r -1- 9.2
X = 0-62/r^
-4.8/^-1-12
The next section illustrates the complete empirical design
procedure with simulation results for a topology
VI. DEVELOPMENT OF EMPIRICAL DESIGN FORMULAS
(Cr = 3.27) for which the empirical equations are not
In view of the observations on dimensional invariance, available in table II.
as presented in table I, the important factor that initiates
the U-slot design is a knowledge of the ^ ratio from an VII. EMPIRICAL DESIGN TECHNIQUE FoRdJ-SLOT
a-priori knowledge of some nominal specifications.
In this sectionra systematic empirical design proceConsequently, it was decided to examine the relation- dure for design of U-Slot patch antennas on microwave
ships between resonant frequency, fr, substrate parame- substrates is presented from [29]. It is shown, from the
ters Br and h, and, the larger dimension, W, of the U-Slot. results in sees. IV, V and VI, that U-Slot patch anTo that end, following the data in table I, it was decided tennas can be realized which are further optimized usto vary ^ between 2 and 7, in increments of 0.5. A typ- ing IE3D CAD software [24]. VSWR and boresi^t
ical substrate was chosen and from the pre-selected ^ (<t> = 0°,e = 0°) Gain results for the unoptimized and
values, the U-Slot dimensions were found with the aid of optimizedXJ-Slot topologies are included to demonstrate
table I. Following this procedure, various U-Slots were the efficacy of the empirical design procedure.
designed for a wide class of practical substrates availThe limitation of this design procedure, as mentioned
able from Rogers Corp. (We must emphasize that at this before, is that the U-Slot, as shown in Fig. 1, is located
stage the resonant frequency, /,, is unknown and was symmetrically w.r.t coordinate axes with the probe is on
determined as described below.)
the y axis. (For rapid automated calculations, the deThese U-Slot topologies were then characterized by sign procedure can easily be adapted within a computer
the full-wave CAD tool IE3D [24]. The resonant fre- design/simulation code.) It is assumed that the U-slot
quency/r, defined by zero reactance on the Smith Chart, antenna should have a 10 dB return loss bandwidth of
and the corresponding fractional 2:1 VSWR bandwidths > 20%, after final optimization. Another limitation is
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
197
that the dimension a = 6 in this design approach, and
the slot width, t, remains uniform.
(1) From the nominal a-priori specifications for resonant frequency fr, select a commercially available
substrate with fr. thickness h to satisfy the criterion 0.13 < ^^ < 0.18. In the experience of the
present authors the upper and lower limiting values
should be used for low- and high-permittivity substrates, respectively. For intermediate permittivity
0.15
4.5) the criterion
substrates ((r
can be used.
(2) Employ the empirical equations from table 11 to
calculate the ^ ratio, and from step (1), one can
subsequently determine the overall width W. (One
may check forthe additional criterion 3.5 < x 5.5, upon calculation of ^ ratio.)
(3) From table I, one uses ^ « 2.57 to determine
(4) From table I, since ^f » 0.777, calculate L, with
the knowledge of W'j'from step (3).
s 0.144 in table I, calculate
(5) From the relation
the slot width t, with the knowledge of W, from
(3). (This assumes a slot oiuniform width.)
(6) Similarly, from table I, via the relation ^ a 4.5,
calculate b with a knowledge of i, from (4).
(7) Assume rt = 6 in Fig. 1, and calculate L = Ls +
a + b = L, + 2a = Ls + 2b.
(8) Locate the coaxial probe exactly at the center, i.e.,
Xp = 0 and Yp = 0, (orequivalentlyF = f).
(9) Simulate the U-Slot geometry, as obtained via
steps 1 to 8, using IE3D (or any other microstrip
CAD package [27]), and check for the 2:1 VSWR
performance.
(10) By examining the nature of the impedance variation on a Smith Chart, from step (9), adjust the
parameters (as mentioned in sec. V) for further
improvement of wideband performance of L'-Slot.
In order to illustrate the application of the preceding design steps, several topologies were modeled - for
which thcjesults arc available in [29^,130]. The results
contained in [29] (and [30])"mainly demonstrate the applicability of the empirical design procedure for those
substrates for which the empirical relations are contained
in table II. In this paper separate results are presented for
a typical case tr = 3.27 (TMM3), for which no information is available in tables I and 11. Since the permittivity
of (TMM3) lies in the range 2.94 < e/= 3.27 < 4.5,
so either of the equations in table II can be used. These
would result in two different dimensions for the U-Slot,
as ^ would be different for the two cases. The results
from IE3D [24] simulation forthe two U-Slot topologies
arc presented here to illustrate any such differences.
For the TMM3 (tr = 3.27) substrate, an operating frequency fr = 2.3 GHz was chosen. From the condition
' i^ ss 0.15, it was found that h - 1.08 cms. Referring to table II the equations for Cr = 2.94 and 4.5,
corresponding to a substrate thickness h = 1.0 cms.
were chosen. It was found that ^ = 4.738(£r =
2.94) and 3.975(er = 4.5) via the two appropriate empirical relations from table II. The remainder of the dimensions were easily found following steps (1) to (8).
At this stage, the (initial) unoptimized design, obtained by following steps (1) to (8) only, was characterized by the IE3Dcode. The VSWR vs. frequency results
were then examined, and the probe location was changed
from its initial/unoptimized value (Xp = 0 and Yp = 0)
by moving the probe along the y axis to Yp = 0.1 cms,
in view of the results in Fig. 5. (Various other optimization options in section V could have been pursued
as well.)The VSWR and boresight gain vs. frequency
for the unoptimized and optimized U-Slot topologies are
compared to demonstrate the quality of the initial (unoptimized) design obtained via steps (1) to (8).
The final dimensions of the two U-Slot patches (e^ =
3.27) obtained from the two empirical equations in table II are given below. (These geometries include the
optimized probe locations obtained for enhanced bandwidths.)
(i) via the empirical equation for Cr = 2.94 in table II:
W= 4.738, L= 3.424, W, = 1.841, L, = 2.37,
h= 1.0, t= 0.2G5, a=b= 0.527, Xp = 0. Yp =
0:i anddprotc = 0.127 cms; substrate tr = 3.27
(ii) via the empirical equation for Cr = 4.5 in table II:
W= 3.975, L= 2.872, W, = 1.545, L. = 1.988,
h= 1.0, t= 0.223, a=b= 0.442, Xp = 0, Yp =
0.1 and dprobe = 0.127 cms; substrate tr = 3.27
For the two topologies, the VSWR, Gain variations
vs. frequency and the radiation patterns in the cardinal
planes (0 = 0° and 90°) were obtained via IE3D code
[24], and are shown in Figs. 10 to 17. The data in Figs.
10 to 13 compare the performances of the unoptimized
and optimized designs. The results are briefly discussed,
next.
i.r
ij
13
Fig. 10. lllusiratiiig difTercnMS in VSWR for unoptimized (* — • — • ;
Xp =Yp = 0.0 cms), and opiimizcd (o — o — o: Xp = 0.0 and
Yp = 0.1 cms) U-Slot geometry as described in (i).
The results in Figs. 10 and 11 suggest that the unoptimized U-Slot design for tr = 3.27, obtained via
Nalatajan and Chatterjee: Design of Wideband U-Slot Microstrip Patch Antennas
198
Fig. 11. Illustrating differences in borcsight Gain ((^ = 0°,$ = 0°)
for unoptimized (« — • — • ; Xp =Yp = 0.0 cms), and optimized
(o - o - o: Xp = 0.0 and Yy = 0.1 cms) U-Slot geometry as
described in (i).
Fig. 13. Illustrating differences in boresight Gain (^ = 0°, 8 = 0°)
for unoptimized (• -• -*:
Xp =Y, = 0.0 cms), and optimized
(o - 0 - o: Xp = 0.0 and Yp = 0.1 cms) U-Slot geometry as
described in (ii).
the £r = 2.94 equation, performs well within expectations. The 2:1 VSWR bandwidths for the unoptimized'
and optimized topologies are w 33% and 30%, respectively. One also notices that for the imoptimized case the
overall VSWR performance is somewhat inferior compared to the optimzied U-Slot topology. This conclusion can be reached by examining the VSWR behavior
in the vicinity of 2 GHz for the two cases. By moving the
probe, however, the overall VSWR behavior is improved
at the expense of some reduction in bandwidth. The gain
behavior in Fig. 11 does not show marked changes.
would not show perceptible changes in gain behavior
if their wideband performance is enhanced by changing the probe location. Similar influence of probe location was observed for other designs, and are contained
in[29],[30].
Fig. 14.
VSWR characteristics for U-Slot jm substrate (^ =
3.27,tan S = 0.002 and fe = 1.0 cm; * - » - *'and o - o - o refer
to designs derived from empirical equations for Cr = 2.94 and 4.5,
respectively, from Tible II. The data refers to optimized designs only.
Results in Figs. 14 to 17 refer to opfi'mfeerf topologies
as defined earlier in (i) and (ii). The VSWR results in
Fig. 12. Illustrating differences in VSWR for unoptimized (• -• -«:
Fig. 14 indicate that both designs offer « 30% bandXp =Yp = 0.0 cms), and optimized (o — o — o: Xp = 0.0 and
Yp = 0.1 cms) U-SIot geometry as described in (ii).
widths corresponding to a return loss of 10 dB. However
the frequency ranges over which such wideband behavResults in Figs. 12 and 13 exhibit similar characteris- ior occurs is different for the two designs. For both detics when compared to Figs. 10 and 11. Comparing the signs, VSWR < 2 at the original design frequency of 2.3
VSWR variations in Figs. 10 and 12, one notices that un- GHz.
optimized (initial) design is a good estimate that couldbe
The gain characteristics shown in Fig. 15 reveal that
optimized using few iterations on the probe location. In- the two U-Slot designs exhibit overall similarities, except
terestingly, the boresight gain variations in Figs. 11 and that the topology derived from the empirical equation for
13 are rather insensitwe to probe locations. This obser- Cr = 4.5 shows improved'gain-flataess'. Howeverat2.3
vation suggests that initial (unoptimized) U-Slot designs. GHz the gain values are 0.6 and 4 dB, respectively, for
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
199
Signs.
A^HIIC MiXM*
Fig. 15. Boresight (fl = 0°,^ = 0°) gain characterislics forU-Slot
on sub.'itralc fr = 3.27, lani5 = 0.002 and /i = 1.0 cm; • — » — »
and o — o — a refer to designs derived from empirical equations for
fr = 2.94 and 4.n, respectively, from Table II. The data refers to
optimized designs only.
U-SIol geometry defiited in (i) and (ii).
Examining the VSWR and gain characteristics for the
two optimized designs (i) (fr = 2.94), and (ii) (e^ =
4.5), one can explain the differences by realizing that
these are essentially two different U-SIot geometries on
the same substrate material (er = 3.27andh = l.Ocms).
When these two were simulated via the IE3D code the
results were, quite predictably, different and is seen in
Figs. Hand 15.
Fig. 17. Radiation pattern vs. polar angle O'm (ft = 90° plane for
U-Slol on substrate fr = 3.27, tan* = 0.002 and h = 1.0 cm;
o - o - o (2.0 GHz) and D - a - a (2.3 GHz) refer to designs derived
from empirical equations for fr = 2.94 and 4.5, respectively, from
Table n. The data refers to optimized designs only.
One may thus conclude, judging qualitatively the results in Figs. 10 to 17, that the empirical formulas in
table IT are reasonably reliable, and have a wider scope
of applicability as compared to [20]. Follow-up investigations, comparing the present design method and [20,
sec. Ill] are planned for future work.
VIII. DISCUSSION AND FUTURE WORK
Fig. 16. Radiation pattern vs. polar angle 0 in ^ = 0° plane for
U-Slol on substrate fr = 3.27, tan* = 0.002 and /i = 1.0 cm.
o - o - o (2.0 GHz) and D - D - a (2.3 GHz) refer lo designs derived
from empirical equations for fr = 2.94 and 4.5, respectively, from
Table 11. The data rcfcn to optimized designs only.
The principal plane radiation patterns in Figs. 16 and
17. The patterns arc shown for 2.0 and 2.3 GHz, for respective maximum boresight gains as in Fig. 15. The
results show almost no difference for the two U-SIot de-
The empirical design technique developed here for Uslot patches in section VII is by no means exhaustive.
The most important part in the design are embodied in
steps 1 to 8 in sec. VII, that yield a basic design. Steps
9 and 10 are essentially equivalent to the optimization
functions available in the IE3D code [24]. The utility
of stepT9 and 10, from a practical point of view, stem
from the fact that not all microstrip anterma CAD tools
contain this optimization facility [27] like the IE3D code
[24]. In those special situations steps 9 and 10'can play
a critical role in the final design. It is however important to assess how the optimizers, currently available in
IE3D, compare with an U-slot design obtained via the
parametric simulation studies (steps 9 & 10 in sec. VII).
This work is currently in progress, and the final (optimized) results shall also be compared against the Ansofl ENSEMBLE microstrip CAD software for the various optimization options available in IE3D [24]. In addition, in view of the recent work reported in [20], additional investigations are necessary to compare the two
different U-Slot design approaches, i.e., in section VII
and [20, sec. III]. This comparison needs to be carried
out with special emphasis on high permittivity substrates
(r > 6.0.
Natarajan and Chatterjee: Design of Wideband U-Slot Microstrip Patch Antennas
As mentioned earlier, this empirical design technique
is limited because of symmetries associated with U-sIot
topologies. While the dimensional invariance and empirical formulas initiate a design, the complete process
doesn't provide any analytical insight into the electromagnetic behavior of the antenna. Recent investigations
[19],[20] on U-SIot design and performance modeling,
suggest the need for analytical development. To eliminate this present limination, efforts to develop an analytical formulation - based on the generalized cavity [1, pp.
97-102] and multipott network [1, pp. 103-108] models
- are under consideration. The main objective of such
a fiiture effort would be to develop/derive semi-rigorous
formulas that are far less empirical in nature than that
presented here. Such effort would also be aimed at providing information on the nature of currents flowing on
the patch surface, and slot resonance(s). There is some
information on slot resonance in [9] and [10], but they
appear valid for low-permittivity substrates.
IX. CONCLUSION
Earlier investigations have shown that introducing an
U-shaped slot on the radiating sur&ce of a probe-fed,
rectangular microstrip patch antenna, on a single layer
substrate, resulted in ultra-wideband topologies with better than 10 dB return loss performance. Since most of
the published results on U-slot were for foam or air substrates, and no systematic design procedure was available, a technique has been proposed in this paper for the
design of such a class of wideband microstrip antennas
on low to high permittivity microwave substrates. By
examining the earlier data for U-slots, the key feature
of dimensional invariance was established, and empirical equations have been derived based on data obtained
from extensive (IE3D) simulations. Subsequently, a systematic design procedure, for rapid parametric simulation and design of U-slot patch antennas, has been proposed incorporating the two above-mentioned features.
Using this proposed technique, VSWR and boresight
(4> = 0° and 9 = 0°)_gain comparisons were performed
for the unoptimized and optimized U-Slot topofogiesr The results suggest that initial designs for U-slots that
can be optimized within a few iterations via parametric simulations or state-of-art microstrip antenna analysis
CAD tools, such as IE3D or ENSEMBLE.
ACKNOWLEDGEMENT
This work was largely supported in parts through several research contracts firom Honeywell Federal Manafacturing and Technologies, KC, MO, and, the University Missouri Research Board (UMRB). In addition, the
authors wish to thank Richard D. Swanson, Senior Staff
Engineer, Honeywdl Federal Manafacturing and Technologies, KC, MO, for his continued active interest in
this work, and to Dr. Jian-X. Zheng of Zeland Corp.,
Fremont, CA, for many usefiil hints and clarifications regarding the IE3D software. Finally, the authors thank
200
Ms. Julie Burchett with Rogers Corp., Chandler, AZ, for
supplying TMM and Ultralam substrate samples.
REFERENCES
[IJ R. Garg, P. Bhaitia, L Bahl and A. Ittipiboon, Micmslrtp Antenna Design Handbook. Norwood, MA, USA: Aitech House,
Inc.,200L
[2] K. Fujimoto and J. R. James, Mobile Antenna Systems Handbook
(2''e6.).ibid.,
[3] K.-L. Wong, Compact and Broadband Microstrip Antennas. NY,
USA: John-Wley & Sons, Inc., 2002.
[4] T. Huynh and K.-F. Lee, "Single-Layer Single Patch Wideband
Microstrip Antenna," £te<;/ronicsie«e«, vol.31,no. 16, pp. 13101312, August, 1995.
[5] K.-F. Lee, K.-M. Luk, K.-F Tong, S.-M. Shum, T. Huynh & R. Q.
Lee, "Experimental and Simulation Studies of the Coaxially Fed
U-Slot Rectangular Patch Antenna,'V£E Proc. Mierow. Antennas
Pmpag., part H, vol. 144, no. 5, pp. 354-358, October 1997.
[6] Y. L. Chow, Z.-N. Chen, K.-F. Lee & K.-M. Luk, "A Design Theoiy on Broadband Patch Antennas with Slot"lS9eIEEE AP Symposium Digest, vol. 2, pp. 1124-1127, Atlanta, GA June 21-26,
1998.
[7] M. Clenet and L. Shafai, "Multiple Resonances and Polarisation
of U-Slot Patch Antenna," Electronics Letters, vol. 35, no. 2, pp.
101-103, January 21,1999.
[8] K.-F. Tong, K.-M. Luk, K.-F. Lee and R. Q. Lee, "A BroadBand (/-Slot Rectangular Patch Antenna on a Microwave Substiate,"/£££ Trans. Antennas Propagat., vol. 48, no. 6, pp. 954960, June 2000.
[9] R. Bhalla and L. Shafai, "Resonance Behavior of Single U-Slot
and Dual U-Slot Antenna,"'2001/£E£y</'-S5jpm/ios/iim Digest,
vol. 2, pp. 700-703, Boston, June 2001..
[10] R. Bhalla and L. Shafai,"Resonance Behavior of Single U-slot
Microstrip Patch Antenna," Micrmvave and Optical Technology
letters, v. 32, no. 5, pp. 333-335, March 5,2002.
[11] A. K. Shackelford, K.-F Lee, K.-M. Luk and R. C. Chair, "U-Slot
Patch Antenna with Shoituig fin"Electronics Letters, vol. 37, no.
12, pp. 729-730, June 7,2001.
[12] A. Shackelford, K.-F. Lee, D. Chatterjee, Y.-X. Guo, K.-M. Luk
and R. C. Chair, "Small Size Wide-Bandwidth Microstrip Patch
PMtmari'i'^'^IEEE AP-SSymposium Digest, vol. 1, pp. 86-89,
Boston, June 2001.
[13] Y.-X. Guo, A. Shackelford, K.-F. Ue and K.-M. Luk, "Broadband Quarter-Wavelength Patch Antennas with a U-shaped
Slot,"Wicn)wave and Optical Technology Letters, v. 28, no. 5, pp.
328-330, Mareh 5,2001.
[14] K.-F. Tong, K.-M. Luk and K.-F. Lee,"Wideband H-shaped
Aperture-Coupled U-Slot Patch Antenna,''A/icrowave and Optical
Technology Letters, v. 28, no. 1, pp. 70-72, Januarys, 2001.
[15] Y. X. Guo, K. M. Luk, K.-F. Lee and Y. L. Chow, "Double USlot Rectangular Patch Antenna,"£fec//onicj Letters, vol. 34, no.
19,pp. 1805-1806, Sept. 17,1998.
[16] M. Clenet, C. B. Ravipati and L. Shafai, "BandwidUi Enhancement of U-Slot Microstrip Antenna Using a Rectangular Stacked
Patch,"AficnJWflve and Optical Technology Letters, v. 21, no. 6, pp
393-395, June 20,1999.
[17] Y.-X. Guo, K-M. Luk, K-F. Lee and R. Chair, "A Quarter-Wave
U-Slot Patch Antenna wifli Two Unequal Anns for Wideband
and Dual-Frequency Operation," IEEE Ti-ans. Antennas Propagat.,
vol. 50, no. 8, pp. 1082-1087, August 2002.
[18] B.-L. Ooi, S. Qin and M.-S. Leong, "Novel Design of BroadBand Stacked Patch Antenna," IEEE Ihu». Antennas Propagat.,
vol. 50, no. 10, pp. 1391-1395, October 2002.
[19] A. K. Shackelford, K.-F. Lee and K. M. Luk, "Design of SmallSize Wide-Bandwidth Microstrip-Patch Antennas,";££E Antennas and Propagation Magazine, vol. 45, no. 1, pp. 75-83, Febniaiy
2003. (See also IEEE Antennas and Propagation Magazine, vol.
45, no. 2, p. 58, April 2003 for collections.)
[20] S. Weigand, O. H. Huff, K. H. Pan and J. T. Bemhard, "Analysis
and Design of Broad-Band Single-Layer Rectangular U-SIot Microstrip Patch Antennas," /£££ Trans. Antennas Propagat., vol.
51, no. 3, pp. 457-468, March 2003.
[21] R, A. Sainati, CAD of Microstrip Antennas for Wireless Applications. Norwood, MA, USA; Artech House, Inc., 1996.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
201
[22] Jian-X. Zheng, "A 3-D Ekclromagnetic Simulator for HighFrcquehcy Applications," IEEE Trans, on Components. Packaging
andManafacluring Technology Parl-B, vol. 18, no. 1. pp. 112-118,
February 1995.
[23]
, "Three Dimensional Electromagnetic Simulation of Electronic Circuits of General Shaperinlernalional Journal of Microwave and Millimcter-Wavc Compuler-Aided Engineering, vol.
4, no. 4, pp. 384-395, 1994.
[24] Zeland Software, Inc., IE3D User's Manual. Re!ease(s) 7.0, 8.0
& 9.0. Fremont, CA, USA.
[25] E. K. Miller, "Characterizalion, Comparison and Validation of
Electromagnetic Modeling So!>vii>n,"ACES Journal, vol. 3, no. 2,
pp. 8-24,1989.
[26]
."Development of EM Modeling Software Performance
Standards,"! 988 IEEE APS Symposium Digest, vol. 3, pp. 13401343, June 6-10, 1988.
[27] D. M. Pozar, S. Duffy, S. Targonski and N. Herscovici, "A Comparison of Commercial Software Packages for Microstrip Antenna
Analysis,"2000/£££/</'-5SynViojmm Digest, vol. I, pp. 152155, Salt Lake City, UT, July 16-21,2000.
[28] K.-F. Lee, "Design and Simulation of Miniature Wideband U-slol
and L-Probc Microstrip Patch Antennas,"; Part I of "Simulation.
Design and Development of Miniature, Broadband, Low-Profile
Antennas for Cellular and PCS Communications"; technical report
ECE-UMKCOO-HN YWL-TROl, prepared under contract for Honeywell Federal Manafacluring and Technologies, KC, MO, ECE
Division, SICE, UMKC, September 2000.
[29] V. N,nlarajan and D. Ch.ittcrjee, "An Empirical Design Technique for U-SIot Micro.strip Patch Antennas on Microwave Substrates; Part-I: Infinite Ground Plane Analysis,"technical report
ECE-UMKC02-HNYWL-TR01, li/rf.. September 2002.
[30] Venkataraman Natarajan, "Some Design and Modeling Aspects
of a Class of Wideband Microstrip Antennas Using Commercially
Available CAD Softwares,"M.S. Thesis, EE Depanment, University of Missouri, Columbia, MO. December 2002.
Venkataraman Natarajan was bom in
Chennai, India on July 2, 1978. He received
his M.S. and B.S. degrees in Electrical Engineering from the University of MissouriColumbia, USA, and University of Madras,
Chennai, India, in 2002 and 2000, respectively. From August 2000 he has been working as a Graduate Research Assistant at the
Computation.il Electromagnetics Laboratory
at the University of Missouri-Kansas City.
He was awarded the Outstanding Graduate
Student award in April 2001. He has been an IEEE Snidcnt Member
since August 1998.
Deb Chatlcrjec received D.E.TeI.E (Bachelor's) degree in Electronics and Telecommunication Engineering from Jadavpur University, Calcutta, India, in 1981, M.Tech^(Master's) degree in Electronics and Electrical Communication Engineering from IIT
Kharagpur, India, in 1983, M.A.Sc (Master
of Apphed Science) degree in Electrical and
Computer Engineering from Concordia University, Montreal, Canada, in 1992, and, the
Ph.D. degree in Electrical Engineering from
the Department of Electrical Engineering and Computer Science, University of Kansas, Lawrence, USA in 1998. From 1983-1986 he
worked with the Antenna Group at Hindustan Aeronautics Limited
(HAL), Hyderabad, India, on monopulse feeds for applications to the
Fire Control Radar. Since Augu.st 1999 he is an Assi-slant Professor
with the Electrical and Computer Engineering. University of Missouri
Kansas City (UMKC). His research interests are in analytical and numerical techniques in electromagnetics, with an emphasis on highfrequency (asymptotic) applications, performance modeling of adaptive (sniiirt) ,ind phased array antennas, and electromagnetic effects in
biological .sy.stcms. Dr. Chatlcrjec is a member of the IEEE Antennas
and Propagation Society, and reviewer for the IEEE Transactions on
Antennas and Propagation, IEEE Antennas and Wireless Propagation
Letters and Radio Science.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
202
UTILIZING PARTICLE SWARM OPTIMIZATION IN THE FIELD COMPUTATION
OF NON-LINEAR MAGNETIC MEDIA
A.A. Adly' and S.K. Abd-El-Hafiz^
'Cairo University, Electrical Power and Machines Department, Giza, 12211, Egypt
^Cairo University, Department of Engineering Mathematics, Giza, 12211, Egypt
Abstract
It is well known that the computation of magnetic fields
in nonlinear magnetic media may be carried out using
various techniques. In the case of problems involving
complex geometries or magnetic media, numerical
approaches become especially more appealing. The
purpose of this paper is to present an automated particle
swarm optimization approach using which field
computations may be carried out, via energy
minimization, in devices involving nonlinear magnetic
media. The approach has been implemented and
simulations were carried out for different device
configurations. It is found that the computations
obtained using the proposed approach are in good
qualitative and quantitative agreement with those
obtained using the finite-element approach. Details of
the proposed approach, simulations and comparisons
with finite element results are given in the paper.
L INTRODUCTION
It is known that magnetic field computation in
nonlinear magnetic media may be carried out using
various techniques (refer, for instance, to [1,2]).
Obviously, for cases involving complex geometries
and/or magnetic media, numerical and artificial
intelligence approaches become especially more
appealing (see, for example, [3,4]).
Irrespective of the adopted approach, geometrical
domain subdivision is usually performed and local
magnetic quantities are considered. One way to obtain
an electroinagnetic field solution is through the
minimization of the problem's energy functional, which
may take complicated non-quadratic forms. The
purpose of this paper is to present an automated particle
swarm optimization (PSO) approach [5,6] using which
2-D field computations may be carried out in devices
involving nonlinear magnetic media. More specifically,
the magnetic energy functional is fu^t formulated in
terms of the unknown magnetic vector potentials
corresponding to the discretization scheme.
A swarm of particles, each designated by a position
vector that represents the unknown potentials, is
initially randomly generated. These position vectors
may be regarded as potential solutions to the energy
minimization problem. The swarm is repeatedly moved
(i.e., modified) by the optimization algorithm and, upon
convergence, the unknown magnetic vector potentials
are found. Consequently, the field distribution is
computed everywhere.
Among the advantages of the proposed approach are
its simplicity, ability to handle complex magnetic media
and computational efficiency. The proposed approach
has been implemented and computations were carried
out for different electromagnetic device configurations.
These computations showed good qualitative and
quantitative agreement with results obtained using the
finite-element (FE) approach. Details of the approach,
computations and comparisons with results obtained
using the FE approach are given in the following
sections.
n. PROBLEM FORMULATION
In general, 2-D electromagnetic field problems may
be reduced into 1-D problems using the following wellknown magnetic vector potentiaMjT^ formulation:
B=fJl=iVxAu^),
(1)
where fi is the permeability, u^is a unit vector
orthogonal to the problem plane, H is the magnetic field
vector and B is the magnetic flux density vector.
For nonlinear media, the magnetostatic energy
functional E may be expressed in the form:
jybdb-JA \dCl,
(2)
n
where £1 represents the problem domain area, y is (he
reciprocal of the magnetic permeability (i.e., y=M~^).
and J is the current density along u^.
Neglecting hysteresis effects, the B-H relation of
most non-linear magnetic materials, especially those
used in electromagnetic power devices, may be
reasonably approximated by:
B=C !j|l[ ef{ and H
1054-4887 ©2003 ACES
m"
«B.
(3)
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
203
where, n is an odd number, C is a constant while e^
and e^arc unit vectors along the field and flux density
directions, respectively.
By subdividing the problem domain into finite
magnetic and nonmagnetic regions, and from
expressions (l)-(3), the magnetostatic energy
formulation becomes:
III
t^^'Mrf^n+o
r"
db-JrA^r
III.
AQr
(4)
where, /i°'' is the average value of the magnetic vector
potential within the r'* subdivision, AQ^ is the area of
the r"' subdivision, //^ is the permeability of free space,
while P,„ and P„,„ represent the number of magnetic
and non-magnetic domain subdivisions.
In the case when triangular domain subdivisions are
adopted, the value of the vector potential within the r'*
subdivision may obviously be expressed by:
Ar (A:, y) = a^o (Ar],Ar2-Ar3 )+ "rx (Arl Ar2-Ar3 )«
+ ary{Ar].Ar2.Ar3)y' (5)
where aro.Orx^ndafyaTC constants within the r'
triangular subdivision and may be formulated in terms
of the vector potential values at its three vertices (i.e.,
ArhAr2.Ar3)Hence, from (1) and (5), corresponding flux density
components can be written in the form:
Bry=-arxiAr]'Ar2^Ar3)
. (6)
Br\ = ^a?x(Arl>Ar2^Ar3)+a^iAr],Ar2, A^s)
Sub.stituting (6) into (4), we obtain:
r=l
It turns out that obtaining a solution of a
magnetostatic field problem involving non-linear
magnetic media should always correspond to a
minimization of expression (7). In this work, obtaining
the solution is achieved through the utilization of the
PSO approach to directly search for the appropriate
vector potential values as explained in the following
section.
(<4 iArJ.Ar2A3 )+ «o> (>^r;.>*r2 Ar3 )/''
in + 2)C"
j iArl+A,2+Ar3)
"
3
Aa^
ENERGY
MINIMIZATION
USING
PARTICLE SWARM OPTIMIZATION
Particle Swarm Optimization (PSO) is an
evolutionary computation technique, which simulates
the social behavior of insect swarms or bird flocks.
Different from traditional search algorithms, such an
evolutionary computation technique works on a
population of potential solutions of the search space.
Through cooperation and competition among the
potential solutions, this technique is suitable for finding
global optima in complex optimization problems [5].
The idea of our PSO implementation is that the
whole swarm proceeds in the direction of the swarm
member with the best fitness in a more or less
stochastic way. A swarm consists of several particles
M, where each particle keeps track of its own attributes.
The most important attribute of each particle is its
current position as given be an W-dimensional vector A*
= (Aj<A2
.A^). corresponding to a potential
solution of the function to be minimized.
Along with the position, each particle has a current
velocity, v'=(»/5^,v|,....,v^), keeping track of the
speed and direction in which the particle is currently
traveling. Each particle also has a current fitness value,
which is obtained by evaluating the objective function
at the particle's current position. Additionally, each
particle remembers its own personal best position, p
= (P^<P2'—-'PN^- ^^ "^* swarm level, the best
overall position among all particles, p' , is also
recorded. The index of the best particle is represented
by the symbol g. Upon termination of the algorithm;/>*'
will serve as the solution.
During each epoch, every particle is accelerated
towards its own personal best as well as in the direction
of the global best position. This is achieved according
to the following two expressions [6].
vf = wxvf +c\Xrand}0'x(.pf-Af)
+ C2^rand20x(.pf-Af)
2^0
iAr] + Ar2 + Ar3)
r=l
*
3
An..
(8)
(9)
(7)
where k = I
M, / = 1
N, C] and cj are two
Adly and Abd-El-Hafiz: Particle Swarm Optimization in Non-Linear Magnetic Media
positive constants, randl( ) and rand2( ) are two
random functions in the range [0, 1], while w is the
inertia weight. Expression (8) is used to calculate the
particle's new velocity according to its previous
velocity and the distances of its current position from
its own best experience (position) and the group's best
experience. If the velocity is higher than a certain Mmit,
Vmaxy this limit will be used as the new velocity for this
particle in this dimension. Thus, keeping the particles
within the search space. Then, the particle flies toward
a new position according to expression (9).
A number of factors affect the performance of the
PSO. First, the number of particles in the swarm, M,
affects the runtime significantly. A balance between
variety (more particles) and speed (fewer particles)
must be considered. Another factor is the maximum
velocity parameter, v^a. which controls the maximum
global exploration ability PSO can have. A very large
value for this parameter can result in oscillation. On the
other hand, a small value can cause the particle to
become trapped in local minima. The inertia weight, w,
is employed to control the impact of the previous
history of velocities on the current velocity. Thus, it
influences the trade-off between global and local
exploration abilities of the particles. A larger inertia
weight facilitates global exploration and the search of
new areas. A smaller inertia weight tends to facilitate
local exploration to fine-tune the current search area.
Suitable selection of the inertia weight can provide a
balance between global and local exploration abilities
[6].
For the implementation used in this paper, the
number of particles, M, has been set equal to 3S. The
maximum velocity, v„um is set to a constant value,
which is equal to half the size of the search domain.
The inertia weight is used to attenuate the magnitude of
the velocity updates over time. This attenuation is a
linear function of the current epoch number. Thus, w
linearly decays from about 0.52 to 0.48. The constants
Ci and C2 are set to the default value of 2.
IV.
NUMERICAL IMPLEMENTATION
SIMULATION RESULTS
AND
204
First, simulations were carried out, using the
proposed PSO approach, for an electromagnet having
length, width, depth and air-gap length of 1.0m, l.Ora,
0.25m and 0.1m, respectively. Magnetic field and
current density distributions were investigated subject
to coil excitations corresponding to current densities of
1.5x10*^ A/m^
and
0.75x10^ A/m^.
These
excitations were chosen to drive the core in the initially
linear and nonlinear magnetization ranges. Beside
giving an idea about the relative dimensions of the coil
in comparison with the core, Figs. 2-5 demonstrate the
simulation results. In these figures a vector is plotted at
the center of every triangular sub-domain. In other
words, the figures give also information about the
discretization scheme.
In order to check the accuracy of the proposed
approach, computations were compared to FE identical
simulations. This comparison revealed good qualitative
and quantitative agreement. For instance, maximum
flux density obtained using the FE approach for the
high and low excitation levels were found to be
1.I3T and 0.59T, respectively. Furthermore, Figs. 67 demonstrate the extent to which results obtained using
the proposed approach match FE computations.
Fig.l. Assumed exact B-H properties
approximation using expression (3).
and
its
iffegg^^
The proposed approach has been implemented and
computations were carried out for three different
electromagnetic device configurations. Throughout the
simulations,
comparisons
were
made
vwth
computational results obtained using the Quickfield FE
package (Version 4.3) and the core nonlinear B-H
relation was assumed as shown in Fig. 1. As can be
seen from the same figure, this relation was reasonably
approximated using expression (3) by choosing n = 3
Fig.2. Flux density distribution obtained using the
proposed PSO approach for the elecfromagnet subject
and C = ^-^y,,
to an excitation of 1.5x10^ A/m^ (Bjasx =1.08 T).
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
205
— FE
£=Z=i3
-fi-PSO
w
TX
,
^A.^:=S=/
(m)
Fig.6. Comparison between PSO and FE computations
for the particular horizontal contour line located at the
center of the electromagnet window.
Fig.3. Magnetic field distribution obtained using the
proposed PSO approach for the electromagnet subject
to an excitation of 1.5x10 A/m .
* -2.0E*04
;\\l___J._il' I
Fig.7. Comparison between PSO and FE computations
for the particular vertical contour line located at the
center of the electromagnet window.
Fig.4. Flux density distribution obtained using the
proposed PSO approach for the electromagnet subject
to an excitation of 0.75x10^ A/m^ (Bn,ax =0.55T).
Fig.5. Magnetic field distribution obtained using the
proposed PSO approach for the electromagnet subject
to an excitation of 0.75x10 A/m .
Second, simulations were carried out for an
electromagnetic actuator having window diniensions,
depth, core opening, and plunger dimensions of
0.8mx0.9m, 0.2m, 0.50m and 1.0mx0.4m, respectively.
Once more, magnetic field and current density
distributions were investigated for; case#l when the
plunger is in contact with the core and, case#2 when it
is 0.1m apart. Excitation was kept constant
5
2
corresponding to a cupent density of 3x10 A/m .
Due to the symmetry of this configuration, results of
these distributions are only given in the right half of the
solution domain as demonstrated by Figs. 8-11. Again,
these figures give some information about the adopted
discretization scheme as well as the relative coil to core
dimensions. Figs. 12-13, on the other hand, reveal the
good agreement between results obtained using the
proposed PSO approach and those obtained through the
FE technique. It should also be mentioned that, using
the FE technique, the maximum flux density values
were found to be 1.16T for case#l and 0.311 for
case#2.
Adly and Abd-El-Hafiz: Particle Swann Optimization in Non-Linear Magnetic Media
206
Hh7
Fig.8. Flux density distribution obtained using the
proposed PSO approach for the electromagnetic
actuator when the plunger is in contact with the core
(Bmax=112T).
Fig.ll. Magnetic field distribution obtained using the
proposed PSO approach for the electromagnetic
actuator when the plunger is 0. Im apart from the core.
0.8
-*-Bx{FE)
By(FE)
o-By (PSO)
E 0.0
0.0
Fig.9. Magnetic field distribution obtained using the
proposed PSO approach for the electromagnetic
actuator when the plunger is in contact with the core.
(m)
Fig. 12. Comparison between PSO and FE computations
for the particular horizontal contour line located at the
center of the lower core yoke.
1.8E.»0S
-2.0E+04
Fig. 10. Flux density distribution obtained using the
proposed PSO approach for the electromagnetic
actuator when the plunger is 0.1m apart from the' core
(5max= 0.296 T).
Fig. 13. Comparison between PSO and FE
computations for the particular vertical contour line
located at the center line of the actuator.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
207
Finally, simulations were carried out for a switched
reluctance motor (SRM) whose data is given in Table I.
Particular attention was given to the flux density and
field distributions at the beginning and end of a
stepping phase. Such distributions, driven by a pair of
coil excitation corresponding to a current density of
3x10 A/m , are given in Figs. 14-17. Obviously,
these figures give some information about the adopted
di.scrcfization scheme. In these figures, only steady state
analysis was considered and no motion effects were
incorporated. In addition, the flux was assumed to be
confined inside the motor (i.e., no flux is assumed to
cross the outer stator and inner rotor diameters.
Additional comparisons with results obtained using the
FE technique are given in Figs. 18 and 19. It is also
worth mentioning that the maximum flux density values
computed using the FE technique at the beginning and
end of the stepping phase were found to be 0.08 T and
1.08 T, respectively.
TABLE I
Data of the simulated Switched Reluctance Motor
Number of Stator Poles
8
Number of Stator Poles
6
15°
Stator Pole Arc
Rotor Pole Arc
15°
Stator Core Outer Diameter
0.100 m
Stator Core Inner Diameter
0.750 m
0.450 m
Diameter at Stator Pole
0.445 m
Diameter at Rotor Pole
Rotor Core Outer Diameter
0.250 m
0.150 m
Rotor Core Inner Diameter
Fig.l4. Flux density distribution obtained using the
proposed PSO approach for the SRM at the beginning
of a stepping phase (Bj^ax = 0.07 T).
Fig. 15. Magnetic field distribution obtained using the
proposed PSO approach for the SRM at the beginning
of a stepping phase.
Fig. 16. Flux density distribution obtained using the
proposed PSO approach for the SRM at the end of a
stepping phase (Bmax =1-13T).
Fig. 17. Magnetic field distribution obtained using the
proposed PSO approach for the SRM at the end of a
stepping phase.
Adiy and Abd-El-Hafiz: Particle Swarm Optimization in Non-Linear Magnetic Media
• -Hx(FE)
stem from the non-identical discretization schemes,
An important feature of the proposed search
approach is its simplicity and ability to handle
adjacent discretizations of dramatically different
dimensions without the risk of running into
numerical difficulties related to matrix inversion.
This further highlights the computational memorywise efficiency of the approach.
More efforts are planned to further investigate the
proposed approach in different time harmonic
problems. It should also be stressed that PSO has been
previously used in device dimension optimization (see,
for example, [7]). Coupling this capability to the
proposed field analysis approach might pave the way
towards the ability to simultaneously optimize the
dimensions and compute fields in devices involving
nonlinear media.
c)
-*-Hy(FE)
Q-Hx(PSO) *-Hy{PSO)
'5.0Et03
208
"V/
f 1
0
(mm)
Fig.18. Comparison between PSO and FE computations
for the contour line passing through the centers of the
two energized stator poles at the beginning of the
stepping phase.
REFERENCES
97.5
Angle (degrees)
187.5
Fig. 19. Comparison between PSO and FE computations
for the contour arc located in the center of the air-gap
and subtending one energized stator pole to the other at
the end of the stepping phase.
IV.
DISCUSSION AND CONCLUSIONS
In light of the proposed PSO approach as well as the
experience gained while performing the presented
simulations, the following few remarks should be
stressed:
a) Depending on the problem nature and number of
unknowns, suitable values for the number of
particles M and inertia weight w should be first
investigated. In our case it took some time to
decide upon using the reported M and w values.
Once this sort of tuning is achieved, the proposed
approach gives reasonable results in computational
time comparable to that of the FE technique.
b) It can be observed from the presented results that
good agreement with FE computations can be
achieved using the proposed approach. It should be
pointed out, however, that some of the
discrepancies between results of both techniques
[1] J.K. Sykulski, Computational Magnetics, Chapman
& Hall, London, 1995.
[2] A.A. Adly and S.K. Abd-El-Hafiz, "Automated
two-dimensional field computation in nonlinear
magnetic media using Hopfield neural networks,"
IEEE Trans. Magn., Vol. 38, pp. 2364-2366,2002.
[3] M.V.K. Chari and S.J. Salon, Numerical Methods
in Electromagnetism, Academic Press, San Diego
CA, 2000.
[4] A.A. Adly and S.K. Abd-El-Hafiz, "Utilizing
Hopfield neural networks in analysis of reluctance
motors," IEEE Trans. Magn., vol. 36, pp. 31473149.2000.
[5] L Kennedy and R. Eberhart, "Particle Swarm
Optimization," Proceedings of the IEEE
International Conference on Neural Networks,
Perth, Australia, pp. 1942-1948,1995.
[6] Y. Shi and R. Eberhart, "A Modified Particle
Swarm Optimizer," Proceedings of the IEEE
International Conference
on Evolutionary
Computation, Anchorage, Alaska, pp. 69-73,1998.
[7] G. Ciuprina, D. loan and L-Munteanu, "Use of
intelligent-Particle Swarm Optimization in
Electromagnetics," IEEE Trans. Magn., vol. 38,
pp. 1037-1040,2002.
Amr A. Adly received the BS
and MS degrees in Electrical
Engineering
from
Cairo
University, Egypt, in 1984 and
1987, respectively, and the PhD
degree in Electrical Engineering
from
the
University
of
Maryland,
College
Park,
Maryland, USA, in 1992. In
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
209
1993 he joined LDJ Electronics, Troy, Michigan, USA,
as a Senior Engineer/Scientist. In 1994, he was
appointed as an Assistant Professor at the Electrical
Power and Machines Department, Faculty of
Engineering, Cairo University, Giza, Egypt. Since
1999, he has been working as an Associate Professor in
the same department. His research interests include
electromagnetics, modeling and simulation of complex
magnetic
materials,
electrical
machines,
superconductivity, magnetic recording,
magnetic
measurement instrumentation, and magneto-hydrodynamics. In 1994, he received the Egyptian State Prize
for achievements in Engineering Sciences. He has
authored more than 60 journal and refereed conference
papers. Dr. Adly is a senior member of the IEEE
Magnetics Society and a member of the ACES Society.
Salwa
K.
Abd-El-Hafiz
received the BS degree in
Electrical Engineering from
Cairo University, Egypt, in 1986
and the MS and PhD degrees in
Computer Science from the
. • '^.L'A^ifP';>/ . University of Maryland, College
(!" i^-^r^<'r^> Park, Maryland, USA, in 1990
and 1994, respectively. In 1994, she was appointed as
an Assistant Professor at the Engineering Mathematics
Department, Faculty of Engineering, Cairo University,
Giza, Egypt. Since 1999, she has been working as an
Associate Professor in the same department. Her
research interests include software analysis and
specification, software measurements, software reverse
engineering, neural networks, evolutionary computation
techniques, knowledge-based systems and numerical
analysis. In 2001, she received the Egyptian State Prize
for achievements in Engineering Sciences. Dr. Abd-ElHafiz is a member of the IEEE Software Computer
Society.
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
210
AN APPROACH TO MULTI-RESOLUTION IN TIME DOMAIN BASED ON THE DISCRETE
WAVELET TRANSFORM
C. Represa(*), C. Pereira(*), A.C.L. Cabeceira ( *), I. Barba (**), J. Represa (* )
(*) Dpt. of Electromechanical Engineering. University of Burgos. 09001 Burgos. Spain.
Phone: +34 947258830
FAX: +34 947258831
E-Mail: crepresa@ubu.es
(• *)Dpt. of Electricity and Electron. University of Valladolid. 47011 Valladolid. Spain.
Phone:+34 983423223
FAX:+34 983423225
Abstract'.- In this paper, an approach to multiresolution in time domain (MRTD) is presented.
Maxwell equations are discretized using finite
differences in time and a derivative matrix in
space that allows any desired level of spatial
resolution. This derivative matrix acts on the
coefficients that represent the expansion of the
-field components. These coefficients are
calculated by means of the Discrete Wavelet
Transforms (DWT). In this work hard (PEC
and PMC) boundary conditions have been
introduced into the algorithm using the method
of images. This approach is valid for any kind of
wavelet functions. Stability and dispersion
properties are also investigated. Some numerical
results, showing multi-resolution properties are
presented.
Keywords.- Multi-resolution in time domain,
MRTD, Daubechies wavelets. Discrete Wavelet
Transform, DWT.
1.- Introduction
The development of electromagnetic fields in scale
and wavelet functions [1], has given rise to the
techniques known as Multi Resolution in Time
Domain (MRTD). These techniques are based on
the possibility of increasing the resolution of a
signal, from a coarse level to a fine one, using lowresolution fiinctions (scale functions), combined
with others of intermediate resolution levels
(wavelet functions). Different functions have been
used in this type of analysis: Battle-Lemarie [2],
Haar [3] or Daubechies [4].
In these techniques, the electric and magnetic field
components of Maxwell's curl equations are
represented by a manjrfold expansion in scale and
wavelet functions with respect to space, and step
functions with respect to time. The method of
' This work was sponsored in part by the Spanish Ministry of
Science and Technology, under its project 110-2000-1612-00302.
E-Mail: ibarbafS.ee.uva.es
moments [5] allows to obtain a set of equations
similar to the one used in FDTD (identical, in the
case of Haar functions, to the Yee's algorithm).
Instead of the field values, the coefficients of this
expansion are used in such scheme.
One of the points that affect the complexity of
these schemes is the type of functions used for the
expansion.
For example, if Battle-Lemari6
functions are used, an infinite number of
coefficients must be used, because of the noncompact support of these functions. Actually, we
can truncate the expansion and use a reduced
number of them, obtaining an approximate
solution. This problem does not appear when, for
example, Haar functions are used, since they are
compactly supported, so it is not necessary to
truncate the expansion.
A second point is the desired resolution level for
the solution of the problem. The above mentioned
approaches [2-4] use scale functions for a first
approximation of thfc fields; in this case simple
functions are obtained, but, if an increasing of the
resolution is desired, adding different levels of
wavelet functions, the solution becomes unfeasible,
as the scheme must be modified for every new
level we add. This problem can be solved if a
Discrete Wavelet Transform (DWT) is used to
derive the coefficients of the expansions, as it h^
already been used -to obtain the time domain
solution of electrical networks [6]
In this work, compactly supported Daubechies
wavelet fiinctions [1] have been used to make an
expansion of the fields. Its coefficients are
computed through a Discrete Wavelet Transform.
A derivative matrix, also obtained by means of a
DWT, allows us to compute the value of the
electric and magnetic fields at desired spatial
positions and with the desired resolution. We have
found the stability criterion of the algorithm for
different wavelet functions and in some cases it is
wider than the FDTD stability criterion for the
1054-4887 ©2003 ACES
ACES JOURNAL, VOL. 18, NO. 3, NOVEMBER 2003
211
same spatial resolution (twice in the case of Haar
wavelet functions).
sampling points with a distance of Az between
them. This is illustrated in figure 1.
B
2.- Formulation
The simplest case of a TEM plane wave
propagating in a homogeneous, linear, isotropic and
non-dispersive media, with fields E, and Hy is
analyzed. Then, the equations to solve are:
5H„
= -H-
et
(1)
5H„
dE,
dz~ ^ di
In the proposed scheme, every field component is
expanded in terms of scaling and wavelet fiinctions,
as shown in equation (2), where the coefficients of
the lower resolution level {at) are given by the
scaling function ((|)) and the coefficients of the
successive resolution level (6/) are given by the
wavelet functions (v|/):
I"a:.t.°(z)
F"{z) = 2:"ai<t>i(z) =
k
+zi:"bi;M/i(z)
k
(2)
1=0
The spatial discretization is achieved by dividing
the domain into N cells of size A, as a starting
point, then getting N sampling points in the center
of each cell. This is what we call resolution level 0
(j = 0). If each cell is now divided into 2' points,
the maximum level of desired resolution J(j = J) is
reached and then the simulation domain has N*2'
^-<D
^
'
Decomposition
Fig. I.- Spatial discretization for J=2.
The coefficients a/ are enough to represent the real
values of the field at each sampling point and they
are the initial and the final step of the
decomposition process involved in the Discrete
Wavelet Transform. The multi-resolution analysis
performed with the wavelet transform can be
understood as a digital filtering process where a
signal is decomposed in two parts, one containing
the low frequencies, the scale functions {f), and
other part containing the high frequencies, the
wavelet functions ((f^. The multi-resolution
representation through the Discrete Wavelet
Transform (DWT) is provided by successive filter
banks stages, each one containing a low-pass and a
high-pass filter, described in terms of the
coefficients of their impulse responses L(m) and
H(m) (meZ) respectively [1]. The filtering process
means successive convolutions of the field and the
filter coefficients followed by a decimation process
that retains only the even indexes. The inverse
transform (IDWT) consists in successive
interpolations followed by convolutions that gives
the field values using the coefficients of the
previous decomposition. This is sketched in figure
2.
^^<I^
-(^t
Wavelet Coetfldenis
l£
"
,J
D.
L
-(]>-t^
L?:-<I>
i-5
Reconsinictlon
IDWT
Fig. 2.- A step of the filtering process: the low and high-pass filters outputs the coefficients d and V
respectively.
Then, eqs. (I) arc discretizcd using centered finite
differences for the time derivatives and a spatial
derivative operator computed through a DWT
matrix for the spatial derivatives. Each coefficient
is, then, obtained from its values at the same point
in a previous time step and from the derivation of
the other field coefficients at an intermediate time
step:
Represa, et al.; Multi-Resolution in Time Domain Based on Wavelet Ttansforai
212
m = t^{'r>}j2m-rr\)
(8)
v{t) = |;H(mW2<t.(2t-m)
(9)
(3)
e A
^ '
where At and A are the time and space steps
respectively, and D"* represents the spatial
derivative operator of level J. This operator [7] acts
on the coefficients a and b of the field transform (3)
and it is denoted by a matrix D [8] which can be
computed previously. In this way, eqs. (3) can be
rewritten as follows for an arbitrary level of
resolution J:
E*'
"Hr
being i/the filter length. In figure 3 three different
aspects of the derivative matrix D for three levels
of resolution are depicted, where only the nonzero
elements have been plotted.
E A
M4)
H*'
^l A
E»'
where y=0,...,y-/. This matrix D is banded with a
limited band width, due to the compact support of
the Daubechies wavelets, and this width increases
as the size of the filters does. Its elements are
integrals of the scaling and wavelet functions and
its discrete spatial translations of the first derivative
(eqs. (6.1)-(6.4)). Because of multi-resolution, this
matrix can be decomposed in a set of submatrices
of lower resolution, so if we have named the
original matrix of highest resolution D', it is
constructed with the one level lower matrices A, B,
r, D, then:
ir
(5)
where A' = {a<}.B' = |p5},r' = {yS}.D' = fl/,}.
The components of every matrix can be calculated
in a recursive way, so:
aj=2'ai.,
(6)
and similar expressions for the other coefficients, p,
y and d. The expression ai is obtained by means of
the following expression:
a| = J\|/(z-l)-—v|»(z)dz
(7)
Similar expressions, combining scale and wavelet
functions can be used to calculate the other
coefficients p (scale.wavelet), y (wavelet and scale)
and d (scale and scale). The scale and wavelet
functions in every case are obtained by means of
the low-pass and the high pass filters, L(m) and
H(m) respectively:
(c)
Fig. 3.- Three different aspects of the derivative
matrix G for three resolution levels: (a) j=0, (b)
J=l. (c)j^2.
The use of this derivative matrix allows us to
achieve any level of resolution without
modifications of the algorithm. Moreover, we can
choose at the beginning of the simulation what type
ACES JOURNAL. VOL. 18, NO. 3, NOVEMBER 2003
213
of wavelet function is going to be used, and what
regions of the simulation domain we want tosolve
at higher resolution. This feature will be shown
further in the examples.
3.- Boundary conditions.
The derivative matrix D is built-up in a cyclic form
in such a way that it treats the coefficients on both
boundaries as contiguous. It results in a cyclic
space, that is, an infinite unbounded space. As the
simulation domain is not cyclic at all, the character
of the matrix must be modified. To do so, we
propose to add adjacent columns with the elements
related to the boundaries, and getting an extended
matrix; i.e.
0
a.,
a,
I
a.
0
On
a ,
a„
Oil/
0
y^,
0 i 0
0
0
a.,-1
0 I 0 "l
0 i 0
0
a
a.
and similar expressions for {P), {y}, and {d). The
matrix D acting on the column vector (b'', a'') at
resolution level J implies that its extension must act
on an extended vector too. The additional
coefficients needed to represent the boundaries can
be obtained using the method of images for PEC or
PMC walls. If di and i'i arc the scale and wavelet
coefficients respectively, belonging to a border cell
collocated at position z=Z.Az, we obtain the
additional coefficients a' and b' this way:
for even symmetry:
aK,.i
=+a,'N-i+l
"N
on the rigth
N+i "
ai-i' :+ai
= -b,
on the left
More complex structures
require
sophisticated techniques not studied here.
4.- Stability analysis
For the algorithm to be stable, the pair At and Az
must be chosen in a correct way. We choose these
values following the derivation given in [9] where
the stability problem is treated as an eigenvalue
problem. This means that plane wave cigenmodes
will be assumed to propagate in the numerical data
space. The spectrum of eigenvalues for these
modes due to the numerical space differentiation
process will be determined and compared to the
stable spectrum of eigenvalues determined by the
numerical time differentiation process. By
requiring the complete spectrum of spatial
eigenvalues to be contained within the stable range,
it is ensured that all possible numerical wave
modes in the grid are stable. In this way, equations
(1) can then be rewritten as follow:
H"^-H"-^
At
E"-'-E"
At
1
H-Az
I
Id,Er,
on the rigth
= +b,
on the left
These boundary conditions allow us to implement
in a easy way infinite electric or magnetic walls.
(10)
Id,Hr;
. e-Az ,
and this equations can be split into two eigenvalue
problems concerning time and space:
H"*--H"= AH|'
At
pn+1 _c "
At
(11)
= AEr
—i—VdiE,", =AHP
(12)
Xd,H:;^=AEr'
e-Az 7
where only .scaling functions of resolution level
J=0 have been used. To avofd having any mode
Increasing without limit during normal timestepping it is found from equation (11) that the
stable spectrum of eigenvalues is:
|lm(A)|<2/At
for odd symmetry:
more
(13)
In order to solve equations (12) we introduce a
typical mode of the spatial spectrum like (14) into
the equations:
(14)
Then we get the following set of equations:
Represa, et al.: Multi-Resolution in Time Domain Based on Wavelet Transfonn
H-AzY
(15)
EA27
Simplifying and using Euler's identity we get that
the spatial eigenvalues are given by:
lm(A) = —XdiSen(-klA2)
(16)
Therefore, the maximum value of this eigenvalues
is:
NA)|<;^I|d,|
(17)
To guarantee numerical stability, the range of
spatial modes must be contained completely within
the stable range of time-stepping eigenvalues set by
(17) and so:
Azr'"
(18)
At
Therefore, the upper bound of the time step AIMRTD
is:
^
2
Az
""^ I|d,|c
(19)
214
propagation, that is, the numerical phase velocity
given by the algori&m differ from the phase
velocity of the wave in the medium [II]. This
fictitious dispersive behaviour must be taken into
account specially in considering large structures of
simulation because significant differences between
real and numerical phase can be obtained. Some
studies about dispersion properties of MRTD
methods based on Galerkin procedures [12], [13]
have been done. Now we proceed to study the
dispersion characteristics of the MRTD algorithm
based on the DWT using Daubechies' wavelets
fimctions. To do that, a plane monochromatic
travelling-wave (22) is introduced into the
discretized Maxwell equations (4) and then we
search for the relationship between the angular
frequency m and the numerical wave number if:
"E'' =E e'^'""°"^')
n u; _ u . j[klAz-fflnu)
(22)
As starting point, we expand the field components
using only scaling functions of resolution level
J=0, so ? = (|) and the derivative matrix D* will be
composed of elements^,}. The discretized
Maxwell equations can be written in this way:
"* £-I ="£
- H -.^.^r^d^^H*,
(23)
n-4i_i* n-l
tH*='-fHf-^^2:4''E^
Defming a stability factor like (20):
s=-
AZ
(20)
equation (19) establishes that the stability factor
must be contained within the range:
"Depending on the iype of wavelet fiinction used,
this range implies that the time step can be set to a
different value for the same spatial resolution. If
Haar wavelet functions are used, this value can be
set double than the time step set in FDTD for the
same spatial resolution.
5.- Dispersion properties
The use of numerical techniques for solving
electromagnetic problems as TLM [10], FDTD or
MRTD require always a discretization process.
Such a process result in a phase error of the field
Substituting (22) into (23), simplifying and using
Euler's identity we get this set of equations:
2E,sen
oAt^At 1
;^d|H„sen(-klAz)
2
E Az
_L,
coAt At_l_
At 1
2H„sen
=
5^d|E(,sen(-klAzj
fi Az
(24)
and from them we obtain the dispersion
relationship where cis the speed of light in our
medium:
2Az
(oAt ^ .
-^sen—= ^d,sen|(-iciAz)
(25)
This expression may be reduced to the ordinary
expression to = ck, that is, the non-dispersive case,
if spatial resolution Az is very small in comparison
with the wavelength. If Haar wavelet functions are
used then we also obtain a non dispersive case at
the stability factor s = l (i.e. the maximum stable
value). In any case, dealing with broadband signals,
different expressions will be obtained depending on
the relation between the wavelength A and the
ACES JOURNAL. VOL. 18, NO. 3, NOVEMBER 2003
215
spatial resolution Az, and on the type of wavelet
function we use. In figure 4 it is plotted the
norrnalized phase velocity versus spatial resolution
given by the exprcsion (25) for two different values
of the stability factor. We can appreciate the feature
mentioned before, that is, the normalized phase
velocity tend to unity as we incresc the resolution.
It can also be seen that the higher the order of the
wavelet
function,
the
better
dispersion
characteristics.
NofmtiizBd phaso veloOylor S .05
101
0.C3
OM
0.06
002
CM
006
008
01
012
Oisiancs
0.14
016
0.18
Db—;::_.
I
Db, /-
Z'
OS)
6b,
09
'
/
'
— Db,
'a
/
Db,
(a)
Nomsiizad phaso yehdty forS^ 1
a.,\
fOl
oi,
0.1
012
Dtstan4:»
DM
0.15
018
(b)
^
~--^—
__
Fig. 5.- Three snapshots of a Gaussian pulse
propagation using (a) Haar wavelet functions, and
(b) Daubechies Db2 functions, showing it
dispersion properties.
^—""
y-^'^"^
Db,
oe)
f
> on
008
-
/
01
!
-j
0»
— Db,
6.- Results
DT
(
1
(b)
Fig. 4.- Normalized numerical phase velocity
versus spatial resolution in lambda units for
different waveletfunctions:
(a) with a stability factor s=0.5, and
(b) with a stability factor s=1.0
In figure 5 it is depicted a gaussian pulse
propagation using different types of wavelet
functions: (a) Haar wavelet functions, and (b)
Daubcchics D2 wavelet functions. The case (b)
exhibit a better dispersive behaviour than the case
(a) as it can be expected from previous figure 4.
To validate our proposed technique some aditional
examples have been performed. First of all we have
evaluated the multi-resolution property of the
algorithm. To do that, we have simulated a
gaussian pulse excitation into a one-dimensional
cavity of 1 m long. The simulation domain is split
into two parts of different resolution, one of high
resolution on the left side and 400 mm long, and
other part of low resolution on the right side and
600 mm long. Figure 6 displays this simulation
where it can be clearly seen the differences
between the two parts. Figure 7 shows in detail the
transition between zones.
Represa, et al.: Multi-Resolution in Time Domain Based on Wavelet Transform
2
"!
1
1
216
T
3Vm« =|250 ps!
Time = 1450 ps
I
1.5
W^
I
1
(
•
1
1
1.5
1
05
1
j". 0
-05
^-0.5
.1
0
0.1
_ir
L.
0.2
03
'
I
0.4
05
0.6
Distance
0.7
I
0.9
1_
09
i
'
High Resolution
JTime -ilSM (»
.A i
\
}.
\
\
• 0
1
....J
i
. ..i
Al [ 1 1 1
\
/ 1 M i I
r 1 i i i
\ J
-0.5
-j t44-¥-
.1
-1.5
1
i
1
0.1
02
0.3
0.4
i
!
0.5 06
Distance
1
0.7
1
OB
1
0.9
i
i
(b)
\
\
i
! [....Li
1 1 :/
1
•
1
0.2
\l
1
i
1
1
V;
;
i
;
1
i i i 1 i/
] M i\ 1
i ; i \l
:
;
0.3
1
1
1
1
-1
•
1
Low Refeolution
1
0.4
0.45
0.5
Distance
Fig. 7.- Details of the transition between two zones
of different resolution.
0.35
We have also evaluated the resonant frequencies of
a one-dimensional rectangular cavity with a 100
nun distance between electric walls (PEC). A
spatial resolution of Az = 0.625 mm and a time
discretization of At = 3.127 ps have been chosen
(that means, that, according to the stability criterion
we choose s = 1.5). The electric and magnetic
fields are specified at ^0 and at f=l/2At
respectively,
Ej=exp|
(26)
\1 M 1 1
......J
J.
\
iHighRKouion
!
01
1 v]
\2"
i Tinwi-MSOJps
\
_J^
-0.5
(a)
*
K
i /
:Lcrw*R?soiiiiiop
i
!
0.4
0.S
Q.S
tXslsnce
I
07
\J\
1
0.8
H^=exp
(27)
where the width of the pulse is ^^ = 10 and its
center kc = 80 in terms of the index of the spatial
mesh. After 4096 time steps with a spectral
resolution of Af = 78.1 MHz, the resonance
spectrum obtained have been plotted in figure 8
i
09
(C)
300
Fig. 6.- Gaussian pulse propagation using Haar
wavelet Junctions. The simulation domain is split
into two zones ofdifferent resolution.
"150
100
iJ uu ..ill
\J
J. j^
Fig. 8.- Resonant frequencies obtained in a cavity
oflength 100 mm using Haar waveletfunctions.
ACES JOURNAL, VOL. 18. NO. i, NOVEMBER 2003
217
7.- Conclusions
A new approach to Multi-Resolution in TimeDomain has been investigated. In this study a
derivative matrix is used to calculate the spatial
derivatives of the electromagnetic fields. This
matrix acts on the coefficients of the wavelet
expansions of the fields obtained from the Discrete
Wavelet Transform. The use of this matrix allows
us to solve the discretized Maxwell equations at
any level of desired spatial resolution and wherever
we may want. It has been shown that different
types of Daubechies wavelet functions exhibit
better dispersion properties for the same spatial
resolution and that the time step can be chosen
bigger than in other time-domain methods with the
same mesh size. In this way, wc can also choose at
the beginning of the simulation different types of
wavelet functions in order to improve our results.
9.- References
[1] I. Daubechies, Ten Lectures on Wavelets,
CBMS-NSF Ser. App. Math. 61, SIAM:
Philadelphia 1992.
[2] M. Krumpholz and L.P.B. Katehi, "MRTD:
New Time-Domain schemes based on
multiresolution analysis," IEEE Trans.
Microwave Theor)- Tech. 1996; 44(4):555571.
[3] M. Fujii and W.J.R. Hoefer, "A 3-D HaarWavelet-Bascd Multiresolution Analysis
Similar to the FDTD Method - Derivation and
Application -, " IEEE Trans. Microwave
Theory' Tech. 1998; 46(12):2463-2475.
[4] Y. W. Cheong, Y. M. Ue, K. H. Ra, and C.
C. Shin, "Wavelet-Galerkin scheme of timedependent inhomogeneous electromagnetic
problems," IEEE Microwave Guided Wave
Lett. 1999,9(8):297-299.
[5] R. F. Harrington, Field Computation by
Moment Methods, IEEE Press: New York
1992.
[6] S. Barmada and M. Raugi, "A general tool for
circuit analysis based on wavelet transform,"
-Int. J. Circuit Theory and Applications. 2000;.
28:461-480.
(7] G. Beylkin, R. Coifman and V. Rokhlin, "Fast
wavelet transform and numerical algorithms,"
I. Comm. PureAppl. Math. 1991, 44:141-183.
[8] G. Beylkin, "On the representation of
operators in bases of compactly supported
wavelets," SIAM J. Numer. Anal. 1992;
6(6): 1716-1740.
[9] A. Taflove, Computational Electrodynamics:
The Finite-Difference Time-Domain Method,
Artcch House 1995.
[10] C. Christopoulos, The Transrvission-Line
Modelling Method: TLM, IEEE MTT: New
York 1995.
[11] J. Represa, C. Pereira, M. Panizo, F. Tadeo,
"A Simple Demonstration of Numerical
Dispersion under FDTD," IEEE Trans.
Education. \991\W{\)m-\()2.
[12] E.M. Tentzeris, R.L. Robertson, J.F. Harvey
and L.P.B. Katehi, "Stability and Dispersion
Analysis of Battle-Lemarie-Based MRTD
Schemes, " IEEE Trans. Microwave Theory
Tec/;. 1999; 47(7):1004-10I3.
[13] M. Fujii and W.J.R. Hoefer, "Dispersion of
Time Domain Wavelet Galerki.n Method
Based on Daubechies' Compactly Supported
Scaling Functions with Three and Four
Vanishing Moments," IEEE Microwave
Guided Wave Lett. 2000; 10(4): 125-127.
Cesar Represa was bom in
Medina del Campo, Spain,
in 1971. He received the
Licenciado degree in Physics
from the University of
Valladolid in 1995, and the
PhD degree in 2001 from the
University of Burgos. He
was Profesor Asociado at
the University of Burgos
from 1997 to 2000. Since 2000 he has been
Assistant Professor at the University of Burgos. His
research interest includes numerical methods in
electromagnetics.
Carmen Pereira was bom in
Zamora, Spain, in 1951. She
received the Licenciado
degree in Physics in 1974,
and the PhD degree in 1983,
-.jvr -J ^ both from the University of
\ ""^"j^
Valladolid, Spain. She was
Assistant Professor at the
\-.
University of Valladolid from
1974 to 1985, and since 1985
has been Profesor Titular in Electromagnetics at
the Universities of Valladolid and -Burgos. Her
current research interest includes numerical
methods in electromagnetics and microwave
devices.
Ana Cristina L. Cabeceira
was bom in Pontevedra,
Spain, in 1969. She received
the Licenciada Degree on
Physics in 1992, and the Ph.
D. Degree in 1996, both from
the University of Valladolid,
Spain. From 1992 to 1999
she was Assistant Professor,
Represa, et al.: Multi-Resolution in Time Domain Based on Wavelet Transfonn
then she become Profesora Titular in
Electromagnetism of the University of Valladolid.
Her main research interests include numerical
methods for Electromagnetism, as well as
characterization of electromagnetics properties of
materials.
Ismael Barba was bom in
Palencia, Spain, in 1970. He
received his Licenciado
degree in Physics from the
University of Valladolid in
1993. He received his MA
degree
in
Electronical
Engineering in 1995, and his
PhD in Physics in 1997, all
from the same University. He
was Assistant Professor from 1994 to 1999, and
since 1999, he is Profesor Titular in
Electromagnetics in the University of Valladolid.
His main research interest includes numerical
methods in electromagnetics, and characterization
of electromagnetics properties of materials.
Jos£ Represa was bom in
Valladolid, Spain, in 1953. He
received the Licenciado degree
in Physics in 1976, and the
PhD degree in 1984, both from
the University of Valladolid,
Spain. He was Assistant
usfrt^-aMHTi
Professor from 1976 to 1985,
and since 1985 Professor in Electromagnetics at the
University of Valladolid. His current research
interest includes
numerical
methods
in
electromagnetics,
characterization
of
electromagnetics properties of materials and
microwave devices.
218
2003 INSTITUTIONAL MEMBERS
AUSTRALIAN DEFENCE LIBRARY
Northcott Drive
Campbell, A.C.T. 2600 AUSTRALIA
EDINBURGH DSTO
PO Box 830673
Birmingham, AL 35283-0673
HRLLABS, RESEARCH LIBRARY
3011 Malibu Canyon
Malibu.CA 90265
BAE SOWERBY RESEARCH CTR.
FPC267PoBox5
Filton. Bristol, BS134 7QW
UNITED KINGDOM
ELEC. COMM. LAB LIBRARY
Hikarinooka Yokosuka Shi
Kanagawa-Ken
239-0847 MZ JAPAN
lEE INSPEC
Acquisitions Section
Michael Faraday House
6 Hills Way
Stevenage, Herts UK SGI 2AY
BAE SYSTEMS
W423A, Warton Aerodrome
Preston, PR41AY
UNITED KINGDOM
EM SOFTWARE & SYSTEMS
PO Box 1354
Stellenbosch, S. AFRICA 7599
IIT RESEARCH INSTITUTE
185 Admiral Cochrane Drive
Annapolis, MD 21401-7396
BAE SYSTEMS, ADV. TECH. CTR
W. Hanningfield Rd, Tech Centre Lib.
Gt. Baddow, Chelms., CM2 8HN UK
EMBRAER SA
Av Brig. Faruna Lima 2170
Sao Paulo, BRAZIL 12227-901
IND CANTABRIA
PO Box 830470
Birmingham, AL 35283
BEIJING BOOK COMPANY, INC
701 E Lindon Ave.
Linden, NJ 07036-2495
ENGINEERING INFORMATION, INC
PO Box 543
Amsterdam, Netherlands 1000 Am
INSTITUTE FOR SCIENTIFIC INFO.
Publication Processing Dept.
3501 Market St.
Philadelphia, PA 19104-3302
CENTER HP COMPUTING
PO Box 830657
Birmingham, AL 35283-0657
ERA TECHNOLOGY, LTD
Cleeve Rd
Leatherhead, Surrey, UK OX5 1JE
INSYS LTD.
Reddings Wood, Ampthill
Bedford, MK45 2HD UK
CULHAM EM
BIdg F5, Culham
Abingdon, Oxford, 0X14 3ED UK
ETSE TELECOMUNICACION
Biblioteca, Campus Lagoas
Vigo, 36200 SPAIN
IPS RADIO & SPACE SVC/LIBRARY
PO Box 5606
W. Chatswood, 2057 AUSTRALIA
DARTMOUTH COLL-FELDBERG LIB
6193 Murdough Center
Hanover, NH 03755-3560
FANFIELD LTD
Braxted Park
WItham, Essex, CM83XB UK
LAND ENGINEERING CO.
661 Bourke St.
Melbourne, AUSTRALIA 3000
DARMSTADT U. of TECHNOLOGY
Schlossgartenstrasse 8
Darmstadt, Hessen
GERMANY D-64289
FLORIDA INTERNATIONAL U.
10555 W.FIagler St.
ECE Dept., EAS-3983
Miami, FL 33174
LEMA-EPFL
ELB-ECUBLENS
Lausanne, VD
Switzerland CH-1015
DEFENCE RESEARCH ESTAB. LIB.
3701 Carling Avenue
Ottowa, ON, K1A 0Z4 CANADA
GEORGIA TECH LIBRARY
225 North Avenue, NW
Atlanta, GA 30332-0001
LIBRARY of CONGRESS
Reg. Of Copyrights
Attn: 40T Deposits
Washington DC, 20559
DEPARTMENT OF PHYSICS
Faculty of Sciences of Tunis
El Manar, TUNISIA 2092
HANGANG UNIVERSITY / IS4.-tB.
17 Haengdang-Dong, Seongdong KU
Seoul, S KOREA 133-791
LINDA HALL LIBRARY
5109 Cherry Street
Kansas City, MO 64110-2498
DTIC-OCP/IBRARY
8725 John J. Kingman Rd. Ste 0944
Ft. Belvoir,VA 22060-6218
HKUST, UNIVERSITY LIBRARY
Clear Water Bay Road
Kowloon, HONG KONG
MAAS
16 Peachfield Road
Great Malvem, Wore.
WR144AP UNITED KINGDOM
ECI/NCHC
12 Fl, No. 82, Sec. 2
Fu Shing South Road
Taipei, TAIWAN 106
HOKKAIDO UNIVERSITY
Kita13,Nishi8, KITA-KU
Sapporo, JAPAN, 060-8628
MISSISSIPPI STATE UNIV LIBRARY
PO Box 9570
Mississippi State, MS 39762
NAKANSAI KINOKUNNA CO.
Attn: M. MIYOSHI
POBox36(NDLAKANSAI)
Hongo, Tokyo, JAPAN 113-8688
QUEEN MARY & WESTFIELD
COLLEGE
Mile End Rd.
London E1 4NS UK
TELSTRA RES. LABS LIBRARY
770 Blackburn Road
Clayton, VIC. 3168 AUSTRALIA
NATIONAL INSTITUTE OF AIST
AIST Tsukuba Central 2-1-1-1
Umezono
Tsukuba-Shi, Ibariki
305-8568 JAPAN
QUEENSLAND CENTER, LIBRARY
2643 Moggill Road
Brisbane, QLD, 4069 AUSTRALIA
TOKYO KOKAUNIV
1404-1 Katakura-Cho
Hachioji, Tokyo, JAPAN 192-0914
NATL GROUND INTELLIGENCE
220 7" Street N.E.
Charlottesville, VA 22902-5307
RENTON TECH LIBRARY/BOEING
PO BOX 3707
SEATTLE, WA 98124-2207
UNIV OF CENTRAL FLORIDA LIB.
PO Box 162666
Orlando. FL 32816-2440
NATL RADIOLOGICAL PROT. BD.
Chilton. Didcot, OXON,
0X11 ORG UK
ABDULSHUKKOOR
Sheikh. Zayed Road, Apt. #1208
PO Box 72409
Dubai. UNITED ARAB EMIRATES
UNIV OF COLORADO LIBRARY
Campus Box 184
Boulder. CO 80309-0184
NAVAL POSTGRADUATE SCHOOL
Attn: J. Rozdal
411 Dyer Rd./Rm 111
Monterey. CA 93943-5101
SOUTHWEST RESEARCH INST.
6220 Culebra Road
San Antonio, TX 78238
UNIV OF MISSOURI-ROLLA LIB.
1870 Miner Circle
Rolla. MO 65409-0001
NAVAL RESEARCH LABORATORY
C. Office. 4555 Overlook Avenue. SW
Washington, DC 20375
SWEDISH DEFENCE RESEARCH
AGENCY
P 0 80x1165
Linkoping, SWEDEN S-58111
US ARMY COLD REGION RES, LIB.
72 Lyme Road
Hanover. NH 03775-1290
OVIEDO LIBRARY
LRC 04440
PO Box 830679
Birmingham. AL 35283
SWETS BLACKWELL
440 Creamery Way, Suite A
Exton, PA 19341
RUBENS VALERIO
2129 Flalbush Ave.. #4000
Brooklyn.NY 11234-4336
PENN STATE UNIVERSITY
126 Patemo Library
University Park. PA 16802-1808
TECHNISCHE INFORMATIONS
BIBLIOTECH.
1-1-2 Zeitschriftenzuzgang
Welfelgarten IB
Hannover, GERMANY 30167
VECTOR FIELDS. LTD
24 Bankside
Kidlington. Oxford 0X51JE
PHILIPS RESEARCH LAB LIBRARY
Cross Oak Lane, Salfords
Redhill. RH1 5HA SURREY UK
TECHNISCHE UNIV. DELFT
Mekelweg4, Delft, Holland. 2628 CD
NETHERLANDS
ACES COPYRIGHT FORM
This form is intended for original, previously unpublished manuscripts submitted to ACES periodicals and conference publications. The signed forni,
appropriately completed, MUST ACCOMPANY any paper in order to be published by ACES. PLEASE READ REVERSE SIDE OF IBIS FORM FOR
FURTHER DETAILS.
TITLE OF PAPER:
RETURN FORMTO:
Dr. Atef Z. Elsherbeni
University of Mississippi
Dept. ofElectrical Engineering
Andeison Hall Box 13
University, MS 38677 USA
AUTHORS(S)
PUBLICATION TITLE/DATE:
PART A - COPYRIGHT TRANSFER FORM
(NOTE: Conqjany or other forms may not be substituted for this form. U.S. Government employees whose work is not subject to copyright may so certify
by signing Part B below. Authas whose work B subject to Crown Copyright may sign Part C overieaf).
The undeisigned, desiring to publish the above paper in a publication of ACES, hereby transfer their copyrights in the above paper to The Applied
Computational Electromagnetcs Society (ACES). The undersigned hereby represents and warrants that the paper is original and that he/she is the author
of the paper or otherwise has the power and aufcority to make and execute thk assignment.
Returned Rights: In return for these rights, ACES hereby grants to the above authors, and the employers for whom the work was performed, royalty-ftee
permission to:
1. Retain all proprietary rights otherthan copyright, such as patent lights.
2. Reuse all or portions of the above paper in other wodcs.
3. Reproduce, or have reproduced, the above paper for the authcr's personal use or for internal company use provided that (a) the source and
ACES copyright are indicated, (b) the copies are not used in a way that implies ACES endorsement of a product or service of an employer, and (c) the
copies per se are not offered for sale.
4. Make limited distribution of all or portions of the above paper prior to publication.
5. In the case of work perfonned under U.S. Government contract, ACES grants the U.S. Government royalty-free peimission to reproduce all
or portions of the above paper, andto authorize others to do so, for U.S. Government puiposes only.
ACES Obligations: In exercising its rights under copyright, ACES will make all reasonable efforts to act in the interests of the authois and employeis as
well as in its own interest. In particular, ACES REQUIRES that:
1. The consent ofthe first-named author be sought asa condition in grantbg re-publication pennission to otters.
2. The consent ofthe undeisigned employer be obtained as a condition in granting permission to others to reuse all orpoitions ofthe paper forpromotion
or mariceting puiposes.
In the eventthe above paper is not accepted and published by ACES or is withdrawn by the author(s) before accqitance by ACES, this agreement becomes
null and void.
AUTHORIZED SIGNATURE
TITLE (IF NOT AUTHOR)
EMPLOYER FOR WHOM WORK WAS PERFORMED
DATE FORM SIGNER
Part B - U.S. GOVERNMENT EMPLOYEE. CERTIFICATION
(NOTE: if your work was performed under Government contract but you are not a Government employee, sign transferfoim above and see item 5 under
Returned Rights).
This certifies that all authors ofthe above paper are employees of the U.S. Government and performed this work as part of their employment and that the
paper is theiefor not subject to U.& copyright protection.
AUTHORIZED SIGNATURE
TITLE (IF NOT AUTHOR)
NAME OF GOVERNMENT ORGANIZATION
DATE FORM SIGNED
PART C - CROWN COPYRIGHT
(NOTE: ACES recognizes and will honor Crown Copyright as it does U.S. Copyright, It is understood that, in asserting Crown Copyright, ACES in no
way diminishes its rights as publisher. Sign only HALL authors are subject to Crown Copyright),
This certifies that all authors of the above Paper are subject to Crown Copyright. (Appropriate documentation and instructions regarding form of Crown
Copyright notice may bcattached).
AUTHORIZED SIGNATURE
TITLE OF SIGNEE
NAME OF GOVERNMENT BRANCH
DATE FORM SIGNED
Information to Authors
ACES POLICY
ACES distributes its technical publications throughout the world, and it may be necessary to translate and abstract its publications, and articles contained
therein, for inclusion in variais compcndiums and similar publications, etc. When an article is submitted for publication by ACES, acceptanceof the
article implies that ACES has the rights to do all of the things it nomially does with such an article.
In connection with its publishingactivities, it is the policy of ACES to own the copyrights in its technical publications.and to the contributions contained
therein, in order to protect the interests of ACES, its authors and their employers, and at the same time to facilitate the appropriate re-use of this material
by others.
The new United States copyright law requires that the transfer of copyrights in each contribution from the author to ACES be confirmed in writing. It is
therefore necessary that you execute either Part A-Copyright Transfer Form or Part B-U.S. Government Employee Certification or Part C-Crown
Copyright on this sheet and return it to the Managing Editor (or person who supplied this sheet) as pronptly as possible.
CLEARANCE OF PAPERS
ACES must of ncces.sity assume that materials presented at its meetings or submitted to its publicatbns is properly available for general dissemination to
the audiences these activities arc organized to serve. It is the responsibility of the authors, not ACES, to determine whether disclosure of their material
requires the prior consent of other parties and if so, to obtain it. Furthemiore, ACES must assume that, if an authoruses within his/herarticlepreviously
published and/or copyrighted material that permission has been obtained for such use and that any required credit lines, copyright notices, etc. are duly
noted.
AUTHOR/COMPANY RIGHTS
If you arc employed and you prepared your paperas a part of your job, the rights to your paper initially rest with your employer. In that case, when you
sign the copyright form, we assume you are authorized to do so by your employer and that your employer has consented to all of the terms and conditions
of this fonu. If not, it should be signed bysomeonc so authorized.
NOTE RE RETURNED RIGHTS; Just as ACES now requires a signed copyright transfer form in order to do "business as usual", it is the
intent of this form to return rights to the author and employer so that they too may do "business as usual". If further clarification is required, please
contact: Tlic Managing Editor, R. W. Adler, Naval Postgraduate School, Code BC/AB, Monterey, CA, 93943, USA (408)556-2352,
Please note that, although authors are permitted to re-use all or portions of their ACES copyrighted material in other wprks, this does not include granting
third party requests for reprinting, republishing or other types of re-use.
JOINT AUTHORSHIP
For jointly authored papers, only one signature is required, but we assiime all authois havebeen advised and have consentedto the terms of this form.
U.S. GOVERNMENT EMPLOYEES
Authors who arc U.S. Government employees arc not required to sign the Copyright Transfer Form (Part A), but any co-authors outside the Government
arc.
Part B of the fonn is to be used insteadof Part A only if all authois are U.S. Government employees and prepared the paper aspart of their job.
NOTE RE GOVERNMENT CONTRACT WORK: Authors whose work was performed under a U.S. Government contract but who are not
Government employees are required so sign Part A-Copyright Transfer Form. However, item 5 of the form returns reproduction rights to the U. S.
Government when required, even ihcugh ACES copyright policy is in effect with respect to the reuse of material by the general public.
January 2002
APPLIED COMPUTATIONAL ELECTROMAGNETICS SOCIETY JOURNAL
http://aces.ee.oIemiss.edu
INFORMATION FOR AUTHORS
PUBLICATICXV CRITERIA
Each paper is required to manifest some relation to applied
computational electromagnetics.
Papers may address
general issues in applied computatianal electromagnetics,
or they may focus on specific applications, techniqcs,
codes, or computational issues. V/hile the following list is
not exhaustive, each paper will generally relate to at least one
of these areas:
1
Code validation. This is done using internal checks or
experimental, analytical or other computational data.
Measured data of potential utility to code validation
efforts will also be consideicd for publication.
2 Code performance analysis.
This usually involves
identification of numerical accuracy or other limitations,
solution convergence, numerical and physical modeling
error, and parameter tradeoffe. However, it is also
peraiissible to address issues such as ease-of-use, set-up
time, run time, special outputs, or other special features.
3
Computational studies of basic physics. This involves
using a code, algorithm, or computational technique to
simulate reality in such a way that better, or new
physical insight or understanding is achieved. -
i New computational techniqcs , or new applications for
existing computational techniques orcodes.
S Tricks of the trade"
and techniques.
6
in selecting and applying codes
New codcs,algorithms,code enhancement,and code
fixes. This category is self-explanato^', but includes
significant changes to existing codes, such as
applicability extensions, algoridim optimization, problem
correction, limitation removal, or other performance
improvement. Note: Code (fr algorithm) capability
descriptions are not acceptable, unless they contain
sufficient technical material to justify consideration.
Applications of interest include, but are not limited to,
antennas (and their electromagnetic environments), networks,
static fields, radar cross section, shielding, radiation hazards,
biological
effects,
electromagnetic
pulse
(BMP),
electromagnetic interference
(EMI),
electromagnetic
compatibility (EMC), power transmission, charge transport,
dielectric, magnetic and nonlinear materials, microwave
components, MEMS technology, MMIC technology, remote
sensing and geometrical and physical optics, radar and
communications systems, fiber optics, plasmas, particle
accelerators, generators and motors, electromagnetic wave
propagation, non-destructive evaluation, eddy currents, and
inverse scatteririg.
Techniques of interest include frequency-domain and timedomain techniques, integral equation and diftrential equation
techniques, diffraction theories, physical optics, moment
methods, finite differences and finite element techniques,
modal expansions, perturbation methods, and hybrid methods.
This list is not exhaustive.
A unique feature of the Journal is the publication of
unsuccessful
efforts
in
applied
computational
electromagnetics. Publication of such material provides a
means to discuss problem areas in electromagnetic modeling.
Material representing an unsuccessful application or negative
results in computational electromgnetics will be considered
for publication only if a reasonable expectation of success
(and a reasonable effort) are reflected. Moreover, such
material must represent a problem area of potential interest to
the ACES membership.
Where possible and appropriate, authors are required to
provide statements of quantitative accuracy for measured
and/or computed data. This issue is discussed in "Accuracy
& Publication: Requiring, quantitative accuracy statements to
accompany data," by E. K. Miller, ACES Newsletter, Vol. 9,
No. 3, pp 23-29,1994, ISBN 1056-9170.
EDITORIAL REVIEW
7
8
Code input/output issues.
This normally involves
innovations in input (such as input geometry
standardization, automatic mesh generation, or
computer-aided design) or in output (whether it be
tabular, graphical, statistical, Fourier-transformed, or
otherwise signal-processed).
Material dealing vnth
input/output database management, output interpretation,
or other input/output issues will also be consideied for
publication.
Computer hardware issues. This is the categoiy for
analysis of hardware capabilities and limitations of
various types of electromagnetics computational
requirements. Vector and parallel computational
techniques and implementation are of particular interest.
In order to ensure an appropriate level of qality control ,
papers are peer reviewed. They are reviewed both for
technical con'ectness and for adherence to the listed
guidelines regarding information content.
JOURNAL CAMERA-READY SUBMISSION DATES
March issue
July issue
November issue
deadline 8 January
deadline 20 May
deadline 20 September
Uploading an acceptable camera-ready article after the
deadlines will result in a del^f in publishing this article.
The ACES Journal is flexible, within reason, in reganJ to
style. However, certain cquirements arc in effect:
responsibility to provide acceptable camera-ready pdf files.
Incompatible or incomplete pdf files will not be processed,
and authore will be requested to re-upload a revised
acceptable version.
1.
SUBMITTAL PROCEDURE
STYLE FOR CAMERA-READY COPY
2.
The paper title shoul'd NOT be placed on a separate page.
The title, author(s), abstract, and (space permitting)
beginning of the paper itself should all be on the first
page. The title, author(s), and author affiliations should
be centered (ccntir-justificd) on the first page.
An abstract is REQUIRED. The abstract should be a
brief summary of the work described in the paper. It
should state the computer codes, computational
lechniques, and applicatfons discussed in the paper (as
applicable) and should othenvise be usable by technical
abstracting and indexing services.
3.
Either British English or American English spellings
may be used, provided that each word is spelled
consistently throughout the paper.
4.
Any commonly-accepted format for referencii^ is
permitted, provided that internal consistency of format is
maintained. As a guideline for authors who have no
other preference, we fecomrmnd that references be given
by author(s) name and year in the body of the paper
(with alphabetical listing of all references at the end of
the paper). Titles of Journals, monographs, and similar
publications should be in italic font or should be
underlined. Titles of papers or articles should be in
quotation marks.
5.
Internal coasistency shall also be maintained for other
elements of style, such as equation numbering. As a
guideline for authors who have no other preference, we
suggest that equation numbers be placed in parentheses
at the right column margin.
6.
The intent and meaning of all text must be clear. For
authors who arc NOT masters of the English language,
the ACES Editorial Staff will provide assistance with
gramiTiar (subject to clarity of intent and ircaning).
7.
Unused space should be minimized. Sections and
subsections should not normally begin on a new page.
PAPER FORMAT
The preferred format for initial submission and camera-ready
manuscripts is 12 point Times Roman font, single line spacing
and double column format, similar to that used here, with top,
bottom, left, and right 1 inch margins. Manuscripts should be
prepared on standard S.Stl I inch paper
Only camera-ready electronic files are accepted for
publication. The term Camcra-rcady"mcans that the
material is ncat,legible,and reproducible. Full details can
be found on ACES site. Journal section.
ACES reserves the right to edit any uploaded material,
however, this is not generally done. It is the author(s)
All submissions should be uploaded to ACES server through
ACES web site (http://aces.ee.olemiss.edu) by using the
upload button, journal section. Only pdf files are accepted for
submission. The file size should not be larger than 5MB,
otherwise permission from the Editor-in-Chief should be
obtained first. The Editor-in-Chief will acknowledge the
electronic submission after the upload process is successfully
completed.
COPYRIGHTS AND RELEASES
Each primary author must sign a copyright form and obtain a
release from his/her organization vesting the copyright with
ACES. Copyright forms are available at ACES, web site
(http://aces.ee.olemBS.edu). To shorten the review process
time, the executed copyright form should be forwarded to the
Editor-in-Chief immediately after the completion of the
upload (electronic submission) process. Both the author and
his/her organization are allowed to use the copyrighted
material freely for their own private purposes.
Permission is granted to quote short passages and reproduce
figures and tables from and ACES Journal issue provided the
source is cited. Copies of ACES Journal articles may be
made in accordance with usage permitted by Sections 107 or
108 of the U.S. Copyright Law. This consent does not extend
to other kinds of copying, such as for general distribution, for
advertising- or promotional purposes, for creating new
collective works, or for resale. The reproduction of multiple
copies and the use of articles or extracts for commercial
puiposes requins the consent of the author and specific
permission from ACES. Institutional members are allowed to
copy any ACES Journal issue for their internal distribution
only.
PUBLICATION CHARGES
ACES members are allowed 12 printed pages per paper
without charge; non-membeis are allowed 8 printed pages per
paper without charge. Mandatory page charges df $75 a page
apply to all pages in excess of 12 fotmembers or 8 for nonmembeis. Voluntary page charges are requested for the free
(12 or 8) pages, but are NOT mandatoiy or required for
publication. A priority couitesy guideline, which favors
members, applies to paper backlog. Authoi^ are entitled to
15 free reprints of their articles and must request these from
the Managing Editor. Additional reprints are available to
authors, and reprints available to non-authois, for a nominal
fee.
ACES Journal is abstracted in INSPEC.in Engineering
Indcx,and in DTIC.