Electromagnetic Simulation Using the FDTD Method with Python. Dennis M. Sullivan

by caption."/>

Figure 1.1 Interleaving of the E and H fields in space and time in the FDTD formulation. To calculate H_y, for instance, the neighboring values of E_x at k and k + 1 are needed. Similarly, to calculate E_x, the values of H_y 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 ε₀ and μ₀ differ by several orders of magnitude, E_x and H_y will differ by several orders of magnitude. This is circumvented by making the following change of variables (2):

(1.5)

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

(1.6a)

(1.6b)

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

(1.7)

      where c₀ is the speed of light in free space. (The reason for this will be explained in Section 1.2.) Therefore, remembering that ε₀μ₀ = 1/(c₀)²,

(1.8)

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

(1.9a)

(1.9b)

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 E_x field is positive in both directions, but the H_y 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/c₀. 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):

      (1.10)

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

      (1.11)

      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