1. No category

5/5/2015
Instructor
Dr. Raymond Rumpf
(915) 747‐6958
rcrumpf@utep.edu
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
222 Grid
44 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 1E5  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 1E5  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 H1  x H0.5 

  
1 0 0 0 H 2   x H1.5 
1 1 0 0 H3   x H 2.5 

  
0 1 1 0 H 4  x H3.5 
0 0 1 1 H5  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 H1  x H0.5 

  
1 0 0 0 H 2   x H1.5 
1 1 0 0 H3   x H 2.5 

  
0 1 1 0 H 4  x H3.5 
0 0 1 1 H5  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 
 1x  1x 0
 0 1 1
0
0  E2 
x
x

Dexe   0
0  1x  1x 0  E3 

 
0
0  1x  1x E4 
 0
 0
0
0
0  1xE5 
0
0  E1 
 1x  1x 0
 0 1  1
0
0  E2 
x
x

Dexe   0
0  1x  1x 0  E3 

 
0
0  1x  1x E4 
 0
 1x 0
0
0  1x E5 
Slide 32
16
5/5/2015
Derivative Operators on a 33 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 100200 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 44 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 . D1 D1
D 2
Dx1 Dx1 
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
Dx2
 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
Deye z  Deze y  μ xxh x
Deze x  Dexe z  μ yy h y
Dexe y  Deye x  μ zz h z
Dhyh z  Dhzh y  ε xxe x
Dhzh x  Dhxh z  ε yy e y
Dhxh y  Dhyh 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
Dexe y  Deye x  μ zz h z
Dhyh z  Dhzh y  ε xx e x
Dhzh x  Dhxh z  ε yy e y
Dhxh y  Dhyh 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
 xyi , j ,k H yi , j ,k   xyi 1, j ,k H yi 1, j ,k   xyi , j 1,k H yi , j 1, k   xyi 1, j 1,k H yi 1, j 1,k  xzi , j ,k H zi , j ,k   xzi , j , k 1 H zi , j ,k 1   xzi 1, j , k 1 H zi 1, j , k 1   xzi 1, j ,k H zi 1, j ,k
E zi , j 1,k  Ezi , j ,k E y

  xxi , 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
 yzi , j ,k H zi , j ,k   yzi , j , k 1 H zi , j ,k 1   yzi , j 1, k 1 H zi , j 1, k 1   yzi , j 1,k H zi , j 1,k
  yxi , 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  zxi , j ,k H xi , j , k   zxi 1, j , k H xi 1, j ,k   zxi 1, j , k 1 H xi 1, j , k 1   zxi , j ,k 1 H xi , j ,k 1  zyi , j , k H yi , j ,k   zyi , j 1,k H yi , j 1,k   zyi , j 1,k 1 H yi , j 1,k 1   zyi , j ,k 1 H yi , j ,k 1

  zzi , j ,k H zi , j , k


4
4
x 
y 
i , j ,k
i , j , k 1
 xyi , j ,k E yi , j , k   xyi , j 1,k E yi , j 1,k   xyi 1, j 1,k E yi 1, j 1,k   xyi 1, j ,k E yi 1, j ,k  xzi , j ,k Ezi , j ,k   xzi , j ,k 1 Ezi , j ,k 1   xzi 1, j ,k 1 Ezi 1, j , k 1   xzi 1, j ,k Ezi 1, j , k
H zi , j ,k  H zi , j 1,k H y  H y

  xxi , 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
 yzi , j ,k Ezi , j ,k   yzi , j , k 1 Ezi , j ,k 1   yzi , j 1,k 1 Ezi , j 1, k 1   yzi , j 1,k Ezi , j 1,k
  yxi 1, j 1, k Exi 1, j 1,k   yxi 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


  yyi , 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  zxi , j ,k Exi , j ,k   zxi 1, j ,k E xi 1, j ,k   zxi 1, j , k 1 Exi 1, j ,k 1   zxi , j ,k 1 Exi , j ,k 1  zyi , j , k E yi , j ,k   zyi , j 1,k E yi , j 1,k   zyi , j 1,k 1 E yi , j 1,k 1   zyi , j ,k 1 E yi , j ,k 1



  zzi , 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
  xxi , j ,k H xi , j ,k
z 
 i , j ,k H yi , j ,k   xyi 1, j ,k H yi 1, j ,k   xyi , j 1,k H yi , j 1,k   xyi 1, j 1,k H yi 1, j 1,k
 xy
4
 xzi , j ,k H zi , j ,k   xzi , j ,k 1 H zi , j ,k 1   xzi 1, j ,k 1 H zi 1, j ,k 1   xzi 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  zxi , j ,k H xi , j ,k   zxi 1, j ,k H xi 1, j ,k   zxi 1, j ,k 1 H xi 1, j ,k 1   zxi , j ,k 1 H xi , j ,k 1


x
y
4
 zyi , j ,k H yi , j ,k   zyi , j 1,k H yi , j 1,k   zyi , j 1,k 1 H yi , j 1,k 1   zyi , 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

  xxi , j ,k Exi , j ,k
y
z 


 xyi , j ,k E yi , j ,k   xyi , j 1,k E yi , j 1,k   xyi 1, j 1,k E yi 1, j 1,k   xyi 1, j ,k E yi 1, j ,k
4
 xzi , j ,k Ezi , j ,k   xzi , j ,k 1 Ezi , j ,k 1   xzi 1, j ,k 1 Ezi 1, j ,k 1   xzi 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  zxi , j , k Exi , j ,k   zxi 1, j ,k Exi 1, j ,k   zxi 1, j ,k 1 Exi 1, j ,k 1   zxi , j ,k 1 Exi , j ,k 1


x
y 
4

 zyi , j ,k E yi , j ,k   zyi , j 1,k E yi , j 1,k   zyi , j 1,k 1 E yi , j 1,k 1   zyi , j ,k 1 E yi , j ,k 1
4
  zzi , j ,k Ezi , j ,k
Lecture 10
Slide 57
Matrix Form
i , j , k 1
 E yi , j ,k
 xyi , j ,k H yi , j ,k   xyi 1, j ,k H yi 1, j ,k   xyi , j 1,k H yi , j 1,k   xyi 1, j 1,k H yi 1, j 1,k  xzi , j ,k H zi , j , k   xzi , j , k 1 H zi , j ,k 1   xzi 1, j ,k 1 H zi 1, j ,k 1   xzi 1, j ,k H zi 1, j ,k
Ezi , j 1,k  Ezi , j , k E y

  xxi , 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
 yzi , j ,k H zi , j , k   yzi , j , k 1 H zi , j ,k 1   yzi , j 1,k 1 H zi , j 1,k 1   yzi , j 1,k H zi , j 1,k
  yxi , j 1,k H xi , j 1, k   yxi 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


  yyi , 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  zxi , j , k H xi , j ,k   zxi 1, j ,k H xi 1, j , k   zxi 1, j , k 1 H xi 1, j ,k 1   zxi , j ,k 1 H xi , j ,k 1  zyi , j ,k H yi , j ,k   zyi , j 1,k H yi , j 1,k   zyi , j 1,k 1 H yi , j 1,k 1   zyi , j ,k 1 H yi , j ,k 1



  zzi , j ,k H zi , j ,k
x 
y 
4
4
Deye z  Deze y  μxx h x  R x R y μxy h y  R x R z μxz h z
Deze x  Dexe z  R x R y μyxh x  μyy h y  R y R z μyz h z
Dexe y  Deye x  R x R z μzxh x  R y R z μzy h y  μzz h z
i , j ,k
i , j , k 1
 xyi , j ,k E yi , j , k   xyi , j 1,k E yi , j 1, k   xyi 1, j 1,k E yi 1, j 1,k   xyi 1, j ,k E yi 1, j ,k  xzi , j ,k Ezi , j ,k   xzi , j ,k 1 Ezi , j ,k 1   xzi 1, j , k 1 Ezi 1, j ,k 1   xzi 1, j ,k Ezi 1, j ,k
H zi , j , k  H zi , j 1,k H y  H y

  xxi , j ,k Exi , j , k 

y 
z 
4
4
i , j ,k i , j ,k
i , j 1, k i , j 1, k
  yxi 1, j 1,k Exi 1, j 1,k   yxi 1, j ,k E xi 1, j ,k
Ex
 yzi , j ,k Ezi , j ,k   yzi , j ,k 1 Ezi , j ,k 1   yzi , j 1, k 1 Ezi , j 1,k 1   yzi , 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


  yyi , 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  zxi , j ,k E xi , j ,k   zxi 1, j , k E xi 1, j ,k   zxi 1, j ,k 1 Exi 1, j ,k 1   zxi , j ,k 1 Exi , j ,k 1  zyi , j , k E yi , j , k   zyi , j 1,k E yi , j 1,k   zyi , j 1,k 1 E yi , j 1, k 1   zyi , j ,k 1 E yi , j ,k 1



  zzi , j ,k E zi , j ,k
x 
y 
4
4
Dhyh z  Dhzh y  εxxe x  R x R y εxy e y  R x R z εxz e z
Dhzh x  Dhxh z  R x R y εyx e x  εyy e y  R y R z εyz e z
Dhxh y  Dhyh 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)
Deye z  Deze y  μxx h x  R x R y μxy h y  R x R z μxz h z
Deze x  Dexe z  R x R y μyxh x  μyy h y  R y R z μyz h z
Dexe y  Deye 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 
Dhyh z  Dhzh y  εxxe x  R x R y εxy e y  R x R z εxz e z
Dhzh x  Dhxh z  R x R y εyxe x  εyy e y  R y R z εyz e z
Dhxh y  Dhyh 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