FDTD update across inhomogeneous media

could anyone guide me as to how to formulate electric field update equations at the bounday between media with different electrical properties?

The easiest answer is, don't worry about it. If you place your fields on a Yee grid, they will implicitly satisfy Maxwell's divergence equations. You only have to worry about updating the fields based on the curl equations.

If you are updating your Ex field, for example, you will perform one update for every cell in your problem grid. Use the epsilon value that exists at the same point as the Ex field you are currently updating. The same can be said for Ey, Ez, Hx, Hy, and Hz, except the H values are associated with a relative permeability. Does this make sense?

For the most part, FDTD works pretty well if the material interface is normal to the x-, y-, or z-axis. If you want to represent curved surfaces there are some other tricks like smoothing or averaging the material properties at the interface. If you are including PEC or metallic structures, you probably want to place tangential fields along the interface, but included in the metal. This stuff is pretty well documented in the literature.

You may want to read Chapter 3 and 4 from my dissertation. Chapter 3 describes a frequency-domain method very similar to FDTD. It also uses a Yee grid and I describe some techniques for representing devices on the Yee grid in section 3.3, pp. 61-64. All of this applies to the time-domain as well. I devote all of Chapter 4 to FDTD. You can download my dissertation here:

https://www.edaboard.com/download.php?id=124081

Hope this helps! Good luck!
-Tip

Thanks very much for the reply. I have not read your dissertation yet, so i might be asking a question already addressed in it. I'll ask anyway...
I am constructing my grid such that the E fields are present on the edges of a cell and the H fields are within the cell. So i should have no problem updating the H fields since the field being updated is WITHIN the cell.
But when updating the E fields, since say Ez could lie on an edge which could be on a boundary between two different media, we could have on value of H used in the Ez update equation lying in medium1 and the other value in medium2, for example. what value of epsilon do i use in this case?

When you use any type o staggered grid, you have to realize the field components are at physically different points in space. You have to resist the temptation to think that all field components in the same cell are seeing the same material. Since they are in physically different locations it is certainly possible, and most likely, that they will see different materials despite being in the same cell. Also, when you implement sources you will need to take this into account as the different field components will be slightly out of phase.

I build my materials on a Yee grid this way. First, I construct what I call the "2X" grid. That is, I create a an array in memory that has twice the grid resolution of the Yee grid I want to implement. Second, I build my device on the 2X grid as normal without considering anything about a staggered grid. Third, I interpolate, or parse, the material values onto my actual grid from the 2X grid.

When you do this, you will notice that each of your field components is assigned its own unique material property. It is even possible, and likely, that electric fields within the same grid cell will be assigned different material properties.

I have found this approach to work very well. Does this make sense?

-Tip

rrumpf, thanks a lot.
it does make sense but would it be possible to provide some code descibing what you mean. I tend to be able to understand programming code better than english language sometimes (lol).

Unfortunately the issue of tangential fields in interface boundaries still remains a elusive to me. since they are tangential to the boundary between two media, in your material property assignment, what value of material property will you assign this field? the value corresponding to medium1 or that corresponding to medium2? thats the crust of what am trying to ask really.

Thanks for your time and hoping for a reply!

Here is the implementation for a two-dimensional grid in MATLAB:

% DEFINE MATERIALS
er1 = 1.0; %relative permittivity outside circle
er2 = 10.0; %relative permittivity of circle
r = 0.75; %radius of circle

% DEFINE GRID
Nx = 20; dx = 0.1;
Ny = 20; dy = 0.1;

% COMPUTE 2X GRID
Nx2 = 2*Nx; dx2 = dx/2;
Ny2 = 2*Ny; dy2 = dy/2;

% COMPUTE 1X GRID AXES
xa = [0:Nx-1]*dx; xa = xa - mean(xa); %x-axis
ya = [0:Ny-1]*dy; ya = ya - mean(ya); %y-axis

% COMPUTE 2X GRID AXES
xa2 = [0:Nx2-1]*dx2; xa2 = xa2 - mean(xa2); %x-axis
ya2 = [0:Ny2-1]*dy2; ya2 = ya2 - mean(ya2); %y-axis

% CREATE CIRCLE ON 2X GRID
[Y,X] = meshgrid(ya2,xa2); %meshgrid
R = (X.^2 + Y.^2) <= r^2; %binary circle
ER2 = er1*(~R) + er2*(R); %ER on 2x grid

% PARSE MATERIALS INTO YEE GRID
ERx = ER2(2:2:Nx2,1:2:Ny2); %ER for Ex fields
ERy = ER2(1:2:Nx2,2:2:Ny2); %ER for Ey fields
ERz = ER2(1:2:Nx2,1:2:Ny2); %ER for Ez fields

% SHOW RESULTS
subplot(232); imagesc(xa2,ya2,ER2');
colorbar; axis equal tight;
title('\epsilon_r on 2X grid');

subplot(234); imagesc(xa,ya,ERx');
colorbar; axis equal tight;
title('\epsilon_r for E_x on 1X grid');

subplot(235); imagesc(xa,ya,ERy');
colorbar; axis equal tight;
title('\epsilon_r for E_y on 1X grid');

subplot(236); imagesc(xa,ya,ERz');
colorbar; axis equal tight;
title('\epsilon_r for E_z on 1X grid');

The code computes ERx, ERy, and ERz. These are 2D arrays that contain the permittivity that you would associate with the Ex, Ey, and Ez fields respectively.

-Tip

I think the simplest method is to average the properties.
In more detail, you average the properties in a cube which has the same size as a Yee-cell
but the field location you want to update exactly in its center.
This works well if the material interfaces all run along the boundaries of the Yee-cells.

There are more complicated methods for the more general case but often the above methods works well.

I'll add that averaging works great for large curved surfaces. I have found it to work poorly for resolving very fine features. For example, if you have to resolve a dielectric structure with size on the order of one grid cell or smaller, averaging will tend to "erase" those structures leading to inaccurate simulations. In this case, the method I summarized above seems to work the best. One could argue that you should use finer grid resolution and averaging in this case. That is true unless you run out of memory doing so. You could also move to an unstructured or non-uniform grid, but that is a bit of work and tends to be rare in FDTD.

Good luck!
-Tip