close

Вход

Забыли?

вход по аккаунту

?

New methods for time -domain modelling of RF/microwave passive structures and active devices

код для вставкиСкачать
New Methods for Time-Domain Modelling of
RF/Microwave
Passive Structures and Active Devices
by
S huiping Luo
Subm itted in partial fulfillm ent of the requirem ents
for the degree of
DOCTOR OF PHILOSOPHY
Major Subject: Electrical and Computer Engineering
at
D alhousie U niversity
H alifax, N ova Scotia
January 2007
© C opyright by Shuiping Luo, 2007
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
Library and
Archives Canada
Bibliotheque et
Archives Canada
Published Heritage
Branch
Direction du
Patrimoine de I'edition
395 W ellington Street
Ottawa ON K1A 0N4
Canada
395, rue W ellington
Ottawa ON K1A 0N4
Canada
Your file Votre reference
ISBN: 978-0-494-27644-0
Our file Notre reference
ISBN: 978-0-494-27644-0
NOTICE:
The author has granted a non­
exclusive license allowing Library
and Archives Canada to reproduce,
publish, archive, preserve, conserve,
communicate to the public by
telecommunication or on the Internet,
loan, distribute and sell theses
worldwide, for commercial or non­
commercial purposes, in microform,
paper, electronic and/or any other
formats.
AVIS:
L'auteur a accorde une licence non exclusive
permettant a la Bibliotheque et Archives
Canada de reproduire, publier, archiver,
sauvegarder, conserver, transmettre au public
par telecommunication ou par I'lnternet, preter,
distribuer et vendre des theses partout dans
le monde, a des fins commerciales ou autres,
sur support microforme, papier, electronique
et/ou autres formats.
The author retains copyright
ownership and moral rights in
this thesis. Neither the thesis
nor substantial extracts from it
may be printed or otherwise
reproduced without the author's
permission.
L'auteur conserve la propriete du droit d'auteur
et des droits moraux qui protege cette these.
Ni la these ni des extraits substantiels de
celle-ci ne doivent etre imprimes ou autrement
reproduits sans son autorisation.
In compliance with the Canadian
Privacy Act some supporting
forms may have been removed
from this thesis.
Conformement a la loi canadienne
sur la protection de la vie privee,
quelques formulaires secondaires
ont ete enleves de cette these.
While these forms may be included
in the document page count,
their removal does not represent
any loss of content from the
thesis.
Bien que ces formulaires
aient inclus dans la pagination,
il n'y aura aucun contenu manquant.
i
*i
Canada
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
D A L H O U S IE U N IV E R S IT Y
To co m p ly w ith the C a n a d ia n P riv acy A ct the N a tio n a l L ibrary o f C anada has requested
th at the fo llow ing p a g e s b e rem o v ed from this c o p y o f the thesis:
P re lim in a ry P ages
E xam iners S ignature P age
D alhousie L ibrary C opyright A greem ent
A p p en d ices
C opyright R eleases ( if applicable)
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
Table of Contents
1
2
3
4
5
6
List of F ig u res.................................................................................................................................vi
List o f T a b le s ............................................................................................................................... viii
List o f A b b rev iatio n s....................................................................................................................ix
D ed ication........................................................................................................................................xi
A ck now ledgem ent....................................................................................................................... xii
A bstract.......................................................................................................................................... xiii
I n tr o d u c ti o n ....................................................................................................................................1
1.1
Frequency-dom ain m e th o d s ........................................................................................ 2
1.1.1
M ethod o f m om ents........................................................................................... 2
1.1.2
Finite elem ent m e th o d .......................................................................................4
1.1.3
Finite difference m eth o d ................................................................................... 5
1.2
T im e-dom ain m eth o d s................................................................................................... 6
1.2.1
F D T D ..................................................................................................................... 7
1.2.2
T L M ........................................................................................................................8
1.3
M otivation and contributions........................................................................................9
C o n c e p t o f F D T D ........................................................................................................................12
2.1
Introduction o f classical F D T D ................................................................................. 12
2.1.1
Y ee’s grid and FD TD fo rm u la....................................................................... 12
2.1.2
N um erical s ta b ility ...........................................................................................15
2.1.3
N um erical d isp ersio n ....................................................................................... 16
2.2
L um ped Elem ents in F D T D ....................................................................................... 17
2.3
C o n c lu sio n ...................................................................................................................... 20
N ew C o m p a c t 2D F D T D M e th o d f o r 3D W a v eg u id e S tr u c tu r e s ..............................21
3.1
In tro d u c tio n ....................................................................................................................21
3.2
Form ula o f the New C om pact 2D F D T D ............................................................... 22
3.3
N um erical V a lid a tio n .................................................................................................. 25
3.4
D iscussion and C onclusion ......................................................................................... 28
C o m p a c t I D F D T D M e th o d f o r W a v eg u id e S t r u c t u r e s ...............................................29
4.1
In tro d u c tio n .................................................................................................................... 29
4.2
Form ula o f the Proposed ID FD TD M eth o d ..........................................................30
4.3
N um erical D ispersion of the Proposed M e th o d ....................................................33
4.4
A pplications o f the Proposed M e th o d ..................................................................... 36
Efficient generation o f an incident w a v e ..................................................................... 36
Absorbing boundary condition.......................................................................................37
4.5
N um erical E xam ples..................................................................................................... 38
4.6
D iscussion and C onclusion......................................................................................... 43
F D T D -B ased M o d a l P M L ........................................................................................................ 44
5.1
In tro d u c tio n .................................................................................................................... 44
5.2
D erivation o f the Proposed ID M odal P M L ..........................................................44
5.3
N um erical E xam ples..................................................................................................... 46
5.4
D iscussion and C onclusion ......................................................................................... 49
A N ew S u b g rid d in g M e th o d fo r 2D F D T D ........................................................................50
6.1
In tro d u c tio n .....................................................................................................................50
6.2
D erivation o f the Stable Subgridding S c h e m e .......................................................50
iv
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
6.3
6.4
7
N um erical E x a m p le ..................................................................................................... 59
C o n clu sio n ......................................................................................................................62
Extraction of Causal Time-Domain Parameters Using Iterative Methods
63
7.1
In tro d u c tio n ................................................................................................................... 63
7.2
The E rror Feedback B ased FFT M e th o d ................................................................65
7.3
The H ilbert T ransform B ased FFT M ethod ............................................................ 68
7.4
N um erical V a lid a tio n ..................................................................................................73
A) Triangular p u lse ...................................................................................................................... 73
B) R ectangular p u ls e ...................................................................................................................75
C) FET a m p lifie r..........................................................................................................................77
7.5
D iscussion and C onclusion .........................................................................................82
8 Extraction of Causal Time-Domain Parameters Using Rational Function
Approximation.................................................................................................................... 83
8.1
8.2
8.3
8.4
9
In tro d u c tio n ....................................................................................................................83
Causal T im e-D om ain Extraction W ith The R ational Fitting T echnology ...83
N um erical V a lid a tio n .................................................................................................. 88
C o n clu sio n ...................................................................................................................... 94
Conclusion and Future Work.................................................................................... 95
9.1
9.2
S um m ary..........................................................................................................................95
Future W ork ....................................................................................................................96
References:.......................................................................................................................... 97
Appendix A: The Numerical Dispersion of Compact ID FDTD for TEmn mode in a
Rectangular Waveguide...................................................................................................110
Appendix B: The Formula of ID Modal CFS-PML.................................................... 114
v
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
List of Figures
Figure 2.1 Positions of the field com ponents in Y ee’s grids..................................................... 15
Figure 2.2 T he position o f a one-port lum p-elem ent device in the FD TD g rid .................... 18
Figure 2.3 A tw o-port lum p-elem ent device (transistor) in the FD TD g rid ...........................18
Figure 2.4 T he equivalent tw o-port netw ork of an active d evice............................................. 19
Figure 3.1 T he electric and m agnetic fields positions in a unit cell o f the com pact 2D
FD TD grid...............................................................................................................................................24
Figure 3.2 T he rectangular w aveguide cavity o f length d=0.04m in the z direction, width
a=0.03m in the x direction, and height b=0.02m in the y direction.........................................25
Figure 3.3 T he rectangular w aveguide cavity w ith a groove in the m iddle o f the z
direction...................................................................................................................................................27
Figure 4.1 T he electric and m agnetic field positions in Y ee’s lattice for a w aveguide
structure, with Z the w ave propagation direction......................................................................... 31
Figure 4.2 T he E y recorded at a point 300Az or 602 aw ay from the source p lan e ..............35
Figure 4.3 T he difference betw een the two E ys produced by the ID FD TD and the
reference 3D F D T D .............................................................................................................................. 36
Figure 4.4 T he proposed absorbing term ination using the ID FD TD m esh [76][77]........ 37
Figure 4.5 T he m ode extraction and com bination diagram at the interface betw een the
waveguide and the ID FD TD absorbing term ination [76][77].................................................38
Figure 4.6 T he Ey values of the first 0.5 ns recorded at a point lAz aw ay from the source.
...................................................................................................................................................................39
Figure 4.7 T he relative errors of Ey at a point lAz away from the source obtained by the
proposed ID F D T D .............................................................................................................................. 40
Figure 4.8 T he reflection coefficient from ID FD TD term ination for T E n m ode............ 41
Figure 4.9 T he reflection coefficient from ID FD TD term ination for T E 84 m ode..............42
Figure 4.10 T he reflection coefficient from ID FD TD term ination for m ulti-m ode case.42
Figure 5.1 T he difference betw een the e obtained by the proposed ID m odal C F S-PM L
and that obtained by the original 3D C F S-PM L for TEio m o d e...............................................47
Figure 5.2 T he difference betw een the e s obtained by the proposed ID m odal CFSPM L and that obtained by the original 3D C FS-PM L for T E u ............................................... 47
Figure 5.3 The w aveguide structure with two strips under study............................................ 48
Figure 5.4 T he reflected Ev obtained by the proposed ID C FS-PM L for the first 9 m odes
and the original C F S-P M L ................................................................................................................. 49
Figure 5.5 T he com puted reflection coefficients by the PM L in the corresponding
frequency-dom ain................................................................................................................................. 49
Figure 6.1 A subgridding structure with a m esh ratio o f 3:1 in a FD TD lattice.................. 58
Figure 6.2 T he parallel plate w aveguide w ith w idth o f 10mm and coarse m esh size of
Ajc=lmm and Ay=lmm............................................................................................................................. 59
Figure 6.3 T he reflection coefficient from the subgridded m eshes in the center
corresponding to the parallel plate w aveguide in Figure 6 .2 .....................................................60
Figure 6.4 The rectangular resonator of 6mmx5mm with a fin o f length 2mm in the m iddle.
61
vi
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
Figure 7.1 T he tim e-dom ain yn (r) obtained with the direct inverse Fourier transform of
the frequency-dom ain Y -param eters............................................................................................... 64
Figure 7.2 T he flow chart o f the error feed-back based FFT m ethod..................................... 67
Figure 7.3 T he flow chart o f the H ilbert transform based FFT m ethod...................................72
Figure 7.4 T he tim e-dom ain pulse and its Fourier transform obtained with the direct cut­
o ff m ethod............................................................................................................................................... 74
Figure 7.5 T he tim e-dom ain w aveform s and its spectra extracted w ith the proposed two
iterative m ethods w hen the know n frequency range is from 0 to 6G H z for the triangular
p u lse ..........................................................................................................................................................74
Figure 7.6 T he tim e-dom ain w aveform and its spectra extracted with the direct cut-off
approach.................................................................................................................................................. 76
Figure 7.7 T he tim e-dom ain w aveform s and their spectra extracted w ith the tw o proposed
iterative m ethod m ethods when the know n frequency range o f interest is from 0H z 6G H z for the rectangular p u lse......................................................................................................... 77
Figure 7.8 T he tim e-dom ain Y -param eters extracted with the error feedback FFT m ethod
for the FET used in the am plifier......................................................................................................78
Figure 7.9 T he tim e-dom ain Y -param eters extracted with the H ilbert transform based
FFT m ethod for the FET used in the am plifier............................................................................. 79
Figure 7.10 T he Y -param eters in frequency-dom ain after using the iteration m ethods....80
Figure 7.11 T he S-param eters after using the iteration m ethods and direct cut-off
81
Figure 7.12 T he Layout o f the FET am plifier circuit.................................................................. 81
Figure 7.13 T he com puted S21 of the am plifier in Figure 7.10 from different m ethods. .82
Figure 8.1 T he Y -param eters in the tim e-dom ain corresponding to T able 8-2 and (8.13)
but w ithout the aS(t) term ................................................................................................................. 90
Figure 8.2 T he Y -param eters in the frequency-dom ain obtained with the pro p o sed
91
Figure 8.3 S-param eters in the frequency-dom ain obtained w ith the proposed m ethod and
direct cut-off...........................................................................................................................................92
Figure 8.4 T he overall S21 of the am plifier...................................................................................93
Figure 8.5 T he com putation tim e versus the num ber o f iterations for the am plifier
94
vii
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
List of Tables
T able 3-1 The com parison betw een num erically com puted and analytical resonant
frequencies..............................................................................................................................................26
T able 3-2 The m em ory and CPU tim e needed for the 2D FD TD and 3D FD TD
corresponding to the first e x a m p le .................................................................................................. 26
Table 3-3 The com parison o f resonant frequencies obtained by 2D FD T D and the 3D
FD TD in the second e x a m p le ............................................................................................................27
Table 3-4 The m em ory and CPU tim e needed for the 2D FD TD and 3D FD T D in the
second exam ple..................................................................................................................................... 27
T able 4-1 The m em ory and CPU tim e used by the proposed ID m ethod and the
referenced 3D FD TD m eth o d ............................................................................................................40
T able 6-1 T he values o f a for subgridding ratios K = 3 :l, 5:1, and 7 :1 ...................................58
T able 6-2 The relative errors o f the com puted resonant frequency o f the fin structure. ...61
T able 8-1 The initial poles (Hz) of Y 1 1, Y12, Y 21, and Y 22 as in (8 .7 ).............................. 89
T able 8-2 T he final poles and residuals o f Y 1 1, Y12, Y21, and Y22 as in (8 .1 2 )..............89
Table 8-3 The constant term o f Y l l , Y 12, Y21, and Y 22 as in (8 .1 2 ).................................. 89
viii
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
List of Abbreviations
ABC
A bsorbing B oundary Condition
CAD
C om puter A ided D esign
ADI
A lternating D irection Im plicit
ADS
A dvanced D esign System
CHE
C om bined Field Integral Equation
CFL
C ourant-Friedrich-Levy
CFS-PM L
C om plex Frequency Shifted Perfectly-M atched Layer
CPU
Central Processing U nit
EFIE
Electric Field Integral Equation
EM
Electrom agnetic
FD TD
Finite-D ifference T im e-D om ain
FDM
Finite D ifference M ethod
FEM
Finite E lem ent M ethod
FET
Field E ffect T ransistor
FFT
Fast Fourier Transform
IE
Integral Equation
M CD
M odified C entral D ifference
M FIE
M agnetic Field Integral Equation
MM
M om ent M ethod
M oM
M ethod o f M om ent
M RTD
M ulti-R esolution Tim e-D om ain
PM L
Perfectly M atched Layers
PSTD
Pseudo Spectral Tim e-D om ain
RA M
Random A ccess M em ory
RCS
R adar Cross Section
RF
R adio Frequency
TE
Transverse Electric
TLM
T ransm ission Line M ethod
TM
Transverse M agnetic
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
VF
V ector Fitting
ID
O ne D im ensional
2D
Tw o D im ensional
3D
Three D im ensional
x
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
Dedication
This w ork is dedicated to
m y parents, my wife, and m y son.
xi
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
Acknowledgement
First and forem ost, I w ould like to thank m y supervisor, Dr. Z hizhang (D avid) C hen for
his very valuable guidance and financial support throughout m y research. I w ish to thank
as well the other m em bers o f the G uiding C om m ittee for their help. I thank also my
colleagues in the M icrow ave and W ireless R esearch L aboratory o f D alhousie U niversity
for their helpful discussion and support. Special thanks go to M r. Yan C hen and Dr.
C hangning M a for their help during m y research.
I extend my deepest thanks and
appreciation to m y wife W eijun W u for her understanding and m oral support during my
Ph.D program .
xii
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
Abstract
This thesis focuses on im proving existing FD TD m ethods and tim e-dom ain m odeling o f
active devices.
In this thesis, four m ethods are proposed to im prove the efficiency o f FD TD . First, a
new com pact tw o dim ensional (2D) FD TD m ethod is proposed to analyze the w aveguide
discontinuities along the direction o f w ave propagation, w hich is uniform in one o f the
transverse directions; it reduces the original three-dim ensional (3D) problem s to twodim ensional (2D) problem s. Second, a com pact one-dim ensional (ID ) FD T D form ula is
proposed for uniform ly filled w aveguide; it has the sam e num erical dispersion as the
original 3D FD TD form ula and can be used as an efficient incident w ave generator and
perfect absorbing boundary condition for the single m ode. Then, a ID m odal Perfectly
M atched Layers (PM L) is developed by applying the proposed ID FD TD to the
traditional PM L form ula. Finally, a new 2D FD TD subgridding m ethod is proposed that
is not only very sim ple w ith controllable low reflections but also has been proven to be
stable.
W hen the FD TD m ethod is used to sim ulate R F/m icrow ave circuits w ith active
devices, the governing voltage-current equations of a device can be incorporated into the
FD TD m arching equations. H ow ever, the param eters o f m ost electronic com ponents are
often given or m easured in the frequency-dom ain and are band-lim ited. F or instance, Sparam eters of a field-effect transistor are usually given or obtained only in the frequencydom ain and over a lim ited frequency range or band o f interest.
O btaining a causal tim e-dom ain m odel from the band-lim ited frequency-dom ain
param eters is a challenge. In this thesis, three m ethods are proposed to solve this problem .
The first two m ethods are iterative m ethods based on Fast Fourier T ransform (FFT). One
applies the FFT in com bination with the error feedback principle and the second applies
the H ilbert transform in conjunction with FFT. The third m ethod uses the rational
function fitting technique. T he extracted tim e-dom ain param eters using the three m ethods
not only are causal but also contain alm ost the sam e frequency-dom ain inform ation as the
original param eters over the given lim ited frequency range.
xiii
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
1
1 Introduction
Com petition in science and technology today is intense, both betw een businesses and
states. Industries and governm ents are dedicating m ore and m ore resources to research
and developm ent. N ew system s are becom ing increasingly com plex at the sam e tim e as
the developm ent period is becom ing ever shorter. B ecause of the advance o f com puter
technology, num erical sim ulation has becom e a very im portant and valuable tool in
reducing the cost and tim e o f prototyping and testing. In the fields of electrom agnetic and
m icrow ave engineering, com putational electrom agnetics is very popular and provides a
practical solution to large com plex electrom agnetic system problem s. For exam ple, the
designers of high speed and very large scale integrated circuits and com plex m icrow ave
integrated circuits m ust use num erical sim ulation tools for functional verification and
signal integrity analysis before com m itm ents are m ade to m ass m anufacture the products.
Com putational
electrom agnetics
also
provides
a
pow erful
tool
for
research
in
electrom agnetic fields and m icrow ave technology theory.
The success o f com putational electrom agnetics depends on tw o factors. One is the
advance of com puter hardw are; the other is the availability of good algorithm s and
com puter softw are. D evelopm ent o f com puter hardw are is very rapid now, and is
characterized by high perform ance and low prices. T he hardw are developed provides a
pow erful
platform
for
software,
and
contributes
to
m aking
the
com putational
electrom agnetics softw are practical and w idespread. Low price and high perform ance
hardw are also enables m ore people to conduct their research on their personal com puters.
H ow ever, com puter hardw are, no m atter how advanced, is inherently lim ited in its
capability. T herefore, the availability o f efficient softw are alw ays plays a role in speeding
up com putations.
M axw ell’s equations and the related boundary conditions are often used to solve
electrom agnetic problem s directly. Com putational electrom agnetics use all kinds of
num erical m ethods to approxim ate the differential or integral M axw ell’s equations and
the related boundary conditions. The electrom agnetic problem s can be described in
frequency-dom ain or tim e-dom ain. T he num erical m ethods to solve them are often called
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
2
frequency-dom ain m ethods and tim e-dom ain m ethods.
Frequency-dom ain form ula states the electrom agnetic problem for a given frequency.
N um erical m ethods for frequency-dom ain problem s do not suffer from stability problem s
and can use non-uniform m eshes easily, so they are very efficient in solving large system
problem s w ith very fine structures and a narrow w orking frequency band. T hey are also
very good at solving system s containing frequency-dependent m aterials.
H ow ever, if nonlinear elem ents are present in the system s being m odeled or the
w orking frequency band is very wide, then frequency-dom ain m ethods are not efficient.
In these situations, tim e-dom ain m ethods provide better solutions because tim e-dom ain
m ethods can include the nonlinear elem ents directly. If using a w ideband exciting source,
one sim ulation o f a tim e-dom ain m ethod can provide w ideband frequency-dom ain results
using the Fast F ourier transform (FFT).
Today, a num ber of efficient num erical m ethods for solving electrom agnetic
problem s are available. Each m ethod has its particular advantages and disadvantages, and
is w ell-suited for a certain type of problem s but not for others. In the follow ing sections,
the frequency-dom ain m ethods are discussed first, after w hich the tim e-dom ain m ethods
are introduced.
1.1
Frequency-domain methods
Frequency-dom ain m ethods solve the electrom agnetic problem s in the frequency-dom ain
and obtain the frequency-dom ain results directly. They can be based on differential
M axw ell’s equations or related integral equations. M any kinds o f frequency-dom ain
m ethods exist. The m ost popular m ethods are M ethod of M om ents (M oM ), Finite
Elem ent M ethods (FEM ), and Finite-D ifference M ethods (FDM ).
1.1.1
Method of moments
M ethod of M om ents (M oM ) is a num erical procedure to solve a linear operator equation
by transform ing it into a system o f sim ultaneous linear algebraic equations. T he linear
operator can be differential, integral, or integro-differential. This m ethod can also be
called M om ent M ethod (M M ). It solves the original operator equations by using w eighted
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
3
residuals. H arrington has played an im portant role in popularizing this m ethod in
electrom agnetic engineering, describing it in a detailed and system atic fashion [1], As a
pow erful num erical m ethod, the M oM has been used extensively for m ore than three
decades. It is very suitable for radiating and scattering problem s [1]-[19], for exam ple,
antenna analysis, w aveguide discontinuity analysis, and R adar C ross Section (RCS)
analysis.
T he general procedure o f the m ethod o f m om ents can be described as follow s.
A ssum e that an electrom agnetic problem can be m odeled w ith the linear operator
equation below.
U = F
(U )
w here L is the linear operator, F is the know n force or source function, and 7 is the
unknow n quantities need to be solved.
T he first step is to expand 7 as a finite sum of basis functions:
M
j
=
( 1.2 )
w here 7, are the expansion coefficients that need to be determ ined and bj are the know n
basis functions pre-selected by a user.
By substituting (1.2) into (1.1), the follow ing equation is obtained:
M
F = L % J J b i = Y dJ iLbi
M
(1.3)
T he second step is to test the error of the above equation with another pre-selected M
linear independent w eighting functions vtf (j=l,2,...M). M ore specifically, the follow ing
equations are obtained by taking the inner product o f each w eighting function on both
sides of equation (1.3).
M
M
(1.4)
The above equations can be used to find the unknow n coefficient J j .
Equation (1.4) can be rew ritten in a m atrix form as:
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
4
(1.5)
where:
Z ;, =< wp L b >
( 1.6 )
Vj =< Wj,F >
(1.7)
( 1.8)
1.1.2
Finite element method
In general, M oM is very good for analyzing unbounded radiation problem s o f perfect
electric conductor configurations and hom ogeneous dielectrics. For problem s with
com plicated geom etries and m any arbitrarily shaped dielectric regions, how ever, M oM is
not very efficient due to the possible fast change o f the m edium properties. T he finite
elem ent m ethod (FEM ), by contrast, is good at m odeling inhom ogeneous dielectric
bodies as it requires the entire volum e of the solution dom ain to be m eshed. T herefore
each m esh elem ent m ay have com pletely different m aterial properties from those o f
neighboring
elem ents
and
this
m akes
FEM
excellent
at
m odeling
com plex
inhom ogeneous configurations. H ow ever, FEM m ay not m odel unbounded radiation
problem s as effectively as M oM because o f the potential errors due to the absorbing
boundaries introduced.
T he FEM has becom e one of the m ost popular frequency-dom ain m ethods in
com putational electrom agnetics [20]-[30]. O ne way to construct the FEM form ulation is
to use variational techniques and w ork by m inim izing or m axim izing an expression that is
know n to be stationary about the true solution. For exam ple, solutions o f M axw ell's
equations alw ays require that the energy w ithin a structure is m inim ized; the finite
elem ent m ethod can solve for the unknow n field quantities by m inim izing the energy
functional.
T he finite elem ent m ethod involves the subdivision o f the problem region into
subdom ains or finite elem ents and approxim ation o f the field in each elem ent in term s of
the linear com binations of basis functions. T he basis functions usually are sim ple
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
5
functions (often linear) defined over each elem ent. The elem ent m odel contains
inform ation about the geom etry, m aterial properties, excitations and boundary constraints.
For a three-dim ensional tim e-harm onic problem , the energy functional can be w ritten as:
F =
+ e\E \2
X
2
(19)
2
2co
The first tw o term s in the integrand are the energy stored in the m agnetic and electric
fields, and the third term is the energy dissipated due to conductivity.
E xpressing H in term s o f E , the functional F becom es a function o f E . T he field
region is then divided into a num ber o f small hom ogeneous elem ent areas. T he elem ents
can be small in regions w ith geom etric details and m uch larger in uniform regions. The
electric field E can be expanded as a linear com bination o f know n basis functions with
unknow n coefficients. F then becom es a functional about these unknow n coefficients.
Setting the derivative o f this functional w ith respect to unknow n coefficients to zero, a
system of equations w ith all the unknow n coefficients is obtained:
'K
j2
I'll
-Vl2
y%\
^22
~K
•
e
•
2
ym . A .
w here J i (i= l,2 ,..,n ) are the sources term , E t (i= l,2 ,..,n ) are the unknow n coefficients, and
y..(i= l,2,..,n; j= l,2 ,...,n ) are functions o f the target geom etry, m aterial properties, and
boundary constrains.
A fter solving equation (1.10), the field distribution is known and other physical
quantities o f interest can be com puted from the know n field distribution.
T he m ajor advantage o f FEM over other num erical m ethods is that the geom etric and
m aterial properties o f each elem ent can be defined independently and flexibly. This
m akes FEM com petitive in m odeling problem s w ith com plicated geom etries and m any
arbitrarily shaped dielectric regions, since FEM can use small elem ents in regions of
com plex geom etry and big elem ents in large uniform regions.
1.1.3
Finite difference method
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
6
Finite difference m ethod was one of the first num erical m ethods to be used in
com putational electrom agnetics because o f its sim plicity in concept [31]-[34], Sim ilar to
the FEM , the finite difference m ethod is based on the differential equations o f the
electrom agnetic problem s. The finite difference m ethod uses the difference operator to
approxim ate the original differential operator and convert the original differential
equations into a system o f algebraic equations. W ith the advance o f com puter technology,
the finite difference m ethod developed extensively, and has becom e one o f the num erical
m ethods used currently [35]-[39]. H ow ever, it is not as popular as the finite elem ent
m ethod for frequency-dom ain problem s. This m ay result from the fact that finite elem ent
m ethod attracted m uch attention in civil and m echanical engineering.
T he general procedure of finite difference m ethod consists of three steps:
a) D ivide the problem dom ain into grids and use the node values o f fields as
unknow ns.
b) R eplace the differential equation by the approxim ate difference equations, which
relate the field value in a node to the field values at neighboring nodes.
This
converts the original differential equation into a system o f linear algebra equations
in the unknow n field values at the discretized nodes.
c) Solve the system of linear algebraic equations with the know n boundary
conditions.
1.2
Time-domain methods
The frequency-dom ain form ulae are not good at handling nonlinear elem ents and
m aterials; they can obtain the result about one frequency point for each sim ulation and
are efficient only for narrow band applications. Tim e-dom ain m ethods, how ever, can
obtain the results of wide band frequencies w ith only one sim ulation and handle the
nonlinear elem ents and m aterials easily. Today, w ideband applications are very popular.
T im e-dom ain m ethods are very good options to obtain the w ideband results efficiently,
and they are becom ing m ore and m ore com petitive. U sually, each frequency-dom ain
m ethod has a tim e-dom ain counterpart. For exam ple, there are tim e-dom ain FE M [40]
and M oM [41] associated w ith the frequency-dom ain FEM and M oM . H ow ever, the m ost
popular tim e-dom ain m ethods are finite-difference tim e-dom ain m ethod (FD TD ) and
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
7
transm ission line m ethod (TLM ). In the FD TD and TLM , the sim ple and explicit
updating form ulations m ake the FD TD and TLM relatively easy to program . H ow ever,
the explicit tim e-dom ain updating m ethods can experience num erical stability problem s.
To ensure the tim e-dom ain com putation to be stable, the tim e step has to be sm aller than
an upper lim it w hich is related to spatial steps. For fine structures that require small
spatial steps, the tim e step has to be sm all. This will reduce the efficiency o f the explicit
tim e-dom ain m arching m ethods.
1.2.1
FDTD
The FDTD is one o f the m ost popular tim e-dom ain m ethods. It was proposed by K. S.
Yee [42] in 1966. Like the FEM , it uses the partial differential equations w hich describe
the target electrom agnetic system . U nlike the FEM , the FD TD m ethod does not use
variational concepts or w eighted residuals; instead, it directly approxim ates the tim e and
space derivatives in the tim e dependent M axw ell curl equations by sim ple 2nd order
central differences. FD TD uses a staggered grid in tim e and space, and the electric and
m agnetic fields are com puted on the staggered grid with a m arching-on-in-tim e technique.
It provides a direct solution o f M axw ell’s curl equations and is very flexible at solving
com plicated electrom agnetic problem s.
Although the FD TD m ethod was proposed in 1966 [42], it began to garner public
attention only after the 1980s. T he early FD TD m ethods needed to discretize the w hole
problem dom ain by uniform grids. For com plex practical problem s, this required a huge
am ount of com puter m em ory and CPU pow er, m ore than com puters could provide at the
tim e. From the 1980s on, tw o factors pushed forw ard the developm ent o f FD TD . One
was the sw ift advance in com puter perform ance, including CPU speed and com puter
m em ory, which provided the enabling technology and hardw are base. T he other factor
was application requirem ents. Practical electrom agnetic engineering problem s m ethods
becam e m ore and m ore com plicated; take, for exam ple, the electrom agnetic com patibility
analysis o f com plex electronic system s and the signal integrity analysis o f large-scale
integrated circuits. T hese com plex electrom agnetic problem s require flexible and
pow erful electrom agnetic sim ulation tools. As a direct solution to M axw ell’s curl
equations, using a sim ple and explicit updating form ula, the FD TD m ethod becam e one
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
o f the m ost popular num erical m ethods [43]-[54], A fter B erenger proposed the perfectly
m atched layer absorbing boundary conditions [50], the m ethod could be m ade accurate
on sm aller m eshes.
W hen using the FD TD m ethod, there are tw o constraints [55]: one is the C ourantFriedich-Levy (CFL) lim it that guarantees stability; the other is to sufficiently resolve the
problem to reduce the num erical dispersion resulting from the discretization. T he first
lim it can be rem oved by the unconditionally stable alternating direction im plicit finitedifference tim e-dom ain (A D I-FD T D ) m ethod [56] [57]; the second lim it can be rem oved
by the application o f high-order schem es such as the M ulti-R esolution T im e-D om ain
(M RTD) m ethod [58] and the Pseudo-Spectral T im e-D om ain (PSTD ) m ethod [59].
1.2.2
TLM
The transm ission line m ethod (TLM ) is another popular tim e-dom ain m ethod proposed
by Johns [60]. TLM is a num erical technique for solving field problem s using circuit
equivalent. It em ploys the analogy betw een the M axw ell’s equations for electric fields
and m agnetic fields and the equations for voltages and currents on a m esh o f continuous
transm ission lines [61]. B y com puting the voltage and current in the netw orks, the
electric and m agnetic fields can be solved. The sym m etrical condensed node form ula
introduced by Johns [62] has becom e the standard for 3D TLM m ethods.
Like other num erical m ethods, the TLM m ethod solves the electrom agnetic problem
by approxim ating the original continuous problem using a discretized system . H ow ever,
the transm ission line m atrix m ethod uses a physical discretization process; the other
num erical m ethods, such as the finite difference m ethod and the finite elem ent m ethod
use the m athem atical discretization process. B ecause o f the sim plicity in form ula and
program m ing, the T L M m ethod obtained m uch attention in m any applications and
becam e one o f the m ost popular tim e-dom ain m ethods [63]-[69].
The general procedure of TLM m ethods includes tw o basic steps:
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
a) R eplace the field region by an equivalent transm ission line netw ork and derive the
equivalence betw een the field and the netw ork quantities.
b) Solve the equivalent transm ission line system through the repeated scattering and
transm ission of voltage waves.
T LM shares both the m ajor advantages and disadvantages of FD TD . O n one hand,
they both used explicit update form ula and easily can handle com plex structures and
nonlinear elem ents and m aterials; on the other hand, they both suffer from the num erical
dispersion problem [70] and a tight tim e-step constraint.
1.3
Motivation and contributions
A lthough the FD TD m ethod is sim ple in concept and program m ing as well as robust and
flexible in applications, it still requires intensive com puter resources, especially for
com plex and electrically large structures. T herefore, efforts in im proving com putational
efficiency and accuracy continue unabated. T here are m any techniques to achieve high
com putational efficiency [55]; for exam ple, sem i-analytical m ethods w ere applied to
reduce the problem dim ensions, incident w ave generator and absorbing boundary
conditions w ere developed with high perform ance and efficiency, and subgridding
m ethod w ere form ulated for com plex structures w ith fine geom etric details.
For w aveguide structure, a sem i-analytical m ethod, called com pact 2D FD TD , was
proposed [71][72]. It can obtain the propagation characteristics of w aveguide structure by
reducing a 3D FD TD problem to a 2D FD TD problem , but it cannot be used to analyze
the discontinuity along the direction o f w ave propagation. T o solve the issue, a new
com pact 2D FD TD m ethod is proposed in the third chapter of this thesis; it can be used to
analyze the discontinuity along the direction o f w ave propagation.
D uring the application o f w aveguide structures, an incident w ave as the source and
absorbing boundary condition as the m atching load are often required in order to obtain
the S-param eters [73]. U sually, the incident w ave is
obtained by sim ulating
uniform w aveguide and it needs a large am ount of com puter m em ory and long run
a long
tim e.
T here are m any kinds o f absorbing boundary conditions in w aveguide applications
[50][74]-[86]; how ever, they can not provide good absorption for frequencies around the
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
10
cut-off frequency. To solve the issue, a com pact ID FD TD form ula is proposed in the
fourth chapter, w hich can provide very good absorption perform ance for the full
frequency spectrum including the cut-off frequency. W hen applied to the C om plex
Frequency Shifted Perfectly-M atched L ayer (C FS-PM L) m ethods [83] [84], an efficient
m odal based PM L is produced in the fifth chapter.
In order to solve com plex problem s w ith fine geom etric details, m any FD TD
subgridding m ethods have been proposed [87]-[98], How ever, these subgridding schem es
have suffered the problem s o f late-tim e instability or uncontrollable reflections from the
interfaces. To overcom e the problem , a new subgridding form ula is proposed in the sixth
chapter, which is robust and has controllable low reflection.
W hen the FD TD m ethod is used to sim ulate R F/m icrow ave circuits w ith active
devices, it is necessary to consider how to com bine the FD TD with active devices. One
approach is to use the tim e-dom ain lum p-elem ent circuit m odel o f the active device
directly [99] [101]; another approach is to use the frequency-dom ain param eters o f the
active device [102][103]. The tim e-dom ain lum p-elem ent circuit m odel o f an active
device is not available often and only the frequency-dom ain param eters are available. In
such a case, the frequency-dom ain param eters should be converted into causal tim edom ain param eters.
How ever, the frequency-dom ain param eters o f m any electronic com ponents are
often given or m easured in a band-lim ited frequency range. D irect application o f the
inverse Fourier transform to the band-lim ited frequency-dom ain param eters usually leads
to non-causal param eters in the tim e-dom ain. Such non-causal tim e-dom ain param eters
are neither physical nor com patible with, or capable of, being incorporated into a
sim ulator based on a causal m athem atical m odel such as the FD TD m ethod. To obtain a
causal tim e-dom ain m odel from the band lim ited frequency-dom ain param eters, three
m ethods have been used in this thesis.
First, two iterative m ethods based on the Fast Fourier T ransform (FFT) are presented
in the seventh chapter: one applies the FFT in com bination with the error feedback
principle; the other applies the H ilbert transform in conjunction w ith the FFT. These
m ethods are conceptually sim ple and easy to im plem ent. The extracted tim e-dom ain
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
11
param eters are not only causal but also contain alm ost the sam e frequency-dom ain
inform ation as the original param eters over the given lim ited frequency range in both
m agnitude and phase.
H ow ever, the tim e-dom ain param eters extracted using the iterative m ethods are only
in num erical form and are com putationally tim e consum ing w hen a convolution is
perform ed w ith them in a tim e-dom ain sim ulation. To solve this problem , a
rational
function fitting technique is introduced in the eighth chapter, w hich extracts the causal
tim e-dom ain param eters o f an active device from their know n band-lim ited frequencydom ain counterparts. One o f the m ajor advantages o f this technique is that the resulting
tim e-dom ain param eters can be expressed in the form o f exponential functions. The
convolution w ith these exponential functions can then be perform ed in a recursive
fashion w ithout requiring a com plete past history o f the tim e-dom ain param eters, and the
CPU tim e for each tim e-m arching step is constant. The total CPU tim e and m em ory o f a
convolution will increase linearly w ith the tim e steps. T herefore, com putational
efficiency is im proved significantly, especially for a sim ulation w ith a large num ber o f
iterations.
Finally, it is w orth to m ention that the m ost w ork o f the thesis has been published in
[104]-[112].
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
12
2 Concept of FDTD
2.1
Introduction of classical FDTD
The FD TD m ethod is a direct approxim ation o f M axw ell's curl equations in the tim edom ain. It discretizes the differential form o f M axw ell’s equations directly by using 2nd
order central finite difference approxim ation. In the FD TD , the electric field E -grid is
offset both spatially and tem porally from the m agnetic field H -grid. T he resulting update
form ulae for the E and H fields are know n as the leap-frog scheme.
2.1.1
Yee’s grid and FDTD formula
In a linear and isotropic m edium , M axw ell's curl equations can be w ritten as:
//— = - V x £ - p 'H
dt
(2 . 1)
e ™ = V xH -aE
dt
( 2 .2 )
w here £ is the electrical perm ittivity, p
is the m agnetic perm eability, a
is the
conductivity, and p ' is the m agnetic conductivity respectively.
In C artesian coordinates, (2.1) and (2.2) can be rew ritten as the follow ing six scalar
equations:
(2.3)
(2.4)
dt
p
dx
dz
(2.5)
( 2 .6 )
dt
£
dy
dz
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
13
(2.7)
dt
£
dt
dz
£
dx
dx
'
dy
( 2 .8 )
1
U sing 2nd order central finite difference to approxim ate the spatial and tem poral
derivatives in (2.3) - (2.8), Y ee’s FD TD form ulae can be obtained as follow s [55]:
12 P i,j+ i,k + ±
H.
*+4
At
1+
2
H. ”
i,j+±,k+i
2 P i , j + ± , k +±
(2.9)
At
n
"
.
, - E
i,j+±,k+1
+ -
At
1+
z
, - E
1,7+1,t + 4
Z
i ,j,k + i
Ay
Az
-
2a
1H.
F
y
n + 2
P '^,
.k+i
At
2Pi+i j , k+±
—
l~
r
1+
J'
2
TJ
’
P \ +u ^ At
i+4 ,j,k + i
i+\j,k+±
(2 . 10)
At
i , j,k + i
+
-
1+
P
At
i +k,j,k+ ±
F
n
X
- E
i+k,j,k+1
Ax
x
i + 2>J ’k
Az
2Pi+\,j,k^
H.
P ' i + ± J +l , k A t
2 P j+ \,j+ \,k
t+\.j+-k,k
1+
P
,+J. i + 1 t
2 ’J
^
At
2'
2 P i+i . i+i , k
( 2 . 11)
At
n
F
Pi+j.i+j.k
,
2
! ^ i+j -J+r k
At
.r
i+\,j+\,k
- E
A y
x
i*\,i,k
E
^ y
i+l,7+f*
4.*.
U ./+-
A x
2P,+a +u
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
14
<7i + i,. j ., k.At
1-
2e.i + i J . k
n+1
f+^. j.k
1+
a t + j,. J . k.At
i+i.j.k
l£
( 2 . 12 )
At
n + 1 /2
+-
<t , . Ar
n+1/2
n+1/2
-H,
i+\,j+±,k
£ i+ \,i.k
Ay
tj
n+1/2
i+4
Az
1+
2e. ,
1-
<xl , J + y,, k, Af
2s.
E,
+ ±.k
l.J+2’kAt
7. . ,
1+
-
■'
i,j+ h k
. , .
t,J+-k,k
(2.13)
At
n + 1/ 2
n + 1 /2
^ > i,j'+4,*+4 - H .
,-+± k
1,7 ”
+
1+
//
n+1/2
-H ,
n+1/2
— (—
(7I , J.+ \,, k.At
Ax
Az
2 f l,J+4*
• • ,*
<7 ., . ,A?
1-
•
' 2
<,J,k+\
n +1
1+
,/t+i
<7u j ■, k + \
2£ i . j , k + - L
(2.14)
Ar
H
+ -
1+
°
2 f.
++iAf
n + 1/2
_
i + \ . J . k +-£
(—
n + 1/ 2
tj
v
Ax
i-fM +l
n+1/2
n + 1 /2
i,j'+4,*+4 -
H .
Ay
i-jX+j
w here U\"J k = U ( i A x , j A y , k A z , n A t ) represents the value o f U ( x , y , z , t ) at the point
( x , y , z ) = ( i Ax , j Ay , k A z ) at tim e t = nAt . Ax, Ay, Az and At are the spatial and
tem poral increm ents respectively.
T he corresponding field com ponents in the Y ee’s grids are show n in Fig. 2.1.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
15
Zi
y
E,
F igure 2.1 Positions o f the field components in Y ee’s grids.
It can be seen from equations (2.9) - (2.14) and Figure 2.1 that the com ponents o f
electric and m agnetic fields are interlaced w ithin the FD TD grid and com puted at
alternate half tim e steps. G iven the initial values and boundary values o f electric and
m agnetic fields, equations (2.9) - (2.14) can be updated explicitly. It is, therefore, an
explicit tim e-dom ain m ethod. T here is no need to solve system s o f linear algebraic
equation, each field com ponent value can be evaluated directly from the neighboring
field com ponents and its ow n value in the last tim e step. T he sim plicity of concept and
ease of program m ing m akes FD TD one o f the m ost popular electrom agnetic num erical
methods.
2.1.2
Numerical stability
FD TD is an explicit update form ula in the tim e-dom ain. U sually, explicit update
form ulae have stability problem s. For the given spatial increm ents Ax, Ay, and Az , in
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
16
order to keep the sim ulation stable, the tim e step size At m ust satisfy the follow ing
stability condition [55]:
1
At <■
1
1
1
+
y[/U£ y A x2
(2.15)
1
r + -
A y2
A z2
This stability condition is called C F L stability condition because the analysis m ethod
was presented by C ourant, Friedrich, and Levy (CFL).
2.1.3
Numerical dispersion
W hen using FD TD to sim ulate electrom agnetic w ave propagation, the num erical phase
velocity of the sim ulated w ave m ode in the FD TD lattice can differ from the actual wave
velocity c ; this phenom enon is called num erical dispersion of FD TD [55]. T he num erical
dispersion causes the num erical phase velocity to differ from the actual w ave velocity,
and it will reduce the accuracy o f com puted results. H ence it is desirable to understand
num erical dispersion’s operation and its effect on accuracy, especially for electrically
large structures.
C onsidering a plane m onochrom atic traveling w ave in a uniform lossless m edium ,
Taflove derived the follow ing num erical dispersion relationship for FD TD [55]:
+
*
Ay
sin( '
2
)
12
1
2
(2.16)
£C/5
. £ Ax
sm(
) = * sin( " )
Ax
2
|_cAf
2 J
1
2
1
1
2
<
coAt
1
1
+
LAz
2
J
w here k x , k y , and k_ are, respectively, the x, y, and z com ponents o f the num erical
w avevector, and CD is the w ave angular frequency.
In contrast to the num erical dispersion relationship of FD TD, the analytical
dispersion relationship o f plane w ave in uniform lossless m edium is m uch sim pler [122]:
c
<2'17)
w here kx , Jk , and k. are, respectively, the x, y, and z com ponents o f the analytical
wavevector.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
17
It can be seen from equations (2.16) and (2.17) that the num erical dispersion is
different from
the analytical dispersion. In contrast to the num erical dispersion
relationship o f FD TD , the analytical dispersion relationship o f plane w ave in a uniform
lossless m edium is m uch sim pler. In order to have an acceptably sm all error from
num erical dispersion, the spatial increm ents Ax, Ay, Az should be less than one tenth of
the sm allest w avelength [55].
2.2
Lumped Elements in FDTD
W hen the FD TD is used to analyze structures including lum ped circuit elem ents, one way
to account for the lum ped circuit elem ent is to add a lum ped electric current density term
in the M axw ell’s equations. T he process can be described as follows.
M axw ell’s curl equations applied to the cells that contain the device are:
£ — - = ' Vx H - J
dt
j
(2.18)
(2.19)
w here J j is the additional current density term resulting from the lum p-elem ent device.
To com pute J <i, for sim plicity, consider an one-port lum p-elem ent device that is
oriented in the z-direction as show n in Figure 2.2. Suppose that the current flow ing
through the device is id . Then the current density J j can be obtained as
J - —^ —
‘
( 2 . 20 )
AxAy
w here Ax and Ay are the space increm ents o f the FD TD mesh that covers the cross
section of the device.
id can be determ ined through the know n device I-V relationship,
(2 .21)
where
is the device voltage oriented in the z-direction. For instance, for a resistor of
re sista n c e R, id = v d / R .
T he device voltage vd can be obtained from its relations to the electric field:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
18
( 2 . 22 )
v,= E M
C om bination o f (2.20)-(2.22) reads
_ f{E M )
d
(2.23)
AxAy
T he above equation can be substituted into (2.18), form ing the M axw ell’s equations that
include the lum p device m odel for FD TD com putations. T he unknow ns are the field
quantities to be solved.
Ax
lump-element
device
a -
Az
Ax
Figure 2.2 The position of a one-port lump-element device in the FDTD grid.
The above concept can be extended to m ulti-port device m odels. Figure 2.3 shows a
tw o-port transistor in the FD TD grid and Figure 2.4 presents the extracted tw o-port
netw ork representation of the device.
lum p-elem ent device
Ax
y
▲
X
Az
Figure 2.3 A two-port lump-element device (transistor) in the FDTD grid.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
19
lump-element
device
F igure 2.4 The equivalent two-port network o f an active device.
v dx(t) and v(/2( r ) , and id] (t) and idl (t) are the voltages
of the netw ork. T hey can be related by the know n device
and currents at the tw o ports
adm ittance param eters (Y-
param eters):
yn (0 ® v di (0 + y u (0 ® v di (0
*rfi(0
.y 2i (0 ® v di (0 +
jd2(t)_
(2.24)
(0 ® ^ d i (0
w here ® represents tim e-dom ain convolution, and yn (f)> y 12(0» T2 i( f) and ^ 2 2 ( 0 are
the tim e-dom ain Y -param eters.
v dx(t) a n d v d2(t), and idx(t) and id2(t) are also related to electric fields and current
densities in a w ay sim ilar to that described before for the case of one-port device. For
instance, voltages v^, and vd2 are related to electric fields as shown below :
v di
=- f E
dl
i = l,2
(2.25)
w here Po is the voltage reference point (for exam ple, the grounded via) and P; is the point
need to be com puted. idx(t) and id2{t) are related to the current densities through (2.20),
or m ore specifically,
J*= -
L
AxAy
i = 1,2
(2.26)
C om bination o f (2.24)-(2.26) reads:
Jdi = f i ( E , y lX, y i2)
i = 1,2
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(2.27)
20
Substitution o f the above equation into (2.18) forms the M axw ell’s equations that include
the lump device model for FDTD computations. Again, the unknowns are the field
quantities to be solved.
Once the M axw ell’s equations including lump-element devices are formed, they can
be solved with details described in [55][99]-[l 03].
2.3
C onclusion
This chapter has reviewed the basic concepts and formulas associated with FDTD and
introduced briefly the stability condition and numerical dispersion o f FDTD. Finally, the
method incorporating the governing voltage-current equations o f a device into the FDTD
frame was discussed.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
21
3 New Compact 2D FDTD Method for 3D
Waveguide Structures
3.1
Introduction
C om putational efficiency still rem ains a challenge with the FD TD m ethod. For a threedim ensional guided-w ave structure, three-dim ensional num erical FD T D
grids are
generally needed. To reduce the m em ory and CPU tim e requirem ents, tw o-dim ensional
com pact FD TD m ethods have been developed [71][72] w here field variations along the
propagation (or longitudinal) direction are assum ed to be com plex exponential functions.
T he m ethods are very useful and effective in analyzing the dispersion characteristics of
guided-w ave structures, such as m icrostrip lines or coplanar w aveguides.
H ow ever, the com pact 2D FD TD m ethods developed so far assum e the structure
uniform ity along the propagation direction and cannot be used for com puting geom etric
or m aterial discontinuities in the propagation direction. Therefore, an efficient m ethod is
desired that can
account for w aveguide
structure discontinuities.
A m ong
m any
w aveguide structures, structure uniform ity does exist in one o f the transverse directions.
For exam ple, m any rectangular w aveguide connectors and filters popular in industry
[1 13]-[115] have such a feature. N avarro et al. studied these types o f structures partially
[116], but only considered the TE10 m ode in the H -plane w aveguides. Field variations in
one of the transverse directions (for exam ple, the y direction in this case) are assum ed to
be non-existent. Such an assum ption is not suitable for m any applications w here higherorder m odes exist. Therefore, a m ore general com pact 2D FD TD m ethod is needed. Since
the structures are uniform in one transverse direction, field distributions in this direction
can be expanded w ith sine or cosine functions in space. B ased on this, we propose a new
com pact 2D FD TD m ethod.
First, w e can classify the traveling m odes based on the num ber o f periods o f field
variations in this uniform direction; then, we can use this fact to reduce the 3D w aveguide
problem s to a 2D problem s and solve them with the new 2D com pact FD TD form ulation.
In contrast to the 2D com pact FD TD m ethods reported before [71] [72], the proposed
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
22
m ethod can be applied to discontinuity problem s in the longitudinal direction.
This chapter is organized in the follow ing way. First, the derivations o f the proposed
2D com pact form ulae are given. Tw o num erical exam ples are then presented to
dem onstrate the validity and effectiveness o f the proposed technique. Finally, discussion
and conclusions are presented.
3.2
Formula of the New Compact 2D FDTD
W ithout losing generality, we assum e that the geom etry o f the w aveguide structures
under study are uniform in the y direction. W e also assum e that the m odes under
consideration have n standing waves in the y direction; w e can call them Tn m odes for
sim plicity.
Since the perfect conductors are placed at the end w alls o f the y-direction, the
norm al m agnetic field com ponent h v and tangent electric field com ponents Ex and e ,
have to be zero at y = 0 and b. As a result, variations o f the field com ponents along ydirection can be expanded in term s o f sin(7ryl b) or cos( n y / b ) . M ore specifically, the six
field com ponents can be expressed as:
Ex (x, y , z , t ) = Ex(x, z, t) s in ( ^ ^ )
b
£ ( x , y , z , t ) = E (x, z,t) cos(”~ ^ “)
b
nTzy
E, (x, y, z, t) = E, (x, z,t) sin(— —)
L
H X( x , y , z , t ) = H x{ x , z , t ) c o s ( ^ y - )
(2U )
H (x, y , z , t ) = H (x, z , t ) sinC-^^")
b
H_{x, y , z , t ) = H . ( x , z , t ) c o s ( r ^ - )
b
w here b is the w idth o f the w aveguide in the y direction, Ex - H_ represent the rem aining
part of each field value function after the term including the variations along the y
direction is extracted.
By inserting (3.1) into M axw ell’s equations in the C artesian coordinates, the
follow ing equations are obtained:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
23
dE
.
at
1,nn
.ri7I y
-
.
dH
nTZy.
+ ) = — (— H , + — ± ) s m (— i )
b
£
b
dz
b
dE
,rwry 1 , d H
3tf_. ,«rry.
— cos(— - ) = - (—-i - — cos(— i )
of
o ^ dz
ox
b
dE, .
,rwry 1 3//
n;r _
nrrv
i r sm(— ) “ <- ^ + T , i - )s,n(- r )
( 3 ' 2 )
at
b
ju
3tf . /wry
— sm(—-i) =
dr
o
b
'
dz
b
1 ,3£ 3£ . . n n y .
(— 2 -- —-) sin(— i )
/i
dz
ox
b
3H .
,nny
1 3£ n x _
,/wry.
—- cos(— - ) = ---- (—- - - - E x) cos(—-i)
3r
6
ii dx
b
b
Removing the common terms on both sides of (3.2), we can reduce (3.2) to 2D
M axw ell’s equations:
d E x( x , z , t )
1. t i n dH (x,z,t)
— ^ ------= — (— H , ( x , z , t ) + — -------)
dt
£
b
3E y ( x , z , t ) _ l
dt
dz
£
dz
dx
3E X x , z , t ) 1 d H y ( x , z , t )
— ------ = - ( ----^
dt
3H x( x , z , t )
----- =
of
£
dx
n7t
,
+ — H x(x,z,t))
b
1/wr-,
N d E y( x , z , t \
(— E , ( x , z , t ) --------)
fi
3H y { x , z , t )
i
dt
fl
b
(3 3 )
oz
^d E x ( x , z , t )
3E r ( x , z , t )
dz
3H , ( x , z , t )
1 dE
------= ---- (—
dt
3H r ( x , z , t ) .
d H x( x , z , t )
dx
(x,z,t)
nn -
^
-----------— E x(x, z , t ) )
jJ.
dx
b
By replacing the derivatives with their 2nd order central finite-difference counterparts,
the discretized equations shown in (3.4) - (3.9) are obtained:
(3.4)
£
£
AZ
Ax
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(3.5)
24
E'r'(i, k + j ) = e : (i, k +$)+— —
£
| At
H
?
h
"+' (i, k +j)
b
(i + ± , k + ± ) -
£
H"+i (i -
4 , A:
(3.6)
+ 4)
Ax
/7 ;+>a, k + v =
(i , k + i )
^
fi b
e:
a ,k + \)
'
(3.7)
At E"y {i,k + \)-E"y (i, k + 1)
A
Az
H ;+" a + i ^ + i ) = /7 ;-= a + i,z )
At E"(i + ± , k + l ) -E"z (i + ± , k )
H
| At £;'(» +
1,
fc +
//
H ' r H i + ^ , k )
(3.8)
Az
4 ) -£ ;'(,- ,J k + X )
Ax
= H ' r H i + \ , k ) + —
-
n
—
E : { i + \ , k )
b
(3.9)
At E ; a + U ) - E ; ( i , k )
LI
Ax
T he grid positions and field com ponents for the 2D FD TD m ethod is show n in
Figure 3.1.
H
E.
E_
E.
E
H.
H
H
XA
y /
E W
(WO
H
E
WE
Figure 3.1 The electric and magnetic fields positions in a unit cell of the compact 2D FDTD grid.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
25
Equations (3.4) - (3.9) are tw o-dim ensional and easy to program . In the follow ing
sections, two num erical exam ples are given to show the efficiency o f this new m ethod.
3.3
Numerical Validation
To validate the proposed m ethod, w e first apply it to a rectangular w aveguide cavity with
length d=0.04m in the z direction, w idth a=0.03m in the x direction, and height b=0.02m
in the y direction (w hich is assum ed to be the geom etry uniform direction). T he analytical
resonant frequencies w ere docum ented in [117]. T he cavity structure is show n in Figure
3.2 and the num erical grid plane is in the x-z plane.
y
b
Figure 3.2 The rectangular waveguide cavity of length d=0.04m in the z direction, width a=0.03m in the x
direction, and height b=0.02m in the y direction.
In order to facilitate the com parisons, a reference 3D FD TD sim ulation is also run
for the same cavity structure. The m esh sizes w ith both the proposed m ethod and the
reference m ethod are taken to be the same: A x= 0.001m , A y=0.001m , and A z= 0.001m . The
tim e step size is At=AiCFL = 1.9245x10“12j , and the AtCFL is the C F L tim e step lim it com puted
for the conventional 3D FD TD. The total num ber of iterations taken is 8192. T he source
for 2D FD TD is a point source w ith E z and H z. The source for 3D FD TD is a line source
w ith Ez that varies as
s m ( x y / b )
along y direction (for T M Z m ode), and w ith H z that varies
as c o s (^ y lb ) along the y direction (for T E Z m ode). The source in the tim e-dom ain is two
pulses, which are equal to 1 at tim e step 1 and equal to -1 at tim e step 2. T here are two
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
26
reasons to select this kind o f source: one is that it is easy to be im plem ented, the other is
that it includes all the possible frequencies needed. H y and E y are recorded and Fouriertransform ed to obtain the resonant frequencies by identifying the am plitude peak points.
Table 3-1 show s the analytical resonant frequency of the first five m odes o f Tj and
the relative errors o f num erical results com puted w ith the proposed 2D FD TD and the
reference 3D FD TD m ethods.
Table 3-1 The comparison between numerically computed and analytical resonant frequencies
o
H
m
The
resonant
modes
TM-no
TE 111 /TM 111
TE 012
TE 112 /TM 112
Analytical
results
8.3853GHz
9.0139GHz
9.7628GHz
10.607GHz
11.726GHz
errors
with the
proposed
method
-0.15%
-0.08%
0.06%
-0.13%
0.07%
errors
with the
reference
FDTD
-0.15%
-0.08%
0.06%
-0.13%
0.07%
Table 3-2 shows the respective m em ory and CPU tim e used w ith the 2D FD TD and
the 3D FD TD for the first exam ple. The running platform is the M atlab on a laptop
Pentium -IV PC w ith 1.8-GHz CPU and 512-M B RAM .
Table 3-2 The memory and CPU time needed for the 2D FDTD and 3D FDTD corresponding to the first
example
Memory
CPU time
The reference 3D
FDTD
1507KB
399s
The Proposed
Method
316KB
9s
It can be seen from Table 3-1 that the proposed m ethod yields alm ost the same
results as those w ith the reference 3D FD TD m ethod. H ow ever, the proposed 2D FD TD
m ethod uses about 1/44 o f the CPU tim e (including FFT) o f the reference 3D FD TD and
about 1/5 o f the m em ory (including FFT) o f the reference 3D FD TD.
The second exam ple is to validate the discontinuity in a w aveguide structure. It is a
rectangular w aveguide cavity o f the sam e size as the previous exam ple but w ith a groove
in the m iddle o f the z direction as shown in Figure 3.3. The w idth o f the groove is
0.006m and the depth is 0.01m .
T he m esh size, the tim e step size, and the source are the sam e as in the first exam ple.
The total num ber o f iterations is 16384. H y and E y are again recorded and Fourier
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
27
transform ed to obtain the resonant frequencies by identifying the am plitude peak points.
Table 3-3 show s the values o f resonant frequency of the first five m odes o f Ti com puted
by the proposed 2D FD TD and the reference 3D FD TD m ethods. T able 3-4 show s the
com putation resources used w ith the tw o m ethods.
O nce again, it can be seen from Table 3-3 that the proposed m ethod gives alm ost the
sam e results as the reference 3D FDTD. H ow ever, the proposed 2D FD TD m ethod uses
about 1/42 of the CPU tim e (including FFT) o f the reference 3D FD TD and about 1/3 of
the m em ory (including FFT) o f the reference 3D FDTD.
Table 3-3 The comparison o f resonant frequencies obtained by 2D FDTD and the 3D FDTD in the second
example
Mode
Mode
Mode
Mode
Mode
1
2
3
4
5
The reference
3D FDTD
8.0873GHz
9.0387GHz
9.4827GHz
10.720GHz
11,988GHz
The Proposed
Method
8.0873GHz
9.0387GHz
9.4827GHz
10.720GHz
11.988GHz
Table 3-4 The memory and CPU time needed for the 2D FDTD and 3D FDTD in the second example
Memory
CPU time
The reference 3D
FDTD
1763KB
786.5s
The Proposed
Method
572KB
18.8s
cl
Figure 3.3 The rectangular waveguide cavity with a groove in the middle of the z direction.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
28
3.4
Discussion and Conclusion
In this chapter, a new com pact 2D FD TD m ethod for 3D w aveguide problem s has been
proposed. N um erical exam ples presented in this chapter show that this m ethod has the
sam e accuracy as the traditional 3D FD TD m ethod, but w ith m uch less m em ory
requirem ents and CPU tim e. N evertheless, it should be noted that the num erical analysis
is not a theoretical one; m ore com prehensive analytical studies are needed on the exact
stability condition and num erical dispersion relationship o f the proposed com pact 2D
FD TD. In addition, further applications to com plicated w aveguide structures are needed.
These will be left for future work.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
29
4 Compact ID FDTD Method for Waveguide
Structures
4.1
Introduction
W hen using the FD TD m ethod to com pute w aveguide structures, tw o things are often
needed: a know n incident w ave for calculating electrical param eters (e.g., scattering
param eters) and effective absorbing boundary conditions for term inating open structures.
To obtain an incident w ave, a separate sim ulation o f a long w aveguide structure is
usually run [73]. H ow ever, for a three-dim ensional (3D) structure, such a sim ulation is
often inefficient because it requires a large am ount of m em ory and CPU tim e. In order to
solve the problem s as w ell as to develop efficient absorbing boundary conditions, m any
one-dim ensional (ID ) m ethods have been proposed [74]-[81]. M ost of them use
analytically or sem i-analytically generated G reen’s functions. H ow ever, these analytical
continuous G reen’s functions often have characteristics quite different from those o f the
3D FD TD form ulations o f a discrete dom ain, for instance, at the cut-off frequencies.
Consequently, they do not offer highly accurate results, for instance, near or below the
cut-off frequency o f a w aveguide m ode.
To solve the problem , w e propose a new sim ple ID FD TD m ethod in this chapter.
U nlike other m ethods developed so far, this m ethod is derived directly from the FD TD
form ulae; therefore, it has the sam e num erical characteristics as that o f the 3D FD TD
m ethod with w hich it interfaces. As a result, it not only allow s efficient com putation o f
an incident w ave due to its ID nature, but it also enables an extrem ely high absorption o f
num erical incident w aves (e.g. better than -200dB even at or below a cu t-o ff frequency).
This chapter is organized as follows. Section 4.2 presents the derivation o f the
proposed ID m ethod, w hile Section 4.3 illustrates the dispersion analysis o f the proposed
m ethod. Section, 4.4 describes the two applications of the proposed m ethod: generation
of the incident w ave and absorption o f w aves. Section 4.5 sets forth the num erical results
that dem onstrate the validity and effectiveness of the proposed technique. Further
discussion and conclusions relating to the proposed m ethod follow in Section 4.6.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
30
4.2
Formula of the Proposed ID FDTD Method
For linear and lossless m edium , the conventional 3D FD TD form ulae can be w ritten as:
E I"*1. = £ 1"^ ..
+^eAy <H=S U
(4.1)
At
A'
eAz
-
' I t *
H
f
,
c
W
(4.2)
T he equations for the other field com ponents can be expressed sim ilarly.
For a hom ogeneously filled w aveguide, field distribution pattern of a given m ode on
a cross-section do not change w ith tim e or frequency or w ith the longitudinal coordinate.
They can be found analytically or num erically (e.g. [71][72]). Supposing that the m ode is
traveling in the z-direction and Y ee’s grid is applied as shown in Figure 4.1, the
discretized form o f equations (4.1) - (4.2) can then be rew ritten as:
M U ,,-M U ,,
+ eAy
H>
(4.3)
k~2' = H >■ '4.2T-H
At
+ / jA x
~//A( 7E , I'Url-k
j t ~E>
* 1
(4.4)
1
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
31
12.k
1 - 1 /2 ,k
F igure 4.1 The electric and magnetic field positions in Y ee’s lattice for a waveguide structure, with Z the
wave propagation direction.
T he equations for other field com ponents can be w ritten sim ilarly as follow s:
p i»+i _ p I"
At f Tt ,11+4
r i l"+4 \
PA7
•*+2
2 2
£Ax
t r'
-1)//. I"**
' ' '-1
^
=
(4.5)
e _ r'
At
'“i-J
~
—
ehy
{ a „
I.
■ ,
*> 'v4
(4.6)
>
- 1 ) H
X
r * , .
,
1 'o-iH
Ar
jlA z
-— ( 1 - 5 . 1,,)£. I" ,
//Ay
^ M ‘ '■■'•‘T
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(4.7)
w here
H T l',+7 . ,
(4.9)
" 2’1 1
I"",1* z x ii ~ i / —L - - T„ j _ lg i l
r' 1 H. 1*7.
(4.10)
H
or...
” '
1/ |,,+2
,
CV*
^ U = 7 7 7 ^
‘ ' -H*-k
(4.11)
(4-12)
£ . I" , i- i
P .k r - fr ^
(4.13)
(4.14)
- iM-J
* H-/.*
y i.j-kyk
(4-15)
(4*16)
Coefficients a and p are ratios o f field quantities on the nodes o f a cross section o f
the w aveguide. They are constant and can be found from the know n unchanged field
distributions o f a given m ode. Each m ode has its ow n set o f coefficients a and p . N ote
that in com puting or and p , one should choose non-zero field points for the
denom inators in equations (4.9) - (4.16) for a given mode.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
33
Careful exam ination of the equations reveals that equations (4.3) - (4.8) are
essentially ID FD T D recursive form ulae. C om putations need to be carried out only along
the z-direction (or the ^-direction) w ith the specified i and j for each m ode. A t any other
i and j , the field quantities can be obtained from the know n field distributions on the
sam e cross section.
4.3
Numerical Dispersion of the Proposed Method
In order to com pare the num erical dispersion o f the proposed m ethod w ith the
conventional FD TD m ethod, we consider the TE™ m ode in a rectangular w aveguide as
an example. F or the T M mn m ode, a sim ilar analysis can be m ade and sim ilar conclusions
reached.
Suppose the rectangular w aveguide has w idth a in x direction and height b in y
direction. The field com ponents for the T E mn m ode along z-direction can be w ritten as:
Ex = E x0 cos(kxx) sin(kyy) eJ{ktZ ffl,)
Ey = Ey0 sin(/;jx)cos(fcv>')e'(M_“ )
E, = 0
H x = H x0 sir\(kxx)cos{kyy ) e 1'-kzZ~‘a)
H y = H y0 cos(kxx)sin(kyy ) e nt'-z' ,a)
/ / , = H , 0 cos(kxx)cos(kyy ) e liklZ~'a)
w here kx = —
a
, k, is the spatial frequency in the z direction, and co is the
, k
’
b
tem poral angular frequency.
Substituting (4.17) into (4.3) - (4.8), we obtain:
(4.18)
(4.19)
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
The above equations form a system of five hom ogeneous equations w ith unknow ns
Exo, Ey0, Hx0, H yo, and H z0. B ecause the solutions of the system m ust be nontrivial, the
determ inant of its coefficient m atrix should be equal to zero. This leads to:
(4.23)
(4.24)
+
Ax2
w here k
a
+
Ay2
(4.25)
Az2
b
E quation (4.23) corresponds to® = 0 , and represents the static solution. Equation
(4.24) will lead to Hz0 = 0 , w hich does not agree w ith the assum ption o f T E m odes.
Therefore, the rem aining equation (4.25) is the num erical dispersion relationship o f the
TE m odes. M ore details on the num erical dispersion relationship can be found in
A ppendix A.
It is obvious that equation (4.25) is the sam e as the num erical dispersion relationship
of the 3D FD T D m ethod [55] for T E mn m ode w ith kx =—
a
and k = — .
b
To verify the above claim num erically, we considered a rectangular w aveguide with
width a= 0.02m in x direction, height £>=0.01m in y direction and z was the w ave traveling
direction. T he m esh size was Ax=0.001m, Ay=0.001m, andAz=0.001m for the 3D FD TD
m esh that discretizes the w aveguide. The tim e step size was taken as Ar=Armax w here
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
35
Armaxis the C F L tim e step lim it o f the 3D FD TD form ulae. TEio m ode w as excited w ith a
m odulated G aussian pulse sin(2;r/0f)e~l'~'o)’/:r' in the center o f the structure. P aram eter T
equaled 20A t ,
t0 equaled 60Ar , / 0 equaled 60.465GFIz
and the corresponding
wavelength A in the w aveguide was about 5Az = 0 .0 0 5 m . The recording point was 300Az
or 60As away from the source plane. Such a long distance betw een the source and the
field recording point allows us to observe effects o f the num erical dispersion on the field
solutions. Figure 4.2 shows the Ey s com puted w ith the 3D FDTD and the proposed ID
FD TD . The difference o f the two Ey s com puted w ith the 3D FD TD and the proposed ID
FD TD is show n in Figure 4.3.
As can be seen, tw o e
s
overlap com pletely. T he m axim um difference betw een the
tw o Evs is less than 2xl(T15(V/m), w hereas the m axim um field value is around 0.5 (V/m ).
Such sm all differences suggest that they are due to com puter rounding-off errors. N ote
that the selection o f the specified i and j in equations (4.3) - (4.8) has very little effect on
the error level. These verify experim entally the claim we had before: the num erical
dispersion relationship of the proposed ID FD TD is the sam e as that o f the original 3D
FDTD.
0.5
proposed 1D FDTD
3D FDTD
0.4
0.3
0.2
E
>
>.
L
U
- 0.1
- 0.2
-0.3
-0.4
-0.5
t(ns)
Figure 4.2 The Ey recorded at a point 300Az or 6 0 2 , away from the source plane.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
36
LU
o
a>
o
c
8?
® -0.5
%
-1.5
t(ns)
Figure 4.3 The difference between the two Eys produced by the ID FDTD and the reference 3D FDTD.
4.4
Applications of the Proposed Method
The proposed m ethod as represented by equations (4.3) - (4.8) have tw o m ajor
applications in sim ulating a w aveguide structure: to efficiently generate num erical
incident waves that are required for com puting electrical param eters such as scattering
param eters and to serve as a w ideband absorbing boundary that is com puted only in one
dimension.
Efficient generation o f an incident wave
B ecause of its ID nature, the proposed m ethod can be used to obtain an incident w ave by
com puting a long w aveguide structure; the w aveguide is long enough that the wave
reflected by any im perfect term ination at the ends cannot return to the m easurem ent point
and contam inate the incident w ave [73]. In the num erical exam ple presented in Section
4.5, the incident w ave obtained w ith the proposed m ethod is found to be fundam entally
the sam e as that obtained using the conventional 3D FD TD m ethod; the differences was
less than -200dB. In other w ords, the proposed m ethod can produce an incident w ave for
all intents and purposes identical to that obtained with the conventional 3D FD TD . This
stem s from the fact that the proposed ID m ethod is derived from the 3D FD TD m ethod
and therefore has the same num erical characteristics as that o f the 3D FD TD m ethod, as
proven before.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
37
Absorbing boundary condition
Since the proposed m ethod can easily sim ulate a long w aveguide structure, it can also be
used to term inate a w aveguide structure as illustrated in Figure 4.4. In Figure 4.4, a
waveguide is connected to a discontinuity w here both of them are m odeled using the
conventional 3D FD TD grid. The w aveguide is then term inated w ith the absorbing
boundary that is m odeled using the proposed ID FD TD m ethod. Field com ponents
Ex I , a n d
Ey I.
are used to pass the field values from the 3D FD T D grid into the
proposed ID FD TD grid, w hile field com ponents Ey |(
and Ex \._^Jk are used to pass
the field values in the proposed ID FD TD grid into the 3D FD TD grid.
Waveguide discontinuity
3D FDTD inesh
i LX
3D waveguide
Proposed ID FDTD
mesh
kJu;: j,k-i
Interface layer between the 3D FDTD and ID FDTD
F igure 4.4 The proposed absorbing termination using the ID FDTD mesh [76][77],
M ulti-m odes generally exist in the w aveguide. H ow ever, equations (4.3) - (4.8) are
valid only for a single m ode. To solve the problem , each m ode has to be extracted at the
interface betw een the 3D FD TD m esh and the proposed ID FD TD m esh. T he m ode
extraction can be perform ed using the orthogonality o f m odes as described in [76] [77],
Figure 4.5 illustrates such extraction operations. For the com putation o f the total field
value in the 3D FD TD from the know n m ode field values obtained by the ID FD TD , it
ju st needs to add the field values o f different m ode together. The 3D field values of each
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
38
m ode can be obtained by the ID FD TD field values and the field distribution pattern of
this m ode.
In Figure 4.5, each m ode corresponds to an independent ID FD TD m esh line. The
positioning o f the ID FD TD m esh lines, i.e. the specific values o f i and j in (4.3) - (4.8),
can be the sam e or can be different. T he requirem ent for the positioning o f each ID
FD TD m esh line is that it should not be at the null field points o f the m ode it sim ulates.
3D FDTD mesh in 3D waveguide
Proposed ID FDTD mesh
for each mode
i
Mode 1
'E.
h
M odel
£ X' ■
E
•a
ModeN-1
■a; 'Ex
E
a
■
Mode \
'E_
Interface layer between the 3D FDTD and ID FDTD
>Y
Figure 4.5 The mode extraction and combination diagram at the interface between the waveguide and the
ID FDTD absorbing termination [76][77].
In the num erical exam ple presented in Section 4.5, it is shown that the proposed
term ination provides an absorption of better than -200dB even at or below the cut-off
frequencies in both the single and m ulti m ode case.
4.5
N u m e ric a l E x a m p le s
W e considered again the sam e w aveguide as that used in Section 4.3: a rectangular
w aveguide with w idth a= 0.02m in x direction, height b=0.01m in the y direction and z as
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
39
the wave propagating direction. T he m esh size isA x = 0 .0 0 1 m , A y= 0.001m , and A z= 0.001m .
The tim e step size is taken as At=Attmx. T he total num ber of the iterations is 4096 (which
am ounts to 7.8808ns o f the real tim e). T he source used in the FD TD sim ulation is the
D irac im pulse function S(n) . M atlab was used to program the m ethods and the sim ulations
w ere run on a laptop Pentium -IV PC w ith 1.8-GHz CPU and 512-M B R A M .
For the first application, w e com puted the incident waves for T E n m ode. The
w aveguide is long enough to ensure there is no reflection from the far end to the
recording points. T he tim e-dom ain signatures o f
e
, obtained with a reference full-w ave
3D FD TD sim ulation and the proposed ID FD TD sim ulation, were recorded at a point
lAz aw ay from the source plane. T he result are show n in Figure 4.6 (for clarity, only first
0.5ns is shown). It can be seen that they overlap com pletely. T he corresponding
frequency-dom ain relative errors o f the incident waves, obtained by the proposed ID
FD TD sim ulation com pared w ith the results obtained by the reference full-w ave 3D
FD TD sim ulation, is shown in Figure 4.7. The m axim um relative error in the frequencydom ain is less than -200dB , w hich is due to com puter rounding-off errors.
0.2
proposed 1D FDTD
3D FDTD
0.15
2T 0.05
-0.05
- 0.1
0.05
0.15
0.2
0.25
t(ns)
0.3
0.35
0.4
0.45
0.5
Figure 4.6 The Ey values o f the first 0.5 ns recorded at a point lAz away from the source.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
40
-240
-250
-260
o
t
-270
s>
S
-290
-300
-310
100
150
200
250
f(GHz)
Figure 4.7 The relative errors o f Ey at a point lAz away from the source obtained by the proposed ID
FDTD.
Table 4-1 show s the com putational expenditures used for a w aveguide o f length
206.8cm , w hich was the w aveguide length used in all the num erical exam ples. It is long
enough to ensure the reflection from the far end w ould not reach the interface o f ID and
3D regions during the sim ulation tim e. As can be seen, the m em ory used by the proposed
ID m ethod is about 0.6% o f that used by the 3D FD TD , w hile CPU tim e is about 0.23% .
The proposed m ethod, therefore, saves significant am ounts o f m em ory and C PU tim e
usage.
T able 4-1 The memory and CPU time used by the proposed ID method and the referenced 3D FDTD
method
Memory
CPU time
The reference 3D
FDTD
2334KB
4271s
The Proposed 1D
FDTD
14KB
9.9s
For the second application, we used the proposed m ethod as the absorbing
boundaries to term inate the rectangular w aveguides at both ends. W e then m easured the
effectiveness of the absorption. The source was placed 2Az away from the absorbing
boundary and the e
was recorded at a point lAz aw ay for the com putation o f the
reflection coefficient. Such placem ents o f the source and recording points allow us to
m easure the absorption o f evanescent m odes by the absorbing boundary. T he reflection
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
41
coefficient r caused by the ID FD TD term ination was calculated using the follow ing
form ula:
Ir 1= 20 log I
1 (dB)
( 4.26 )
V
w here E rf is the incident w ave and obtained with a separate reference sim ulation w here a
long w aveguide had been com puted w ith the conventional full-w ave 3D FD TD m ethod.
Figure 4.8 illustrates the calculated reflection coefficient w hen only the single
dom inant m ode
w as excited in the rectangular w aveguide, w hile Figure 4.9 shows
TEn
the calculated reflection coefficient w hen only the m ode TEg4 was excited. Figure 4.10
shows the calculated reflection coefficient w hen m ulti-m odes o f
T E n , T E 21, T E 31,
and
TE41
T E io , T E 20, T E 30, T E 40,
w ere excited sim ultaneously w ith equal m agnitude (the w orst
case in w hich m ultiple m odes exist).
-230
-240
-250
-260
re -270
-290
-300
-310
100
150
200
250
f(GHz)
Figure 4.8 The reflection coefficient from ID FDTD termination for TEn mode.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
42
•240
-250
■270
CD
|
-280
(0
(3
-290
-300
-310
-320
100
150
200
250
f{GHz)
F igure 4.9 The reflection coefficient from ID FDTD termination for TE84 mode.
-250
-260
-270
-280
S'
I
-290
a
O
-300
-310
-320
-330
0
50
100
150
200
250
f(GHz)
Figure 4.10 The reflection coefficient from ID FDTD termination for multi-mode case.
It can be seen from Figure 4.8 to Figure 4.10 that the proposed ID FD TD m ethod
provides alm ost perfect absorbing term inating conditions in the entire frequency
spectrum from D C to 250G H z. In all the cases, the absorptions are better than -200dB
even at or below the cu t-off frequencies (the cu t-off frequency o f T E 10 is 7.5G H z and the
cut-off frequency of TEg 4 is 84.85GHz).
It should be noted that in the above num erical experim ents the num bers o f m odes
excited were chosen arbitrarily to test the perform ance of the proposed m ethod. In
solving a realistic structure, the num ber o f m odes to be considered can be decided in the
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
43
same m anner as that em ployed in the m ode m atching techniques o r in the techniques
described in [74] [75] [77] [81 ][73] [ 118] [ 119].
4.6
Discussion and Conclusion
In this chapter, a new com pact ID FD TD m ethod for uniform ly filled w aveguide
structures has been proposed. It has the sam e num erical characteristics as the
conventional 3D FD TD m ethod. Therefore, it can be used to either efficiently generate an
incident w ave or to effectively serve as a perfect absorbing term ination for specified
m ode. The errors or the reflections w ith the proposed m ethod w ere found to be extrem ely
sm all, less than -200dB even at or below the cut-off frequencies. In addition, despite such
a high degree o f effectiveness, the program m ing o f the proposed m ethod is very easy and
little analytical pre-processing is required.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
44
5 FDTD-Based Modal PML
5.1
Introduction
For an unbounded structure, an absorbing boundary condition is necessary to reduce the
size of the sim ulation dom ain and thus, the C PU tim e. PM L is a very efficient and
flexible absorbing boundary condition, and it can be used for full open problem s or
w aveguide problem s. T he PM L schem e was initially proposed by B erenger in 1994 as an
electrom agnetic absorbing boundary condition [50], Since then, other variations o f the
PM L have been developed and proven to be very effective [82]-[84], In particular, the
com plex-frequency-shifted (CFS) PM L was show n to be effective for arbitrary m edia
[84], In term inating w aveguide structures, m odal PM Ls w ere proposed, w hich reduce 3D
PM L operations or 2D PM L operations to ID PM L operations [85][86]. They save
significant am ounts of com putational m em ory and C PU tim e. H ow ever, these m odal
PM L form ulae are not derived directly from conventional FD TD grids and therefore have
num erical dispersion characteristics different from those o f the FDTD. A s a result, when
connected w ith a FD TD grid, they do not perform as w ell as the original 3D PM L. B ased
on the ID FD T D form ulae in chapter 4, a new ID m odal PM L schem e is proposed in this
chapter that has alm ost the sam e absorption perform ance as the original PM L for single
m ode situations.
5.2
Derivation of the Proposed ID Modal PML
For a hom ogeneously filled w aveguide, field distributions (m odal solutions) o f a given
m ode on the cross-section do not change and can be found analytically or num erically.
U sing this inform ation, a ID FD TD schem e is obtained for m odeling the w aveguide in
chapter 4. This idea can be sim ilarly applied to the P M L form ulae.
Now, consider a w aveguide w ith z-direction as the w ave traveling direction and the
CFS-PM L schem e described in [84] w ith s =1, S =1, and s z
+ — ^ ------- H ers 5 , S
1 a z + jcoeo
k
and s are the stretched coordinate and a z is a dam ping factor that can be optim ized for
better absorption perform ance for a specific m ode. For sim plicity, E x is considered here.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
45
The updating form ula o f Es in C FS-PM L is described in equation (5.1) and (5.2):
__ p
«+•
i~ .j,k
I I , —it Z h+1/2
11+ 1 / 2
+— (H,
*
x
' r 1*
eAy
\
(5.1)
At
m+1/2
- — ( f f v y1 L / A+J.
eAz
|n + l/2 \
t j
n v / - I i *_L '
■(' 2'-' 2
a,k^ + <jr
£0kr
h+i
TJ
= At______ 2___
£„k. a .k . + cr.
At
t-i.M
2
(5.2)
£o+^
Ar
£n&,
2
+ a.
or
At
2
£&,
a ,k , + cr,
°-^+ _ l.-S
2Ar
2
'
2
At
w here Px is an auxiliary or interm ediate variable.
C om paring (5.1) w ith (4.1), it can be readily seen that the sam e ID form ula in (4.3)
can be applied to (5.1) in the sam e m anner. This results the ID m odal C F S-P M L equation
for E r ;
i - j -t
= P. ", t +
At
eAy
( a U . , ■, - 1 )H. 1 7 . ,
(5.3)
L/ iM'
eAz
e 0k,
a jc ^ r o ,
~AT~
/—
i.J.*
£ 0C
2
o r ,/:,
+ cr,
H.j.*
Ar
—
Ar
£nfc_
- 8- ^ +
Ar
(5.4)
— i
2
"
2
Ar
p
GT.& + (7 ,
— =-
*
e0k,
4- +
2_
a ,k , + cr.
HJJt
-
Ar
w here the definition a z is the sam e as (4.9).
The PM L in a w aveguide acts like a w aveguide filled w ith a lossy m edium . T he
m edium is uniform on the cross section o f the w aveguide and therefore it does not change
the cross-sectional field patterns o f a m ode. In other w ords, the cross-sectional field
distribution pattern in the PM L region is the sam e as that in the w aveguide for a given
m ode.
The equations for other field com ponents can be found in a sim ilar m anner and are
listed in appendix B. They form the new ID m odal PM L. Since it is derived directly from
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
46
the FD TD grid, it should have num erical characteristics sim ilar to the original FD T D grid.
It should also perform at least the sam e level as the original 3D PM L for a single m ode.
N um erical results presented in the next section serve to validate these claim s.
N ote that like other m odal PM Ls, the proposed m odal PM L can be applied to one
m ode only. In a m ultim ode situation, the interesting m odes can be extracted using the
m ethod described in chapter 4 and then term inated w ith the m odal PM L associated w ith it.
5.3
N u m e ric a l E x a m p le s
T o validate the proposed ID m odal PM L, a rectangular w aveguide was considered. The
w aveguide had a rectangular cross section w ith a width o f a=0.04m in the x-direction, a
height of £>=0.02m in the y-direction, and the w ave propagation direction was in the zdirection. T he m esh sizes for FD TD sim ulation w ere Ax=0.001m , Ay =0.001 m , and
A z=0.001m . T he tim e step w as taken as Af=Armaj = 1.9245p s w here Armax w as the C F L tim e
step lim it o f the FD TD form ulae. T he total num ber o f the iterations was 8192. O ne end of
the w aveguide was term inated with a 10-layer PM L and another end was connected to a
very long w aveguide. T he conductivity profile for the PM L w as a z(z) = crzm( l - z / d ) 4 ,
r z(z) = 1, and a . ( z ) = a „ „ ( \ - z / d ) i for 0 < z < d ( d was the total length o f the PM L). The
source plane was placed 2 A z s away from the PM L. E v was recorded at a point 1 Az
from the PM L. T he source in the tim e-dom ain was the D irac im pulse function S(n).
a
= 2 tz e J c
w ith
fc
being the cut-off frequency o f the considered m ode.
Tw o m odes, TEio and T E n , w ere excited in the w aveguide, respectively. Tw o
separate sim ulations w ere run: one w ith the proposed ID m odal C F S -P M L and the other
w ith the original 3D CFS-PM L. The relative differences betw een the
e
s o f the two
sim ulations are show n in Figure 5.1 and Figure 5.2. Figure 5.1 is for the TEio m ode,
w hile Figure 5.2 is for the T E n mode.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
47
-320 •
.330*______ 1________1_______ 1_______ i_______
0
50
100
150
200
250
f(GHz)
Figure 5.1 The difference between the E y obtained by the proposed ID modal CFS-PML and that
obtained by the original 3D CFS-PML for TE10 mode.
-220
•230
-240
— -250
-260
•270
-280
-290
-300
-310
-320
100
150
200
250
f(GHz)
Figure 5.2 The difference between the E s obtained by the proposed ID modal CFS-PML and that
obtained by the original 3D CFS-PML for TEn.
As can be seen from Figure 5.1 and Figure 5.2, the differences betw een the proposed
ID PM L and the 3D PM L are extrem ely sm all, less than -220dB across the w hole
frequency spectrum , including those frequencies at and below cutoff. W e therefore
conclude that the proposed ID m odal PM L has num erical properties very close to those
o f the original 3D PM L o f FD TD grid.
To assess the absorption perform ance of the proposed m odal PM L in a general
m ultim ode situation, inductive fins w ere inserted into the w aveguide (see Figure 5.3).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
48
The TE 10 m ode was excited with a m odulated G aussian pulse sin(2^/0f)e'<'“',,)J/r2 in the
center o f the structure. Param eter T was equal to l50Ar and tQ to 6 0 0 A r . T he 10-layer
C FS-PM L was 5 Az s aw ay from the closest fins. E y field was recorded at a point 1 Az
before the P M L and 2 Ax aw ay from the side wall. In this case, the presence of the fins
caused the Ey field to consist o f m any higher order m odes. For conventional C FS-PM L,
f c is the cut-off frequency o f the
T E io
m ode; for
ID
m odal
C F S -P M L ,
f c is the cut-off
frequency o f each m ode considered. In this exam ple, the m odes considered in the
m odal C F S-PM L include
T E i0
to
ID
T E 90.
S ource line for E,,
o
Fins
CFS-PML
Long
waveguide
o
10m m
F igure 5.3 The waveguide structure with two strips under study.
Again, tw o separate sim ulations w ere run: one w ith the proposed ID m odal P M L and
another with the original 3D CFS-PM L. Figure 5.4 shows the reflected E y s obtained w ith
the proposed
ID
m odal PM L and with the original 3D CFS-PM L. Figure 5.5 show s the
corresponding reflection coefficients in the frequency-dom ain with the total num ber o f
the iterations in the tim e-dom ain being 4096.
As can be seen from Figure 5.4 and Figure 5.5, the
ID
C FS-PM L m ethod actually
perform s better than the original 3D PM L w ith sm aller num erical reflections.
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
49
x 10'3
Original CFS-PML
1D Modal CFS-PML
0.4
I
-0.2
-0.4
-0.6
0.5
2.5
3.5
t(ns)
Figure 5.4 The reflected E obtained by the proposed ID CFS-PML for the first 9 modes and the original
CFS-PML.
-10
" ■ Original CFS-PML
1D Modal CFS-PML
-20
e
-30
-40
-50
-60
f(GHz)
F igure 5.5 The computed reflection coefficients by the PML in the corresponding frequency-domain.
5.4
Discussion and Conclusion
In this chapter, a new efficient FD TD -based ID m odal PM L m ethod was proposed. It
was derived directly from the original FD TD form ulae and therefore has num erical
characteristics very close to those o f the original 3D PM L o f FD TD. N um erical results
show that:
a) T he relative differences betw een the results o f the proposed m ethod and the
original PM Ls are extrem ely sm all (less than -220dB) for specified m odes.
b) T he proposed m ethod perform s at least as well as the original PM Ls, if not better.
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
50
6 A New Subgridding Method for 2D FDTD
6.1
Introduction
For com plex problem s w ith a large solution dom ain and fine geom etry, FD TD
sim ulations still require a large am ount of m em ory and a long C PU tim e. O ne o f the
principal reasons is that fine geom etry requires small m esh sizes to resolve the fields
around them. In order to solve the problem , m any subgridding schem es have been
developed for the FD TD m ethod [87]-[98], Fine m eshes are em ployed only in the regions
that contain fine geom etry w hile coarse m eshes are used as m uch as possible. Schem es
are then developed to interface the fine m esh regions and the coarse m esh regions.
H ow ever, these subgridding schem es suffer from tw o problem s: A) long-term or late­
tim e instability and B) uncontrollable reflections from the interfaces. In certain cases,
long-term stability is achieved but at the cost o f high com plexity of the subgridding
algorithm s [89][90][92] [96]. The schem es presented in [94] [97] are sim ple and o f low
reflections, but they do not guarantee long term stability.
In this chapter, a sim ple and stable subgridding schem e w ith low reflection is
proposed for 2D FD TD sim ulations. It is not only very sim ple w ith controllable low
reflections but also proven to be stable.
The chapter is organized as follow s. First, we present the proposed sim ple and stable
subgridding schem e and discuss how low reflections are achieved. T w o num erical
exam ples are then com puted to dem onstrate the validity and effectiveness o f the proposed
schem es. Finally, discussions and conclusions concerning this m ethod are presented.
6.2
Derivation of the Stable Subgridding Scheme
In a linear and lossless m edium , the spatially discretized M axw ell’s equations can be
expressed as follows:
ij+±,k+±
M
i,
dt
F
y i,j+y,k+1
-
A z
E
v
i,j+f . *
F
z /,;+ l,* + -L
-
F
z
i.j.k+l
A y
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
(6 . 1)
51
z
i+\,j,k+j
i+k.j.k+±
Q J
M
F
i+k,j,k+k
F
- FZ
U .k+ k
X
H-k,j,k+l
- F
x
i+i,j,k
( 6 .2 )
- Fv
i,j+i,k
(6.3)
1
dHy
dHz
i+k.j+i.k
Ax
F
/+ j,7 + j.*
x
i+-k,j+\,k
dt
Az
F
- Fx
y
/+ |,7 ,*
i+ lj+ j,*
Ay
Ax
dE.
e
(6.4)
H
i+ k j~ k ,k
i+i,j,k+±
H y
Ay
dE„
i.M .*
H
i+-l- i fc— L
i+i.j.k
Az
i.M.k
dt
(6.5)
l,J+k,k+-k
X
i,j+ \,k-j
Z
,+ 2
■J y
Ax
Az
£
Hy
dE, ;,M+i
i J M
i
( 6 .6 )
H
Hy i-\,j,k+k
i+k,j,k+k
H x i,i+\,k+\
Ax
Hx
i.j-
i.*+i
Ay
Equations (6.1)-(6.6) can be rew ritten as:
dB r
AyAz-
i.j+ j.k + j
dt
= Ay ( E y
i,j+ k,k + \
- Fy
i,7+-UJ
Az ( E z i . y + 1- , Fz
dB,.
AzAx-
i+i,j,k+i
= Az ( E z
dt
-A x (£
i+ \,j,k+ \
i,j,k+i
- Fz
i+ k j,k + l
i,j,k+ \
x
(6.7)
)
)
i+ \,j.k^
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
(6 .8)
52
dB,
AxAy-
i+ \,j+ k ,k
- Ax(E
-
i+±j+l,k
dt
-A y(Ey
dD
Ay Az-
i+ i.j .k
dt
■- AyAzJx
i+ \,j+ \,k
(6.9)
- FV
i+ k .j.k
&z(Hz i+jj+kk
( 6 . 10)
Ay ( f f ,
dDA
AzAx— ^ 7 2’— Az Ax Jy
dt
M+l
dt
Hy
i,j+ \,k
A x ( H x i,M ,t4 “ Hx
dD,
AxAy-
Fx
( 6 . 11 )
) - Az(H
- A x A y J z\.J k , =
Ay ( H y
,,k+k ) - ^ ( H x
(6 . 12)
i,j+ \,k + k
- H .x
i.j-U +P
where
B.
B.
B.
D,
D,
i.j+ k,k+ -k
i+ k ,j,k + k
=n
i+ \,j+ \,k
=v
i + i,j,k
i,j+ \,k
D, i.l.k+y
I i , j + 2 'k + 2
^ jl,J+ ^ ,* + -2
= £
= £
= £
i + i ,j .k + \
i+\ J + \ , k
i+ i.j,k
i.j+ i.k
F
z
i+ \,j+ \,k
(6.14)
(6.15)
x i+ i,j .k
(6.16)
v
i,j+ \,k
(6.17)
z
:J,k+ \
F
i.j.kH
H
(6.13)
Equations (6.7)-(6.18) can be rew ritten as:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(6.18)
53
a*,
(ey 1,7'+^,fc+1
ez i,M+p
/,7+U+j
>’
i.M .fc
i ,j+±,k+-k
) = -
(6.19)
dt
dby
i+lj.k')
(«x
(6 .20 )
) = -
i+l,7,fc+4
dt
dbz
(gv
;j+±k)
i +l,j+i:k
y
(*z
i+±,j+-k,k
—h.
(K
i,j+\,k+\
^
i+4,7+1,fc ^
i+^J+^.k
( 6 .21 )
dt
d d r i+^j.k
-i
2-7
i+kr.i.k-k')
i,j ,k
2
(6 .22)
i+ i.j.k
dt
J+i,k
i , j + - k , k —k
by
(K
b.
-h.
i h ±’ k) ~ ( h •’ i+\,i,k+\
1
^
z
i-k;,j+\,k')
ddz
i'k+7 ^x i-H-'+v ~ “a T
i - i j , k + ^ ~ ( b x
dt
i,j,k+-k
(6.23)
~Jy
i, j+ i,k
Jz
j>k+i
_ AyAz
i,j+i,k+±
n ij+ u + ih*
fa
i,j+\,k+\
AzAx
i+k,j,k+k
i+\,j,k+\ hy
Ay
i+ k j,k + i
AxAy
i+kj+-k,k
i+ j ,j + i ,k
Az
AyAz ^
fa
i+±j,k
_
1 ,
‘■J+r k
= -------- £
z i+ i,j+ i,k
p
i+^,j,k
*
i+ i,j ,k
A z A x
■-
Ay
i.j+ \,k e y i.j+i.k
AxAy
<-J-k+1z
A z
£
i,j,k+\ e z
,i,k+ \
where
R ep ro d u ced with p erm ission o f th e copyright ow ner. Further reproduction prohibited w ithout perm ission.
(6.24)
(6.25)
(6.26)
(6.27)
(6.28)
(6.29)
(6.30)
54
= AyAzB
(6.31)
i+x.j.t+x
- AxAyB
i + i,j+ i,k
i+ 4 ,j'+ 4 ,t
(6.33)
i+i , k = * y a z o , i+jjX
(6.34)
= AzAxD
(6.35)
i.jM i
i.j+ i.k
- AxAyD
(6.36)
= A >'AZi^ t+-i, /,fe
Jv
(6.32)
. ..
.
J y i t j + ± tk
= AzAx/,)
- AxAyJ
= AxE,
= AyE
u .>H= A zE-
i,j,k + j
i +4 r , j , k
i, i + ± , k
iJMi
i-s;
A x//,
i-c '
(+4j,i+i
'-s;
i+l M , k
i,j+j,k
(6.37)
(6.38)
(6.39)
(6.40)
(6.41)
(6.42)
(6.43)
i,j+ \,k + \
(6.44)
Ay H y
i + i,j+ i,k
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(6.45)
55
H ere e is denoted as an electric voltage along the edge of the m esh and b is the m agnetic
flux through a m esh cell facet. T he definitions o f h and d are the sam e as e and b except
they refer to m agnetic voltage and electric flux. ] is the current through a m esh cell facet.
Equations (6.19)-(6.44) can be w ritten as a m atrix form:
Ce =
db
dt
ru -
m
~
J
(6 -46)
d = D ee
b = D Mh
w here e is a vector containing all electric voltages, h is a vector containing all m agnetic
voltages, d is a vector containing all electric fluxes, b is a vector containing all m agnetic
fluxes, 7 is a vector containing all the currents through m esh cell facets. T he m atrices C
and C contain the signs for the sum m ation o f the field quantities. The diagonal m atrices
D£ and d
contain m aterial properties and m esh dim ensions.
E quation (6.46) can also be obtained by discretizing the follow ing integral form o f
the M axw ell’s equations directly, especially for non-uniform m eshes:
3
rr~
,77
= - ( j E*dl
(6.47)
= (§H »dl -
(6.48)
dt
—
dt
It is shown in [120] that the follow ing equation is a sufficient condition to ensure the
stability of (6.46):
C = CT
(6.49)
This stability condition can only ensure that the spatially dicretized system is stable.
W hen the tem poral discretization is also applied, an additional stability condition is
needed.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
56
The above concept is now applied to the construction o f a new 2D FD TD
subgridding schem e. For a subgridding structure w ith a ratio o f 3:1 show n in Figure 6.1,
by discretizing (6.47) directly and using equation (6.49) (with the help o f M aple), the
follow ing sim ple subgridding form ulations are obtained:
H ' ^ (2,1) = H"^ (2,1) +
fjAY
[£" (2,1) - E" (2,0)]
x
(6.50)
Af [e ;(2 , d -
fjA X
H
e ;(U )]
(1,2) = H ' ^ (1,2 ) + - ^ —[E: (1,2) - E"x (1,1)]
fjAY
(6.51)
_^
[£"(1’2)" £ "(0’2)]
C * (1,1) = hn" (1,1) + A . [< (1,1) - E: (2,1)]
flAy
(6.52)
- ^ - [ < ( l , l ) - £ ; ' (1,2)]
fjAx
hT" (2,i) = c * ( 2 , i ) + [<: ( 2 ,i)- £ ; (2, i)]
fjA y
(6.53)
fjAx
d , 2 ) = ( i , 2)+^
fjAy
’
[< a, 2) - < a, l)]
(6.54)
—
le” (1>2) —E" (1,2)]
fjAx
>
£ ;+1(2,i) = £;'(2,i)+
At [ ^ ( 1 , 1 ) + ^ ( 2 , 1 ) + ^ ( 3 , 1 )
eAy
3
j r H (2 1 )]
(6 -5 5 )
1
£"+1(1, 2) = £" (1, 2) Af X % l ) + C " (1 .2 ) + A : % 3 )
eAx
3
(6 -5 6 )
— - H -----------------■— i ----------*-------------- (1 ,Z)J
w here
Ex, £ v, and H , are field com p onents in coarse m esh es, ex, e , and hz are field com p onents
in subgridded m esh es, AX = 3A x, and AY = 3 A y .
A lthough the subgridding form ulations as derived above are stable w hen
the C FL
stability condition for fine m eshes is satisfied, num erical experim ents have show n that
the reflection from the interface betw een the coarse m eshes and the fine m eshes is too
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
57
large for general applications. In order to reduce the reflection from the interface w hile
retain the sim plicity and stability, a param eter a is introduced as follows:
At
[a E xn (2,1) - E" (2,0)]
juA Y 1
At
/zAX
x K’ J
n
(6.57)
[£ ;(2 ,i)-£ ;a ,D ]
+^
1 £ ;iU ) - £;<U)1
—
[«£; (1,2) -
piAx
(6.58)
(0,2)1
2(1,1)
=
+ A . K ( u ) - aE: ( i m
//Ar 1
*
(659)
(1, 2)]
h"** (2,1) = h'r* (2,1)
-<2.1)1
(6.60)
/C+=(l,2) = ^ ( l , 2 )
+^
[e:(1’2)“ e"(1,1)1
-
le"(1.2) - aE" (1,2)]
(6-61)
jj A x
En
; \ 2 , \ ) = E"x {2,\) +
At
+ hT*(2,1) + C * ( 3 ,l)
fAy
(6 6 2 )
3
-//;^ (2 ,1 )]
E"+l (1,2) = E" (1,2) -
A? rh
(1,1) + h
eAx
(1,2) + kn+*(1,3)
3
-fl!"4 (1,2)]
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(6.63)
---- P— -----►
— — •—
■ j= 6
:
■ *
•
t
*
.> •
»
t. •
i> «
> ♦
>
j'"-5
■ .H
j=3
•
E..( 1,2)1
i
'
.> •
•
-■ •
>- •
‘> 0
.
t «
>. •
i~ ^
H /1 .2 )
V
E x( l , l )
,
E v i l,! } I
H .O . I )
1=1
, > ♦
i ♦
■.
j= i
E 5( 2 , t )
•
0
H ,(2 , }
1=2
1=3
F igure 6.1 A subgridding structure with a mesh ratio of 3:1 in a FDTD lattice.
It can be verified using the M aple that the new subgridding form ulations (6.57)-(6.63)
still satisfy the stability condition (6.49). If the tim e step size Ar is determ ined by the
C FL stability condition o f fine m eshes, the tim e-dom ain m arching should be stable. This
has been verified by num erical exam ples. T herefore, the focus now is to m inim ize the
reflections due to the subgridding. T he trial-and-error technique is used to obtain a by
m inim izing the reflection from the interface based on the geom etry show n in Figure 6.2.
Table 6-1 gives the optim ized values o f a w e found for subgridding ratios 3:1, 5:1 and
7:1.
Table 6-1 The values o f a for subgridding ratios K =3:l, 5:1, and 7:1.
subgridding ratio K
The value of a
3:1
0.718
5:1
0.5875
7:1
0.509
It should be noted that only a single a w as introduced for the reflection
m inim ization. M ore param eters m ay be introduced in (6.57)-(6.63) w ith increased
com plexity of the optim ization process.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
59
6.3
N u m e ric al E x a m p le
To verify the proposed m ethod, w e first apply the proposed m ethod to a perfect parallel
plate w aveguide for the TM m ode that propagates in the x direction (see Figure 6.2). The
w idth betw een tw o plates is 10mm.
Source point
for Hj,
v |
Observation
line for
Subgridded
region
X
-rrr
- —*
v,
T
Source point for H
Figure
6.2
The
parallel
plate
waveguide
with
width
of
10mm
and
coarse
mesh
size
of
Ax= 1mm and A y =1 mm •
The coarse m esh size is A x = lm m and A y = lm m . 2 x 2 = 4 coarse m eshes in the center
are then replaced by a subgridded fine m esh o f Ax=-^ mm and A p j m m . K is the
subgridding ratio. T he tim e step for both the coarse m esh and the fine m esh is the same,
and At=AtCFL AtCFLis the m axim um tim e step determ ined for the fine m esh region based
on the CFL stability condition. The source is a G aussian pulse in tim e-dom ain. In order to
contruct propagating m ode before the subgridded fine region, the source points are
located in tw o sym m etric points ju st beside the tw o plates and are 20 m esh sizes away
from the subgridded fine region. The observation points are a line 5 m esh sizes away
from the fine region. T he coarse m esh size is equal to
at the m axim um frequency
considered (15G H z in this case); the m axim um reflections produced by the subgridded
fine m eshes are recorded and show n in Figure 6.3. As can be seen, the reflections are
below -63dB for the three given subgridding ratios, K = 3 : 1 , 5 : 1 , 7 : l . This is a strong
indication that the proposed m ethod is effective.
The second exam ple com puted is a 2D resonator o f 6mm x 5mm w ith a fin o f length
2m m in the m iddle (see Figure 6.4). The coarse m esh size is A x=lm m and A y = lm m . 2 x 2 = 4
coarse mesh around the end of the fin are replaced by a subgridded m esh o f A x = { m m
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
60
and Ay=jmm • T he tim e step for both coarse m esh region and fine m esh region is the same,
At=AtCFL and AtCFLis the m axim um tim e step determ ined for the subgridded fine m esh
region based on the C F L stability condition. The point source used in the FD TD
sim ulation equals 1 at n = l and equals 0 w hen n > l in tim e-dom ain. T able 6-2 shows the
relative errors o f the com puted resonant frequency o f the first m ode in com parisons with
the results o f the subgridding m ethod presented in [94], It can be seen from T able 6-2 that
the proposed m ethod gives alm ost the sam e results com pared w ith the low reflection
subgridding m ethod.
-60
-70
-80
-90
subgridding ratio=3:1
~
-100
subgridding ratio—5:1
subgridding ratio=7:1
<5 -120
-130
-140
-150
-160
f(GHz)
Figure 6.3 The reflection coefficient from the subgridded meshes in the center corresponding to the
parallel plate waveguide in Figure 6.2.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 6.4 The rectangular resonator o f 6mm x 5mm with a fin o f length 2mm in the middle.
It can also be found that the relative error o f the com puted resonant frequency
decreases as the subgridding ratio increases. H ow ever, this decrease w ill level o ff w hen
the subgridding ratio increases beyond a certain value. This is due to the fact that there
are tw o types o f errors in the subgridding schem e: the error due to the interface betw een
the coarse m esh and fine m esh and the error due to the approxim ation errors o f the FD TD
m ethod in com puting fine geom etric structures. W hen the error due to the FD TD
approxim ation is prom inent, increase in the subgridding ratio can reduce the overall error.
W hen the error due to the interface is prom inent, increase in the subgridding ratio will
have no effect in reducing the overall error.
To check the long tim e stability, we have tried different positions o f the fin and
different sizes o f the subgridded m eshes around the end o f the fin w ith iterations larger
than 500,000. In all cases, the m ethod is found to be stable. This has verified the stability
o f the proposed m ethod num erically.
Table 6-2 The relative errors o f the computed resonant frequency of the fin structure.
subgridding
ratio K
3:1
5:1
7:1
Errors of the method
in this chapter
-2.05%
-1.26%
Errors of the method
in [94]
-2.05%
-1.26%
The total time step
numbers used
1x8192
2x8192
-0.47%
-0.47%
2x8192
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
62
6.4
Conclusion
In this chapter, a sim ple and stable 2D FD TD subgridding schem e is proposed. T he term
stable doesn’t m ean the m ethod is absolutely stable for any tim e step size A t. It ju st
m eans it is stable w hen the tim e step size At is determ ined by the C F L stability condition
o f the fine m esh. This schem e adapts the stability condition proposed in [120] and
introduces a controlling param eter that m inim izes the reflections due to the subgridding.
N um erical experim ents w ere perform ed to confirm the stability and effectiveness of the
proposed subgridding schem e. H ow ever, further study is required to study the sensitivity
o f a to different structures. B ecause of its low reflection and stability, this new 2D
subgridding schem e is a very good option for the optim ization o f 2D photonic crystal
structures and H -plane filters.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
63
7 Extraction of Causal Time-Domain Parameters
Using Iterative Methods
7.1
Introduction
R ecent interest in the tim e-dom ain m odeling m ethods has been largely m otivated by the
dem ands for broadband electronic system s and high-speed digital circuits. W hen FD TD
is used to analyze circuits including active devices, the m ost convenient w ay is to solve
the tim e-dom ain lum p-elem ent circuit m odel of the active device directly. H ow ever, the
tim e-dom ain lum p-elem ent circuit m odel of an active device is often unavailable and
only the frequency-dom ain netw ork param eters are available. The frequency-dom ain
param eters need to be converted to tim e-dom ain param eters and the param eters in the
tim e-dom ain have to be causal in order to be com patible w ith the FD TD sim ulator.
U nder norm al circum stances, netw ork param eters are given or m easured only in the
frequency-dom ain and over a lim ited frequency range. For instance, m ost m anufacturers
provide the S-param eters for a FET transistor over only a lim ited frequency range under
certain biasing conditions. To convert these frequency-dom ain param eters to their tim edom ain counterparts, transform ation techniques such as the inverse F ourier T ransform
can be applied. H ow ever, direct application o f these transform ation techniques usually
leads to non-causal tim e-dom ain param eters.
For instance,
the
S-param eters
of
NE425S01 w ere given by m anufactures over the frequency range o f 0.5G H z - 18GHz
[121]. One can supposedly obtain the frequency-dom ain Y -param eters from the given
frequency-dom ain S-param eters [122] and then apply the inverse F ourier transform
directly to obtain the tim e-dom ain Y -param eters. Figure 7.1 show s the tim e-dom ain
y,,(r) obtained as such for N E425S01. As can be seen, y u (t) has som e significant values
at t<0 and therefore is non-causal. In other w ords, the sim ple inverse F ourier transform is
not adequate. Therefore, schem es need to be developed that can extract the causal tim edom ain param eters from the band-lim ited frequency-dom ain param eters.
It should be noted that the proposed m ethods do not guarantee the obtained
frequency-dom ain param eters outside the given frequency range are close to actual
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
64
values; how ever, it is not o f concern to a user.
0.5
>.
-0.5
1
-0.5
0
0.5
t(ns)
1
1.5
2
Figure 7.1 The time-domain y n ( 0 obtained with the direct inverse Fourier transform of the frequencydomain Y-parameters.
Perry et al. [123] and C hen [124] studied such phenom ena and proposed the use of
the H ilbert transform to tackle the problem . The m ethod proposed by Perry et. al.
achieves the causality by m aintaining the m agnitudes o f the original param eters but
m odifying their phases for the m inim um -phase system s [123]. T he technique proposed by
Chen has difficulties in handling the singularity of the H ilbert transform ation integrals
[124], The errors w ere found to be som ew hat large because o f the sim ple direct
extrapolation.
To obtain causal tim e-dom ain param eters from the band-lim ited frequency-dom ain
param eters w ith good accuracy in the given frequency range, two iterative techniques to
extract causal tim e-dom ain param eters are explored. They are conceptually sim ple and
easy to im plem ent. T he extracted tim e-dom ain param eters not only are causal but also
contain the sam e frequency-dom ain inform ation as the original param eters over the given
lim ited frequency range in both m agnitude and phase. C om prehensive num erical studies
on the effectiveness and validity o f the proposed approach are presented.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
65
7.2
The Error Feedback Based FFT Method
Suppose that Yon( / ) = Yn n r ( / ) + jY on j ( f ) is the given frequency param eter that is
know n over a finite frequency range o f interest. W e propose the follow ing procedure to
obtain its causal tim e-dom ain counterpart:
Step 1:
Let the values o f the target frequency-dom ain Y -param eters, denoted
as YT( f ) = YT r( f ) + jYT J ( f ) , be equal to Yon( f ) w ithin the given
frequency range and equal to zero outside the given frequency range.
Transform YT( f ) to its tim e-dom ain counterpart, denoted as y ( t ) , w ith a
straightforw ard inverse fast F ourier transform (IFFT). y(t) is usually
non-causal.
Step 2:
Force the causality into the obtained y(t) by sim ply setting all the values
to zero for t < 0 . That is, let y(t) = 0 for t < 0 .
Step 3:
Fourier Transform the tim e-dom ain param eters obtained in Step 2 back
into frequency-dom ain and obtain the corresponding frequency-dom ain
param eter,
denoted
as
Tmod( / ) = Fmod_r ( / ) + jYmod_ ,( / )
.
C om pute
differences or errors betw een Ymod( f ) and the original Yori( f ) over the
frequency range o f interest. T he difference is defined as
* Y ( f ) = Ymoi( f ) - Y ori( f )
(7-1)
T he error functions are defined as
m ax(A7r ( / ) )
erro r1 = ------- :---------- r-
n 2)
„
m a x (A T (/))
error2 = ------- — ------ :-
(7 3)
error = m ax {error\,error!)
(7.4)
H ere m ax m eans the m axim um value w ithin the frequency range o f
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
66
interest. If error is acceptably sm all, the tim e-dom ain param eters
obtained in Step 2 are the causal tim e-dom ain param eters that w e can
accept and the com putations term inate here; otherw ise, go to the next
step.
Step 4:
M odify
the
target
frequency-dom ain
Y -param eters
YT( f ) = YT_r( f ) + j YT i ( f ) in the frequency range o f interest as follows:
YT( f ) = YT( f ) —A Y ( / )
(7.5)
A Y ( f ) = AY r( f ) + j ( a A Y . ( f ) )
(7.6)
w here a is a coefficient used to balance the convergence speed betw een
the real part and im aginary part, a is selected em pirically in the
follow ing m anner: the initial a
is set to be
1;
in the subsequent
iterations, it is given by
a = a + (error2 —errorl)
E xtrapolate
the
target
(7.7)
frequency-dom ain
Y -param eters
YT{ f ) - YT r( f ) + j YT . ( f ) outside the frequency range o f interest w ith
the follow ing form ulas:
YT( f n) = YT( f n) - A Y ( f n)
(7.8)
n = m + 1, ..., N . f n is the FFT frequency point, n is the index in the
FFT. T he first frequency point f m is the border frequency point inside
and outside of the frequency range o f interest. N is the FFT index o f the
highest frequency point used in FFT. AY ( f n) is obtained w ith the
follow ing recursive form ulas:
^ Y ( f „ ) - 0.9* A F ( / n_j),
m +l<n< N
(7.9)
The factor 0.9 is used to sm ooth the change of the param eter values
outside the given frequency range. It is selected em pirically.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
67
Step 5:
Take the inverse Fourier Transform o f the frequency-dom ain param eters
obtained in Step 4 and go back to Step 2.
In the procedure described above, the error AY ( f ) is introduced into
the iterative loop.
Therefore, w e call the m ethod the error feedback based m ethod. The flow chart in Figure
7.2
illustrates the iterative procedure.
O ur num erical experience shows that the above iterations norm ally lead to the
convergence. H ow ever,
if the
know n param eters
are not the frequency-dom ain
param eters o f causal param eters, the iterations m ay not converge. To
deal w ith the
problem , like in any other iterative m ethods, a prescribed num ber can b e introduced into
the iterations such that the com putation will be term inated once the num ber o f iterations
exceeds the prescribed num ber. In the num erical exam ples later, the prescribed num ber is
set to be 500 w ith errors in (7.4) set to be less than 1%.
to Yoji(f) w ith in th e
cy ra n g e
In v e rs e F o u r ie r - tr a n s f o r m \^ ( f ) a n d o b ta in
th e tim e d o m a in yT(t)
F o rc e c a u s a lity b y s e ttin g y T(t)= 0 fo r t<0
F o u r ie i- tr a n s f o r m yT(t) to o b ta in its
fre q u e n c y d o m a in c o u n te r p a r t ’5(nod (f)
C o m p u te th e e r r o r s b etw een
;f) a n d
Y j(f) w ith in th e fre q u e n c y ra n g e .
A re th e e r r o r s sm a ll e n o u g h ?
C o r r e c t th e YT(f)
w ith th e e r r o r s
NO
Y ES
y T(t) is th e c a u sa l tim e -d o m a in p a r a m e te r
END
Figure 7.2 The flow chart o f the error feed-back based FFT method.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
68
7.3
The Hilbert Transform Based FFT Method
H ilbert transform can also be used w ith FFT to form another iterative m ethod for
extracting the causal tim e-dom ain param eters. For a real causal signal y ( t) , let Y ( f ) be
its Fourier transform or the frequency-dom ain counterpart:
Y ( f ) = F{ y ( t ) }
(7.10)
w here F{ } represents
the Fourier transform . Then the
im aginary part of Y ( f )
should satisfy the follow ing H ilbert transform relation [125]:
W ) = - Jt
x -L f-f
=
w here
real part o f Y ( f )
and the
=3r® W >
*¥
4T' = ~ ®
x -L f-f
(7-11)
W )
#
Yr( f ) = r e a l ( Y ( f ) ) ,
(7-12)
Yi ( f ) = i m a g ( X ( f ) )
and Y ( f ) = Yr( f ) + j Y i ( f )
.
®
represents the convolution in the frequency-dom ain.
T he inverse Fourier transform o f — — is a sign function:
¥
t >o
f l
s g n (,) =
,< o
<7-13)
- - ^ = F {s g n (0 }
(7 .1 4 )
-1
or sim ply,
Consequently, the H ilbert transform (7.11) and (7.12) can be rew ritten as:
Yr( f ) = F { F - ' { ^ - ® j Y i ( f ) } }
(7 .1 5 )
xf
j Yi( f ) = F { F - i[ ^ ® Y r( f ) } }
KJ
or,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
(7 .1 6 )
69
yr ( / ) = F{sgn(OF-1{ ^ ( / ) } }
(7.17)
y W ) = -H sgn(r)F-'{yr( / ) } }
(7.18)
W ith the application o f FFT, equation (7.17) and (7.18) can be expressed as
Yr( f ) = F F T ^ g n i t W F n j Y t f ) } }
(7.19)
j Y , ( f ) = FFT{sgn(t)IFFT{Yr( f ) } }
(7.20)
T he above relations are used as the basis for the iteration m ethod; therefore, the
m ethod is nam ed the H ilbert transform -based FF T m ethod.
In the follow ing paragraphs, the procedure for the H ilbert transform -based FFT
m ethod is described. To do so, let Yori r{ f ) and Yori , ( / ) be denoted again as the real and
im aginary parts o f the original frequency-dom ain Yon( f ) = Yon r ( / ) + jY ori , ( / ) that is
know n over a finite frequency range of interest. T he procedure involves five steps:
Step 1:
Set Yr{ f ) to be equal to Yori r ( f ) w ithin the given frequency range and
equal to zero outside the range. U se equation (7.20) to com pute the
im aginary part ! ) ( / ) of the frequency-dom ain Y param eter
Step 2 : Force Yt ( f ) to be equal to YoriJ( f ) w ithin the frequency range o f interest.
For the values outside the range, an extrapolation technique sim ilar to that
used in the previously discussed error feedback FFT m ethod is applied.
T hat is, Y j ( f ) outside the interest frequency range are updated w ith the
recursive formula:
(7.21)
w here n = m + l, ..., N . f n is a FFT frequency point that lies outside the
frequency range o f interest, n is the FFT index. The first frequency point
f m is the border frequency point o f the frequency range o f interest. N is the
FFT index o f the highest frequency point used in FFT. A! ( ( /„ ) is updated
w ith the follow ing recursive form ulas:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
70
AYi( f J = 0.9*(Yi( f n) - Y i (/„_,))
(7.22)
N ote that ^ ( / „ ) in the right-hand side o f (7.21) and (7.22) is the value
obtained in the previous iteration. The factor 0.9 is used to sm ooth the
change o f the param eter values outside the given frequency range. It is
selected em pirically.
Step 3 : C om pute the causal tim e-dom ain param eter y (t) and check the errors of
its frequency-dom ain counterpart:
A) Store the values o f F ( / ) in Yr o[d( f ) as tem porary values and use
equation (7.19) to com pute the new values Yr new( f ).
B) Set Yr( f ) = (Yr M ( f ) + Yrnew( f ) ) / 2 and take the inverse Fourier
T ransform of Y ( f ) = Yr( f ) + j Yi( f ) to obtain the corresponding tim edom ain param eter y { t ) .
C)
Force the causality into the obtained y(t) by setting all the values
to zero for t <
0
and then transform it back into the frequency-dom ain and
obtain the corresponding frequency-dom ain param eter that becom es the
new Y ( / ) = Yr( f ) + j Yi( f ) .
D)
C om pute
differences
or
errors
betw een
the
frequency-dom ain
param eter Y ( f ) and the originally specified or given frequency-dom ain
param eters Yori( f ) w ithin the frequency range o f interest. If the differences
or errors are acceptably sm all, the tim e-dom ain param eters y(t) are the
causal tim e-dom ain param eters w e can accept and the com putations
terminate here; otherwise, go to the next step.
Step 4:
Force Yr( f ) to be equal to Yari r( f ) w ithin the frequency range of interest.
To update the values outside the range, an extrapolation technique sim ilar
to that used in Step 2 is applied. That is, Yr( f ) outside the interest
frequency range are updated w ith the recursive form ula:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
71
W
J = W j - A r r( / n)
w here n = m + 1,
(7.23)
N . f n is a FFT frequency point that lies outside the
frequency range o f interest at the FFT index o f n. The first frequency point
f m is the border frequency point o f the frequency range o f interest. N is the
FFT index of the highest frequency point used in FFT. AYr( f n) is obtained
w ith the follow ing recursive form ulas:
A Yr( / J = 0.9*(Y r( / J - Y r( / B_1))
(7.24)
N ote that Yr{ f n) in the right-hand side o f (7.23) and (7.24) is the value
obtained in the previous iteration. T he factor 0.9 is used to sm ooth the
change o f the param eter values outside the given frequency range. It is
selected em pirically.
Step 5 : U se equation (7.20) to com pute the new Yt( f ) and go to Step 2.
T he factor 0.9 in step 2 and step 4 is used to sm ooth the change o f the param eter
values outside the given frequency range. It is a trial and em pirical selection.
The error functions m ust be defined in Step 3 to assess the differences betw een the
com puted frequency-dom ain Y param eter and the prior-know n or -specified frequency Y
param eter w ithin the frequency range of interest. In our case, we take follow ing error
functions, with the note that the other error functions can be used depending on a u ser’s
preference:
maxi
errorl
(7.25)
maxi
errorl
error = m zx ierro rl, e r ro r l)
(7.26)
(7.27)
T he iterative procedure above is illustrated by the flow chart in Figure 7.3. As
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
72
described above for the error feedback m ethod, the above iterations should be halted if
the com putations cannot converge after the num ber of iterations becom es too large or
exceeds a prescribed num ber. In our com putations described in the follow ing section, the
pre-set num ber is 500 and the error in (7.27) is set to be less than 1%.
Let Yr(f) be equal to
r(f) within
the frequency range
and use (7.21) to compute Y,(f)
Force Y,(f) to be equal to Y0_ yft within
the frequency range and ”
extrapolarion technique to obtain Y.(f)
outside the range
Store Yr{l) in Yr iW(f)
Use (7.20) to obtain Y r ww(f)
Use (7.21) to compute ¥
r 5ld<f>+V
Inverse Fourier-transform
Y(f)=Yr(f)-t-jYj(|) and obtain y(t)
Force the causality by setting
v(t)= 0 for t< 0
Fourier-transform y(t) and
obtain the new Y(f)
Compute the errors between Y(f) and Y ^ f )
within the frequency range,
Are the errors small enough?______
Force Yt(f) to be equal to
T(f) within
die frequency range and apply an
extrapolation technique to obtain Yr(f)
outside the range
NO
YES
v(t)is the causal
time-do main parameter
F igure 7.3 T he flow chart of the Hilbert transform based FFT method.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
73
7.4
Numerical Validation
To validate the tw o proposed m ethods, w e first apply them to tw o theoretically know n
cases: a tem poral triangular pulse and a tem poral rectangular pulse w ith their spectra only
partially know n. O nce the m ethods are proven valid, they are applied to a real FET
am plifier.
A) Triangular pulse
The triangular pulse considered is expressed as follows:
y(t)-- *
M
j
0
f <T
\t\
(7.28)
otherwise
w here T = 0.1 (ns) is the duration of the pulse in tim e. B y taking the F ourier transform , a
com plete spectrum o f the triangular pulse is obtained as
j l T f r , s in ^ T J )
Y ( / ) = (TVf •
'
/
/* \ 7
CT f ) 2
(7.29)
The pulse has a finite duration in the tim e-dom ain but extends to infinity in
frequen cy-dom ain.
Suppose now that Y ( f ) is only given or know n over OFlz - 6 GHz. T hat is, w e now
have
(7.30)
unknown
f > 6 GHz
To extract the pulse in the tim e-dom ain with the above lim ited inform ation, the
sim plest approach is to directly take the inverse Fourier transform o f Y ( f )
assum ption Of Y { f ) being zero beyond
6
with
GHz. The resulting time-domain signal is found
to be not causal. To ensure the causality, one can sim ply cut off the tim e-dom ain values
for t<0 (i.e. set the values for t<0 to be zero). W e call the approach the direct c ut-off
method. The causal tim e-dom ain signal obtained with the direct cut-off m ethod can be
converted to the frequency-dom ain to check its spectrum . T he results are show n in Figure
7 .4
, w hich dem onstrates that there is a big difference betw een the com puted spectrum
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
74
and the original value w ithin the frequency range o f 0-6G H z. T he m axim um relative
error is nearly
1 0 0 %.
A.
/V
------ ideal w aveform
------d ire c t cut-off
s"r
.
t(ns)
x 10
6
4
2
0
■2
th eo retical
d ire c t cut-off
-4
0
1
3
2
4
5
6
f(GHz)
x 10
-
------ th eo retical
------d ire c t cut-off
.
~
f(GHz)
Figure 7.4 The time-domain pulse and its Fourier transform obtained with the direct cut-off method.
ideal w aveform
fe e d b a c k FFT
Hilbert tran sfo rm b a s e d FFT
- -
—
theoretical
fe e d b a c k FFT
Hilbert tran sfo rm b a s e d FFT
-
theoretical
fe e d b a c k FFT
Hilbert tran sfo rm b a s e d FFT
Figure 7.5 The time-domain waveforms and its spectra extracted with the proposed two iterative methods
when the known frequency range is from 0 to 6GHz for the triangular pulse.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
75
Now the tw o m ethods proposed in this chapter are applied. Figure 7.5 show s the
results. B ecause the extracted tim e-dom ain signals are fully causal, their values for
negative tim e are equal to zero and are not plotted in the figures for space lim itation.
As can be seen, there is big difference betw een the extracted tim e-dom ain pulse and
the original triangular pulse. The reason is that the proposed m ethods only w arrant the
frequency-dom ain param eters are close to the actual values w ithin the given frequency
range but not outside the given frequency range. T he differences in the frequency-dom ain
w ithin the range o f 0 to 6 G H z are less than 1%.
B) Rectangular pulse
The second exam ple o f a rectangular pulse is defined as
1
0 <t<T
)
otherwise
(7.31)
where T = 0.1 (ns) is the duration of the pulse. Fourier transform o f the rectangular pulse
gives
(7.32)
Now suppose that Y ( f ) is know n only from 0 to
6 GH z,
application o f the direct
cut-off approach as described in the last section leads to the results show n in Figure 7.6.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
76
1.5 r
1
ideal w aveform
direct cut-off
f 0-5 :
0
-0.5
0
0.2
0.4
0.6
0.8
t(ns)
x 10
10
5
"to
theoretical
direct cut-off
0
■5
0
1
2
3
4
5
6
f(GHz)
theoretical
direct cut-off
f(GHz)
Figure 7.6 The time-domain waveform and its spectra extracted with the direct cut-off approach.
It can be seen from Figure 7.6 that the difference betw een the com puted spectrum
and the original value is large, especially the real part. H ow ever, the results from the
proposed m ethods are m uch closer to the original values w ithin the concerned frequency
range o f 0 to 6 GH z. They are show n in Figure 7.7.
From the tw o exam ples relating to triangular and rectangular pulses, w e can see that
although the extracted tim e-dom ain pulses are different from the original shapes, the
corresponding frequency spectra are very close to the original values w ithin the
concerned frequency range. Therefore, the extracted tim e-dom ain pulses show n in Figure
7.6 and Figure 7.7 can be used to adequately represent the original triangular and
rectangular pulses in term s o f the frequency range o f interest.
In sum m arizing the above results, we can conclude that the proposed iterative
m ethods are effective and useful in extracting tim e-dom ain signals from the given
frequency-dom ain inform ation in a limited frequency range o f interest. T he tim e-dom ain
signals extracted as such are causal and can be used to represent the original pulses in
light of the frequency range o f interest.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
77
1.5
1
ideal waveform
feed b a ck FFT
Hilbert transform b a s e d FFT
f 0-5
0
-0.5
(
10
e
Tao>
5
nu
theoretical
feed b a ck FFT
Hilbert transform b a s e d FFT
-5
(
0
theoretical
feed b a ck FFT
Hilbert transform b a s e d FFT
fe-0-5
<o
E
-1
0
1
2
3
f(QHz)
4
5
6
Figure 7.7 The time-domain waveforms and their spectra extracted with the two proposed iterative method
methods when the known frequency range o f interest is from 0Hz -6 G H z for the rectangular pulse.
In the follow ing section, we apply the proposed iterative m ethods to a linear FET
am plifier.
C) F E T a m p lifie r
The FET am plifier used NE425S01 as its active device [121], As indicated before,
the S-param eters of the active device are given by the m anufacturer only over the
frequency range o f 0.5G H z to 18GHz. They can be easily converted to the Y -param eters
in the frequency-dom ain. Since the given frequency data is quite sparse and not suitable
for FFT directly, the spline interpolation was applied.
N ow the proposed m ethods are applied to extract the causal tim e-dom ain param eters
o f the device. Figure 7.8 and Figure 7.9 show that the extracted tim e-dom ain param eters
are causal. T he spectra o f the extracted tim e-dom ain param eters and their com parisons
with the original param eters are shown in Figure 7.10. T he param eters extracted w ith the
two m ethods are sim ilar in general shapes in the tim e-dom ain but different in details. For
exam ple, y ,, (t) extracted with the H ilbert transform based FFT iterative m ethod has one
m ore sm all ripple at the initial tim e than the y u (t) extracted with the error feedback FFT
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
78
m ethod. These differences can be attributed to the fact that the different values o f Yn ( f )
outside the frequency range with the two m ethods.
x 109
2
cT“m
>
x 10
0
■2
0
5
CM
0.6
1.2
1.8
0.6
1.2
1.8
x 10
0
>.
•5
0
X 10
t(n s)
Figure 7.8 The time-domain Y-parameters extracted with the error feedback FFT method for the FET used
in the amplifier.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
Figure 7.9 The time-domain Y-parameters extracted with the Hilbert transform based FFT method for the
FET used in the amplifier.
Figure 7.11 shows the results o f the S-param eters converted from the Y -param eters.
As can be seen, the S-param eters extracted from the direct cut-off approach are vastly
different from the original data but the results from the proposed m ethods are very close
to the original data w ithin the frequency range of interest w ith their tim e-dom ain
counterparts being causal and com patible w ith the FD TD m odeling.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
80
0.2
0.1
0.1
—
original data
—-
error feedback FFT
—
Hilbert transform based FFT
—
original data
- -
error feedback FFT
—
Hilbert transform based FFT
- 0.1
10
0.5
15
18
10
0.5
15
18
x 10'
0.01
’
/
\
\
•10
—
original data
- -
error feedback FFT
—
Hilbert transform based FFT
/
h
~
cn
<0
E
/
V
10
15
original data
- -
error feedback FFT
—
Hilbert transform based FFT
0
-------------------------- 1
---------------------------1--------------- 1
-0.01 ------------------------1
0.5
15
18
-20 ---------------------------------------------------------------- 1------------- 1--------------0.5
—
10
18
0.5
3.5
—
~
0
: -0.5
—
—
original data
- -
error feedback FFT
—
Hilbert transform based FFT
original data
-
—
error feedback FFT
Hilbert transform based FFT
10
0.5
15
-0.5
0.5
18
0.1
10
15
18
15
18
0.05
—
—
—
original data
-
error feedback FFT
CM
CM
Hilbert transform based FFT
> ^ 0 .0 5
en
co
£
0
—
—
—
0.5
10
f (GHz)
15
18
-0.05
0.5
original data
-
error feedback FFT
Hilbert transform based FFT
10
f (GHz)
F igure 7.10 The Y-parameters in frequency-domain after using the iteration methods.
A fter the causal tim e-dom ain Y -param eters o f the FET are obtained, they are
included into an in-house tim e-dom ain sim ulator that is based on the m odified central
difference m ethod (M C D ) [126] for m odeling the overall am plifier. Figure 7.12 shows
the layout o f the am plifier. The characteristic im pedances o f m ain transm ission line TL1
and open-circuited stub line TL2 are both 5 0 Q . The phase velocities on TL1 and T L 2
are both 2 .1 2 5 3 5 x l0 8 (m/s). T he length of TL1 is 6.4(m m ) and the length o f TL2 is
3.6(mm).
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
81
2
4
rea)(S11)
2
co
0
0
2
0 .5
10
5
15
2
18
0 .5
0.2
5
10
15
18
5
10
15
18
5
10
15
18
10
15
18
real(s12)
0.2
CO
-0.2
-
0 .5
10
5
15
18
10
0.2
0 .5
reaKs21)
10
1—
CM
<0
I
-10
-20
0 .5
10
5
15
-10
0 .5
18
1
------------- original d a ta
------------- d ire c t c u t-o ff
— - — e r ro r f e e d b a c k
real(s22)
1
</>
0
— H ilbert tr a n s f o r m
0
1
0 .5
5
10
f (G H z)
15
18
0 .5
5
f (G H z)
Figure 7.11 The S-parameters after using the iteration methods and direct cut-off.
T
L
2
5------►
PO RT
1
TL1
N E425S01
PORT
2
o-------Figure 7.12 The Layout of the FET amplifier circuit.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
82
Figure 7.13 is the com puted overall S21 of the w hole am plifier circuit. For reference,
the A gilent A dvanced System D esign (a frequency-dom ain sim ulator) was also em ployed
to sim ulate the circuit. The results are also plotted in Figure 7.13. A s can be seen, the
results from M C D that em ploys the extracted causal tim e-dom ain Y -param eters are very
close to the results from ADS.
C\l
CO
$o
-10
-15
result from ADS
error feedback FFT
Hilbert transform based FFT
-20
-25
f ( GHz)
Figure 7.13 The computed S21 o f the amplifier in Figure 7.10 from different methods.
7.5
Discussion and Conclusion
In this chapter, tw o iterative approaches w ere proposed for extracting causal tim e-dom ain
param eters from their frequency-dom ain counterparts that are know n only over a
frequency range o f interest. Both m ethods are show n to be effective and useful w ith a
good degree o f accuracy of less than 1%. T he tim e-dom ain param eters extracted as such
can be included in a tim e-dom ain sim ulator that has a stringent requirem ent for tim e
causality.
It should be noted that the m ethod described in this chapter can only be applied to
small signal linear param eters. Extension to large signal nonlinear situations is the subject
o f future research.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
83
8 Extraction of Causal Time-Domain Parameters
Using Rational Function Approximation
8.1
Introduction
N etw ork param eters o f a lum p-elem ent device are usually given in a lim ited frequency
band o f interest, or an operation frequency range.
To include them in tim e-dom ain
sim ulations, they need to be converted to their causal tim e-dom ain counterparts.
In
chapter 7, tw o iterative m ethods w ere discussed. H ow ever, the resulting tim e-dom ain
param eters are in num erical data form at. T he convolution w ith them in the tim e-dom ain
is very tim e consum ing for long sim ulations. In [127], various rational function fitting
techniques except the vector fitting (VF) w ere discussed. In [128], the vector fitting was
developed w ith the aim o f better curve fitting o f frequency-dom ain data, but only
considered passive structures in frequency-dom ain.
In this chapter, w e propose the use o f the rational fitting technique to extract the
causal tim e-dom ain param eters o f a lum ped device, in particular an active device, from
their known band-lim ited frequency-dom ain counterparts. O ne o f the m ajor advantages
o f this approach is that the resulting tim e-dom ain param eters can be expressed in the
form of exponential functions. The convolution with these exponential functions can then
be perform ed in a recursive fashion w ithout requiring a com plete past history o f the tim edom ain param eters. T he C PU tim e for each tim e-m arching step is constant, and the CPU
tim e and m em ory can therefore be significantly reduced, especially for a sim ulation with
a large num ber o f iterations.
In the follow ing sections, the proposed rational fitting technique is described and the
recursive convolution is shown.
Finally, a num erical exam ple is presented to
dem onstrate its validity and effectiveness.
8.2
Causal Time-Domain Extraction With The Rational Fitting Technology
The objective in this section is to find rational functions in the frequency-dom ain that
m atch the original param eters w ithin the frequency range of interest. Suppose that
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
84
is an adm ittance param eter know n at frequency points
con( n
= 1,2,..., N F ) over a
specified frequency range o f interest. A rational function Y is used to approxim ate Y . It
should differ little o r not at all from the original param eter f ( s ) |J=;(B on the know n
frequency points.
The rational function can then be expressed as [127]:
(8 . 1 )
N ( s ) = a 0 + a ] s + a 2 s 2 + ... + a N s N
( 8 .2 )
+ ... + b N s N
(8.3)
D ( s ) = b0 + b ]s + b2s 2
w here s = jco = j
. a ,
j=012
N and bt,
i=012
N are the coefficients to be determ ined; N
is the order o f approxim ation and (N + l ) is not bigger than the num ber NF.
To find the coefficients at and bt , one m ay force Y( s ) to be equal to the original
param eter Y ( s )| s=jw at the given or know n N F frequency points f n. T hat is,
(8.4)
or
(8.5)
with sn - jcon = j ' l n f n and n= l,2,..., NF.
E quation (8.5) contains N F com plex equations, or equivalently, 2N F real equations.
How ever, equation (8.5) is hom ogeneous.
(8.5)
An additional condition is needed to m ake
inhom ogeneous so the solutions w ill not be trivial. To do so, w e can set bN = 1 [127].
As long as (N + l) < N F is set, the rem aining (2N + l) coefficients can then be solved by
using the least squares technique.
A fter the coefficients are found, the rational function is determ ined. To obtain the
corresponding tim e-dom ain expression, the poles o f the rational function need to be
decided. M any softw are packages can be used. In our case, M A TLA B was em ployed.
O nce all the poles are found, Y ( 5 ) can be rew ritten as:
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
85
^ ) = In =-\ ■»T 7I n + fco
( 8 -6 )
w here yn is the pole and Rn is the residue o f y (j) at s = yn .
In order to ensure the stability o f the tim e-dom ain counterpart o f 7 ( s ) and to
express it in an exponential form , any poles that are in the right h a lf s-plane need to be
rem oved or flipped to the left half plane, follow ed by a refitting o f the residues.
A lthough
the above procedure appears theoretically feasible, there are two
draw backs in practice: first, the process m ay suffer from num erical problem s such as illconditioning o f the system m atrix for high-order approxim ations, and secondly,
m ultiplication w ith the denom inator in (8.5) m ay result in a frequency-dom ain w eighting
that induces large errors or lim its the m ethod to a low -order approxim ation w hen the
frequency range is wide. To overcom e these problem s, the vector fitting proposed in [128]
is then applied. T he procedure is sum m arized below .
Suppose an,
n=12
M are the poles in the left h a lf s-plane after the poles o f ( 8 .6 ) in the
right h a lf s-plane are discarded. T he resulting Y ( 5 ) becom es:
M D'
f(f> = Z — z r + K
(8.7)
» -l S ~ a n
w here R 'n is the residue o f Y(s) at s - a n .
Let an,
i]=12
Mbe the starting poles. Then the follow ing three steps are taken:
Step # 1 : an unknow n auxiliary function cr(s) is introduced and o { s ) and cr(s)Y(s ) are
approxim ated with rational functions w ith the poles an,
n=12
M:
M
m
c
c x ( s ) Y ( s ) ^ — ^ + d = d ^ ----------n=i s - a
M cI K '- * - )
cr(5) = S - ^ = r + l = - ^ -----------
(8 .8 )
(8-9)
n -\
w here cn, cn, and d are the unknow n coefficients to be found. z„ and zn are the zeros of
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
86
crO)y(s) and c n » , respectively. They can be found once cn, cn, and d are determ ined.
Equation (8.9) is then m ultiplied w ith the original Y( s) and the resulting equation is
set equal to ( 8 . 8 ):
M c
Z ^ ^ + d = [ Z - ^ + 1] ^ )
^n
n=l
n=l $
(8 . 1 0 )
&n
Equation (8.10) is enforced at the know n frequency points and solved for
cn, c„ and d w ith the least squares techniques.
Step # 2 : ( 8 . 8 ) is divided with (8.9) :
M
n < '- o
Y(s) = d - %-----------
(8.11)
ri(j -z j
n—\
Equation (8.11) shows that the poles o f
7 (5 )
are equal to the zeros o f a ( s ) in (8.9),
w hich can be calculated from the solution o f ( 8 . 1 0 ) by solving an eigenvalue problem
[128],
Step # 3 : If E (s )|J=J(a obtained in (8.11) is not sufficiently close to the original param eters
w ithin the frequency range o f interest, an is replaced w ith znand another iteration
can be
started by going back to Step #1.
The final approxim ating function can be rew ritten in the follow ing form:
M r
7(5) = 7(5) = X ^ - + «
i=i s Pi
(8 . 1 2 )
w here r is the residue at pole p t . a is a constant.
U sing the inverse L aplace transform [129], the tim e-dom ain form o f (8.12) can be
expressed as:
M
y(r) = y(t) = aS{t) + u{t)'YJ riep>’
(8.13)
i= l
w here S(t) is the D irac D elta function and u(t) is the unit step function that ensures the
causality.
y( t )
is
the
causal
tim e-dom ain
representation
of
the
param eter 7 (j)| s^j01.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
original
87
Since (8.13) is in exponential form in the tim e-dom ain, its convolution in tim edom ain can be com puted very efficiently in a recursive fashion [130][131]. M ore
specifically, a current i(t) can be found as:
* (0
U Af = [ y ( 0 ® v ( 0 ] U Af
(gi4)
=«v(ouA
, +Zw(ou
w here v(r) is a voltage , ® is the convolution in tim e-dom ain, At is the tim e step used in
a tim e-dom ain sim ulator and k is the iteration num ber, and
= e p,Ay
/( 0
+ f A'
U(t_1)Af
(8.15)
e ^ kA,-T)v{T)dT
B y using the trapezoidal rule, the integration in the above equation can be num erically
evaluated as:
r
J ( k - l) A r
e p'aA,- r)v(T)dT
(8.16)
= Y [ ^ ( O U „ , + v ( r ) U Af]
The above recursive convolution takes a total o f 3 M k additions and (4M +1 )-k
m ultiplications to reach the k-th tim e step (i.e. the com putation tim e for each tim e step is
constant), w hile the num erical convolution involving a w hole past history o f a param eter
needs 0.5-fc (A:-l) additions and 0 . 5 k - ( k + 3) m ultiplications (i.e. the com putation tim e
for each tim e step increases w ith the iteration num ber). In other w ords, the com putation
tim e with the recursive convolution is proportional to the num ber of iteration, w hile the
regular convolution is proportional to the square o f the num ber of iterations. T hat is, for a
sim ulation w ith a large num ber o f iterations, the com putation tim e w ith the recursive
convolution will be m uch sm aller than that w ith the regular convolution.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
8.3
Numerical Validation
In order to verify the efficiency o f the proposed m ethod, it is tested w ith the sam e
am plifier show n in Figure 7.12 in chapter 7. As in chapter 7, the frequency-dom ain Y
param eters o f the FET transistor need to be converted into causal tim e-dom ain Y
param eters. In this case, N in (8.1) is taken to be 12 for Y l l ( f ) , Y 12(f), Y21(f), and
Y22(f). The num ber of frequency points provided by the m anufacturer is N F =21. The
com puted poles for the rational approxim ations o f Y l l ( f ) , Y 12(f), Y 21(f), and Y 22(f),
before the application o f the vector fitting as expressed by (8.7), are show n in T able 8-1.
The poles, residues, and the constant term c r, after application o f the vector fitting as
shown in (8.12), are listed in Table 8-2 and T able 8-3.
T he frequency-dom ain differences betw een the Y -param eters obtained w ith the
proposed m ethod and the original data are less than 1.3% w ithin the frequency range of
interest. T he difference or error is defined by:
maxi
errorl
maxi
(8.17)
maxi
error2
maxi
error = m ax(em ?rl, error 2 )
(8.18)
(8.19)
w here m ax m eans the m axim um value w ithin the frequency range o f interest. Yori_r( f )
and Yori , ( / ) are the real part and im aginary part o f the original frequency-dom ain data.
Yr ( / ) and Yt ( / ) are the real part and im aginary part of the rational approxim ation data.
Figure 8.1 shows the causal tim e-dom ain Y -param eters obtained w ith the proposed
technique ( a S ( t ) term is not plotted for clarity), w hile Figure 8.2 shows their frequencydom ain correspondents that includes the a S ( t) term . For com parison, the original data
are also plotted on the sam e figure. For further validations, the corresponding Sparam eters are also com pared in Figure 8.3
As can be seen, the param eters obtained w ith the proposed m ethod are basically
overlapping the m anufacturer’s data, w hile the results obtained w ith the direct cut-off
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
89
m ethod show significant differences w ithin the frequency range o f 0.5G H z to 18GHz.
Therefore, w e conclude that the causal tim e-dom ain responses represented in an
exponential form by (8.13) are valid and can be used to m odel the F E T in the tim edom ain, w ith little frequency-dom ain inform ation inaccuracy.
O nce the causal representation o f the Y -param eters o f the FET are obtained w ith the
proposed m ethod, they can be put into a FD TD based sim ulator for m odeling the overall
am plifier as show n in Figure 7.12. In our case, the m odified central difference m ethod
(M C D ) [126] was em ployed in constructing the sim ulator. Figure 8.4 is the com puted
overall S21 o f the am plifier. For reference, the A gilent A dvanced System D esign was
em ployed to sim ulate the circuit. T he results are also plotted in Figure 8.4. A s can be
seen, the two curves are very.
Table 8-1 The initial poles (Hz) o f Y l l , Y12, Y21, and Y 22 as in (8.7)
Yu( f )
-1.6390e+08
-2.637 le+08
-7.8294e+08
-2.3385e+08
-2.1594e+08
-2.841 le+08
-6.4923e+08
-3.4879e+08
-8.9298e+07
Y22( f )
Y2l( f )
Yn ( f )
±jl.0627e+ll
±j9.6443e+10
± j8.5994e+10
±i7.8650e+10
±jl.0652e+ll
±j9.7359e+10
±j8.6824e+10
± j7.9242e+10
±j6.5207e+10
-2.2991e+06
-1.7657e+06
-9.5890e+07
-5.5625e+08
-3.3949e+08
-6.9017e+06
± jl.l302e+ll
±jl.0639e+ll
±j9.6994e+10
±j8.6497e+10
± j7.9156e+10
±j6.5179e+10
-9.7073e+07
-1.5727e+08
-5.8508e+08
-3.4192e+08
±jl.0656e+ll
±j9.7509e+10
±j8.6887e+10
± j7.9307e+10
T able 8-2 The final poles and residuals of Y 1 1, Y12, Y21, and Y22 as in (8.12)
Poles (Hz)
-3.9005e+10
-2.7216e+11
-6.1162e+09
-6.5318e+09
-5.1854e+09
Y11(f)
Residues (Hz)
7.2486e+07
-8.8473e+09
-5.1316e+06 +j3.0844e+06
± j 6.2023e+10
+ j 7.4647e+10
-3.8004e+06 ±j2.9213e+06
6.7926e+08 +j4.1983e+07
+ j 8.2233e+10
Y21(f)
Residues (Hz)
6.7096e+06 +jl.3006e+06
±j4.2532e+10
3.7669e+07 +j8.6430e+06
±j6.3089e+10
Poles (Hz)
-3.4487e+09
-7.3245e+09
-1.3942e+09 ±j6.9142e+10
4.9240e+06 ± j2.301 le+04
-1 .8 9 4 7 e + 1 0 ± j7 .5 7 5 9 e + 1 0
1 .5 9 7 4 e + 0 8 ± |1 .7 4 8 1 c + 0 8
-5.1568e+09 ±j8.2266e+10
-6.9454e+!0 ± j 1.0174e+11
-2.8295e+09 +j8.0159e+08
2.0019e+09 +jl.3756e+09
Y12(f)
Residues (Hz)
±j6.3163e+10
1.9162e+05 ±j4.8670e+05
1,6009e+06 +j2.2171e+06
±j6.9695e+10
±j7.4970e+10
1.6074e+05 ±j4.1186e+04
-7.9590e+07 +2.1019e+06
±j8.2227e+10
-4.1183e+06 +jl,1193e+07
±jl.3625e+ll
Y22(f)
Poles (Hz)
Residues (Hz)
-2.5779e+10 ±j4.5645e+10
4.9684e+06 ±j5.2463e+07
-1.0385e+10 ±j6.7829e+10
-8.8947e+06 ±j4.6101e+05
-5.2475e+09 ±j8.2188e+10
3.2774e+08 ±jl.3077e+08
- 6 .6 4 5 0 e + 0 9 ± j l . l 9 1 5 e + l l
- 2 .4 2 1 0 e + 0 7 ± j 2 . 12 7 0 e + 0 7
Poles (Hz)
-4.6629e+09
-1.2964e+10
-2.2396e+09
-5.1509e+09
-6.4961e+09
Table 8-3 The constant term o f Y 1 1, Y12, Y21, and Y22 as in (8.12)
a
Y ll (f)
Y 12(f)
Y2l(f)
Y22(f)
2.6856e-002
-1.3213e-004
8.8978e-003
1.2232e-002
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
90
x 10
2
0
■2
0
0.1
0.2
0.3
0.4
0.5
t(ns)
0.6
0.7
0.8
0.9
1
1
0.2
0.3
0.4
0.5
t(ns)
0.6
0.7
0.8
0.9
1
0.2
0.3
0.4
0.5
t(ns)
0.6
0.7
0.8
0.9
1
x 10
5
c\j
>.
0
■5
x 10
1
C\l
0
1
Figure 8.1 The Y-parameters in the time-domain corresponding to Table 8-2 and (8.13) but without the
aS(t) term.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
91
0.2
0.1
original data
>
1c
P
0.1
vector fitting
o>
co
E
0
original data
vector fitting
-
0.5
10
5
15
18
0.1
0.5
0.01
10
5
15
18
15
18
15
18
15
18
0.01
original data
vector fitting
CM
CM
>
-
0.01
0
O)
«
original data
E
vector fitting
-0 .0 2 L-
0.5
-
5
10
15
0.01
0.5
18
0.5
5
10
0.5
original data
5
0
3 -0.5
0.5
vector fitting
CM
>
5
original data
CO
vector fitting
E
10
f (GHz)
CO
15
18
0.5
5
10
0.05
original data
b
CM
CM
vector fitting
CM
CM
>
CO
cc
0.05
co
original data
E
0.5
5
10
f (GHz)
15
18
vector fitting
-0.05
0.5
5
10
f (GHz)
F igure 8.2 The Y-parameters in the frequency-domain obtained with the proposed.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
92
1 1 i
onginal data
original data
direct cut-off
direct cut-off
W
---------- vector fitting
vector fitting
S'
co
g
10
0 .5
15
0
-2
0 .5
18
0.2
CM
cc
CO
original oaia
y
CD
-
0.2
15
18
15
18
0 .2 r -
CM
"co
10
original data
'5>
CO
direct cut-off
-
E
0.2
-
---------- vector fitting
10
0 .5
15
direct cut-off
vector fitting
-0.4
0.5
18
10
CM
- original data
original data
-10
direct cut-off
—
-10
vector fitting
-20
direct cut-off
- vector fitting
0 .5
0.5
3
- original data
original data
direct cut-off
—
- direct cut-off
CM
CM
- vector fitting
vector fitting
0 .5
0.5
f (GHz)
f (GHz)
F igure 8.3 S-parameters in the frequency-domain obtained with the proposed method and direct cut-off.
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
93
CM
result using MCD
result using ADS
-15
-20
-25
-30
0.5
f(GHz)
F igure 8.4 The overall S21 o f the amplifier.
To show the advantage o f the proposed m ethod in convolution, Figure 8.5
dem onstrates the com putation tim e versus the num ber o f iterations using recursive
convolution form ulas. For com parison, the C PU tim e using non-recursive convolutions
is also plotted.
As can be seen, the total com putation tim e with recursive convolution
using (8.14) and (8.15) is proportional to the num ber o f iterations (i.e., the C PU tim e for
every tim e-m arching step o f the M CD com putation is the same).
H ow ever, the
com putation tim e with non-recursive convolution is proportional to the square o f the
num ber of tim e-m arching steps (i.e., the CPU tim e for an M CD tim e-m arching step is
increased as the num ber o f tim e-m arching steps increases).
In other w ords, as the
num ber of tim e-m arching steps becom es larger and larger, the com putation w ith non­
recursive convolutions becomes slower and slower and will eventually come to stand-still.
W hen the num ber o f M C D tim e-m arching steps reaches the order o f 105, the saving in
C PU tim e w ith the recursive convolutions is approxim ately thirty tim es.
N evertheless, from observing Figure 8.5, one can see that there is a threshold of the
num ber of tim e-m arching steps when the recursive convolution starts to use less
com putation tim e than the non-recursive convolutions. This is because at the beginning,
Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.
94
non-recursive convolutions only involve a short history o f data w hile the recursive
convolutions alw ays involve M term s (see (8 .1 4 )). In this case, the threshold num ber is
approxim ately 1500.
900
l
<
1
I
I
1
I
"l '
'" 1
800
700 -
/
c 600 O
o
CD
O 500 E
(
_
/
----- MCD using recursivs convolution
—
MCD using numerical convolution
/
/
s
400 -
—
/
3
§ 300 o
/
/
//
-
s
200 -
-
100 -
-
0
D
1
2
3
4
5
6
iteration number
7
8
9
1
x 104
Figure 8.5 The computation time versus the number o f iterations for the amplifier.
8.4
C o n clu sio n
In this chapter, a m ethod using rational approxim ation was proposed for extracting casual
tim e-dom ain netw ork param eters. It results in a tim e-dom ain param eter in an exponential
form , which leads to an efficient recursive tim e-dom ain convolution com putation. The
saving in the C PU tim e in com parison w ith the non-recursive convolutions can be tens,
hundreds and even thousands o f tim es, depending on the num ber o f tim e-m arching steps.
N um erical exam ples have been provided to num erically verify the effectiveness and
accuracy o f the proposed m ethod.
It should be noted that the m ethod described in this chapter can only be applied to
sm all signal linear param eters. Extension to large signal nonlinear situations is the subject
o f future research.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
95
9 Conclusion and Future Work
9.1
Summary
The FD TD m ethod is one o f the m ost flexible and pow erful num erical tim e-dom ain
techniques to be w idely used to sim ulate com plex electrom agnetic system s and
R F/m icrow ave circuits, including passive structures and active devices. T he focus o f this
w ork can be divided into tw o parts:
the first relates to bringing about som e
im provem ents to the FD T D m ethod and the second relates to tim e-dom ain m odeling o f
active device.
In order to im prove the efficiency of the FD TD m ethod, four m ethods w ere proposed.
First, in chapter 3 a new com pact 2D FD TD m ethod was proposed, w hich can reduce a
3D w aveguide problem into a 2D problem that is uniform in one o f the transverse
directions. Second, a com pact ID FDTD m ethod w as proposed in chapter 4 that can be
used as incident w ave generators and m odal absorbing boundary conditions in uniform ly
filled w aveguides. Third, a new ID m odal PM L was proposed in chapter 5 that utilizes
the proposed ID FD T D form ula to reduce the 3D P M L into a ID m odal PM L. Finally, a
new 2D FD TD subgridding m ethod was proposed in chapter 6 that is not only stable and
sim ple, but also has low reflections.
W hen the FD TD m ethod is used to analyze R F/m icrow ave circuits w ith elem ent
active devices, it needs to convert the frequency-dom ain param eters into tim e-dom ain
param eters because the netw ork param eters o f m any elem ent active devices are given by
the m anufactures or m easured in frequency-dom ain and are band-lim ited. In order to
extract the causal tim e-dom ain param eters from the band-lim ited frequency-dom ain
param eters, three m ethods w ere proposed.
First tw o iterative m ethods based on FFT w ere proposed in chapter 7: one uses the
negative feedback and another uses the H ilbert transform technique. Then, the rational
function fitting technique was used to obtain a tim e-dom ain m odel from the band-lim ited
frequency-dom ain in chapter 8. The tim e-dom ain m odel obtained from the iterative
m ethods is in num erical data form at, while the tim e-dom ain m odel extracted from the
rational function fitting is in exponential function form .
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
96
W hen the extracted tim e-dom ain m odel are used in the tim e-dom ain sim ulator, for
exam ple, FD TD , there exists tim e-consum ing num erical convolution w ith these extracted
tim e-dom ain param eters. B ecause the convolution w ith the exponential function can be
realized very efficiently by a recursive form , the tim e-dom ain m odel extracted from
rational function fitting is m ore useful. H ow ever, the rational function fitting technique
has som e lim itation on applicable functions, the shape o f target function m ust be suitable
for the rational function approxim ation; w hile the iterative m ethods do not have this
lim itation.
9.2
Future Work
There are som e aspects for m ethods proposed in this thesis that are w orthy o f further
research:
a) For the com pact 2D FD TD m ethod, only a few sim ple exam ples w ere used. The
next step can use this m ethod in the optim um design of w aveguide devices. A
com prehensive analytical study of its stability condition and num erical dispersion
relationship is also needed.
b) For the com pact ID FD TD m ethod and the ID m odal PM L, hybrid absorbing
boundary conditions can be obtained by using the proposed ID m ethods to absorb
a few of the low er order m odes and the traditional P M L to absorb the other higher
order m odes.
c) For the 2D FD TD subgridding m ethod, further study on the extension to 3D
FD TD and on the sensitivity o f a on different structures is needed.
d) The three m ethods in chapter 7 and 8 can only be used for tim e-dom ain m odeling
o f active devices in sm all signal situations. H ow ever, active devices can also w ork
in a large signal m ode such as a transistor in a pow er am plifier or m ixer. Further
study is needed to extend these m ethods to tim e-dom ain large signal m odeling.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
97
References:
[1]
H arrington, R. F., Field Computation by Mo ment Methods, T he M acm illan Co.,
N ew York, 1968.
[2]
Am itay, N. and G alindo, V., "On Energy Conservation and the M ethod o f
M om ents in Scattering Problem s," I E E E Trans. Antennas and Propagat., vol.
17, pp. 747-7 5 1 , Nov. 1969.
[3]
Rahm at-Sam ii, Y. and M ittra, R., “Integral equation
solution
and RCS
com putation o f a thin rectangular plate,” IEEE Trans. Antennas Propagat., vol.
22, pp. 608-610, July 1974.
[4]
D avis, W. A. and M ittra, R., “A new approach to the thin scatterer problem using
the hybrid equations,” IE EE Trans. Antennas and Propagat., vol. 25, pp.402406, M ay 1977.
[5]
N ew m an, E. H. and Pozar, D. M ., "Electrom agnetic M odeling o f C om posite
W ire and Surface G eom etries," IE EE Trans. Antennas and Propagat., vol. 26,
pp. 784-7 8 9 , Nov. 1978.
[6]
Sarkar, T. K., "A N ote on the C hoice W eighting Functions in the M ethod of
M om ents," IEEE Trans. A ntennas and Propagat., vol. 33, pp. 436-441, April
1985.
[7]
Ney, M. M ., "M ethod o f M om ents as A pplied to E lectrom agnetic Problem s,"
IEEE Trans. Microwave Theory Tech., vol. 33, pp. 972-980, Oct. 1985.
[8]
N ew m an, E. H., "TM and T E Scattering by a D ielectric/Ferrite C ylinder in the
Presence o f a H alf-plane," IEEE Trans. A ntennas Propagat., vol. 34, pp. 804813, June 1986.
[9]
Leviatan, Y., H udis, E. and Einziger, P. D., "A M ethod o f M om ents A nalysis of
E lectrom agnetic C oupling T hrough Slots U sing a G uassian B eam Expansion,"
IEEE Trans. A ntennas Propagat., vol. 37, pp. 1537-1544, D ec. 1989.
[10]
G oggans, P. M. and Shum pert, T. H., "CFTE M M Solution for TE and TM
Incidence on a 2-D C onducting B oday w ith D ielectric Filled Cavity," IEEE
Trans. A ntennas Propagat., vol. 38, pp. 1645-1649, Oct. 1990.
[11]
Peters, M. E. and N ew m an, E. H ., "M ethod o f M om ents A nalysis o f A nisotropic
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
98
A rtificial M edia C om posed of D ieletric W ire O bjects," IE E E Trans. M icrow ave
Theory Tech., vol. 43, pp. 2023-2027, Sept. 1995.
[12]
Khalil, A. I., Y akovlev, A. B. and Steer, M. B., "Efficient M ethod-of-M om ents
Form ulation for the M odeling of Planar C onductive Layers in a Shielded
G uided-W ave Structure," IEEE Trans.
M icrow ave T heory T ech., vol. 47, pp.
1730-1736, Sept. 1999.
[13]
Li, J. and Li, L., "Electrom agnetic Scattering by a M ixture o f C onducting and
D ielectric
O bjects:
A nalysis U sing
M ethod
o f M om ents,"
IEEE
Trans.
M icrow ave T heory Tech., vol. 53, pp. 514-520, M ar. 2004.
[14]
Y uceer, M ., M autz, J. R. and Arvas, E., "M ethod o f M om ents Solution for the
R adar Cross Section o f a Chiral B ody o f Revolution," IEEE Trans. A ntennas
Propagat., vol. 53, pp. 1163-1167, Mar. 2005.
[15]
U beda, E. and Rius, J. M ., "Novel M onopolar M FIE M oM -D iscretization for the
Scattering A nalysis o f Sm all Objects," IEEE Trans. A ntennas Propagat., vol. 54,
pp. 50-57, Jan. 2006.
[16]
Rao, S., W ilton, D. and G lisson, A., "Electrom agnetic scattering by surfaces o f
arbitrary shape," IE E E Trans. A ntennas Propagat., vol. 30, N o. 3, pp. 409-418,
M ay 1982.
[17]
W eim in, S. and Balanis, C. A.,
"M FIE analysis and design
o f ridged
w aveguides," IEEE Trans. M icrow ave T heory Tech., vol. 41, No. 11, pp. 19651971, Nov. 1993.
[18]
H uddleston, P., M edgyesi-M itschang, L. and Putnam , J., "C om bined field
integral equation form ulation for scattering by dielectrically coated conducting
bodies," IE E E Trans. A ntennas Propagat., vol. 34, No. 4, pp. 510-520, April
1986.
[19]
Hubing, T. H ., Survey o f N um erical Electrom agnetic M odeling Techniques,
R eport N um ber TR 91-001.3, U niversity o f M issouri-R olla, EM C L aboratory,
1991.
[20]
W inslow , A. M ., "M agnetic field calculations in an irregular triangular m esh,"
Law rence R adiation Laboratory, Liverm ore, California, U C R L -7784-T , Rev. 1,
1965.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
99
[21]
Silvester,
P.,
"A
G eneral
H igh-O rder
Finite-Elem ent
A nalysis
Program
W aveguide," IE E E Trans. M icrow ave T heory Tech., vol. 17, No. 4, pp. 204-210,
A pr 1969.
[22]
Chari, M. V. K., B edrosian, G., D 'A ngelo, J. and Konrad, A., "Finite elem ent
applications in electrical engineering," IE E E Transactions on M agnetics, vol. 29,
No. 2, pp. 1306-1314, M ar 1993.
[23]
D aly, P., "H ybrid-M ode A nalysis o f M icrostrip by Finite-Elem ent M ethods,"
IE E E Trans. M icrow ave T heory T ech., vol. 19, No. 1, pp. 19-25, Jan. 1971.
[24]
Hara, M ., W ada, T., Fukasawa, T. and K ikuchi, F., "A three dim ensional
analysis o f R F electrom agnetic fields by the finite elem ent m ethod," IEEE
Transactions on M agnetics, vol. 19, N o. 6, pp. 2417-2420, Nov. 1983.
[25]
H ano, M ., "Finite-Elem ent A nalysis o f D ielectric-Loaded W aveguides," IEEE
Trans. M icrow ave Theory Tech., vol. 32, N o .10, pp. 1275-1279, O ct. 1984.
[26]
A ngkaew , T., M atsuhara, M. and K um agai, N., "Finite-Elem ent A nalysis of
W aveguide M odes: A N ovel A pproach that Elim inates Spurious M odes," IEEE
Trans. M icrow ave Theory Tech., vol. 35, No. 2, pp. 117-123, Feb. 1987.
[27]
Jin, J. M ., V olakis, J. L. and Liepa, V. V., "Fictitious absorber for truncating
finite elem ent m eshes in scattering," IE E Proc H M icrow aves A ntenna Propag.
vol. 139, No. 5, pp. 472-476, Oct. 1992.
[28]
Riblet, G. P., "Treatm ent o f tw o-dim ensional E-plane bandstop filters using the
finite elem ent m ethod," IEE Proc H M icrow aves A ntenna Propag. vol. 144, No.
6, pp. 411-414, Dec. 1997.
[29]
K oning, J., Rieben, R. N. and Rodriguez, G. H., "V ector finite-elem ent m odeling
o f the full-w ave M axw ell equations to evaluate pow er loss in bent optical
fibers," Journal o f Lightw ave Technology, vol. 23, No. 12, pp. 4147-4154, Dec.
2005.
[30]
G iannacopoulos, D. D., Kak, K. F. and M irican, B., "Efficient load balancing for
parallel adaptive finite-elem ent electrom agnetics with vector tetrahedra," IEEE
T ransactions on M agnetics, vol. 42, No. 4, pp. 555-558, April 2006.
[31]
Schneider, M. V., "C om putation o f Im pedance and A ttenuation o f TEM -Lines
by Finite D ifference M ethods," IEEE Trans. M icrow ave T heory T ech., vol. 13,
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
100
No. 6, pp. 793-800, Nov. 1965.
[32]
Beaubien, M. J. and W exler, A., "An A ccurate Finite-D ifference M ethod for
H igher O rder W aveguide M odes," IEEE Trans. M icrow ave T heory T ech., vol.
16, No. 12, pp. 1007-1017, Dec. 1968.
[33]
H ornsby, J. S. and G opinath, A., "Num erical A nalysis o f a D ielectric-L oaded
W aveguide w ith a M icrostrip L ine-Finite-D ifference M ethod," IE E E Trans.
M icrow ave Theory Tech., vol. 17, No. 9, pp. 684-690, Sep. 1969.
[34]
Beaubien, M. J. and W exler, A., "U nequal-A rm Finite-D ifference O perators in
the Positive-D efinite Successive O verrelaxation (PDSO R) A lgorithm ," IEEE
Trans. M icrow ave T heory T ech., vol. 18, No. 12, pp. 1132-1149, D ec. 1970.
[35]
C hrist, A. and H artnagel, H. L., "Three-D im ensional Finite-D ifference M ethod
for the A nalysis of M icrow ave-D evice Em bedding," IEEE Trans. M icrow ave
T heory Tech., vol. 35, No. 8, pp. 688-696, Aug. 1987.
[36]
Kim, C. M. and R am asw am y, R. V., "M odeling of graded-index channel
w aveguides using nonuniform finite difference m ethod," Journal o f Lightw ave
Technology, vol. 7, No. 10, pp. 1581-1589, Oct. 1989.
[37]
B eilenhoff, K., H einrich, W . and H artnagel, H. L., "Im proved finite-difference
form ulation in frequency dom ain for three-dim ensional scattering problem s,"
IEEE Trans. M icrow ave T heory Tech., vol. 40, No. 3, pp. 540-546, M arch 1992.
[38]
Rew ienski, M . and M rozow ski, M ., "Iterative application o f boundary conditions
in the parallel im plem entation of the FD FD m ethod," IEEE M icrow ave and
G uided W ave Letters, vol. 10, No. 9, pp. 362-364, Sept. 2000.
[39]
Zscheile, H., Schm uckles, F. J. and H einrich, W ., "Finite-difference form ulation
accounting for field singularities," IEEE Trans. M icrow ave T heory T ech., vol.
54, N o. 5, pp. 2000-2010, M ay 2006.
[40]
Lou, Z. and Jin, J. M ., "M odeling and sim ulation o f broad-band antennas using
the tim e-dom ain finite elem ent m ethod," IEEE Trans. A ntennas Propagat., vol.
53, No. 12, pp. 4099-4110, Dec. 2005.
[41]
Lem diasov, R. A., Obi, A. A. and Ludw ig, R., "Time dom ain form ulation o f the
m ethod o f m om ents for inhom ogeneous conductive bodies at low frequencies,"
IEEE Trans. M icrow ave T heory Tech., vo. 54, No. 2, Part 2, pp. 706-714, Feb.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
101
2006.
[42]
Yee, K. S., “N um erical solution o f initial boundary value problem s involving
M axw ell’s equations in isotropic m edia,” IEEE Trans. A ntennas Propagat., vol.
14, No. 3, pp. 302-307, M ay 1966.
[43]
Taflove, A. and U m ashankar, K., "A hybrid m om ent m ethod/finite-difference
tim e-dom ain approach to electrom agnetic coupling and aperture penetration into
com plex geom etries," IEEE Trans. A ntennas Propagat., vol. 30, N o. 4, pp. 617627, Jul. 1982.
[44]
Choi, D. H. and H oefer, W. J. R., "The Finite-D ifference-T im e-D om ain M ethod
and its A pplication to Eigenvalue Problem s," IEEE Trans. M icrow ave Theory
Tech., vol. 34, No. 12, pp. 1464-1470, Dec. 1986.
[45]
G w arek, W . K., "A nalysis of arbitrarily shaped tw o-dim ensional m icrow ave
circuits by finite-difference tim e-dom ain m ethod," IEEE Trans. M icrow ave
Theory T ech., vol 36, N o. 4, pp. 738-744, April 1988.
[46]
Chu, S. T. and Chaudhuri, S. K., "A finite-difference tim e-dom ain m ethod for
the design and analysis o f guided-w ave optical structures," Journal of L ightw ave
Technology, vol. 7, No. 12, pp. 2033-2038, Dec. 1989.
[47]
R eineix, A. and Jecko, B., "Analysis o f m icrostrip patch antennas using finite
difference tim e dom ain m ethod," IEEE Trans. A ntennas Propagat., vol. 37, No.
11, pp. 1361-1369, Nov. 1989.
[48]
Kim, I. S. and H oefer, W . J. R., "A local m esh refinem ent algorithm for the tim e
dom ain-finite difference m ethod using M axw ell's curl equations," IE E E Trans.
M icrow ave T heory Tech., vol. 38, No. 6, pp. 812-815, June 1990.
[49]
Sullivan,
D.
M .,
"A
frequency-dependent FD TD
m ethod
for
biological
applications," IEEE Trans. M icrow ave T heory Tech., vol. 40, No. 3, pp. 532539, M arch 1992.
[50]
Berenger, J. P., “A perfectly m atched layer for the absorption o f electrom agnetic
w aves,” Journal o f C om putational Physics, vol. 114, pp. 185-200, 1994.
[51]
Kapoor, S., "Sub-cellular technique for finite-difference tim e-dom ain m ethod,"
IEEE Trans. M icrow ave T heory Tech., vol. 45, No. 5, Part 1, pp. 673-677, M ay
1997.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
102
[52]
Noda, T. and Y okoyam a, S., "Thin w ire representation in finite difference tim e
dom ain surge sim ulation," IEEE T ransactions on P ow er D elivery, vol. 17, N o. 3,
pp. 840-847, July 2002.
[53]
K arkkainen, M. K., "FDTD surface im pedance m odel for coated conductors,"
IEEE Transactions on E lectrom agnetic C om patibility, vol. 46, No. 2, pp. 222233, M ay 2004.
[54]
K okkinos, T., Sarris, C. D. and E leftheriades, G. V., "Periodic FD TD analysis o f
leaky-w ave structures and applications to the analysis o f negative-refractiveindex leaky-w ave antennas," IEEE Trans. M icrow ave Theory Tech., vol. 54, No.
4, Part 1, pp. 1619-1630, June 2006.
[55]
Taflove, A., C om putational Electrodynam ics, N orw ood, M A: A rtech H ouse Inc.,
1995.
[56]
Zheng, F., Chen, Z. and Zhang, J., “T ow ard the developm ent o f a threedim ensional unconditionally stable finite - difference tim e - dom ain m ethod,”
IEEE Trans. M icrow ave T heory Tech., vol. 48, N o.9, pp. 1550-1558, Sept. 2000.
[57]
N am iki, T. and Ito, K., "Investigation o f num erical errors o f the tw o-dim ensional
A D I-FD TD m ethod," IEEE Trans. M icrow ave T heory Tech., vol. 48, N o. 11,
Part 1, pp. 1950-1956, Nov. 2000.
[58]
K rum pholz, M . and K atehi, L. P. B., "M RTD: new tim e-dom ain schem es based
on m ultiresolution analysis," IEEE Trans. M icrow ave Theory Tech., vol. 44, No.
4, pp. 555-571, A pril 1996.
[59]
Liu, Q. H., “T he pseudospectral tim e dom ain (PSTD): A new algorithm for
solution o f M axw ell’s equations,” IEEE A ntennas and Propagation Society
International Sym posium , vol. 1, pp. 122-125, July 1997.
[60]
Johns, P. B. and Beurle, R. L., “N um erical solution o f tw o-dim ensional
scattering problem s using a transm ission-line m atrix,” Proc. Inst. Elect. Eng.,
vol. 118, No. 9, pp. 1203-1208, 1971.
[61]
Sadiku, M . N. O. and Agba, L. C., "A sim ple introduction to the transm issionline m odeling," IEEE Transactions on C ircuits and System s, vol. 37, N o. 8, pp.
991-999, Aug. 1990.
[62]
Johns, P. B., "A Sym m etrical C ondensed N ode for the TLM M ethod," IEEE
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
103
Trans. M icrow ave T heory Tech., vol. 35, No. 4, pp. 370-377, A pril 1987.
[63]
Chen, Z., Ney, M . M. and H oefer, W . J. R., "A bsorbing and connecting
boundary conditions for the T L M m ethod," IEEE Trans. M icrow ave T heory
Tech., vol. 41, No. 11, pp. 2016-2024, Nov. 1993.
[64]
Righi, M ., H oefer, W. J. R., M ongiardo, M. and Sorrentino, R., "Efficient T L M
diakoptics for separable structures," IEEE Trans. M icrow ave T heory T ech., vol.
43, N o. 4, pp. 854-859, Part 1-2, A pril 1995.
[65]
Trenkic, V., C hristopoulos, C. and Benson, T. M., "D evelopm ent o f a general
sym m etrical condensed node for the TLM m ethod," IEEE Trans. M icrow ave
T heory Tech., vol. 44, No. 12, Part 1, pp. 2129-2135, Dec. 1996.
[66]
Pena, N. and Ney, M. M ., "A bsorbing-boundary conditions using perfectly
m atched-layer (PM L) technique for three-dim ensional TLM sim ulations," IEEE
Trans. M icrow ave T heory Tech., vol. 45, N o. 10, Part 1, pp. 1749-1755, Oct.
1997.
[67]
Pierantoni, L., Tom assoni, C. and R ozzi, T., "A new term ination condition for
the application o f the TLM m ethod to discontinuity problem s in closed
hom ogeneous w aveguide," IEEE Trans. M icrow ave Theory T ech., vol. 50, No.
11, pp. 2513-2518, Nov. 2002.
[68]
So, P. P. M., D u, H. and H oefer, W . J. R., "M odeling o f m etam aterials w ith
negative refractive index using 2-D shunt and 3-D SCN TLM netw orks," IEEE
Trans. M icrow ave T heory T ech., vol. 53, No. 4, Part 2, pp. 1496-1505, A pril
2005.
[69]
V alderas, D., Legarda, J., G utierrez, I. and Sancho, J. I., "D esign o f UW B
folded-plate m onopole antennas based on TLM ," IEEE Trans. A ntennas
Propagat., vol. 54, No. 6, pp. 1676-1687, June 2006.
[70]
K rum pholz, M . and Russer, P., "On the dispersion in TLM and FD TD ," IEEE
Trans. M icrow ave T heory Tech., vol. 42, No. 7, Part 1-2, pp. 1275-1279, July
1994.
[71]
Xiao, S., V ahldieck, R. and Jin, H., "Full-w ave analysis o f guided w ave
structures using a novel 2-D FD TD ," IEEE M icrow ave and G uided W ave
L etters, vol. 2, No. 5, p p .165-167, M ay 1992.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
104
[72]
Xiao, S. and V ahldieck, R., "An efficient 2-D FD TD algorithm using real
variables," IEEE M icrow ave and G uided W ave Letters, Vol. 3, N o. 5, p p .127129, M ay 1993.
[73]
H uang, T. W ., H oushm and, B. and Itoh, T., “Efficient m odes extraction and
num erically exact m atched sources for a hom ogeneous w aveguide cross-section
in a FD TD sim ulation,” IEEE M TT-S International M icrow ave Sym posium
D igest, Vol. 1, pp. 31-34, M ay, 1994.
[74]
Esw arappa, C., So, P. P. M . and H oefer, W . J. R., “New procedures for 2-D and
3-D m icrow ave
circuit analysis w ith the TLM
m ethod,” IE E E
M TT-S
International M icrow ave Sym posium D igest, vol. 2, pp. 661-664, M ay 1990.
[75]
M oglie, F., R ozzi, T., M arcozzi, P. and Schiavoni, A., “A new term ination
condition for the application of FD TD techniques to discontinuity problem s in
close hom ogeneous w aveguide,” IEEE M icrow ave and G uided W ave Letters,
vol. 2, No. 12, pp. 475-477, D e c ,1992.
[76]
H uang, T. W ., H oushm and, B. and Itoh, T., “The im plem entation o f tim edom ain diakoptics in the FD TD m ethod,” IE E E Trans. M icrow ave T heory Tech.,
vol. 42, No. 11, pp. 2149-2155, Nov. 1994.
[77]
Righi, M ., H oefer, W. J. R., M ongiardo, M . and Sorrentino, R., “E fficient TLM
diakoptics for separable structures,” IEEE Trans. M icrow ave T heory T ech., vol.
43, No. 4, Part 1-2, pp. 854-859, A pril 1995.
[78]
M rozow ski, M ., N iedz'w iecki, M. and Suchom ski, P., “A fast recursive highly
dispersive absorbing boundary condition using tim e dom ain diakoptics and
Laguerre polynom ials,” IEEE M icrow ave and G uided W ave Letters, vol. 5, No.
6, p p .183-185, June 1995.
[79]
Alim enti, F., M ezzanotte, P., Roselli, L. and Sorrentino, R., “M odal absorption
in the FD TD m ethod: A critical review ,” International Journal o f N um erical
M odeling, vol. 10, pp. 245-264, 1997.
[80]
Lord, J. A., H enderson, R. I. and Pirollo, B. P., “FD TD analysis o f m odes in
arbitrarily shaped w aveguides,” IEE N ational Conference on A ntennas and
Propagation, pp.221-224, 31 M arch-1 April 1999.
[81]
Alim enti, F., M ezzanotte, P., Roselli, L. and Sorrentino, R., “A revised
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
105
form ulation o f m odal absorbing and m atched m odal source boundary conditions
for the efficient FD TD
analysis o f w aveguide structures,” IE E E Trans.
M icrow ave T heory Tech., vol. 48, No. 1, pp. 50-59, Jan. 2000.
[82]
G edney, S. D., “An anisotropic perfectly m atched layer-absorbing m edium for
the truncation o f FD TD lattices,” IEEE Trans. A ntennas Propagat., vol. 44, no.
12, pp. 1630-1639, Dec. 1996.
[83]
K uzuoglu, M. and M ittra, R., “Frequency dependence o f the constitutive
param eters o f causal perfectly anisotropic absorbers,” IE E E M icrow ave and
G uided Letters, vol. 6, N o. 12, pp. 447-449, Dec. 1996.
[84]
Roden, J. A. and G edney, S. D., “C onvolution PM L (CPM L): an efficient FD TD
im plem entation o f the C F S-PM L for arbitrary m edia,” M icrow ave and O ptical
T echnology Letters, John W iley and Sons, vol. 27, No. 5, pp. 334-339, Dec.
2000 .
[85]
O koniew ski, M., Stuchly, M. A., M rozow ski, M. and D eM oerloose, J., “M odal
PM L,” IE E E M icrow ave and W ireless Com ponents Letters, vol. 7, No. 2, pp.3335, Feb. 1997.
[86]
Jung, K. and Kim, H ., “An efficient form ulation of a 1-D m odal PM L for
w aveguide structures,” M icrow ave and O ptical T echnology Letters, vol. 21, No.
l ,p p . 48-51, A pril, 1999.
[87]
Kim, I. S. and H oefer, W . J. R . , " A local m esh refinem ent algorithm for the tim e
dom ain-finite difference m ethod using M axw ell’s curl equations,” IEEE Trans.
M icrow ave T heory T ech., vol. 38, No. 6, pp. 812-815, June 1990.
[88]
Zivanovic, S. S., Yee, K. S., and Mei, K. K., “A subgridding m ethod for the
tim e-dom ain finite-difference m ethod to solve M axw ell’s equations,” IEEE
Trans. M icrow ave Theory Tech., vol. 39, No. 3, pp. 471-479, M arch 1991.
[89]
Thom a, P. and W eiland, T., “A consistent subgridding schem e for the finite
difference tim e dom ain m ethod,” Int. J. Num . M odeling, V ol. 9, No. 5, pp. 359374, 1996.
[90]
K rishnaiah, K. M. and Railton, C. J., “Passive equivalent circuit o f FD TD: An
application to subgridding,” Electronics Letters, vol. 33, N o. 15, pp. 1277-1278,
1997.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
106
[91]
O koniew ski, M ., O koniew ska, E. and Stuchly, M. A., “T hree-dim ensional
subgridding algorithm for FD TD ,” IEEE Trans. A ntennas Propagat., vol. 45,
N o.3, pp. 422-429, M arch 1997.
[92]
K rishnaiah, K. M. and Railton, C. J., “A stable subgridding algorithm and its
application to eigenvalue problem s,” IEEE Trans. M icrow ave T heory T ech., vol.
47, No. 5, pp. 620-628, M ay 1999.
[93]
W hite, M. J., Y un, Z. and Iskander, M. F., “A new 3-D FD TD m ultigrid
technique w ith dielectric traverse capabilities,” IEEE Trans. M icrow ave T heory
Tech., vol. 49, No. 3, pp. 422-430, M arch 2001.
[94]
K ulas, L. and M rozow ski, M., “A sim ple high-accuracy subgridding schem e,”
33rd E uropean M icrow ave Conference, M unich, G erm any, pp. 347-350, Oct.
2003.
[95]
Podebrad, O., C lem ens, M., and W eiland, T., “N ew flexible subgridding schem e
for the finite integration technique,” IEEE T ransactions on M agnetics, vol. 39,
No. 5, pp. 1662-1665, M ay 2003.
[96]
M arrone, M. and M ittra, R., “A new stable hybrid three-dim ensional generalized
finite difference tim e dom ain algorithm for analyzing com plex structures,” IEEE
Trans. A ntennas Propagat., vol. 53, No. 5, p p .1729-1737, M ay 2005.
[97]
K ulas, L. and M rozow ski, M., "Low -reflection subgridding," IE E E Trans.
M icrow ave T heory Tech., vol. 53, No. 5, pp. 1587-1592, M ay 2005.
[98]
D onderici, B. and Teixeira, F. L., "Im proved FD TD subgridding algorithm s via
digital filtering and dom ain overriding," IEEE Trans. A ntennas Propagat., vol.
53, N o.9, pp. 2 9 3 8-2951, Sept. 2005.
[99]
Sui, W ., C hristensen, D. and D urney, C., “E xtending the tw o-dim ensional FD TD
to hybrid electrom agnetic system s w ith active and passive lum ped elem ents,”
IEEE Trans. M icrow ave Theory Tech., vol. 40, pp. 724-730, 1992.
[100] Piket-M ay, M . J., Taflove, A. and Baron, J., “FD -TD m odeling o f digital signal
propagation in 3-D circuits w ith passive and active loads,” IEEE Trans.
M icrow ave T heory T ech., vol. 42, pp. 1514-1523, 1994.
[101] Kuo, C., W u, R., H oushm and, B. and Itoh, T., "M odeling o f m icrow ave active
devices using the FD TD analysis based on the voltage-source approach," IEEE
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
107
M icrow ave and G uided W ave Letters, vol. 6, No. 5, pp. 199-201, M ay 1996.
[102]
Zhang, J. Z. and W ang, Y. Y., “FD TD analysis o f active circuit based on the sparam eters,” Proceeding o f 1997 A sia-Pacific M icrow ave C onference, v o l.3,
p p .1049-1052, H ong Kong, D ec. 1997.
[103]
Ye, X. and D rew niak, J. L., "Incorporating tw o-port netw orks w ith S-param eters
into FD TD," IE E E M icrow ave and W ireless C om ponents Letters, vol. 11, No.
2, pp. 77-79, Feb. 2001.
[104] Luo, S. and Chen, Z., "A N ew 2D FD TD M ethod for Solving 3D G uided-w ave
Structures," IEEE M TT-S International M icrow ave Sym posium , pp. 1473-1476,
June 11-16, 2006.
[105]
Luo, S. and Chen, Z., “Efficient one-dim ensional FD TD m odeling o f w aveguide
structures,” Proceeding of Frontiers in A pplied C om putational E lectrom agnetics
2006, V ictoria, BC, Canada, pp. 146-149, June 19-20, 2006.
[106]
Luo, S. and Chen, Z., “An efficient m odal FD TD for absorbing boundary
conditions and incident w ave generator in w aveguide structures," has been
scheduled for publication in Progress In Electrom agnetics R esearch, vol. 68, pp.
229-246, 2007.
[107]
Luo, S. and Chen, Z., “A FD TD -based M odal PM L," IE E E M icrow ave and
W ireless C om ponents Letters, vol. 16, No. 10, pp.528-530, Oct. 2006.
[108] Luo, S. and Chen, Z., “H igh-perform ance one-dim ensional m odal absorbing
boundary conditions for FD TD m odelling o f w aveguide structures," T he 12th
International
Sym posium
on
A ntenna
Technology
and
A pplied
E lectrom agnetics, (A N T EM /U R SI 2006), M ontreal, QC, Canada, July 16 - 19,
2006.
[109] Luo, S. and Chen, Z., "Extracting Causal Tim e D om ain Param eters," T he 10th
International
Sym posium
on
A ntenna
Technology
and
A pplied
Electrom agnetics, (A N TEM 2004/U R SI), O ttaw a, ON, C anada, July 20 - 23,
2004.
[110]
Luo, S. and Chen, Z., "Iterative m ethods for extracting causal tim e-dom ain
netw ork param eters," IEEE Trans, on M icrow ave Theory and T echniques, vol.
53, no. 3, pp. 969-976, M ar. 2005.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
108
[111]
Luo, S. and Chen, Z., "Extraction o f causal tim e-dom ain netw ork param eters
using rational functions," IEEE T ransactions on C ircuits and System s I, V ol. 52,
No. 6, p p .1205-1210, June 2005.
[112]
Luo, S. and Chen, Z., "Causal param eter extractions by vector fitting for use in
tim e-dom ain num erical m odeling," IEEE A ntennas and P ropagation Society
International Sym posium , vol. 3A, pp. 325-328, 3-8 July 2005.
[113]
Arcioni, P., Bozzi, M ., Bressan, M ., Conciauro, G. and Perregrini, L., "Fast
optim ization, tolerance analysis, and yield estim ation of H -/E -plane w aveguide
com ponents w ith irregular shapes," IEEE Trans. M icrow ave T heory T ech., vol.
52, No. 1, pp. 319-328, Jan. 2004.
[114]
Ros, J. V. M ., Pacheco, P. S., G onzalez, H. E., Esbert, V. E. B., M artin, C. B.,
Calduch, M. T., B orras, S. C. and M artinez, B. G., "Fast autom ated design of
waveguide filters using aggressive space m apping with a new segm entation
strategy and a hybrid optim ization algorithm ," IEEE Trans. M icrow ave Theory
Tech., vol. 53, No. 4, pp. 1130-1142, A pril 2005.
[115]
Bandler, J. W ., M oham ed, A. S. and Bakr, M. H., "TLM -based m odeling and
design exploiting space m apping," IEEE Trans. M icrow ave T heory Tech., vol.
53, No. 9, pp. 2801-2811, Sept. 2005.
[116]
Navarro, E. A., G im eno, B., Cruz, J. L., "Analysis o f H -plane w aveguide
discontinuities with an im proved finite-difference tim e dom ain algorithm ," IEEE
Proceedings H, M icrow aves, A ntennas and Propagation, vol. 139, No. 2, pp. 183185, A pril 1992.
[117]
Ishii, T. K., M icrow ave Engineering, N ew York: T he R onald Press Com pany,
1966.
[118]
Loh, T. H. and M ias, C., "Im plem entation o f an exact m odal absorbing boundary
term ination condition for the application of the finite-elem ent tim e-dom ain
technique to discontinuity problem s in closed hom ogeneous w aveguides," IEEE
Trans. M icrow ave T heory Tech., vol. 52, No. 3, pp.882-888, M ar. 2004.
[119]
Lou, Z. and Jin, J. M ., "An accurate w aveguide port boundary condition for the
tim e-dom ain finite-elem ent m ethod," IEEE Trans. M icrow ave T heory Tech.,
vol. 53, No. 9, pp.3014-3023, Sept. 2005.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
109
[120]
Thom a, P. and W eiland, T., "N um erical Stability of Finite D ifference Tim e
D om ain M ethods," IEEE Transactions on M agnetics, vol. 34, No. 5, Part 1, pp.
2740-2743, Sept. 1998.
[121] C alifornia Eastern Labs, D ata sheet o f low noise am plifier N E425S01 [Online],
A vailable: w w w .cel.com /pdf/datasheets/ne425sO 1.pdf [2007, 23 January].
[122] Pozar, D. M ., M icrow ave Engineering, Don M ills, ON: A ddison-W esley, 1990.
[123] Perry, P.
and B razil, T., "Forcing causality on S-param eter data using the
H ilbert transform ," IE E E M icrow ave and G uided Letters, Vol. 8, No. 11, pp.
378-380, Nov. 1998.
[124] Chen, Y., "Design and Sim ulation o f A ctive Integrated A ntennas," M .A.Sc.
Thesis, D alhousie U niversity, 2003.
[125] Papoulis, A., T he Fourier Integral and Its A pplications, N ew York: M cG raw Hill, 1962.
[126] M urakam i, K., H ontsu, S. and Ishii, J., "Transient analysis o f a class o f m ixed
lum ped and distributed constant circuits," Proceedings o f ISC A S'88, pp28432846, June 1988.
[127] M iller, E. K., "M odel-B ased Param eter E stim ation in Electrom agnetics: part I.
B ackground and T heoretical D evelopm ent," IEEE A ntennas and Propagation
M agazine, vol. 40, pp. 42-52, Feb. 1998.
[128] G ustavsen, B. and Sem iyen, A., "Rational A pproxim ation o f Frequency D om ain
R esponses by V ector Fitting," IEEE Trans, on Pow er D elivery, vol. 14, pp.
1052-1061, July 1999.
[129] K uhfitting, P. K. F., Introduction to the Laplace Transform , N ew York, NY:
Plenum Press, 1978.
[130] Sem iyen, A. and D abuleanu, A., "Fast and A ccurate Sw itching T ransient
C alculations on Transm ission Lines w ith G round Return U sing R ecursive
Convolutions," IEEE Trans, on Pow er App. A nd Syst., vol.PA S-94, N o .2, pp.
561-571, M arch/A pril 1975.
[131]
Lin, S. and Kuh, E. S., "Transient Sim ulation o f Lossy Interconnects B ased on
the R ecursive C onvolution Form ulation," IE E E Trans, on C ircuits and System sI: Fundam ental T heory and A pplications, vol.39, p p .879-892, N ov. 1992.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
110
Appendix A: The
Numerical
Dispersion
of
Compact ID FDTD for TEmn mode
in a Rectangular Waveguide
Suppose the rectangular w aveguide has w idth a in x direction and height b in y direction.
The field com ponents for the TEmn m ode along z-direction can be written as:
E x = E x0 cos(kxx ) s m ( k yy ) e ja’-z~w,)
E v = E y0 sin(kIx ) c o s ( k yy ) e nk*-a,)
Ez = 0
,
(A .l)
H x = H x0 s i n ^ x ) co s (k yy ) e ] z
H y = H y0 cos(kxx ) s m ( k yy ) e j{kzZ~ea)
H z = H z0 cos(kxx) cos(kyy ) e j(kzZ~M)
/Tl/T"
wit
where kx = — , k = — , k_ is the spatial frequency in the z direction, and ca is the
' a
'
b
tem poral angular frequency.
Substitution o f (A .l) into (4.3) reads
Ex0 cos[kx( i - ^ ) A x ] s i n ( k J A y ) e j[kzkAz-‘0{n+l)A']
- Ex0 cos[kx (i - ~ ) ^ ] sin(&vj Ay)ejlkzkAz~a*A,]
A t
1
j[k,kAz-w{n+2)At]
1
+ - ( a zy\ ^ ^ - \ ) H z0cos[kx( i - - ) A x ] c o s [ k y( j - - ) A y ] e
sAy
' rJ 2
2
2
‘
2
At
1
jlk,(*+—)Az-o(n+-l)Al]
— — { H y0cos[kx( i - - ) A x ] s i n ( k j A y ) e
2
2
£A z
2
1
- H y0cos[kx( i - - ) A x ] s i n ( k yj A y )e
j[k , ( * — )A z -fi> (n + 2 )A z ]
2
2
}
w here
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
( a .2 )
Ill
1
1
H :o cos[kx( i- - ) & x ] c o s [ k y( j + -)A y] e
a
I
u zr ' i - L j _L
-
=
1
j[l,lA z-(a(n-A )A r]
'
2
__________________________ 2 __________________ :__________2 _________________________________
,
.
1
1
H z0cos[kx( i - - ) A x ) c o s [ k y( j - - ) A y ] e
I
jlk.k&z-aKn^At]
'
2
(A.3)
cos[ky( j + ^)Ay]
cos[ k y( j - ± ) A y ]
D ivision o f both sides o f equation (A .2) by cos[^((-^)A x]e;lMi' oKn+2 w leads to
i,
. i.
jtu-At
jio-At
E x0s in ( k yj A y ) ( e
2 -e 2 )
=
(oczy l(H
At
- l ) / / z0c o s [ ^ v( y - ^ - ) A y ]
(A.4)
jk-—Az
-jk.—Az
■Hy0s in{ kyj A y ) ( e ‘2 - e
2 )
eAz
By substitution o f (A.3) into (A.4), (A.4) can be rew ritten as:
. i
. i4
jo>-At
jco-At
E x0s i n ( k J A y ) ( e
2 -e 2 )
-
A - {c o s [ky ( j + i ) Ay] - co s
s Ay
2
At
{ j - I ) A y ] } H z0
2
j k . —Az
H y0s in ( k yj A y ) ( e ' 2
(A.5)
- j k . —Az
-e
2 )
eAz
B y applying the follow ing identities:
e jd = cos(B) + j s in (# )
..
_ . A + B . A —B
cos(A ) - c o s( B) - - 2 s i n — - — s i n — - —
(A.6)
E quation (A.5) can be sim plified to:
j .
_sm(—
1 . k Ay
j
— sm(— ) H I0 - _
. k zAz
sin( - ^ - ) ^ = 0
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
{A.l )
112
Equations (4.4)(4.8) can be treated sim ilarly. They becom e:
-/-s in A f i . , +
At
2
eAz
1
• A ^m r
liA z
-sin (-^ —
1
2 vu
1 -s in ^ )£ ,0
2
//Ay
1
•/ C O A t
sin(—— ) H
At
2
—
x0
n
= 0
m
juAx
2
■
(A.8)
(A.9)
1 . ,coAt^TJ
- _s i n( — )//* =0
A: Az
.
—
) E y0 +
+ -i-s in (i^ )H ,0 = 0
eAx
2
2
(A -10)
At
sin(^ r 0 A o = 0
2
(A.l 1)
The above equations form a system o f five hom ogeneous equations w ith unknow ns
E xo , E y o ,
H xo, H yo,
and
H zq .
B ecause the solutions o f the system m ust not be trivial, the
determ inant of its coefficient m atrix should be equal to zero. This leads to:
s in (-^ -) = 0
(A. 12a)
. 2//c AZx
sin (—— )
. 2,o At .
//e s in (------)
A z2
A r2
. 2 k Ax
. 2 k Ay .
k Az
.
sin (—— ) sin (-^— ) sin (—5— )
Lie sin (------)
2
■
2
.
2
_
2
Ax2
Ay2
Az2
At2
.
,
w here k
m il
..
(A. 12b)
,. , _ N
(A. 12c)
n7t
= — and k v - — .
a
b
Equation (A. 12a) corresponds to/y = 0, and represents the static solution.
From equation (A. 12b), it can be obtained that
Ar2 sin2( - ^ ) = //£Az2sin2( ^ )
From equation (A. 10), it can be obtained that
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
(A. 13)
113
//A z sin f
(A. 14)
Substitution of (A. 14) into (A.7) results in:
(A. 15)
eAzAt2 sinC^ ^ )
2
B ecause o f equation (A. 13), the first term in equation (A. 15) equals zero. This will
lead to
H ,a
= 0 , w hich does not agree w ith the assum ption o f T E m odes.
It can be seen that equation (A. 12c) approaches the analytical dispersion relationship,
co1/U£ - k 2x + k 2 + k 2, w hen Ax, Ay, Az and At approach zero. Therefore, equation (A. 12c)
is the num erical dispersion relationship o f the TE m odes for the com pact ID FDTD.
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
114
Appendix B: The Formula of ID Modal CFS-PML
The com plete form ulae o f ID m odal C F S-P M L schem e w ith s =1 , 5 = 1 , and
S z = tcz +
can be w ritten as:
a + jm 0
p
n+i
i - ± J 'k
At
"
-(«7 I
p A v
£A>>
' 2 O ’*
x
, -l)Hz
1 1’J 2
■
z
17
1 2’J
1(
2
(B .la )
At
n+1/2
-H ,
£A z (^ V
£0k z
n+l/2
)
a zk z + a z
At
-
i - i ,i,k
At
— —
At
H
(B .lb )
or.
2
i
B+l
At
'-iO ’*
gnfcT | Qrzfcz + cxz
At
At
n+1
_ D
n+-
(H
,
(B.2a)
At
-(«
I.
fA x
z
e0K z
£U
0K zZ _j cci._ K z7______
+ a z_
At
or.
At
j
-
,
- l)//z 17.
'-TO-f*
a zK z + a z
At
■ H*
,
y
'J-iX
2
(B.2b)
At
2
2
At
At
"+1
_ D
U X -j
z
(B.3a)
At
eAy
I- •
xy :’J-
■ l' ) H 1x Il .7J - ,2 ' k. 2,
R e p r o d u c e d with p e r m i s s io n of t h e co p y rig h t o w n e r. F u r th e r r e p ro d u c tio n prohibited w ith o u t p e rm is s io n .
At
2
(B.3b)
CCzKz + a z
M z
M z [ g / , +f f .
+ At
D. n+1
A?
D
£l + ^
Ar
2
A/
= G, " 2
Gx
Ar
+-■■
CE
//A z
2
I"
v f-H-*
-E
I"
v ''.z- 2’*-*
'I
(B .4a)
Ar
8
(i —
i. . ) £
r
My
a 7K 7+ o r
At
H,
//
■+
Ar
(B .4b)
£b
+
Ar
enK.
2
a K + a . -Qx ’H H
Ar
2
'
'..l)£ . L - *
£„A", + or. A', + cr, a
At
At
At
i + ~^A x
'
(B .5a)
At
M z
M
_
n+2
2
-J'k i
M
l
z+
<£ . 'L . , . - £ > r w
^z
At
enK
"O z
|
a z, Kz + a z,
Ar
fo _ M
+
Ar
n+^
2
Ar
■Gy M H ~ M
M l + M lM l
Ar
2z "+T
o
z
2
. M
Ar
2
=
(B .5b)
£g__<£_
-G,
z+ ^ z
^
*—
,Ar—
2
At
(B .6a)
At
fi A x
0
-
4
. W * ,
R e p r o d u c e d with p e r m is s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e r m is s io n .
w here C oefficients a and p are defined in (4.9) - (4.16).
R e p r o d u c e d with p e r m i s s io n of t h e cop y rig h t o w n e r. F u r th e r r e p r o d u c tio n prohibited w ith o u t p e rm is s io n .
Документ
Категория
Без категории
Просмотров
0
Размер файла
2 843 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа