# 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 .

1/--страниц