Full report

Original page 1

NASA/TM 2002 211663 A Two-Dimensional Linear Bicharacteristic Scheme for Electromagnetics John H. Beggs Langley Research Center, Hautpton, Virginia May 2002

Original page 2

The NASA STI Program Office... in Profile Since its founding, NASA has been dedicated to the advancement of aeronautics and space science. The NASA Scientific and Technical Information (STI) Program Office plays a key part in helping NASA maintain this important role. The NASA STI Program Office is operated by Langley Research Center, the lead center for NASA's scientific and technical information. The NASA STI Program Office provides access to the NASA STI Database, the largest collection of aeronautical and space science STI in the world. The Program Office is also NASA's institutional mechanism for disseminating the results of its research and development activities. These results are published by NASA in the NASA STI Report Series, which includes the following report types: TECHNICAL PUBLICATION. Reports of completed research or a major significant phase of research that present the results of NASA programs and include extensive data or theoretical analysis. Includes compilations of significant scientific and technical data and information dccmcd to be of continuing reference value. NASA counterpart of peer-reviewed formal professional papers, but having less stringent limitations on manuscript length and extent of graphic presentations. TECHNICAL MEMORANDUM. Scientific and technical findings that arc preliminary or of specialized interest, e.g., quick release reports, working papers, and bibliographies that contain minimal annotation. Does not contain extensive analysis. CONTRACTOR REPORT. Scientific and technical findings by NASA-sponsored contractors and grantees. • CONFERENCE PUBLICATION. Collected papers fi'om scientific and technical conferences, symposia, seminars, or other meetings sponsored or co-sponsored by NASA. • SPECIAL PUBLICATION. Scientific, technical, or historical information from NASA programs, projects, and missions, often concerned with subjects having substantial public interest. • TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and technical material pertinent to NASA's mission. Specialized services that complement the STI Program Office's diverse offerings include creating custom thesauri, building customized databases, organizing and publishing research results.., even providing videos. For more information about the NASA STI Program Office, see the following: • Access the NASA STI Program Home Page at http://www.sti.nasa.gov • E-mail your question via the Internet to hclp_sti.nasa.gov • Fax your question to the NASA STI Help Desk at (301) 621 0134 • Phone the NASA STI Help Desk at (301) 621 0390 • Write to: NASA STI Help Desk NASA Center for AeroSpace Information 7121 Standard Drive Hanover, MD 21076 1320

Original page 3

NASA/TM 2002 211663 A Two-Dimensional Linear Bicharacteristic Scheme for Electromagnetics John H. Beggs Langley Research Center, Hautpton, Virginia National Aeronautics and Space Administration Langley Research Center Hampton, Virginia 23681-2199 May 2002

Original page 4

Available from: NASA Center for AeroSpace Information (CASI) 7121 Standard Drive Hanover, MD 21076 1320 (301) 621 0390 National Technical Information Service (NTIS) 5285 Port Royal Road Springfield, VA 22161 2171 (703) 605 6000

Original page 5

Abstract The upwind leapJ?og or Linear Bicharacteristic Scheme (LBS) has previously been extended to treat lossy dielectric and lossy magnetic materials. This report extends the Linear Bicharacteristic Scheme Jot computational electromagnetics to the two-dimensional case, which includes treatment of lossy dielectric and magnetic materials and perJect electrical conductors. This is accomplished by implementing the LBS JOt homogeneous lossy dielectric and magnetic media and Jor perJect electrical conductors. Heterogeneous media are modeled by applying surJace boundary tions at dielectric material boundaries are required. conditions, and no special extrapolations or interpola- The PerJectly Matched Layer (PML) outer boundary concept is also developed Jot this scheme. Results are presented Jot two-dimensional model problems on uniJorm grids, and the FDTD algorithm is chosen as a convenient reJerence algorithm Jot comparison. The results demonstrate that the explicit LBS is a dissipation-J?ee, second-order accurate algorithm which uses an upwind computational stencil rather than a central difference stencil, and yet it has approximately onethird the phase velocity error. Computational requirements are also discussed. 1 Introduction Numerical solutions of the Euler equations in Computational Fluid Dynamics (CFD) have illustrated the importance of treating a hyperbolic system of partial differential equations with the theory of characteristics and in an upwind manner (as opposed to symmetrically in space). These two features provide the motivation to use the Linear Bicharacteristic Scheme (LBS), also called the upwind leapfrog (UL) method, for the construction of many practical wave propagation algorithms. The upwind leapfrog (UL) method is based upon the Method of Characteristics, which is a widely used numerical solution concept in CFD [1] [16]. In a hyperbolic system, the solutions (i.e. waves) propagate in preferred directions called characteristics. A characteristic can be defined as a propagation path along which a physical disturbance is propagated [17]. The relevance to Maxwell's equations is intuitively obvious because electromagnetic waves have preferred directions of propagation and finite propagation speeds. Characteristic-based methods have also been successfully implemented and demonstrated primarily electromagnetic problems [18] [31 ]. for free space and perfect electrical conductor (PEC) This report extends the LBS to the two-dimensional case to model both homogeneous and heterogeneous lossy dielectric and magnetic materials and perfect electrical conductors (PECs). The LBS was originally developed to improve unsteady solutions in computational acoustics and aeroacoustics [32]-[38]. It is a classical leapfrog algorithm, but it uses a one-sided (or upwind) stencil for the spatial derivatives, which follows the wave characteristic more closely when compared with a classical leapfrog method. This approach preserves the time-reversibility of the leapfrog algorithm, which results in no dissipation, and it permits more flexibility by the ability to adopt a characteristic based method. Clustering the stencil around the characteristic enables high accuracy to be achieved with a low operation count in a fully discrete way [33]. The use of characteristic variables allows the LBS to treat the outer computational boundaries naturally using the exact compatibility equations. The LBS treats the outer boundary condition naturally without nonreflecting approximations. The interior point algorithm predicts the outgoing characteristic variables at the domain boundaries. For multidimensional applications, in principle, through knowledge of the wave propagation angle, the local coordinates can be rotated to align with the characteristics, at which the boundary condition

Original page 6

becomesalmostexact.Therefore,noextraneousboundaryconditionis required.Inthecaseswherethis coordinatetransformationisnotimplemented,thecharacteristic-basedalgorithmprovidesonlyanapproximationattheoutergridboundaries.However,thePerfectlyMatchedLayer(PML)outerboundaryconcept canbeappliedtothisscheme,whichisdiscussedlaterinthisreport.TheLBSalsooffersanaturaltreatment ofdielectricinterfaces,withoutanyextrapolationorinterpolationoffieldsormaterialpropertiesnearmaterialdiscontinuities.Exactboundaryconditionsonthetangentialfieldcomponentsaredirectlyenforcedat materialinterfaces.TheLBSoffersacentralstorageapproachwithlowerdispersionthantheYeealgorithm [39].It haspreviouslybeenappliedtotwoandthree-dimensionalfree-spaceelectromagneticpropagation andscatteringproblems[34],[37],andit wasrecentlyextendedtotreatlossydielectricandmagneticmaterialsfortheone-dimensionalcase[40]. TheobjectiveofthisreportistopresenttheextensionoftheLBStothetwo-dimensionalcase,which includeslossydielectricandmagneticmaterials.Resultsarepresentedforseveraltwo-dimensionalmodel problems,andtheFDTDalgorithmis chosenasaconvenientreferenceforcomparison.Theprinciplesto extendthisproceduretothethree-dimensionalcasearestraightforward.Sections3 and4 presenttheLBS implementationfortheTMandTEpolarizations,respectively.Section5 outlinesthedielectricmaterial surfaceboundaryconditionandSection6discussestheouterradiationandPMLboundaryconditions.Section8reviewstheFourieranalysisandcomputationalrequirements.Finally,Section9presentsresultsfor two-dimensionalmodelproblemsandSection10providesconcludingremarks. 2 Abbreviation List The following table provides a list of abbreviations Abbreviation CFD FDTD LBS PEC PML TE TM UL 2D 3 TM Polarization and acronyms used throughout this report. Description Computational Fluid Dynamics Finite Difference Time Domain Linear Bicharacteristic Scheme Perfect Electrical Conductor Perfectly Matched Layer Transverse Electric Transverse Magnetic Upwind Leapfrog Two-dimensional Maxwell's equations for linear, homogeneous and lossy media in the two-dimensional TM case (taking O/Oz 0) are OEzot 7_1(OH v OHxoy E ) (1) OH.Ot #1( OEzO---y crH_) (2) Ot # ---_z crHv (3)

Original page 7

where_ and_* aretheelectricandmagneticconductivities,respectively.Usingtheelectricdisplacement D cE and making the substitution c 1/v/-fi-g gives ODz ( OH x OHy _ cr O--7--+ \ -y -j ] + - D zc 0 (4) 1 OH, ODz (r* c-7 0---i-+ -aT + Tj Hx 0 (S) 10H: ODz or* c20t Ox + --SHv 0 (6) The procedure for the LBS is to transform the dependent variables Dz, Hx and H:_ to characteristic variables. The algorithm developed here is the simplest leapfrog scheme described by Iserles [41] combined with upwind bias, or simply, the Linear Bicharacteristic Scheme (LBS). To transform (4) (6) into characteristic form, we multiply (5) and (6) by c and then add and subtract from (4) to give --O (D z c1H )+c O (D z c1H Y) + - * + OHx (7) 1 H Z or* OHx --Ot 1H c 8 (Dz + c y) + + (8) 0-7 D+ H +c-5_yo(1)Dz+ H. +ZDz+---H.e#c OHvOx 0 (9) • (0 Dz _ H )0(+c Dz 1),Hx +-Dz+Cr * H. 0H_ 0 (10) Ot c # c Ox Note that these equations are almost identical to the equations for the one-dimensional case [40], except for the addition of the cross-derivative magnetic field terms. The characteristic variables are defined as P Q S to represent the ix and --y propagating solutions, rewritten as OP OP 1 (_ ) 1 Dz -H v (11) C 1 Dz + - H v (12) C D + 1H (13) C 1 Dz - Hz (14) C respectively. Using these definitions, (7) (10) can be 1 ( __) OHz 0-7 c +2 P+g + Q + --y 0 (16) Ot _ + R + _ S Ox 0 (17) Ot c-y + _ R + g + S Ox 0 (18)

Original page 8

It is convenienttodefineandstorethefollowingcoefficientsbeforetime-steppingbegins a -b Equations (15) (18) can be rewritten more concisely as OP OP a__p 0---7-+ c-_-7-x+ 2 oo oo bp Ot c-z + 2 Cr (9"* +-- (19) c # cr (9* (20) c # b OHx + Q + ----y 0 (21) _ OHx + 2Q + O----y- 0 (22) oR oR a_R bs OH:. 0- + c57 + 2 + 2 ox 0 (23) OS OS b a OH:, Ot C-_y + _lg + -S Ox 0 (24) To develop the discretized algorithm for a two-dimensional system, the stencils of Figures 1 and 2 are proposed for the LBS. We discretize time and space as t nAt, x lax, y jAy. To solve the wave (a) (i,j,n) ? (i+i/2, j,n+l) t,lq (i-i/2,j ,n-l) (b) (i-i/2,j,n+l) ? (i,j,n) ! /a, s Y,3 (i+i/2,j,n-l)O x,i Figure 1: Two-dimensional upwind leapfrog computational stencils for right-going (a) and left-going (b) x propagating characteristics. propagation problem without introducing dissipation, it is necessary that the stencil have central symmetry so the scheme employed is reversible in time [33]. The stencil in Figure 1 a is used for +x propagating waves and the stencil in Figure lb is used for x propagating waves. The upwind bias nature of these stencils is clearly evident. Figures 2a and 2b show the stencils for +y propagating waves, respectively. References [32], [33], [36], [37], [38] clearly show that the LBS is second-order accurate. Note that the third and fourth terms in (21) (24) represent the electric and magnetic loss (or source) terms. A key element in developing an accurate LBS scheme is proper treatment of these source terms. The method used here indexes the self source term in (21) (i.e. P) at time level n + 1 and it indexes the coupled source term Q at time level n. This avoids a matrix solution at each grid point, and the formulation

Original page 9

9 (i,j+i/2,n+l) (i,j ,n) '22 < R,S (a) O (i, j-q/2,n-q) (i,j ,n) (i,j-i/2,n+l) Q I y,J x,i // (b) (i, j+i/2,n-l) Figure 2: Two-dimensional upwind leapfrog computational stencils for right-going (a) and left-going (b) y propagating characteristics. easily limits to the perfect conductor condition as G _ co. An identical application is made for equations (22) (24). Using the stencils shown in Figures 1 and 2 and the source term indexing scheme described above, the resulting finite difference equations for (21) (24) are i+l/2,j + pn pn pn \ pn÷l P_Zll/2,j) . ( i 1/2,j.. pni 1/2,j]1 "_ @ C ( i+l/2,j i 1/2,.,_ a pn+l 2At b Q.n 1 (H_!(i,j + 1/2) H_!(i,j 2At b pr_ 1 (Hr_(i.j + 1/2) H_(i,j -2 i 1/ 2,.j @ -'y , (/U++Il/,J+l/2 /_z,j+l/2). @ (i,j12_n 1/2 /ni,j 1 1/2J 2At b S.r 1 -2 i,j+1/2 _ (Hy*(i + 1/2,j) H:;Z(i i,j 1/2 @ ( i,j+1/2 Sn 1 2At " b 1/2 1 (H;(i + 1/2,5) \ -A-Tx ] + 2 i+l/2,j + 1/2)) 0 (25) c " ' + _'i 1/2,j + 1/2)) 0 (26) Ri,j+I/2_Z[i,j 1/2 a Rrt+ 1 + C < Ay ] @ 2 i,.j+1/2 @ 1/2,j)) 0 (27) _ sn sn C " -7, @ 2 L'i,j 1/2 @ 1/2,j)) 0 (28) where P_j denotes the value for P at grid point (i, j) and time level n. Note that the differences are taken with respect to the cell center, i.e. the coordinate (i, j) is located at the center of the cell. Since we know that H - c (R S)/2 and H: - c (Q P) /2, these equations can be rearranged in the form (l+a/kt) Prt+li+l/2,.j prtli 1/2,j @ (1 2lZx) (P'rzi+l/2,j Pftl/2,.j) b/ktQi+l/2,jn

Original page 10

j_'rti,j 1/2) @ b'y (Si,.j+I/2 S[:j . 1/2) (29) b'y( Rn i,.j+1/2 (1 @ a/t) /)'z+1 Oi+l/2,jn 1 (1 "{ 1/2,j j_'rti,j 1/2) @ b'y (i,.j+1/2 [:j . 1/2) (30) b'y( Rn i,.j+1/2 (1 + aAt)/)n+l r+i,j 11/2 @ (1 • vi,j+l/2 Pft 1/2,j) + /"x (Qi+i/2,j Q_I/2,.j) (31) x (P,[I/2,j 2b'x) ( Qni+1/2j on i 1/2,j) bAtpr_l/2,j 2uu) ( Rrt i,j+1/2 j_ni,j 1/2) bat S_ij+l/2 (1 + aAt) S n÷l sn li,j+l/2 (1 2/@)(sni,j+l/2 S_j 1/2) bAtR_Z.i,. 1/2 i,j 1/2 Pft 1/2,j) + /"x (Qi+i/2,j Q_I/2,.j) (32) "x (P[I/2,.j where Ux cAt/Ax and uv cAt/Ay are the x and y Courant numbers. We now rewrite equations (29)-(32) as pn+l i+1/2,j Qn+ I i 1/2,j Rn+l i,j+1/2 i n+ l ,j 1/2 where R[ _ R) are the residuals defined by /[t- P;Z 1}2,j @ ( 1 2l/x)(Pf;1/2,j /2y (:,.j+1/2 [:zj 1/2) /['/ (1 @ a At) (33) (34) /'/ (1 @ a At) R_aZ/(1 + a At) (35) /E_/ (1 @ a At) (36) P(_ 1/2,.j) bAtQri_l/2,j + /2y (i,j+l/2 ,:.j 1/2) (37) R_2t Qi+l/2,jn1 (1 2Ux)( Qni+l/2,j Qnl/2,j) b/ktPf_l/2,j /2y (:t.j+l/2 [:zj 1/2) + /2y (Si,j+l/2 ,:.j 1/2) (38) /l_t n li,j 1/2 @ (1 2/@)(]2_ni,j+1/2 J_j 1/2) bz2xtS::j+l/2,. /2. (Pf;1/2,j P(t 1/2,.j) @ l/x (Qi+l/2,j Q 1/2,.j) (39) R S_-j+I/2r_1 (1 2v',V) (i,j+1/2Sn Sir,j 1/2) b At/', 1/2 /2. (Pf;1/2,j P(t 1/2,.j) . + P'z (Oi+l/2,j O_1/2,.j) (40) Equations (33) (36) are the update equations for the 2D TM LBS scheme at cell (i,j) which can contain as cr that P_+_ lossy dielectric and magnetic materials. Note that Qn+l Rn+l sn+l i 1/2,j' i,j+1/2' and i,j 1/2 0 as required. 4 TE Polarization oc, then we have the PEC condition i+1/2,j, Maxwell's equations for linear, homogeneous and lossy media in the two-dimensional TE case (taking O/Oz 0) are OE. 1 (OHz ) Ot 7 --_y GExj (41)

Original page 11

OEOtv 7_ ( Ot # \ _ OH_Ox E:u) (42) _ ] a---Hz# (43) Using the electric displacement D -- eE and making the substitution c - 1/v/-fi-g gives ODx Ot ODy Ot + _ + 7 Dy 0 (45) OH + adz - 0 (44) Ot e OHz a 10Hz OD. OD v (7" c20t Oy + _ + _-Hz 0 (46) The procedure for the LBS is to transform the dependent variables D., D v and Hz to characteristic variables. To transform (44) (46) into characteristic form, we multiply (46) by c and then add and subtract from (44) and (45) to give 1H a O("y+• 71H z) +c O(Dy+7)• z cODx * ot _ O(D, O(D, Ot v Ox O(Dx :1 H z) 1H O (Dx 70 + c Ot Oy O (Dx + c: H z) 0 (Dx + C Ot Oy The characteristic variables are defined as P Q R S +TDy + 0 (47) _ H 0,. < + -Dy + c Hz 0 (48) e Oy pc . a cODv + -Dx a Hz 0 (49) c Ox pc + aDx + c + a*H 0 (50) c #c 1 Dy + - Hz (51) C Dv 1 Hz (52) C 1 D - Hz (53) C 1 Dx + - Hz (54) C to represent the 4-x and 4-y right and left propagating solutions, respectively. Using these definitions, (47) (50) can be rewritten as OP OP (55) at + c-b7 + : + P+: Q c--_v o P + _) ODx (56) ot _ + 2 + 0+_77 o OR OR + + R + (57) o-T+ o--- : ot _N+: R+: -5 a*# ) S cODy-_f 0 + _--_f) S+--_fc ODy 0 (ss)

Original page 12

Using the a and b coefficients defined in (19) and (20) we can rewrite equations (55)-(58) more simply as OP OP ap b ODx 0 (59) 0--/-+%-7+2 +3 0 Co---7 OQ OO b__p a ODx 0 (60) ot -g-21+2 + 0 + _ o--7 OR OR a b cODv O-T + c-y + R + S _ 0 (61) ON OS b a ODy ot _ + JR + Ts + c 0 (62) To develop the discretized algorithm for a two-dimensional TE system, we use the same stencils as for the TM case, which are shown in Figures 1 and 2. We also employ the same indexing scheme for the self and coupled source terms in (59) (62) and we also use a central difference approximation at the appropriate half-integer indexed cell to evaluate the cross derivative terms. To derive the finite difference equations for (59) (62) we use the same stencils shown in Figures 1 and 2. Since we also know that Dx (R + S)/2 and Dy (P + Q)/2, the TE finite difference equations are (1 @a/kf;)pn+li+l/2,j pni 1/2,jl @ (1 2 Z"x) (Dni+l/2,j P? 1/2,j) b/ktQr_l/2,j @ Rn i,j 1/2) + % Sn i j 1/2) (63) (1 + aAt) on+l Qi+l/2,.jn1 (1 2..) ( Qrt i+1/2.j (T i 1/2,j) bAtpr_l/2,j "i 1/2,j R'Z i,j 1/2) % Sn ij 1/2) (64) (1 @ az2Xt)t.11/2 j_n i,j 1 1/2 @ (1 2 Uy) ( I_n i,j+l/2 R'ni,j 1/2) b/kt X_Ij+I/2 @ (1 + aAt) S n+l snli,j+l/2 (1 2//y) (sni,j+l/2 SPj 1/2) bAtR_zj,. 1/2 i,j 1/2 We now rewrite equations (63)-(66) as p/n+l +i/2,j Qn+ l i 1/2,j iRn+ l i,j+l/2 Sin+ l ,j 1/2 where R[ _ R) are the residuals defined by R[/(1 + a At) (67) n;'/(a + _ At) (68) R/(1 + a At) (69) R2/(1 + a At) (70) R _ pn i 1/2j+(1i 2t,) ( pn i+1/2.j P_i 1/2,j) bat Oi+l/2,j + //Y (i,j+1/2j_n R"i,j 1/2) @ l@ (ij+1/2Sn sn i,j 1/2) (71) R Qi+l/2,j, 1 (1 2u) (Qi+l/2,j' Q_' 1/2,_). bAtP[' 1/2,j .

Original page 13

(72) /r_rtli,j 1/2 + (1 2/@) ( i_r_i,j+l/2 lc_t,j,. 1/2) b_tSi,j÷l/2@n (73) Si,j+l/2n1 (1 2 u_) ( Sni j+1/2 Sni,j 1/2) b AtR_j 1/2 (74) Equations (67) (70) are the update equations for the 2D TE LBS scheme at cell (i,j) which can contain that/+1 lossy dielectric and magnetic materials. Note that as cr _ oc, then we have the PEC condition i+1/2,j, Qr+li1/2,j' Rr_,+ll/2, and Sri+11/2 0 as required. Note that the update equations are identical to the TM case, the differences being in the definition of the characteristic variables and in evaluation of the cross derivative terms. 5 Heterogeneous Materials One of the difficulties with the conventional FDTD algorithm is the error in treatment of material discontinuities. Recent research efforts have attempted to reduce this error source by suitable averaging of material properties across the interface or by interpolation or extrapolation of the electromagnetic fields near these material boundaries [42], [43]. The advantage of the LBS is that the characteristic based nature of the algorithm leads to a very natural treatment of dieletric interfaces. Since the LBS works with characteristic variables, the slope of characteristic curves in each material will be different, and the physical boundary conditions permit an elegant and efficient implementation of a dielectric interface boundary condition. This numerical boundary condition implements the physics exactly, with no averaging, interpolation or extrapolation required. To implement the dielectric material interface boundary condition, consider a portion of a two-dimensional grid shown in Figure 3, which contains material discontinuities in both the x and y directions. We can see that the characteristic variables P and Q are co-located at the center of the cell edges along the y axis. Similarly, variables R and S are co-located at the center of the cell edges along the w axis. Thus, the LBS has a staggered storage scheme, similar to the conventional FDTD method. Spatial derivatives are taken with respect to the cell center, which is where the cell coordinates (i, j) are defined. The characteristic variables at each grid point (i, j) on the interface are split into two components each: PI,j, QI,j, P2,j and Q2,j for interfaces perpendicular to the x axis and/},1, Ri,2, S,z,2 and Si,2 for interfaces perpendicular to the y axis. The terms PI,j, QI,j,/i,1 and Si,1 exist just to the left and bottom of the material interface, respectively, as shown in Figure 3. The remaining terms/3,j, Q2,j,/i,2 and S,i,2 exist just to the right and top of the material interface. Note that the i and j subscripts have been omitted from the dielectric boundary split field components in Figure 3 for clarity. For material 1, equation (33) is used to predict the value for p+l1,j at the boundary and for material 2, equation (34) is used to predict the value for Qr,+l Similarly, equation (35) is used to predict the value of/.+1 and (36) predicts the value for 59 +1 2,j • The procedure for the TE polarization is identical. , i,2 • For example, in the TE case, the characteristic variable P uses field components D:_ and Hz, which both are tangential to material interfaces that are perpendicular

Original page 14

y (i,j+l) P2'Q2 R S X R2,S 2 Figure 3" Section of a two-dimensional computational grid for the LBS showing characteristic variables, dielectric interfaces and corresponding field components and characteristic variables used for the surface boundary condition. to the z axis. To complete the implementation, the Qr_+ll,._ and p,.+s,,terms must be updated. These terms are updated by enforcing the physical boundary conditions on the electromagnetic field at the material boundary. We can then solve for Qr,+11,.and p,d,,+l in terms of the "known" characteristic variables/+11,. and "2,.r+l"To develop this procedure, the electromagnetic boundary conditions on the tangential field components are given by Ezl,j Ez2,. H;I1,j -- H:I2,j Dzl.j Dz2,f _ ' " (75) 1 2 (76) For the right-going wave, substituting (75) and (76) into (11) gives pr+l rv+l lHr+l 1,j zl,j @ Cl :ql,j (77) ¢1 {p,r+l Qn+l c2 {p,r_+l Qn+l_ (78) 2¢2 t, 2,._ + Similarly, substituting (75) and (76) into (12) yields 20 ] +t, 2,. 2,j ] @rz+l /3_z+1 ]--H rt+l (79) 2,j z2,.f c2 v2,j ¢2 {pr+l Qn+l'_ cl {pr_+l Qn+l'_ (80) 2el \ 1,j @ pr_+l and nr_+l 1,j ] _ \ 1,j 1,j ] Since 1,j "2 j are determined at boundary point (i, j) from the usual update equations (we treat them ./7 + 1 g/r_+l as "known" variables), it is necessary to express 2,._ and ,1,. in terms of these variables. Rearranging 10

Original page 15

(78)and(80)gives p,n÷l2,j T1 1,j +F1 2,j pn+l @n+l (81) Qn+l1,j r 2 pn+l1,j + T2 Qn+l2,j (82) where F12 and T1,2 are reflection and transmission coefficients given by F 1 ( C2_2 c1_1) (83) \c2_2 -Cl_l/ rl C2_2 @ Cl_l 2e2cl (84) F2 (c1_1 c2_2 (85) \C2_2 TCl_I/ Ti c2_2 @ Cl_l 2elc2 (86) From (81), it is clear that a right-going wave in material 2 is a sum of a transmitted portion of a right-going wave in material 1 plus a reflected portion of a left-going wave in material 2. A similar argument can be made for the left-going wave in material 1. In fact, the reflection coefficients F1,2 can be shown to be identical to the classical Fresnel reflection coefficients. The transmission coefficients also have the same form as the Fresnel transmission coefficients. Special care needs to be taken when the LBS calculates the solution at grid points near a material discontinuity. For example, for the x interface at grid point (i, j) as in Figure 3, care must be exercised to update the solutions at grid points (i 1,j) and (i + 1,j). At grid point (i 1,j), the term (l,j in (38) becomes Q[t,.i. At grid point (i,j), the terms P?.,.and Qi,j in (37) and (38) become Fll,j and Q_,j, respectively. At grid point (i + 1, j), the term P[_ 1,j in (37) becomes P.,j. Rearranging equations (29) and (30) for grid point i we have (1 + al At)pit,;1 pni 1,j +1 (1 2vl)(P',.j P: 1,j) bl Z tQl,j'rt (87) ,. Q. ,j Q ,j) b AtP, .j (88) where _'l clAt/Ax, and 2 c2At/Ax. The terms al, a2, bl, b2 refer to the a and b coefficients in (19) and (20) for materials 1 and 2, respectively. These equations are now easily solved for PI+1,. and Qn+12,jand then (81) and (82) are applied to obtain/_2+12,j and Qn+12,.i• A similar analysis can be made for the boundary perpendicular to the y axis involving the R and S field components. 6 Outer Boundary Condition The outer radiation boundary condition is used to terminate the computational lattice and permit outgoing waves to pass unreflected through the lattice boundaries [44]. The FDTD algorithm uses a spatial central difference operator where it uses field values from neighboring cells to update solution variables. Thus it cannot be used at the terminating faces of the problem domain. For example, the solution for a wave propagating left to right will eventually require a grid point outside the domain. To terminate the computational lattice, an additional equation (boundary condition) is needed to solve the system and this introduces 11

Original page 16

informationintothesolutionthatis notrequiredbyMaxell'sequations.ThePMLboundarycondition[45] hasrecentlybeenintroduced,whichhasgreatlyincreasedtheaccuracyofFDTDsimulations.However,the PMLcomeswithamoderateincreasein complexityforanFDTDcodedueto additionalvariablestorage andupdateequations. Onthecontrary,theLBSrequiresnoextraneousboundarycondition,andit includesthePMLboundary conditionwithnoextrarequiredstorageorupdateequations.ForthepresentLBSimplementation,likethe MethodofCharacteristics[31], theinteriorpointalgorithmcalculatestheleft-goingcharacteristicattheleft boundary(i.e. i - 0) and the right-going characteristic at the right boundary (i.e. i - imax). Thus for the LBS, at grid point i - 0, equation (34) calculates Q(0, j) and the incoming right-going characteristic, P(0, j), is specified as a boundary condition. This same analysis applies at the right boundary where (33) calculates P(imax,j) and the incoming left-going characteristic, Q(imax,j), is specified as a boundary condition. Shang [20] has noted for characteristic based multidimensional and nonuniform grid problems, in principle, the local coordinate system can be rotated to align with the characteristics, and the compatibility equations provide an exact boundary condition. This transformation has not been implemented in the present work, and will likely be the subject of future studies. A simple, yet effective approximation for multidimensional characteristic based approaches is to set the incoming flux or characteristic variables at the outer boundaries to zero and let the interior point algorithm predict the outgoing variables. When the wave motion is aligned with a coordinate axis, this boundary condition is exact. But this approximation may not be necessary since the LBS automatically includes the PML boundary condition without additional storage or update equations. The linear bicharacteristic form of Maxwell's equations for the 2D TM polarization in free space are OP OP OH:_ (89) 0-7 + c-07x+ Oy o 0(2 OQ OHz (90) Ot c-07 + Oy o OR OR 0-7 + Tv Os Os c Ot Oy In the frequency domain using complex coordinates, P OH: (91) O. o OH:_ 0 (92) Ox we have OHz jw P + C-x + 02 0 (93) O jc_Q c+ R jw R + c_O S jco S c-O_ OHz _ 0 (94) OH_ 02 0 (95) OH_ 02 0 (96) To show how the LBS automatically includes the PML boundary condition, we derive the appropriate update equations using the complex coordinate transformation approach proposed by Chew and Weedon [46]. 12

Original page 17

Specifically,weuse O 02 0 O_ Sx -sy Substitutingtheseinto(93)-(96)gives jwP+ (7:p + 1 O (97) sx 8x 1 0 (98) % Oy (7 x 1 + . (99) 2wco 1 + .(7Y (100) 36060 OP OBx co c-0-x ÷ 0y 0 (101) jc_O+ c + OBxOy 0 (102) o Q jwR+ :vR + OR OBy co jwS+ (7:qS co c Oy Ox 0 (103) OS OBy c 0 (104) Oy Ox where Bx (Sx/%) Hz and B:q (s:q/sx) H:. In typical fashion with a PML FDTD implementation, we let (7z (7y (7, then we have that Bx H,/_v OP OP H u and (101)-(104) become OHx 0----t-+ c -0-_x+ _--P + 0 (105) co Oy OQ OQ (7 OHx Ot C-x ÷ -Qeo ÷ Oy 0 (106) OR OR (7 R OHy O + C_y + co Ox" 0 (107) OS OS ---S OHy Ot c--ffy + co Ox" 0 (108) Furthermore, if we let c -- co, # -- #o and (7*/# 0 - (7/% as required by the PML boundary condition, then the normal LBS update equations given by (21)-(24) can easily be shown to be identical to the LBS PML update equations (105)-(108). This analysis shows how the LBS inherently incorporates the PML boundary condition within the standard update equations. The PML conductivity (7 is still specified using the conventional profiles: linear, quadratic or geometric [43]. 7 Computational Requirements It is instructive to examine the computational requirements of the LBS and the FDTD method. We can use this analysis to determine if the LBS can provide equivalent or better accuracy than FDTD for the same amount of computational resources. Let us assume a 2D grid with N x N cells. The FDTD method requires SF 12N 2 + 16N + 4 (109) 13

Original page 18

total bytes to store the field component arrays, and the LBS requires Sc -- 32N 2 +24N (110) total bytes. Note that this storage calculation does not account for any extra terms such as arrays for boundary conditions, far-field transformations, etc. We can define a storage ratio S_. between the LBS and FDTD as SL S,. SF 32N 2 + 24N (111) 12N 2 + 16N + 4 If the LBS is more accurate than FDTD, we should be able to increase the cell size by a certain factor and still maintain the same accuracy as FDTD. Increasing the cell size decreases the total number of cells required in the grid. Thus, we define a grid reduction factor %, which can be used to determine the breakeven point in storage and accuracy. The grid size for the LBS will be reduced in each dimension by N, giving a new ratio 5,I. 32 (N/Nr) 2 + 24 (N/Nr.) 1 12N2 + 16N +4 Nr?S,. (112) The percentage reduction in grid storage ratio from the FDTD method is then given by (113) Pr 100 (Sr S_) 100 ( 1 _rr1 2Sr ) To determine the breakeven point, we solve Pr - 0 for Nr in terms of N to yield (114) Nr -- v/2N3N 2+4N+1(4N + 3) Taking the limit of the positive root as N _ oc gives % _ 1.63. Thus, the LBS must be at least 1.63 times more accurate than FDTD to achieve equivalent storage for the same accuracy. Factors above 1.63 means the LBS requires less storage than FDTD for the same accuracy. Figure 4 shows a plot of the breakeven 1.6- 1.4-- ] ._ol.2- _0.8 ._0.6 _0.4- 0.2 2 4 6 8 10 12 14 16 18 20 Number of grid cells Figure 4: Breakeven ratio versus number of grid cells. ratio versus the number of grid cells. 14

Original page 19

8 Fourier Analysis Various Fourier analyses of the two-dimensional therefore, only the important results and conclusions LBS have already been completed [36], [37], [38]; from these previous analyses will be reviewed in this report. Most of the information presented is summarized from [36]. The stability condition for the 2D LBS is L,,, L,:v< 1/2, where u_, t,:v are the Courant numbers L,_ - cAt/Am and % - cAt/Ay. Although this stability limit is more restrictive than the standard FDTD method, it is not particularly troublesome because many FDTD simulations use a Courant number of 1/2 for improved accuracy. The complete Fourier analysis will not be outlined here for the sake of brevity. Rather, we present an overview of the procedure followed by a discussion numerical results. The procedure for the Fourier analysis is straightforward. Start with the LBS free space update equations (21)-(24) with a - b - 0 and substitute a solution of the form p_.j _ Poe.(.r, io,_ .jo_) (115) into these expressions. After some algebra, we have the system of equations Tn+l _ V1Tn + V2 Tn 1 (116) which represents the three time-level LBS scheme with T z [przi,j, z,.}, -oz,.,, -,,.j • To complete the Fourier analysis, we make the substitution W +1 - T r_to give IT]n+lgg IV1V21[TIn140 • (117) where/4 is the 4 × 4 identity matrix. The stability matrix G is then given by G [V1_72114 0 (118) which is an 8 × 8 matrix. The stability analysis is completed by calculating the eigenvalues of the stability matrix G for various grid resolutions and grid propagation angles. To that end, we define 0x - 0: -o - 05 -where N is the grid resolution in cells/wavelength 0 cosc (119) 0 sinc (120) 27r/N (121) u0 (122) and c_ is the grid propagation angle. To simplify the analysis, we also set L, - t,, - t,:v. The dispersion relation can be obtained by solving the equation dot [cj¢ G] - 0 (123) for 05. In comparison, the dispersion relation for the FDTD method is sin 2 05 -- r,2 sin 2 (0,/2) ÷ u 2 sin 2 (0;_/2) (124) 15

Original page 20

Fortheone-dimensionalLBS[47],it wasshowntheLBShadlessnumericaldispersionthantheFDTD method.Extensivethree-parameterstudiesofnumericaldispersionforthe2DLBSwereperformedusing thegridresolution(N),Courantnumber(L,)andgridpropagationangle(c_)asparameters.Thesestudies revealedthattheoptimumCourantnumberis L,- 1/2 sincedispersionis minimizedforallpropagation angleswhencomparedtoFDTD. Fora Courantnumber_ - 0.4andpropagationangleof 45°, thenumericaldispersiondecreases smoothlywithincreasinggridresolutionasshownin Figure5.Fromthisfigure,weseethattheLBShas I00 I ,, a, ': i i i i",, i i i r i ",, i i i c_ C ', _"_' ', ', 1 ........ i........ __< .... ]......... [......... ',, ,,, ,".... ,....... ,, -r-I , , ,'. . , 2 4 I I FDTD -- LBS ...... i i i i i i i i i i i i i i i ', ', ', ', ', i......... ] ......... i ......... L......... i........ ,, ,, ,, ,, ,, , , , , , ) 0.1 ........:..................;..................:..............":-"4---:-- '.........:........ &o 0.01 I I I I 0 5 10 15 20 Grid Resolution I I I I I 25 30 35 40 45 50 Figure 5: Phase speed error versus grid resolution N for FDTD method and LBS with , - 0.4 and c - 45°. approximately 1/2 the phase error as FDTD. Generally, the dispersion error for the LBS grows as _ _ 0. When _, - 1/2, numerical dispersion is zero along the coordinate axes and is maximum at 45° as shown in Figure 6 for a grid resolution N -- 10 cells/A. When L, < 1/2, dispersion for the LBS remains substantially less than for FDTD as shown in Figure 7 for N -- 20. From Figures 6 and 7, it is clear that as the grid resolution is doubled, the numerical dispersion decreased by a factor of four; as expected for a second order method. Finally, as shown in Figure 8 for N -- 10 cells/A, numerical dispersion decreases linearly as _, _ 1/2; except for grid propagation angles along 45 ° vectors, where the LBS dispersion is very close to that of FDTD. For propagation along 45 ° vectors, LBS numerical dispersion is minimized around L, - 0.3 and then approaches the FDTD value for , -- 1/2 as shown in Figure 9 for N -- 20. To summarize, the optimal Courant number for the LBS is 1/2. This Courant number offers much lower dispersion for most all propagation angles except those near a 4 vector. For _, < 1/2, numerical dispersion decreases as both grid resolution and Courant number are increased. Typically, LBS dispersion is at least 1/2 that of FDTD, and can be much lower in many instances. 16

Original page 21

1s I '_ '", 1 "' 0 ...... 1 iy o/i 1.5 1 0.5 Figure 6: Phase speed error versus grid propagation N- 10. 0.4 0.3 f.-,x 0.2 0.I 0 0.I 0.2 0.3 I I _ 0.4 0.4 0.30.20.1 Figure 7: Phase speed error versus grid propagation N - 20. I ----_ ' ' h}_ ...... "' 0 0.5 1 1.5 angle c_ for FDTD method and LBS with r, 1/2 and .:. .... i _ I I 0 0.10.20.30.4 angle c for FDTD method and LBS with _, 0.4 and 17

Original page 22

2.5 2 @ @ m @ m 1.5 c_ 1 _q 0 _q q @ 0.5 0 0.I Courant Number Figure 8: Phase speed error versus Courant number L, for FDTD method and LBS with N 10 and oe 0 0.25 I _\ : \ \ : \ : \ \ : I I iFDTD -- : LBS ...... : : 0.2 ........................................................................................... @ v, ', ' i\ i i m 0.15 ............. i- ..... ............................................. c \ ,, ', . '," ,, \ ! \ ', ', ', ', ', : i i i ........... .......... , ', ,, , '," '," ! , 0.1 ...................... _-- ........ ......................... ,,............ ,........... ; ....... o i \i i / '' o\O 0.05 .................................. ........................ ,,................................. : \ i i i 0 / "i'" i i i O. 1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Courant Number Figure 9: Phase speed error versus Courant number t, for FDTD method and LBS with N 20 and oe 4b 18

Original page 23

9 Results To demonstrate the 2D LBS, we consider various canonical problems using the TM polarization. First, we inject an incoming plane wave on the outer boundaries using the LBS, and let the algorithm propagate the signal through the grid using a total field formulation. This is done by specifying the incoming characteristic variable (P, Q, R or S) on the appropriate outer boundary. For example, on the left x boundary, P is specified for all j coordinates at i - 1. We use a 71 x 71 free space grid, with a Ax - Ay - 1 cm, which has a time step of At - 16.67 ps and the incident wave is a Gaussian pulse with FWHM of 35 time steps (or 0.58 ns). We specify the incidence angle as 180 °, and the electric field after 160 time steps is shown in Figure 10. Similar results can be obtained with other incidence angles. It is clear that the LBS easily allows 1 0.8 O6 O4 02 0 7O 6O I0 30 i40 50 60 70 Figure 10: Propagating plane wave injected on outer grid boundaries at 18ff incidence. specification of incoming plane waves in its fundamental algorithm. Next we move on to radiation from a point source in free space. This problem demonstrates that the algorithm can easily treat spherical waves and it also tests the PML boundary condition. Two concurrent grids are used in this problem, each having a cell size of 1 mm. The first is a small test grid of size 101 x 101 cells with an additional 10 cell PML boundary condition. This grid is centered within a large 501 × 501 grid, and the point source is located at the center of both computational grids. The time step is 3.3 ps, and an electric field point source is located at the center of both grids and the total number of time steps is truncated at 512, to allow no reflection from the large grid outer boundaries to reach the field sampling points. The inner grid is terminated with PML for both FDTD and the LBS, and the large grid is terminated with a second-order Liao boundary condition for FDTD and a characteristic based boundary condition for the LBS. The electric field is sampled at the same two locations in both grids, which are located 30 cells in the +x direction from the point source and then +30 cells in the y direction in the smaller grid. The point source is located in the smaller grid at grid point (61, 61) and the two sample points are (61, 91) and (61, 31). Figure 11 shows the electric field at the upper sample point in the large grid for point source radiation in free space. Note the agreement is excellent, and there are no reflections from the outer boundary due to the Liao boundary condition. Similar results were observed at the lower sample point. Figure 12 shows the electric field at the upper sample point (61, 91) in the small test grid using the PML boundary 19

Original page 24

0 0035 l FDTD -- 0.003 ....................................................... LBS 0 0025 0.002 > 0 0015 _J (1) 0.001 t i]iiiiiiiiiiiliiiiiiiiiiiiiiiiiiiiiiiiiliiiiiiiiiiiiiiiiiiii 0 0005 O 4 0 .o O f ; (1) -0 0005 -0.001 ............ . ............... -0 0015 i i -0.002 200 400 Time steps :............... ,,.............. _ ............. i i i 600 800 I000 1200 Figure 11" Electric field versus time sampled at upper grid point in the large grid. condition. Note again the agreement is excellent. Furthermore, we computed the global error in the small test grid with the expression GE - (Elar._e(i,j) Esmc_tt(<.j)) 2 (125) i,j using the difference between the electric fields in the large and small grids. Figure 13 shows this global error using the PML boundary condition for both methods and we see that the PML works very well. The error for the LBS is in the -80 to -100 dB range, which is excellent. Figure 14 shows the time-domain results for the LBS with and without the PML boundary condition. Note the reflections from the outer boundary are clearly visible for the no PML case. 10 Conclusions This report has extended the Linear Bicharacteristic Scheme for computational electromagnetics to the two-dimensional case. Treatment of lossy dielectric and magnetic materials was discussed, and implementation of the PML boundary condition was outlined. It was demonstrated that the LBS has several distinct advantages over conventional FDTD algorithms. First, the LBS is a second-order accurate algorithm which is about 2-3 times as economical. The LBS can also be made to have zero dispersion error in certain instances. Second, the LBS provides a more natural and flexible way to implement surface boundary conditions and outer radiation boundary conditions by using characteristics and an upwind bias technique popular in fluid dynamics. Third, the LBS can provide more flexibility to implement subgridding algorithms because of the compact nature of the computational stencil. A dielectric surface boundary condition was also implemented and results were provided for two-dimensional free space radiation problems. Due to project and time limitations, validation for lossy dielectric materials and heterogeneous materials was not explored 2O

Original page 25

0 0035 0.003 ................................................................ 0 0025 0.002 > 0 0015 (11 0.001 i I FDTD -- LBS ...... i]iiiiiiiiii]iiiiiiiiiiiiiiiiiiiiiiiii[iiiiiiiiiiiiiiiiiiiii 0 0005 o 4 0 J r.) , (11 -0 0005 ........... <............... M i : : : :............... ,.............. _ ............. i i i -0.001 ............ ............... !............... !.............. i ............. i i -0 0015 ............. ............... i i -0.002 200 400 Time steps i i i i i i ',............... :- .............. ; ............. i i i 600 800 1000 1200 Figure 12: Electric field versus time sampled at grid point (61, 91) for point source at grid point (61, 61) in the small grid. I I 0.0001 ........ 7--kX .................................................... FDTD / _ LBS - ..... le-06 __.J ............... _ ............................................... :.............. / le-08 i -;,..........................i.......................................................... g4 o g4 -"........... .........i........................../ ............ le-10 g4 (1) ,--I i........../,".................................i.........-1.............. c6 le-12 0 ,--I © le-14 J........./................................... / le-16 +........t--i...............!..............................!..............._.............. :'.......1.....!i...............!i...............!i...............i...............i.............. le-18 I I i i le-20 0 200 400 Time i i i 600 800 1000 1200 steps Figure 13: Global error versus time in small grid for FDTD method and LBS. 21

Original page 26

0 0035 0.003 ....................................... 0 0025 0.002 > 0 0015 _J :Outer 0.001 i 0 0005 ........ i O _4 0 4J O -0 0005 M -0.001 -0 0015 I I -0.002 200 400 Time steps I I FDTD -- LBS with PML ...... LBS no PML ....... b@undary reflecti:ons iiiiiiii I I I 600 800 I000 1200 Figure 14: Electric field versus time sampled at upper sample point in small grid for LBS with and without PML boundary condition. in the present work. It is anticipated this will be the subject of future reports and articles. The results indicate that the LBS is a very promising alternative to a conventional FDTD algorithm for many applications. Higher-order extensions are available for the 2D case, but were not explored presently [36]. Extensions to three-dimensional problems should be straightforward. References [1] D. S. Butler, "The numerical solution of hyperbolic systems of partial differential equations in three independent variables," Proc. ofthe Royal Soc. of London, vol. 255A, pp. 232 252, 1960. [2] M. B. Abbott, An Introduction to the Method ofCharacteristics, American Elsevier, New York, 1966. [3] J. D. Hoffman V. H. Ransom and H. D. Thompson, "A second-order bicharacteristics method for three-dimensional, steady, supersonic flow," AIAA Journal, vol. 10, no. 12, pp. 1573 1581, Dec. 1972. [4] M. C. Cline and J. D. Hoffman, "The analysis of nonequilibrium, chemically reacting, supersonic flow in three-dimensions using a bicharacteristic method," Journal oJComp. Phys., vol. 12, pp. 1 23, 1973. [5] M. J. Zucrow and J. D. Hoffman, Gas Dynamics, John Wiley and Sons, New York, 1975. [6] R. A. Delaney and R Kavanagh, "Transonic flow analysis in axial-flow turbomachinery cascades by a time-dependent method of characteristics," Transactions of ASME, Journal ofEng. Power, vol. 107, pp. 356 364, Jan. 1983. 22

Original page 27

[7] Y. W. Shin and R. A. Palentin, "Numerical analysis of fluid-hammer waves by the method of characteristics," Journal of Comp. Phys., vol. 20, pp. 220 237, 1976. [8] J. D. Hoffman, "The method of characteristics applied to unsteady one-, two- and three-dimensional flows," Tech. Rep. TR-80-07, Thermal Sciences and Propulsion Center, School of Mechanical Engineering, Purdue Univ., 1980. [9] J. Padyak and J. D. Hoffman, "Flow computations in inlets at incidence using a shock fitting bicharacteristics method," AIAA Journal, vol. 18, pp. 1495 1502, Dec. 1980. [10] J. Vadyak and J. D. Hoffman, "Shock-fitting bicharacteristic algorithm for three-dimensional scarfed nozzle flowfields," AIAA Journal, vol. 21, pp. 23 30, Jan. 1983. [11] B. N. Wang, On the Method of Characteristics and Its' Application to the Calculations of Annular Nozzle Flowfields, Ph.D. thesis, Purdue University, West Lafayette, IN, 1984. [12] D. L. Marcum and J. D. Hoffman, "Calculation of unsteady three-dimensional subsonic/transonic inviscid flowfields by the method of characteristics," in AIAA 22nd Aerospace Sciences Meeting, Reno, NV, Jan. 1984, vol. AIAA 84-0440. [13] D. L. Marcum and J. D. Hoffman, "Calculation of three-dimensional flowfields by the unsteady method of characteristics," AIAA Journal, vol. 23, no. 10, pp. 1497 1505, Oct. 1985. [14] D. L. Marcum and J. D. Hoffman, "Calculation of viscous nozzle flows by the unsteady method of characteristics," inAIAA 23rdAerospace Sciences Meeting, Reno, NV, Jan. 1985, vol. AIAA 85-0131. [15] D. L. Marcum and J. D. Hoffman, "Subsonic/transonic/supersonic nozzle flows and nozzle integration," in Numerical MethodsJbr Engine-AirJ?ame Integration, S. N. B. Murthy and Gerald C. Paynter, Eds., pp. 350498. American Institute of Aeronautics and Astronautics, 1986. [16] C. R Kentzner I. H. Parpia and M. H. Williams, "Multidimensional time dependent method of characteristics," Computers andFluids, vol. 16, no. 1, pp. 105 117, 1988. [17] J. D. Hoffman, Numerical Methods Jot Engineers and Scientists, McGraw-Hill, New-York, 1992. [18] J. S. Shang, "Characteristic based methods for the time-domain Maxwell equations," in AIAA 29th Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1991, vol. AIAA 91-0606. [19] J. S. Shang, "A characteristic-based algorithm for solving 3-d time-domain Maxwell equations," in AIAA 30th Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1992, vol. AIAA 92-0452. [20] J. S. Shang, "A fractional-step method for solving 3D time-domain Maxwell equations," in AIAA 31st Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1993, vol. AIAA 93-0461. [21] J. S. Shang and D. Gaitonde, "Characteristic-based, time-dependent Maxwell equations solvers on a general curvilinear frame," in AIAA 24th Plasmadynamics & Lasers Conference, Orlando, FL, July 1993, vol. AIAA 93-3178. [22] K. C. Hill J. S. Shang and D. Calahan, "Performance of a characteristic-based, 3-d time-domain Maxwell equations solvers on a massively parallel computer," in AIAA 24th Plasmadynamics & Lasers Conference, Orlando, FL, July 1993, vol. AIAA 93-3179. 23

Original page 28

[23]J.S.ShangandR.M.Fithen,"Acomparativestudyofnumericalalgorithmsforcomputationalelectromagnetics,"inAIAA 25th Plasmadynamics & Lasers Conference, Colorado Springs, CO, June 1994, vol. AIAA 94-2410. [24] J. S. Shang, "Characteristic-based algorithms for solving the Maxwell equations in the time domain," IEEEAntennas and Propagation Magazine, vol. 37, no. 3, pp. 15 25, June 1995. [25] J. S. Shang, "A fractional-step method for solving 3d, time-domain Maxwell equations," Journal of Comp. Phys.,vol. 118, pp. 109 119, 1995. [26] J. S. Shang, "Characteristic-based algorithms for solving the maxwell equations in the time domain," IEEEAntennas and Propagation Magazine, vol. 37, no. 3, pp. 15 25, June 1995. [27] J. S. Shang and D. Gaitonde, "On high resolution schemes for time-dependent maxwell equations," in AIAA 34th Aerospace Sciences Meeting & Exhibit, Reno, NV, Jan. 1996, vol. AIAA 96-0832. [28] D. Gaitonde and J. S. Shang, "High-order finite-volume schemes in wave propagation phenomena," in AIAA 27th Plasmadynamics & Lasers Conference, New Orleans, LA, June 1996, vol. AIAA 96-2335. [29] D. C. Blake and J. S. Shang, "A procedure for rapid prediction of electromagnetic scattering from complex objects," in AIAA 29th Plasmadynamics vol. AIAA 98-2925. & Lasers Conference, Albuquerque, NM, June 1998, [30] John H. Beggs and W. Roger Briley, "An implicit characteristic based method for computational electromagnetics," Tech. Rep. MSSU-EIRS-ERC-98-11, Mississippi State University, August 1998. [31] J. H. Beggs, D. L. Marcum and S. L. Chan, "The numerical method of characteristics for electromagnetics," Applied Computational Electromagnetics Society Journal, vol. 14, no. 2, pp. 25 36, July 1999. [32] J. R Thomas and R L. Roe, "Development of non-dissipative numerical schemes for computational aeroacoustics," AIAA, 1993, paper number 93-3382-CR [33] R Roe, "Linear bicharacteristic schemes without dissipation," Tech. Report 94-65, ICASE, NASA/Langley Research Center, Hampton, VA, 1994. [34] B. Nguyen and R Roe, "Application of an upwind leap-frog method for electromagnetics," in Proc. l Oth Annual Review of Progress in Applied Computational Electromagnetic's, Monterey, CA, March 1994, Applied Computational Electromagnetics Society, pp. 446_458. [35] J. R Thomas, C. Kim and P. Roe, "Progress toward a new computational scheme for aeroacoustics," in AIAA 12th Computational Fluid Dynamics Conference. AIAA, 1995. [36] J. R Thomas, An Investigation of the Upwind Leapfivg Method for Scalar Advection and Acoustic/Aeroacoustic Wave Propagation Problems, 1996. Ph.D. thesis, University of Michigan, Ann Arbor, MI, [37] B. Nguyen, Investigation of Three-Level Finite-Difference Time-Domain Methods for Multidimensional Acoustics and Electromagnetic's, Ph.D. thesis, University of Michigan, Ann Arbor, MI, 1996. 24

Original page 29

[38] C. Kim, Multidimensional Upwind Leapfkog Schemes and Their Applications, Ph.D. thesis, University of Michigan, Ann Arbor, MI, 1997. [39] K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302_07, Mar. 1966. [40] J. H. Beggs and S. L. Chan, "The linear bicharacteristic scheme for computational electromagnetics," IEEE Trans. Antennas Propagat., 2001, submitted. [41] A. Iserles, "Generalized leapfrog methods," IMA Journal ofNumerical Analysis, vol. 6, pp. 381 392, 1986. [42] A. Yefet and R Petropoulous, "A non-dissipative staggered fourth-order accurate explicit finitedifference scheme for the time-domain Maxwell's equations," Tech. Report 99-30, ICASE, NASA/Langley Research Center, Hampton, VA, 1999. [43] A. Taflove, Ed., Advances in Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston, MA, 1998. [44] A. Taflove, Computational Electrodynamics: The Finite-Difference Time-Domain Method, Artech House, Boston, MA, 1995. [45] J.-R Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," Journal of ComputationaIPhysics, vol. 114, no. 1, pp. 185 200, 1994. [46] W. C. Chew and W. H. Weedon, "A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates," Microwave and Optical Technologies Letters, vol. 7, no. 13, pp. 599 604, Sept. 1994. [47] John H. Beggs, The Linear Bicharacteristic Scheme fbr Electromagnetics, NASA/Langley Research Center, Hampton, VA, May 2001, NASA-TM-2001-210861. 25

Original page 30

REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704-0188 Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704_)133), Washington, DC 20503. 1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE May 2002 4. TITLE AND SUBTITLE 3. REPORT TYPE AND DATES COVERED Technical Memorandum 5. FUNDING NUMBERS A Two-Dimensional Linear Bicharacteristic Scheme for Electromagnetics 706-31-41-01 6. AUTHOR(S) John H. Beggs 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) NASA Langley Research Center Hampton, VA 23681-2199 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) National Aeronautics and Space Administration Washington, DC 20546-0001 11. SUPPLEMENTARY NOTES 12a. DISTRIBUTION/AVAILABILITY STATEMENT Unclassified-Unlimited Subject Category 33 Distribution: Standard Availability: NASA CASI (301) 621-0390 13. ABSTRACT (Maximum 200 words) 8. PERFORMING ORGANIZATION REPORT NUMBER L-18154 10. SPONSORING/MONITORING AGENCY REPORT NUMBER NASA/TM-2002-211663 12b. DISTRIBUTION CODE The upwind leapfrog or Linear Bicharacteristic Scheme (LBS) has previously been implemented and demonstrated on one-dimensional electromagnetic wave propagation problems. This memorandum extends the Linear Bicharacteristic Scheme for computational electromagnetics to model Iossy dielectric and magnetic materials and perfect electrical conductors in two dimensions. This is accomplished by proper implementation of the LBS for homogeneous Iossy dielectric and magnetic media and for perfect electrical conductors. Both the Transverse Electric and Transverse Magnetic polarizations are considered. Computational requirements and a Fourier analysis are also discussed. Heterogeneous media are modeled through implementation of surface boundary conditions and no special extrapolations or interpolations at dielectric material boundaries are required. Results are presented for two-dimensional model problems on uniform grids, and the FDTD algorithm is chosen as a convenient reference algorithm for comparison. The results demonstrate that the two-dimensional explicit LBS is a dissipation-free, second-order accurate algorithm which uses a smaller stencil than the FDTD algorithm, yet it has less phase velocity error. 14. SUBJECT TERMS computational electromagnetics, FDTD methods 17. SECURITY CLASSIFICATION 18. SECURITY CLASSIFICATION OF REPORT OF THIS PAGE U nclassified Unclassified NSN 7540-01-280-5500 15. NUMBER OF PAGES 30 16. PRICE CODE A03 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT OF ABSTRACT Unclassified Standard Form 298 (Rev. 2-89) Prescribed by ANSI Std. Z39-18 293-102