Unformatted text preview: A 3D Perfectly Matched Medium from Modi ed Maxwell's Equations with Stretched Coordinatesy
Weng Cho Chew and William H. Weedon Electromagnetics Laboratory Department of Electrical and Computer Engineering University of Illinois, Urbana, IL 61801 Key Terms: Maxwell's equations, coordinate stretching, perfectly matched
layer, nite di erence time domain, massively parallel computer. Abstract
A modi ed set of Maxwell's equations is presented that includes complex coordinate stretching along the three cartesian coordinates. The added degrees of freedom in the modi ed Maxwell's equations allow the speci cation of absorbing boundaries with zero re ection at all angles of incidence and all frequencies. The modi ed equations are also related to the perfectly matched layer that was presented recently for 2D wave propagation. Absorbing material boundary conditions are of particular interest for nite di erence time domain (FDTD) computations on a singleinstruction multipledata (SIMD) massively parallel supercomputer. A 3D FDTD algorithm has been developed on a Connection Machine CM5 based on the modi ed Maxwell's equations and simulation results are presented to validate the approach. 1. Introduction
The nite di erence time domain method 1, 2] is widely regarded as one of the most popular computational electromagnetics algorithms. Although FDTD is conceptually very simple and relatively easy to program, the method is actually quite e cient since it involves O(N 1 5) computational complexity in 2D and O(N 1 33) computational complexity in 3D 3]. In fact, FDTD can be considered an optimal algorithm since O(N ) numbers are
: : This work is supported by the O ce of Naval Research under grant N0001489J1286, the Army Research O ce under contract DAAL0391G0339, the National Science Foundation under grant NSFECS9224466, and NASA under grant NASANAG2871. Computer time is provided by the National Center for Supercomputer Applications at the University of Illinois, UrbanaChampaign. File: pml.tex, January 12, 1995 Appeared in Micro. Opt. Tech. Lett. , vol. 7, no. 13, pp. 599604, September 1994. y 1 produced in O(N ) operations. FDTD is also ideally suited for implementation on a singleinstruction multipledata (SIMD) massively parallel computer. The reason is that the stencil operations that must be computed at each node of the space grid involve only nearestneighbor interactions and may be implemented at a minimum communication cost 4]. A major challenge, however, is in implementing absorbing boundary conditions (ABCs) at the edges of the FDTD grid. On scalar and vector computers, these boundary conditions are typically computed using methods such as the EngquistMajda 5], Mur 6], Liao 7] or Higdon 8] ABC. However, these methods are not ideal for parallel supercomputers since they all involve communication with many elements normal to the grid boundary. Such communication can easily surpass the time spent computing core FDTD operations in the grid interior, especially for higherorder boundary conditions, and hence can become a bottleneck in the FDTD code. Also, they do not allow for SIMD operation on a parallel machine without the use of masking. An alternate method of implementing an ABC is to use a conventional absorbing material boundary 4, 9{14] . For SIMD parallel computation, these methods have the advantage that the ABC may be implemented with the same FDTD stencil operation as the interior nodes by modifying the conductivity material parameter at the edge of the FDTD grid. The disadvantage is that the re ection coe cient at the absorbing border is zero only at normal incidence and is both angle and frequency dependent. Consequently, the absorbing material border region must be made quite largetypically 20{100 grid points along each edge in order to minimize re ections. Recently, Berenger 15] suggested a more general method of implementing an absorbing material boundary condition. Berenger proposed a procedure for 2D wave propagation whereby Maxwell's equations are generalized and added degrees of freedom are introduced. The added degrees of freedom allow the speci cation of absorbing borders with zero reection coe cient at all angles of incidence and all frequencies. Moreover, the generalized 2 Maxwell's equations reduce to the familiar Maxwell's equations as a special case and hence the same generalized equations can be used to propagate elds in both the interior region and absorbing region. Although the interface between the interior region and the absorbing boundary is re ectionless, there is still a re ection from the edge of the grid. The advantage of using Berenger's procedure is that much larger conductivity values may be speci ed in the absorbing region, leading to a drastic reduction in the number of grid points required for the absorbing boundary. In the present paper, a formulation similar to the Berenger idea is derived for 3D wave propagation from rst principles using a coordinate stretching approach. The advantage of the new method for SIMD parallel computation is stressed. The method is validated with 3D FDTD numerical computations on a Thinking Machines Corporation Connection Machine CM5. 2. Modi ed Maxwell's Equations
For a general medium, we de ne the modi ed Maxwell's equations in the frequency domain, assuming e; time dependence, as
i!t r r e E = i! H H = ;i! E E= H=0 (1) (2) (3) (4) (5) (6) h r r
where h e r = x e1 @x + y e1 @y + z e1 @z ^ @ ^ @ ^ @
e x y z r = x h1 @x + y h1 @y + z h @z : ^ @ ^ @ ^1 @
h x y z i i e h i i In the above, e , h , i = x y z are coordinatestretching variables that stretch the x y z coordinates for r and r . It shall be shown later that when e and h are complex 3 numbers, the medium can be lossy. Note that (3) and (4) are derivable from (1) and (2). A general plane wave solution to Equations (1) { (4) has the form E = E0 e k r
i (7) (8) and H = H0 e k r
i x y z where k = xk + yk + z k . Substituting Equations (7) and (8) into Equations (1) and (2) ^ ^ ^ above gives k k
e kx ex ky ey kz ez h h e E=! H H = ;! E
kx ky kz hx hy hz (9) (10) where k = x + y + z and k = x + y + z . Combining the above, we have ^ ^ ^ ^ ^ ^ ;! 2 H = k
e H = k (k H ) ; (k k )H :
e h h e e h k (11) But from Equation (9), k H = 0 for a homogeneous medium. This gives the dispersion relation
!2 =k k
e h (12) (13) or
2 = e 1 k2 + e 1 k2 + e 1 k2 h h h
x x x y y y z z z where 2 = !2 . Equation (13) is the equation of an ellipsoid in 3D and is satis ed by
k =
x q e h sin cos
x x (14) (15) (16) k =
y q e h sin sin
y y and
k =
z q e h cos :
z z 4 Note that when e , h , i = x y z are complex, the wave in the x, y, and z directions are attenuative and can be independently controlled. Under the matching condition, e = h , e = h , and e = h , we have jk j2 = jk j2 = 2 . The wave impedance is then given by
i i x x y y z z e h jE j = jk j = r = jH j !
h (17) irrespective of the values for e i = x y z and the direction of propagation.
i 3. Single Interface Problem
Assume that a plane wave is obliquely incident on the interface z = 0 in Figure 1. Furthermore, we may assume that the plane wave is of arbitrary polarization. The incident eld may be decomposed into a sum of two components, one with electric eld transverse to z (TE ) and the other with magnetic eld transverse to z (TM ). We will examine these two components individually.
z z z reflected wave y incident wave x transmitted wave z=0 plane region 1 region 2 Figure 1:
z = 0.
z Plane wave with arbitrary polarization incident on the plane In the (TE ) case, we let the incident eld in region 1 be given as E = E0 e k
i hi i i r: (18) In the above, k E0 = 0, and E0 is in the xy plane. Similarly, we de ne the re ected eld in region 1 as E = RTE E0 e k r (19) 5
r r i r and the transmitted eld in region 2 as E = T TE E0
t ix rx t e kt r :
i iy ry ty (20) Phase matching requires that k = k = k and k = k = k . Hence, we can de ne E0 = E0 = E0 since they all point in the same direction. Applying the boundary condition that the tangential electric eld must be continuous across the plane z = 0*, we have 1 + RTE = T TE: (21)
tx r t The magnetic eld may be determined using Equation (9) for regions 1 and 2 as H1 = k ! E0 e k r + RTE k
ie i i re 1 ! 1 E0 e k i r r (22) (23)
z iz z and
kix ex kiy ey kiz ez where k = x +^ +^ and similarly for k and k . We also de ne k1 = k , k2 = ^ y z k and note that k = ;k1 . Then equating the tangential components of Equations (22) and (23), we have h i k1 e2 2 1 ; RTE = T TE k2 e1 1 : (24)
ie re te tz rz z z z z z H2 = T TE k ! E0 e k
te i t r 2 Combining Equations (21) and (24), we have
RTE = k1 e2 k1 e2
z z z 2 2 z + k2 e1
z z ; k2 e1
z z 1 1 z (25) (26) and 2k1 e2 2 : k1 e2 2 + k2 e1 1 Applying a similar procedure to the TM component, we have
T TE =
z z z z z z RTM = k1 h2 k1 h2
z z z z 2 2 ; k2 h1
z z + k2 h1 z z 1 1 (27) * This boundary condition follows from the modi ed Maxwell's equation (1). 6 and
T TM = 4. A Perfectly Matched Interface
q 2k1 h2 2 k1 h2 2 + k2 h1
z z z z z x z 1 : (28) The phase matching condition requires that k1 = k2 and k1 = k2 , or
x y y 1 e1 h1 sin
x x 1 cos sin 1 = = q 2 e2 h2 sin
x x 2 cos sin 2 (29) (30) and
1 q e1 h1 sin
y y q 1 1 2 e2 h2 sin
y y 2 2 where 1 = !p 1 1 and 2 = !p 2 2. For a perfectly matched medium, we choose 1 = 2 , 1 = 2 , e = h and e = h . Equations (29) and (30) become
x x y y e1 sin
x 1 cos sin 1 = e2 sin 2 cos
x 2 (31) (32) and
e1 sin
y x x y 1 1
y = e2 sin 2 sin 2 :
y If we now choose e1 = e2 and e1 = e2 , then 1 = 2, 1 = 2 and we can show that both RTE = 0 and RTM = 0 for all angles of incidence and all frequencies. If region 1 is a vacuum, then = 0, = 0, and (e1 e1 e1 h1 h1 h1 ) = (1 1 1 1 1 1):
x y z x y z (33) In order to have a lossy region 2 with no re ections at the region 1/region 2 interface, we choose (e2 e2 e2 h2 h2 h2 ) = (1 1 s2 1 1 s2) (34)
x y z x y z where s2 is a complex number. In this case,
k1 = k2 =
x x 0 0 sin cos sin sin (35) (36) k1 = k2 =
y y 7 k1 =
z 0 cos cos (37) (38) where 0 = !p 0 0 . If s2 = s02 + is00 , the wave will attenuate in the z direction. This kind 2 of interface is useful for building material ABCs in a FDTD simulation. k2 =
z 0 s2 5. Modi ed Equations in the Time Domain
For the general case of a matched medium, we let e = h = s , e = h = s and e = h = s . Then, r = r = x 1 + y 1 + z 1 . In Equation (1), we write the ^ ^ ^ curl as @ @ @ r E = s1 @x x E + s1 @y y E + s1 @z z E : ^ ^ ^ (39)
x x x y y y z z z e h @ @ @ sx @x sy @y sz @z e x y z sx sy sz Then, de ning H , H , and H in terms of the components of Equation (39), we let
i! i! @ ^ H = s1 @x x E
sx x (40) (41) (42) @ H = s1 @y y E ^
sy y and
i!
sx sy sz where H = H + H + H . Similarly, we can write Equation (2) as
@ ^ ;i! E = s1 @x x H
sx x @ H = s1 @z z E ^
sz z (43) (44) (45) @ ;i! E = s1 @y y H ^
sy y and where E = E + E + E . Note that H , E , i = x y z are twocomponent vectors. We now let s = 1 + i =! , s = 1 + i =! and s = 1 + i =! . Writing Equations (40) { (42) and (43) { (45) in the time domain, we have
sx sy sz si si x x y y z z @ ;i! E = s1 @z ^ H : z
sz z @H x + @t
s x @ H = ; @x x E ^
sx (46) 8 @H y + @t @H z + @t
s s y @ H = ; @y y E ^
sy (47) (48) (49) (50) (51) z @ H = ; @z z E ^
sz and @E x + @t @E y + @t @E z + @t
s s s x @ E = @x x H ^
sx y @ E = @y y H ^
sy sz z @ E = @z ^ H : z Equations (46) { (51) described 3D wave propagation in a perfectly matched medium. The wave propagation phenomenon described by these equations is very similar to that described by Maxwell's equations with the exception that attenuation may be controlled through the , and variables. The FDTD implementation of these equations on a Yee FDTD grid is straightforward. Absorbing boundaries at the edges of the simulation region may be created by choosing appropriate values of , and . Equations (46) { (51) may be seen to include Berenger's equations 15] as a subset for the 2D TE or TM case. The above equations involve 12 components of electromagnetic elds. For a freespace/lossy medium interface, a scheme may be devised using only 10 eld components for the 3D case, and only 3 components for the 2D case. However, this is achieved at the loss of SIMD operation on a parallel machine.
x y z x y z 6. Computer Simulation Results
In order to demonstrate the new method, a 3D orthogonal grid FDTD algorithm was developed based on Equations (46) { (51). The FDTD algorithm was implemented as a SIMD code on the Thinking Machines Corporation Connection Machine CM5. The algorithm operates very e ciently on the CM5 because the FDTD stencil operations that need to be computed at each node involve only nearestneighbor interactions. The communication operations resulting from the nearestneighbor interactions are at a minimum 9 cost since the neighboring processors are for the most part at the bottom of the fattree communication network, where communication bandwidth is maximum. To validate our 3D FDTD algorithm, we solved a simple problem of computing the eld radiated from an in nitesimal electric dipole in free space. An analytic solution was also computed in the frequency domain for many excitation frequencies. The frequency domain solution was then multiplied by the spectrum of FDTD source pulse and inverse Fourier transformed to yield a timedomain analytic solution for comparison with the FDTD solution. The FDTD solution was solved in a cubic region of dimension (N N N ) = (128 128 32) grid points. The grid parameters chosen were x = y = z = 2:5 mm, t = 4:5 ps and N = 512 time steps were computed. The in nitesimal electric dipole was simulated by exciting the E eld in a single grid cell with the source pulse 1 h4(t= )3 ; (t= )4i e; (52) J (t) =
x y z t y x y z where = 1=4 f0 and a value of f0 = 1:0 GHz was chosen. The dipole source was located
y x y z x y t= at grid location (n n n ) = (91 64 16). The E and E elds were obtained by sampling the elds at grid location (n n n ) = (37 91 16). The absorbing boundaries used for the FDTD simulation consisted of planar layers of thickness 8 grid points on all surfaces. Along the borders parallel to x axis, the value of was speci ed, while and were speci ed on the borders parallel to the y and z axis, respectively. The conductivity values were chosen with a parabolic taper decreasing from the maximum value towards the center of the grid such that the re ection coe cient at normal incidence was R0 = :0001, The E eld computed using both the analytic formulation and the FDTD algorithm are overlaid in Figure 2. The curves due to the analytic and numerical solutions are barely distinguishable, indicating excellent agreement. Similarly, the E eld due to the analytic and numerical solutions are overlaid in Figure 3. Again, we see excellent agreement. 10
x y z x y z x y Any di erence between the analytic and numerical solutions in Figures 2 and 3 may be attributed to modeling errors such as the nite size of the dipole source and the discrete approximation of Maxwell's equations in addition to re ections due to imperfections in the absorbing boundaries.
4000 2000 Amplitude of Ex field 0 2000 4000 6000 8000 10000 0 100 200 300 400 500 600 Time step Figure 2: Analytic and numerical FDTD solution overlaid for the E eld resulting from an in nitesimal electric dipole.
x The CM5 machine used to solve the FDTD problem is located at the National Center for Supercomputing Applications (NCSA) at the University of Illinois. The program was written in CM Fortran and compiled using CMF version 2.1. The CM5 at the NCSA has 512 nodes with vector units. CPU times were determined by running the problem on 32, 64, 128 and 256 node partitions. For this problem, a total of 0.5 million unknown eld quantities (128 128 32 grid) were determined for 512 time steps. The CPU times are shown in Table 1. 7. Conclusions
11 1 x 10 4 0.5 Amplitude of Ey field 0 0.5 1 1.5 0 100 200 300 400 500 600 Time step Figure 3: Analytic and numerical FDTD solution overlaid for the E eld resulting from an in nitesimal electric dipole.
y Table 1: CPU times for FDTD Problem on CM5
Nodes 32 64 128 256 CPU sec (Run 1, Run 2, Run 3, Avg.) 50.5, 50.2, 50.6; 50.4 29.9, 30.0, 30.0; 30.0 17.9, 18.4, 18.4; 18.2 12.4, 13.2, 12.7; 12.8 A modi ed set of Maxwell's equations have been introduced using complex coordinate stretching factors along the three cartesian coordinate axis. This modi cation introduces additional degrees of freedom in Maxwell's equations such that absorbing boundaries may be speci ed with zero re ection coe cient at all frequencies and all angles of incidence. 12 The formulation was shown to be related to the perfectly matched layer that was recently derived by Berenger for 2D wave propagation. A 3D FDTD algorithm was developed from the modi ed Maxwell's equations that uses the re ectionless absorbing interface property to implement radiation boundary conditions at the edges of the FDTD grid. The accuracy of the algorithm was validated by computing the eld radiated from an in nitesimal electric dipole and comparing against a known analytical expression. The FDTD algorithm was implemented on the Connection Machine CM5 and timing results were presented. This breakthrough in absorbing material boundary conditions allows EM scattering to be computed very e ciently on SIMD parallel computers. Acknowledgement
The authors wish to thank J. P. Berenger for sending us a preprint of his work and to A. Ta ove for bringing to our attention J. P. Berenger's work. 13 References
1] K. S. Yee, \Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. Antennas Propagat., vol. AP14, pp. 302{ 307, 1966. 2] A. Ta ove, \Review of the formulation and applications of the nitedi erence timedomain method for numerical modeling of electromagnetic wave interactions with arbitrary structures," Wave Motion, vol. 10, pp. 547{582, 1988. 3] W. C. Chew, Waves and Fields in Inhomogeneous Media. New York: Van Nostrand, 1990. 4] W. H. Weedon, W. C. Chew, and C. M. Rappaport, \Computationally e cient FDTD simulation of openregion scattering problems on the connection machine CM5," in IEEE Antennas and Propagation Society International Symposium Digest, (Seattle, WA), June 19{24, 1994. 5] B. Engquist and A. Majda, \Absorbing boundary conditions for the numerical simulation of waves," Math. Computation, vol. 31, pp. 629{651, 1977. 6] G. Mur, \Absorbing boundary conditions for the nitedi erence approximation of the timedomain electromagnetic eld equations," IEEE Trans. Electromag. Compat., vol. EMC23, pp. 377{382, 1981. 7] Z. P. Liao, H. L. Wong, B. P. Yang, and Y. F. Yuan, \A transmitting boundary for transient wave analysis," Scientia Sinica. (Series A), vol. 27, no. 10, pp. 1063{1076, 1984. 8] R. L. Higdon, \Numerical absorbing boundary conditions for the wave equation," Math. Comput., vol. 49, pp. 65{90, 1987. 14 9] I. Katz, D. Parks, A. Wilson, M. Rotenberg, and J. Harren, \Nonre ective free space boundary conditions for SGEMP codes," Systems, Science and Software, vol. SSSR762934, 1976. 10] R. Holland and J. W. Williams, \Total eld versus scattered eld nitedi erence codes: A comparative assessment," IEEE Trans. Nuclear Sci., vol. NS30, pp. 4583{ 4588, 1983. 11] J.P. Berenger in Actes du Colloque CEM, (Tregastel, France), 1983. 12] C. Cerjan, D. Koslo , R. Koslo , and M. Reshef, \A nonre ecting boundary condition for discrete acoustic and elastic wave equations," Geophysics, vol. 50, pp. 705{708, 1985. 13] R. Koslo and D. Koslo , \Absorbing boundaries for wave propagation problems," J. Computational Physics, vol. 63, pp. 363{376, 1986. 14] C. M. Rappaport and L. Bahrmasel, \An absorbing boundary condition based on anechoic absorber for EM scattering computation," J. Electromag. Waves Appl., vol. 6, no. 12, pp. 1621{1634, 1992. 15] J.P. Berenger, \A perfectly matched layer for the absorption of electromagnetic waves," J. Computational Physics, accepted for publication, 1994. 15 ...
View
Full Document
 Fall '09
 Electromagnet, Electromagnetic wave equation, fdtd algorithm, 3D FDTD algorithm

Click to edit the document details