5/5/2015 Instructor Dr. Raymond Rumpf (915) 747‐6958 [email protected] EE 4395/5390 – Special Topics Computational Electromagnetics (CEM) Lecture #10 Maxwell’s Equations on a Yee Grid These notes may contain copyrighted material obtained under fair use rules. Distribution of these materials is strictly prohibited Lecture 10 Slide 1 Outline • Yee Grid • Maxwell’s equations on a Yee grid • Finite‐difference approximations of Maxwell’s equations on a Yee grid • Matrix form of Maxwell’s equations • Numerical dispersion • Generalization to fully anisotropic materials • Alternative grids Lecture 10 Bonus Slide 2 1 5/5/2015 Yee Grid Lecture 10 Slide 3 Kane S. Yee Kane S. Yee was born in Canton, China on March 26, 1934. He received a B.S.E.E., M.S.E.E and Ph.D. in Applied Mathematics from the University of California at Berkely in 1957, 1958, and 1963, respectively. He did research on electromagnetic diffraction while employed by Lockheed Missile and Space Co. (1959‐1961). He has been associated with the Lawrence Livermore Laboratory since 1963. At present he is a professor in mathematics at Kansas State University. His main areas of interest are electromagnetics, hydrodynamics and numerical solution to partial differential equations. K. S. Yee, “A Closed‐Form Expression for the Energy Dissipation in a Low‐Loss Transmission Line,” IEEE Trans. Nuclear Science, vol. 21, no. 1, pp. 1006‐1008, 1974. Kane S. Yee Lecture 10 Slide 4 2 5/5/2015 3D Grids A three‐dimensional grid looks like this: One unit cell from the grid looks like this: z y x N x 10, N y 10, N z 15 x , y , z grid resolution parameters Lecture 10 Slide 5 Collocated Grid Within the unit cell, where should we place Ex, Ey, Ez, Hx, Hy, and Hz? A straightforward approach would be to locate all of the field components at a common point within in a grid cell; perhaps at the origin. z Ez Hy Ex Hx x Lecture 10 Hz Ey y Slide 6 3 5/5/2015 Yee Grid Instead, we are going to stagger the position of each field component within the grid cells. z Ez Hy Hx Ey x Ex Hz y K. S. Yee, “Numerical solution of the initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Microwave Theory and Techniques, vol. 44, pp. 61–69, 1966. Lecture 10 Slide 7 Stereo Image of Yee Cell To view the Yee cell if full 3D, look past the image above so that they appear double. When the double images overlap so that you see three Yee cells, the middle image will be three‐dimensional. Lecture 10 Slide 8 4 5/5/2015 Yee Cell for 1D, 2D, and 3D Grids 1D Yee Grids Ex 2D Yee Grids z Hy Ez Ex Mode x Hx Ey Ey Mode Hx Hy z Ez Mode z * These are the same for isotropic materials. 3D Yee Grid y Hz Ex x Ey Ez Hy y Hx Ey Hz Mode x Ex Hz Lecture 10 y Slide 9 Reasons to Use the Yee Grid 1. Divergence‐free E 0 H 0 3. Elegant arrangement to approximate curl equations 2. Physical boundary conditions are naturally satisfied Lecture 10 Slide 10 5 5/5/2015 Consequences of the Yee Grid The field components are at physically different positions. This has the following consequences… • Field components within the same grid cell can reside in different materials. – xx multiplies Ex so these functions should be made to exist at the same positions across the grid. – yy multiplies Ey so these functions should be made to exist at the same positions across the grid. – zz multiplies Ez so these functions should be made to exist at the same positions across the grid. – xx, yy and zz are unique arrays and must be constructed separately. • Field components within the same grid cell will be slightly out of phase. – This must be accounted for when constructing sources and when post‐ processing the field data. • The grid causes numerical dispersion where waves propagate slower than a physical wave would. Lecture 10 Slide 11 Visualizing Extended Yee Grids 222 Grid 44 Grid (Ez Mode) j i y x Lecture 10 Slide 12 6 5/5/2015 Finite‐Difference Approximations of Maxwell’s Equations on a Yee Grid Lecture 10 Slide 13 Starting Point We start with Maxwell’s equations in the following form. Ez E y k0 xx H x y z Ex Ez k0 yy H y z x E y Ex k0 zz H z x y H z H y k0 xx Ex y z H x H z k0 yy E y z x H y H x k0 zz Ez x y Lecture 10 Here we have retained diagonally anisotropic material tensors. This will be needed to incorporate a perfectly matched layer boundary condition. Slide 14 7 5/5/2015 Normalize the Grid Coordinates The grid is normalized according to x k 0 x y k0 y z k0 z This “absorbs” the k0 term into the spatial derivatives and simplifies Maxwell’s equations to Ez E y k0 xx H x y z Ex Ez k0 yy H y z x E y Ex k0 zz H z x y Ez E y xx H x y z Ex Ez yy H y z x E y Ex zz H z x y H z H y xx Ex y z H x H z yy E y z x H y H x zz Ez x y H z H y k0 xx Ex y z H x H z k0 yy E y z x H y H x k0 zz Ez x y Lecture 10 Slide 15 Finite‐Difference Equation for Hx z E yi , j ,k 1 Ezi , j ,k Ez E y xx H x y z H xi , j ,k Hy x Ex Hz E yi , j ,k Ezi , j 1,k y i , j , k 1 Ezi , j 1,k Ezi , j ,k E y y Lecture 10 E yi , j ,k xi ,xj ,k H xi , j ,k z Slide 16 8 5/5/2015 Finite‐Difference Equation for Hy z Exi , j ,k 1 Ex Ez yy H y z x Ezi , j ,k Ezi 1, j ,k x Hx H yi , j ,k Ey Exi , j ,k H z y Exi , j ,k 1 Exi , j ,k Ezi 1, j ,k Ezi , j ,k yi ,yj ,k H yi , j ,k z x Lecture 10 Slide 17 Finite‐Difference Equation for Hz z Ez E y Ex zz H z x y Hx Hy E yi , j ,k x Exi , j ,k y H zi , j ,k E yi 1, j ,k E i , j 1, k x E yi 1, j ,k E yi , j ,k Exi , j 1,k Exi , j ,k zi z, j ,k H zi , j ,k x y Lecture 10 Slide 18 9 5/5/2015 Finite‐Difference Equation for Ex i , j ,k i , j , k 1 H zi , j ,k H zi , j 1,k H y H y xxi , j ,k Exi , j ,k y z z H z H y xx Ex y z Ez H zi , j 1,k Hx H yi , j ,k Exi , j ,k H zi , j ,k Ey x y H yi , j ,k 1 Lecture 10 Slide 19 Finite‐Difference Equation for Ey H xi , j ,k H xi , j ,k 1 H zi , j ,k H zi 1, j ,k yi ,yj ,k E yi , j ,k z x z Ez H x H z yy E y z x Hy Ex H zi , j ,k H xi , j ,k Ey x Lecture 10 H xi , j ,k 1 H zi 1, j ,k y Slide 20 10 5/5/2015 Finite‐Difference Equation for Ez H yi , j ,k H yi 1, j ,k H xi , j ,k H xi , j 1,k zzi , j ,k Ezi , j ,k x y z H y H x zz Ez x y H yi 1, j ,k H xi , j 1,k Ezi , j ,k H xi , j ,k H yi , j ,k Ex Hz Ey y x Lecture 10 Slide 21 Summary of Finite‐Difference Approximations of Maxwell’s Equations Ez E y xx H x y z Ex Ez yy H y z x E y Ex zz H z x y H z H y xx Ex y z H x H z yy E y z x H y H x zz Ez x y Lecture 10 i , j , k 1 Ezi , j 1,k Ezi , j ,k E y y E yi , j ,k xxi , j ,k H xi , j ,k z Exi , j ,k 1 Exi , j ,k Ezi 1, j ,k Ezi , j ,k i , j ,k i , j ,k Hy yy z x E yi 1, j ,k E yi , j ,k Exi , j 1,k Exi , j ,k ziz, j ,k H zi , j ,k x y i , j ,k i , j , k 1 H zi , j ,k H zi , j 1,k H y H y xxi , j ,k Exi , j ,k z y H xi , j ,k H xi , j ,k 1 H zi , j ,k H zi 1, j ,k i , j , k i , j ,k Ey yy z x H yi , j ,k H yi 1, j ,k H xi , j ,k H xi , j 1,k ziz, j ,k Ezi , j ,k x y Slide 22 11 5/5/2015 2D Problems In many cases, physical problems can be simplified by assuming the structure is uniform in the z‐direction and wave propagation is restricted to the x‐y plane. For this case, and Maxwell’s equations split into two sets of three z 0 coupled equations. Ezi 1, j , k Ezi , j ,k iyy, j , k H yi , j ,k x E yi 1, j , k E yi , j , k Exi , j 1, k Exi , j ,k zzi , j , k H zi , j ,k x y H yi , j ,k H yi 1, j , k H xi , j ,k H xi , j 1,k zzi , j ,k Ezi , j ,k x y Ezi , j 1, k Ezi , j ,k xi ,xj , k H xi , j ,k y Ezi 1, j , k Ezi , j ,k iy,yj , k H yi , j ,k x Lecture 10 E Mode H yi , j ,k H yi 1, j , k x H zi , j ,k H zi , j 1,k xxi , j ,k Exi , j ,k y H i , j ,k H zi 1, j ,k i , j ,k i , j ,k z yy Ey x i j k i j k , , , 1, H H x x zzi , j ,k Ezi , j ,k y E yi 1, j , k E yi , j , k Exi , j 1, k Exi , j ,k zzi , j , k H zi , j ,k x y H zi , j ,k H zi , j 1,k xxi , j ,k Exi , j ,k y H i , j ,k H zi 1, j ,k z yi ,yj ,k E yi , j ,k x H Mode Ezi , j 1, k Ezi , j ,k xxi , j , k H xi , j ,k y Slide 23 Matrix Form of Maxwell’s Equations Lecture 10 Slide 24 12 5/5/2015 Recall That Fields are Stored in Column Vectors 1‐D Systems E1 2‐D Systems E1 E 2 e E3 E4 E5 E2 E3 E4 E5 E1 E5 E9 E2 E6 E10 E14 E3 E7 E11 E15 E4 E8 E12 E16 Lecture 10 E13 E1 E 2 E3 E4 E5 E6 E 7 E e 8 E 9 E10 E11 E12 E13 E 14 E15 E 16 Slide 25 Matrix Representation of Point‐by‐Point Multiplication (1 of 2) E1 E2 E3 E4 E5 r ,i Ei r1 r 2 r 3 r 4 r 5 r1 E1 r 2 E2 r 3 E3 r 4 E4 r 5 E5 Lecture 10 r1 0 εre 0 0 0 0 E1 r1E1 E r2 0 E 2 r2 2 0 r3 0 E3 r3E3 0 0 r4 0 E4 r4E4 0 0 0 r5 E5 r5E5 0 0 0 0 0 0 Slide 26 13 5/5/2015 Matrix Representation of Point‐by‐Point Multiplication (2 of 2) E1 E2 E3 E4 E5 r ,i Ei r1 r 2 r 3 r 4 r 5 r1 E1 r 2 E2 r 3 E3 r 4 E4 r 5 E5 r1 0 εre 0 0 0 0 E1 r1E1 E r2 0 E 2 r2 2 0 r3 0 E3 r3E3 0 0 r4 0 E4 r4E4 0 0 0 r5 E5 r5E5 0 0 0 0 0 0 Lecture 10 Slide 27 Derivative Operators for Electric Fields (1 of 2) E1 E2 E3 E4 E5 E E Ei i 1 x i 1 x x 2 E2 E1 x E3 E2 x E4 E3 x E5 E4 x E6 E5 x Lecture 10 1 0 1 e 0 De x x 0 0 1 0 0 0 E1 x E1.5 1 1 0 0 E2 x E2.5 0 1 1 0 E3 x E3.5 0 0 1 1 E4 x E4.5 0 0 0 1E5 x E5.5 Slide 28 14 5/5/2015 Derivative Operators for Electric Fields (2 of 2) E1 E2 E3 E4 E5 E E Ei i 1 x i 1 x x 2 E2 E1 x E3 E2 x E4 E3 x E5 E4 x E6 E5 x 1 0 1 e 0 De x x 0 0 1 0 0 0 E1 x E1.5 1 1 0 0 E2 x E2.5 0 1 1 0 E3 x E3.5 0 0 1 1 E4 x E4.5 0 0 0 1E5 x E5.5 Lecture 10 Slide 29 Derivative Operators for Magnetic Fields (2 of 2) H 1 H 2 H 3 H 4 H 5 H H i H i 1 x i 1 x 2 x H 1 H 0 x H H 2 1 x H 3 H 2 x H 4 H 3 x H 5 H 4 x Lecture 10 1 1 1 Dhxh 0 x 0 0 0 0 0 0 H1 x H0.5 1 0 0 0 H 2 x H1.5 1 1 0 0 H3 x H 2.5 0 1 1 0 H 4 x H3.5 0 0 1 1 H5 x H 4.5 Slide 30 15 5/5/2015 Derivative Operators for Magnetic Fields (2 of 2) H 1 H 2 H 3 H 4 H 5 H H i H i 1 x i 1 x 2 x H 1 H 0 x H H 2 1 x H 3 H 2 x H 4 H 3 x H 5 H 4 x 1 1 1 h Dxh 0 x 0 0 0 0 0 0 H1 x H0.5 1 0 0 0 H 2 x H1.5 1 1 0 0 H3 x H 2.5 0 1 1 0 H 4 x H3.5 0 0 1 1 H5 x H 4.5 Lecture 10 Slide 31 Simplest Boundary Conditions Dirichlet Boundary Conditions Assume E6 0 E1 E2 E3 E4 E5 E6 x Periodic Boundary Conditions Assume E6 E1 E1 E2 E3 E4 E5 E1 x E1 E5 x Lecture 10 0 0 E1 1x 1x 0 0 1 1 0 0 E2 x x Dexe 0 0 1x 1x 0 E3 0 0 1x 1x E4 0 0 0 0 0 1xE5 0 0 E1 1x 1x 0 0 1 1 0 0 E2 x x Dexe 0 0 1x 1x 0 E3 0 0 1x 1x E4 0 1x 0 0 0 1x E5 Slide 32 16 5/5/2015 Derivative Operators on a 33 Grid Using Dirichlet Boundary Conditions Nx 1 1 Dirichlet boundary 1 1 conditions 1 0 1 1 1 e Dx 1 1 x 1 0 1 1 1 1 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 e Dy 1 0 0 1 y 1 0 0 1 1 0 0 1 0 1 1 1 1 1 1 0 1 1 Dhx 1 1 x 1 1 0 1 1 1 Dirichlet boundary conditions 1 1 1 0 1 Nx 0 0 1 1 0 0 1 1 Dhy 1 0 0 1 y 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 Note: These matrices have only two diagonals so they are very easy to construct! Lecture 10 Slide 33 2D Derivative Operators for 1D Grids When Nx=1 and Ny>1 Dex Dhx Z zero matrix Dey and Dhy is standard for 1D grid When Nx>1 and Ny=1 Dex and Dhx is standard for 1D grid Dey Dhy Z zero matrix Note: We will do something different when we account for oblique angle of incidence. Lecture 10 0 0 0 e h Dx Dx 0 0 0 0 0 e h Dy Dy 0 0 Slide 34 17 5/5/2015 Size of Derivative Operators 1D Grids If your grid has N x points, your matrices will be N x N x with a total of N x2 elements. 2D Grids If your grid has N x N y points, your matrices will be N x N y N x N y with a total of N x N y elements. 2 3D Grids If your grid has N x N y N z points, your matrices will be N x N y N z N x N y N z with a total of N x N y N z elements. 2 Lecture 10 Slide 35 USE SPARSE MATRICES!!!!!!! WARNING !! The derivative operators will be EXTREMELY large matrices. For a small grid that is just 100200 points: Total Number of Points: Size of Derivate Operators: Total Elements in Matrices: Memory to Store One Full Matrix: Memory to Store One Sparse Matrix: 20,000 20,000 20,000 400,000,000 6 Gb 1 Mb NEVER AT ANY POINT should you use FULL MATRICES in the finite‐difference method. Not even for intermediate steps. NEVER! Lecture 10 Slide 36 18 5/5/2015 Derivative Operators on a 44 Grid Using Dirichlet Boundary Conditions Dex 1 x 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 Dhx 1 x 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 Dey 1 y 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 Dhy 1 y 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 Note: These matrices have only two diagonals so they are very easy to construct! Lecture 10 Slide 37 Dex and on a 3×3×3 Yee Grid Dhx Dex 1 x Dhx Lecture 10 1 x Slide 38 19 5/5/2015 Dey and on a 3×3×3 Yee Grid Dhy Dey 1 y Dhy 1 y Lecture 10 Slide 39 Dez and on a 3×3×3 Yee Grid Dhz Dez 1 z Dhz Lecture 10 1 z Slide 40 20 5/5/2015 Relationship Between the Derivative Operators Transpose Operation A T A j ,i i, j 1 A 1 2 , A T 2 AT = transpose(A); a A 11 a21 1 2 T 1 3 A , A 2 4 3 4 a12 T a11 ,A a22 a22 a21 a22 AT = A.’; Hermitian (Conjugate) Transpose Operation A H A j ,i i, j * a12 a* , A H 1*1 a22 a12 a A 11 a21 AT = ctranspose(A); 1 4 j A 3 2 j * a21 * a22 2 3 j , 4 j 1 4 j AH 2 3 j 3 2 j 4 j AT = A’; Relationship Between the Derivative Operators Dhx Dex H Dhy Dey DHX = -DEX’; H DHY = -DEY’; This means you only have to construct derivative operators for the electric field. The derivative operators for the magnetic field can be computed directly from the electric field derivative operators. This relation does not hold for some boundary conditions such as Neumann. Lecture 10 Slide 41 What About the Second‐Order Derivatives? Recall from Lecture 5, Slide 27 that was a poor approximation of . D1 D1 D 2 Dx1 Dx1 1 2 x 2 1 0 1 0 0 0 2 0 1 0 1 0 2 0 1 0 1 0 2 0 0 0 1 0 1 Dx2 2 1 0 0 0 1 2 1 0 0 1 2 0 1 2 1 0 x 0 0 1 2 1 0 0 0 1 2 What about the derivative operators derived from the Yee grid? 1 1 0 0 0 0 1 1 0 0 1 e 0 0 1 1 0 Dx x 0 0 0 1 1 0 0 0 0 1 Dex Dhx Lecture 10 1 x 2 1 0 0 0 1 1 0 0 1 h 0 1 1 0 Dx x 0 0 1 1 0 0 0 1 2 1 0 0 0 1 2 1 0 0 0 1 2 1 0 0 0 1 2 1 0 0 0 1 1 0 0 0 0 1 The numbers is this matrix may differ slightly from the “ideal” 2nd‐ order derivate operator due to the boundary conditions. Slide 42 21 5/5/2015 Maxwell’s Equations in Matrix Form i , j , k 1 E yi , j , k Ezi , j 1, k Ezi , j , k E y xxi , j , k H xi , j , k y z Ez E y xx H x y z Ex Ez yy H y z x Ey Ex zz H z x y Exi , j , k 1 Exi , j , k Ezi 1, j , k Ezi , j , k i, j,k i, j,k Hy yy z x i 1, j , k i, j,k Ey Ey E i , j 1, k Exi , j , k x zzi , j , k H zi , j , k x y Hz H y xx Ex y z Hx H z yy E y z x Hy H x zz Ez x y i, j,k i , j , k 1 H zi , j , k H zi , j 1, k H y H y xxi , j , k Exi , j , k y z H xi , j , k H xi , j , k 1 H zi , j , k H zi 1, j , k i, j,k i, j,k yy Ey z x H yi , j , k H yi 1, j , k H xi , j , k H xi , j 1, k zzi , j , k Ezi , j , k x y Deye z Deze y μ xxh x Deze x Dexe z μ yy h y Dexe y Deye x μ zz h z Dhyh z Dhzh y ε xxe x Dhzh x Dhxh z ε yy e y Dhxh y Dhyh x ε zz e z Lecture 10 Slide 43 Summary of What We Did No charges E j H H j E Normalized grid & diagonal tensors Ez E y xx H x y z Ex Ez yy H y z x E y Ex zz H z x y H z H y xx Ex y z H x H z yy E y z x H y H x zz Ez x y Lecture 10 Normalized H E k0 r H H k0 r E Finite‐difference approximation Matrix form of Maxwell’s equations i , j , k 1 E yi , j , k Ezi , j 1, k Ezi , j ,k E y De e De e μ h xxi , j , k H xi , j , k y z Exi , j , k 1 Exi , j ,k Ezi 1, j ,k Ezi , j , k i , j ,k i , j ,k yy Hy z x i 1, j , k i , j ,k i , j 1, k i , j ,k Ey Ey E Ex x zzi , j , k H zi , j , k x y i , j ,k i , j , k 1 H zi , j ,k H zi , j 1,k H y H y xxi , j , k Exi , j ,k y z H xi , j ,k H xi , j ,k 1 H zi , j ,k H zi 1, j , k iyy, j , k E yi , j ,k z x i j k i j k , , 1, , H y H y H i , j ,k H xi , j 1, k x zzi , j , k Ezi , j ,k x y y z z y xx x D e D e μ yy h y e z x e x z Dexe y Deye x μ zz h z Dhyh z Dhzh y ε xx e x Dhzh x Dhxh z ε yy e y Dhxh y Dhyh x ε zz e z Slide 44 22 5/5/2015 Numerical Dispersion Lecture 10 Slide 45 Dispersion on a Yee Grid Recall the dispersion relation for an isotropic material with parameters r and r. 2 2 2 2 r r k x k y k z c0 The analogous dispersion relation on a frequency‐domain Yee grid filled with r and r is 2 2 kx x r r sin v 2 x 2 ky y 2 sin y 2 2 2 kz z sin 2 z 2 In this equation, the speed of light c0 is written as v because the velocity changes due to the dispersion of the grid. Lecture 10 46 23 5/5/2015 Drawbacks of Structured Grids (2 of 2) Structured grids exhibit high anisotropic dispersion. Anisotropic Dispersion FDM Lecture 10 Slide 47 Compensation Factor (1 of 2) The numerical dispersion equation is solved for velocity v. 2 2 2 2 k y y 2 kx x 2 kz z v r r sin sin sin 2 2 y 2 z x 1 2 In the absence of grid dispersion, v should be exactly the speed of light c0. Due to the Yee grid, waves slow down by a factor . v c0 We can calculate this factor by combining the above equations. c0 r r Lecture 10 2 kx x sin 2 x 2 ky y 2 sin y 2 2 2 kz z sin 2 z 2 48 24 5/5/2015 Compensation Factor (2 of 2) We can write a simpler and more useful expression for . 2 2 k y y 2 1 2 kx x 2 kz z sin sin sin k0 n x 2 2 y 2 z 2 k0 2 0 n r r Lecture 10 49 Compensating for Numerical Dispersion Given that the wave slows down by factor in the direction of , it k follows that we can compensate for the dispersion by artificially “speeding up” the wave. We do this by decreasing the values of r and r across the entire grid by a factor of . r r r r Notes: 1. We can only compensate for dispersion for one direction k. 2. We can only compensate for dispersion in one set of material values r and r. 3. It is best to choose average or dominant values for these parameters. 4. Choose = 22.5° if nothing else is known. Lecture 10 k 22.5 50 25 5/5/2015 Generalization to Fully Anisotropic Materials Lecture 10 Slide 51 Retain Anisotropic Terms Our analytical equations with just diagonally anisotropic materials were… Ez E y xx H x y z Ex Ez yy H y z x E y Ex zz H z x y H z H y xx Ex y z H x H z yy E y z x H y H x zz Ez x y Lecture 10 For fully anisotropic materials, these are now… Ez E y xx H x xy H y xz H z y z Ex Ez yx H x yy H y yz H z z x E y Ex zx H x zy H y zz H z x y H z H y xx Ex xy E y xz Ez y z H x H z yx Ex yy E y yz Ez z x H y H x zx Ex zy E y zz Ez x y Slide 52 26 5/5/2015 First Guess at Finite‐Difference Approximations E yi , j ,k ? i , j ,k i , j ,k xx H x xyi , j ,k H yi , j ,k xzi , j ,k H zi , j ,k z i , j , k 1 Ezi , j 1,k Ezi , j ,k E y y Exi , j ,k 1 Exi , j ,k Ezi 1, j ,k Ezi , j ,k ? i , j ,k i , j ,k i , j ,k i , j ,k H y yi ,zj ,k H zi , j ,k yx H x yy z x E yi 1, j ,k E yi , j ,k Exi , j 1,k Exi , j ,k ? i , j ,k i , j ,k zx H x zyi , j ,k H yi , j ,k zi z, j ,k H zi , j ,k x y i , j ,k i , j , k 1 ? i , j ,k i , j ,k i , j ,k i , j ,k i , j ,k i , j ,k H zi , j ,k H zi , j 1,k H y H y xx Ex xy E y xz Ez y z H xi , j ,k H xi , j ,k 1 H zi , j ,k H zi 1, j ,k ? i , j ,k i , j ,k yx Ex yi ,yj ,k E yi , j ,k yi ,zj ,k Ezi , j ,k z x H yi , j ,k H yi 1, j ,k H xi , j ,k H xi , j 1,k ? i , j ,k i , j ,k zx Ex zyi , j ,k E yi , j ,k zzi , j ,k Ezi , j ,k x y Lecture 10 Slide 53 The Problem i , j , k 1 E yi , j ,k z xxi , j ,k H xi , j ,k xyi , j ,k H yi , j ,k xzi , j ,k H zi , j ,k Exi , j ,k 1 Exi , j ,k Ezi 1, j ,k Ezi , j ,k z x E yi 1, j ,k E yi , j ,k Exi , j 1,k Exi , j ,k x y yxi , j ,k H xi , j ,k yyi , j ,k H yi , j ,k yi ,zj ,k H zi , j ,k zxi , j ,k H xi , j ,k zyi , j ,k H yi , j ,k zzi , j ,k H zi , j ,k xxi , j ,k Exi , j ,k xyi , j ,k E yi , j ,k xzi , j ,k Ezi , j ,k yxi , j ,k Exi , j ,k yyi , j ,k E yi , j ,k yzi , j ,k Ezi , j ,k zix, j ,k Exi , j ,k zyi , j ,k E yi , j ,k zzi , j ,k Ezi , j ,k Ezi , j 1,k Ezi , j ,k E y y z Ez Hy Hx Ey x Ex Hz y i , j ,k i , j , k 1 H zi , j ,k H zi , j 1,k H y H y y z i , j ,k i , j , k 1 i , j ,k Hx Hx H H zi 1, j ,k z z x H yi , j ,k H yi 1, j ,k H xi , j ,k H xi , j 1,k x y Each term in a finite‐difference approximation MUST exist at the same position in space. The boxed terms exist at physically different locations than the other terms. Lecture 10 Slide 54 27 5/5/2015 The Correction i , j , k 1 E yi , j ,k xyi , j ,k H yi , j ,k xyi 1, j ,k H yi 1, j ,k xyi , j 1,k H yi , j 1, k xyi 1, j 1,k H yi 1, j 1,k xzi , j ,k H zi , j ,k xzi , j , k 1 H zi , j ,k 1 xzi 1, j , k 1 H zi 1, j , k 1 xzi 1, j ,k H zi 1, j ,k E zi , j 1,k Ezi , j ,k E y xxi , j ,k H xi , j ,k 4 4 y z i , j ,k i 1, j , k i 1, j , k i, j,k i 1, j 1,k H xi 1, j 1,k Hx yzi , j ,k H zi , j ,k yzi , j , k 1 H zi , j ,k 1 yzi , j 1, k 1 H zi , j 1, k 1 yzi , j 1,k H zi , j 1,k yxi , j 1,k H xi , j 1,k yx E xi , j , k 1 Exi , j ,k Ezi 1, j ,k Ezi , j ,k yx H x yx i , j ,k i , j ,k yy H y 4 4 z x E yi 1, j ,k E yi , j ,k Exi , j 1,k Exi , j ,k zxi , j ,k H xi , j , k zxi 1, j , k H xi 1, j ,k zxi 1, j , k 1 H xi 1, j , k 1 zxi , j ,k 1 H xi , j ,k 1 zyi , j , k H yi , j ,k zyi , j 1,k H yi , j 1,k zyi , j 1,k 1 H yi , j 1,k 1 zyi , j ,k 1 H yi , j ,k 1 zzi , j ,k H zi , j , k 4 4 x y i , j ,k i , j , k 1 xyi , j ,k E yi , j , k xyi , j 1,k E yi , j 1,k xyi 1, j 1,k E yi 1, j 1,k xyi 1, j ,k E yi 1, j ,k xzi , j ,k Ezi , j ,k xzi , j ,k 1 Ezi , j ,k 1 xzi 1, j ,k 1 Ezi 1, j , k 1 xzi 1, j ,k Ezi 1, j , k H zi , j ,k H zi , j 1,k H y H y xxi , j ,k E xi , j , k 4 4 y z i , j ,k i , j ,k i , j 1, k i , j 1, k Ex yzi , j ,k Ezi , j ,k yzi , j , k 1 Ezi , j ,k 1 yzi , j 1,k 1 Ezi , j 1, k 1 yzi , j 1,k Ezi , j 1,k yxi 1, j 1, k Exi 1, j 1,k yxi 1, j ,k Exi 1, j ,k H xi , j ,k H xi , j , k 1 H zi , j ,k H zi 1, j ,k yx Ex yx yyi , j ,k E yi , j ,k 4 4 z x H yi , j ,k H yi 1, j ,k H xi , j ,k H xi , j 1,k zxi , j ,k Exi , j ,k zxi 1, j ,k E xi 1, j ,k zxi 1, j , k 1 Exi 1, j ,k 1 zxi , j ,k 1 Exi , j ,k 1 zyi , j , k E yi , j ,k zyi , j 1,k E yi , j 1,k zyi , j 1,k 1 E yi , j 1,k 1 zyi , j ,k 1 E yi , j ,k 1 zzi , j ,k Ezi , j ,k x y 4 4 We are forced to interpolate the problem terms so they exist at the same positions as the other terms in the finite‐difference equations. We interpolate the products E and H so that the field and material value being multiplied are at the same points. Lecture 10 Slide 55 Close Up of E H i , j , k 1 Ezi , j 1,k Ezi , j ,k E y y E yi , j ,k xxi , j ,k H xi , j ,k z i , j ,k H yi , j ,k xyi 1, j ,k H yi 1, j ,k xyi , j 1,k H yi , j 1,k xyi 1, j 1,k H yi 1, j 1,k xy 4 xzi , j ,k H zi , j ,k xzi , j ,k 1 H zi , j ,k 1 xzi 1, j ,k 1 H zi 1, j ,k 1 xzi 1, j ,k H zi 1, j ,k 4 i 1, j , k i 1, j , k i , j ,k i , j ,k Hx yxi , j 1,k H xi , j 1,k yxi 1, j 1, k H xi 1, j 1,k Exi , j ,k 1 Exi , j ,k Ezi 1, j ,k Ezi , j ,k yx H x yx z x 4 yyi , j ,k H yi , j ,k yzi , j ,k H zi , j ,k yzi , j ,k 1 H zi , j ,k 1 yzi , j 1,k 1 H zi , j 1,k 1 yzi , j 1,k H zi , j 1,k 4 E yi 1, j ,k E yi , j ,k Exi , j 1,k Exi , j ,k zxi , j ,k H xi , j ,k zxi 1, j ,k H xi 1, j ,k zxi 1, j ,k 1 H xi 1, j ,k 1 zxi , j ,k 1 H xi , j ,k 1 x y 4 zyi , j ,k H yi , j ,k zyi , j 1,k H yi , j 1,k zyi , j 1,k 1 H yi , j 1,k 1 zyi , j ,k 1 H yi , j ,k 1 4 i , j ,k H i , j ,k zz Lecture 10 z Slide 56 28 5/5/2015 Close Up of H E i , j ,k i , j , k 1 H zi , j ,k H zi , j 1,k H y H y xxi , j ,k Exi , j ,k y z xyi , j ,k E yi , j ,k xyi , j 1,k E yi , j 1,k xyi 1, j 1,k E yi 1, j 1,k xyi 1, j ,k E yi 1, j ,k 4 xzi , j ,k Ezi , j ,k xzi , j ,k 1 Ezi , j ,k 1 xzi 1, j ,k 1 Ezi 1, j ,k 1 xzi 1, j ,k Ezi 1, j ,k 4 i , j 1, k i , j 1, k i , j ,k i , j ,k yxi 1, j 1,k Exi 1, j 1,k yxi 1, j ,k Exi 1, j ,k Ex H xi , j ,k H xi , j ,k 1 H zi , j ,k H zi 1, j ,k yx Ex yx z x 4 yyi , j ,k E yi , j ,k yiz, j ,k Ezi , j ,k yzi , j ,k 1 Ezi , j ,k 1 yzi , j 1,k 1 Ezi , j 1,k 1 yzi , j 1,k Ezi , j 1,k 4 H yi , j ,k H yi 1, j ,k H xi , j ,k H xi , j 1,k zxi , j , k Exi , j ,k zxi 1, j ,k Exi 1, j ,k zxi 1, j ,k 1 Exi 1, j ,k 1 zxi , j ,k 1 Exi , j ,k 1 x y 4 zyi , j ,k E yi , j ,k zyi , j 1,k E yi , j 1,k zyi , j 1,k 1 E yi , j 1,k 1 zyi , j ,k 1 E yi , j ,k 1 4 zzi , j ,k Ezi , j ,k Lecture 10 Slide 57 Matrix Form i , j , k 1 E yi , j ,k xyi , j ,k H yi , j ,k xyi 1, j ,k H yi 1, j ,k xyi , j 1,k H yi , j 1,k xyi 1, j 1,k H yi 1, j 1,k xzi , j ,k H zi , j , k xzi , j , k 1 H zi , j ,k 1 xzi 1, j ,k 1 H zi 1, j ,k 1 xzi 1, j ,k H zi 1, j ,k Ezi , j 1,k Ezi , j , k E y xxi , j ,k H xi , j ,k y z 4 4 i, j,k i , j ,k i 1, j , k i 1, j , k Hx yzi , j ,k H zi , j , k yzi , j , k 1 H zi , j ,k 1 yzi , j 1,k 1 H zi , j 1,k 1 yzi , j 1,k H zi , j 1,k yxi , j 1,k H xi , j 1, k yxi 1, j 1,k H xi 1, j 1, k Exi , j ,k 1 Exi , j , k Ezi 1, j ,k Ezi , j , k yx H x yx yyi , j ,k H yi , j , k z x 4 4 E yi 1, j ,k E yi , j , k Exi , j 1,k Exi , j ,k zxi , j , k H xi , j ,k zxi 1, j ,k H xi 1, j , k zxi 1, j , k 1 H xi 1, j ,k 1 zxi , j ,k 1 H xi , j ,k 1 zyi , j ,k H yi , j ,k zyi , j 1,k H yi , j 1,k zyi , j 1,k 1 H yi , j 1,k 1 zyi , j ,k 1 H yi , j ,k 1 zzi , j ,k H zi , j ,k x y 4 4 Deye z Deze y μxx h x R x R y μxy h y R x R z μxz h z Deze x Dexe z R x R y μyxh x μyy h y R y R z μyz h z Dexe y Deye x R x R z μzxh x R y R z μzy h y μzz h z i , j ,k i , j , k 1 xyi , j ,k E yi , j , k xyi , j 1,k E yi , j 1, k xyi 1, j 1,k E yi 1, j 1,k xyi 1, j ,k E yi 1, j ,k xzi , j ,k Ezi , j ,k xzi , j ,k 1 Ezi , j ,k 1 xzi 1, j , k 1 Ezi 1, j ,k 1 xzi 1, j ,k Ezi 1, j ,k H zi , j , k H zi , j 1,k H y H y xxi , j ,k Exi , j , k y z 4 4 i , j ,k i , j ,k i , j 1, k i , j 1, k yxi 1, j 1,k Exi 1, j 1,k yxi 1, j ,k E xi 1, j ,k Ex yzi , j ,k Ezi , j ,k yzi , j ,k 1 Ezi , j ,k 1 yzi , j 1, k 1 Ezi , j 1,k 1 yzi , j 1, k Ezi , j 1,k H xi , j , k H xi , j ,k 1 H zi , j ,k H zi 1, j , k yx E x yx yyi , j ,k E yi , j ,k z x 4 4 H yi , j , k H yi 1, j ,k H xi , j ,k H xi , j 1, k zxi , j ,k E xi , j ,k zxi 1, j , k E xi 1, j ,k zxi 1, j ,k 1 Exi 1, j ,k 1 zxi , j ,k 1 Exi , j ,k 1 zyi , j , k E yi , j , k zyi , j 1,k E yi , j 1,k zyi , j 1,k 1 E yi , j 1, k 1 zyi , j ,k 1 E yi , j ,k 1 zzi , j ,k E zi , j ,k x y 4 4 Dhyh z Dhzh y εxxe x R x R y εxy e y R x R z εxz e z Dhzh x Dhxh z R x R y εyx e x εyy e y R y R z εyz e z Dhxh y Dhyh x R x R z εzxe x R y R z εzy e y εzz e z Lecture 10 Slide 58 29 5/5/2015 Block Matrix Form (1 of 2) Deye z Deze y μxx h x R x R y μxy h y R x R z μxz h z Deze x Dexe z R x R y μyxh x μyy h y R y R z μyz h z Dexe y Deye x R x R z μzxh x R y R z μzy h y μzz h z 0 e D z e D y Dez 0 Dex Dey e x μxx Dex e y R y R x μyx 0 e z R z R x μzx R x R y μxy μyy R z R y μzy R x R z μxz h x R y R z μyz h y μzz h z R x R y εxy εyy R z R y εzy R x R z εxz e x R y R z εyz e y εzz e z Dhyh z Dhzh y εxxe x R x R y εxy e y R x R z εxz e z Dhzh x Dhxh z R x R y εyxe x εyy e y R y R z εyz e z Dhxh y Dhyh x R x R z εzxe x R y R z εzy e y εzz e z 0 h D z Dhy Dhz 0 D Hx Dhy h x εxx Dhx h y R y R x εyx 0 h z R z R x εzx Lecture 10 Slide 59 Block Matrix Form (2 of 2) We can write our two block matrix equations as Ce e μr h h x h h y h z 0 H C Dhz Dhy e x e e y e z Dhz 0 Dhx εxx εr R y R x εyx R z R x εzx Lecture 10 Ch h εr e Dhy Dhx 0 R x R y εxy εyy R z R y εzy Notes: 1. We can handle anisotropic materials just by modifying the material matrices. 2. We do this by incorporating interpolation matrices. 3. Ch Ce H 0 E C Dez Dey R x R z εxz R y R z εyz εzz Dez 0 Dex μxx μr R y R x μyx R z R x μzx Dey Dex 0 R x R y μxy μyy R z R y μzy R x R z μxz R y R z μyz μzz Slide 60 30 5/5/2015 Interpolation Matrices (1 of 2) The derivative operators were constructed from a simple finite‐ difference approximation of the form f1.5 f 2 f1 x x The interpolation matrices are constructed exactly the same way, but uses the following equation for interpolation: f1.5 f 2 f1 2 The interpolation matrices have the following interpretations: R i interpolates along the i -axis using a value from the next cell along i R i interpolates along the i -axis using a value from the previous cell along i Lecture 10 Slide 61 Interpolation Matrices (2 of 2) Think of forming the interpolation this way R x x e D x 2 R y y e D y 2 R z z e D z 2 Strictly speaking, this will not work because the | | operation breaks some boundary conditions. This calculation approach does work for Dirichlet boundary conditions and for periodic boundary conditions that do not include phase. The interpolation matrices are related through: R x R x Lecture 10 H R y R y H R z R z H Slide 62 31 5/5/2015 Alternative Grids Lecture 10 Slide 63 Drawbacks of Uniform Grids Uniform grids are the easiest to implement, but do not conform well to arbitrary structures and exhibit high anisotropic dispersion. Anisotropic Dispersion (see Lecture 10) Lecture 10 Staircase Approximation (see Lecture 18) Slide 64 32 5/5/2015 Drawbacks of Structured Grids (1 of 2) Structured grids are the easiest to implement, but do not conform well to arbitrary geometries. Structured Grid Unstructured Grid Lecture 10 Slide 65 Hexagonal Grids Hexagonal grids are good for minimizing anisotropic dispersion suffered on Cartesian grids. This is very useful when extracting phase information. Phase Velocity as a Function of Propagation Angle x Yee‐FDTD 10 Hex‐FDTD 0° 57° 115° 172° 229° 286° 344° See Text, pp. 101‐103. Lecture 10 Slide 66 33 5/5/2015 Nonuniform Orthogonal Grids (1 of 2) Nonuniform orthogonal grids are still relatively simple to implement and provide some ability to refine the grid at localized regions. See Text, pp. 464‐471. Lecture 10 Slide 67 Nonuniform Orthogonal Grids (2 of 2) Uniform Grid Simulation • 80×110×16 cells • 140,800 cells Nonuniform Grid Simulation • 64×76×16 cells • 77,824 cells Conclusion: Roughly 50% memory and time savings. Lecture 10 Slide 68 34 5/5/2015 Curvilinear Coordinates Maxwell’s equations can be transformed from curvilinear coordinates to Cartesian coordinates to conform to curved boundaries of a device. See Text, pp. 484‐492. M. Fusco, “FDTD Algorithm in Curvilinear Coordinates,” IEEE Trans. Ant. and Prop., vol. 38, no. 1, pp. 76‐89, 1990. Lecture 10 Slide 69 Structured Nonorthogonal Grids This is a particularly powerful approach for simulating periodic structures with oblique symmetries. M. Fusco, “FDTD Algorithm in Curvilinear Coordinates,” IEEE Trans. Ant. and Prop., vol. 38, no. 1, pp. 76‐89, 1990. Lecture 10 Slide 70 35 5/5/2015 Irregular Nonorthogonal Unstructured Grids Unstructured grids are more tedious to implement, but can conform to highly complex shapes while maintaining good cell aspect ratios and global uniformity. Comparison of convergence rates ln x Lecture 10 P. Harms, J. Lee, R. Mittra, “A Study of the Nonorthogonal FDTD Method Versus the Conventional FDTD Technique for Computing Resonant Frequencies of Cylindrical Cavities,” IEEE Trans. Microwave Theory and Techniq., vol. 40, no. 4, pp. 741‐746 , 1992. Slide 71 Bodies of Revolution (Cylindrical Symmetry) Three‐dimensional devices with cylindrical symmetry can be very efficiently modeled using cylindrical coordinates. Devices with cylindrical symmetry have fields that are periodic around their axis. Therefore, the fields can be expanded into a Fourier series in . E , , e even , cos m eodd , sin m m 0 H , , h even , cos m hodd , sin m m 0 Due to a singularity at r=0, update equations for fields on the z axis are derived differently. See Text, Chapter 12 Lecture 10 Slide 72 36 5/5/2015 Some Devices with Cylindrical Symmetry Bent Waveguides Conical Horn Antenna Cylindrical Waveguides Dipole Antennas Diffractive Lenses Focusing Antennas Lecture 10 Slide 73 37
© Copyright 2024