1
ONE-DIMENSIONAL SIMULATION WITH THE FDTD METHOD

This chapter provides a step-by-step introduction to the finite-difference time-domain (FDTD) method, beginning with the simplest possible problem, the simulation of a pulse propagating in free space in one dimension. This example is used to illustrate the FDTD formulation. Subsequent sections lead to formulations for more complicated media.

1.1 ONE-DIMENSIONAL FREE-SPACE SIMULATION

The time-dependent Maxwell's curl equations for free space are

E and H are vectors in three dimensions, so, in general, Eq. (1.1a) and (1.1b) represent three equations each. We will start with a simple one-dimensional case using only Ex and Hy, so Eq. (1.1a) and (1.1b) become

These are the equations of a plane wave traveling in the z direction with the electric field oriented in the x direction and the magnetic field oriented in the y direction.

Taking the central difference approximations for both the temporal and spatial derivatives gives

In these two equations, time is specified by the superscripts, that is, n represents a time step, and the time t is t = ?t?·?n. Remember, we have to discretize everything for formulation into the computer. The term n?+?1 means one time step later. The terms in parentheses represent distance, that is, k is used to calculate the distance z = ?x?·?k. (It might seem more sensible to use ?z as the incremental step because in this case we are going in the z direction. However, ?x is so commonly used for a spatial increment that we will use ?x.) The formulation of Eq. (1.3a) and (1.3b) assume that the E and H fields are interleaved in both space and time. H uses the arguments k?+?1/2 and k?-?1/2 to indicate that the H field values are assumed to be located between the E field values. This is illustrated in Fig. 1.1. Similarly, the n?+ 1/2 or n?-?1/2 superscript indicates that it occurs slightly after or before n, respectively. Equations (1.3a) and (1.3b) can be rearranged in an iterative algorithm:

Notice that the calculations are interleaved in space and time. In Eq. (1.4a), for example, the new value of Ex is calculated from the previous value of Ex and the most recent values of Hy.

Figure 1.1 Interleaving of the E and H fields in space and time in the FDTD formulation. To calculate Hy, for instance, the neighboring values of Ex at k and k?+?1 are needed. Similarly, to calculate Ex, the values of Hy at k?+?1/2 and are needed.

This is the fundamental paradigm of the FDTD method (1).

Equations (1.4a) and (1.4b) are very similar, but because e0 and µ0 differ by several orders of magnitude, Ex and Hy will differ by several orders of magnitude. This is circumvented by making the following change of variables (2):

Substituting this into Eq. (1.4a) and (1.4b) gives

Once the cell size ?x is chosen, then the time step ?t is determined by

where c0 is the speed of light in free space. (The reason for this will be explained in Section 1.2.) Therefore, remembering that e0µ0 = 1/(c0)2,

Rewriting Eq. (1.6a) and (1.6b) in Python gives the following:

Note that the n, n?+?1/2, or n?-?1/2 in the superscripts is gone. Time is implicit in the FDTD method. In Eq. (1.9a), the ex on the right side of the equal sign is the previous value at n?-?1/2, and the ex on the left side is the new value n?+?1/2, which is being calculated. Position, however, is explicit. The only difference is that k?+?1/2 and k?-?1/2 are rounded to k and k?-?1 in order to specify a position in an array in the program.

The program fd1d_1_1.py at the end of this chapter is a simple one-dimensional FDTD program. It generates a Gaussian pulse in the center of the problem space, and the pulse propagates away in both directions as seen in Fig. 1.2. The Ex field is positive in both directions, but the Hy field is negative in the negative direction. The following points are worth noting about the program:

1. The Ex and Hy values are calculated by separate loops, and they employ the interleaving described above.
2. After the Ex values are calculated, the source is calculated. This is done by simply specifying a value of Ex at the point k = kc and overriding what was previously calculated. This is referred to as a hard source because a specific value is imposed on the FDTD grid.

Figure 1.2 FDTD simulation of a pulse in free space after 100 time steps. The pulse originated in the center and travels outward.

PROBLEM SET 1.1

1. Get the program fd1d_1_1.py running. What happens when the pulse hits the end of the array? Why?
2. Modify the program so it has two sources, one at kc - 20 and one at kc + 20. (Notice that kc is the center of the problem space.) What happens when the pulses meet? Explain this from basic electromagnetic (EM) theory.
3. Instead of Ex as the source, use Hy at k = kc as the source. What difference does it make? Try a two-point magnetic source at kc - 1 and kc such that hy[kc - 1] = - hy[kc]. What does this look like? To what does it correspond physically?

1.2 STABILITY AND THE FDTD METHOD

Let us return to the discussion of how to determine the time step. An EM wave propagating in free space cannot go faster than the speed of light. To propagate a distance of one cell requires a minimum time of ?t = ?x/c0. With a two-dimensional simulation, we must allow for the propagation in the diagonal direction, which brings the requirement to . Obviously, a three-dimensional simulation requires . This is summarized by the well-known Courant Condition (3, 4):

where n is the dimension of the simulation. Unless otherwise specified, throughout this book we will determine ?t by

This is not necessarily the best formula; we will use it for simplicity to avoid using square roots.

PROBLEM SET 1.2

1. In fd1d_1_1.py, go to the governing equations, Eq. (1.9a) and (1.9b), and change the factor 0.5 to 1.0. What happens? Change it to 1.1. Now what happens? Change it to 0.25 and see what happens.

1.3 THE ABSORBING BOUNDARY CONDITION IN ONE DIMENSION

Absorbing boundary conditions are necessary to keep outgoing E and H fields from being reflected back into the problem space. Normally, in calculating the E field, we need to know the surrounding H values. This is a fundamental assumption of the FDTD method. At the edge of the problem space we will not have the value of one side. However, we have an advantage because we know that the fields at the edge must be propagating outward. We will use this fact to estimate the value at the end by using the value next to it (5).

Suppose we are looking for a boundary condition at the end where k = 0. If a wave is going toward a boundary in free space, it is traveling at c0, the speed of light. So, in one time step of the FDTD algorithm, it travels

This equation shows that it takes two time steps for the field to cross one cell. A commonsense approach tells us that an acceptable boundary condition might be

The implementation is relatively easy. Simply store a value of Ex(1) for two time steps and then assign it to Ex(0). Boundary conditions such as these have been implemented at both ends of the Ex array in the program fd1d_1_2.py. Figure 1.3 shows the results of a simulation using fd1d_1_2.py. A pulse that originates in the center propagates outward and is absorbed without reflecting...