Why FDTD?

With the continued growth of computing power, modeling and numerical simulation has grown immensely as a tool for understanding and analyzing just about any problem in science. Where in the mid-twentieth century, detailed analyses were required to get any meaningful insight out of complex problems, today we can simply plug the governing differential equations into a computer, the results of which can provide an immense amount of information, which is of course complementary to theoretical analyses. The growth of computing power has brought with it a smorgasbord of modeling methods, applicable in any number of fields. The problem, then, is knowing when to use which method.

In electromagnetic problems, which are of interest to us in this book, there are quite a number of useful numerical methods, including the Method of Moments, Finite Volume methods, Finite Element methods, and Spectral methods, just to name a few. The FDTD method, however, grew to become the method of choice in the 1990s, for a number of reasons. First, it has always had the advantage of being a very simple method; we shall see in Chapter 3 that the derivation of difference equations is very straightforward. However, before the 1990s, the FDTD method was hindered by the need to discretize the simulation space on sub-wavelength scales, with time steps commensurately small. Hence, any reasonable problem would require a large amount of computer memory and time.

At this point, the reader should have all the necessary tools to create a useful FDTD model, in any number of dimensions, in almost any material. In the remainder of this book we cover a number of topics of general interest that have been found to be very useful in a number of FDTD applications. In this chapter, we will discuss modeling of periodic structures, modeling of physical features smaller than a grid cell size, the method known as Bodies of Revolution (BOR) for modeling cylindrical structures in 3D, and finally, the near-to-far field transformation, which is used to extrapolate the far-field radiation or scattering pattern from a confined FDTD simulation.

Modeling periodic structures

Quite a number of problems in numerical modeling involve structures that are periodic in one or more dimensions. Examples include photonic bandgap structures, which involve periodic arrays of dielectric structures; or arrays of antennas, in cases where the array is large enough that it can be analyzed as a periodic structure. These problems can be reduced to the modeling of a single period of the structure through the methods described in this section. Figure 12.1 shows a simple example of a periodic structure, consisting of rows of circular dielectric “balls” that repeat in the y-direction.

The main issue of complexity in modeling periodic structures, and which requires some ingenuity in the modeling effort, arises from the fact that future values of fields are involved in the FDTD update equations.

The von Neumann stability analysis [1, 2] discussed in the previous chapter is widely applicable and enables the assessment of the stability of any finite difference scheme in a relatively simple manner. However, such an analysis reveals little about the detailed properties of the difference scheme, and especially the important properties of dispersion and dissipation. These two metrics together yield information about the accuracy of the finite difference algorithm.

In the continuous world, dispersion refers to the variation of the phase velocity vp as a function of frequency or wavelength. Dispersion is present in all materials, although in many cases it can be neglected over a frequency band of interest. Numerical dispersion refers to dispersion that arises due to the discretization process, rather than the physical medium of interest. In addition, the discretization process can lead to anisotropy, where the phase velocity varies with propagation direction; this numerical anisotropy can also be separate from any real anisotropy of the medium.

In addition, as we saw in Figure 5.3, the finite difference scheme can lead to dissipation, which is nonphysical attenuation of the propagating wave. Recall that the example illustrated in Figure 5.3 was in a lossless physical medium, so any dissipation is due only to the discretization process.

In this chapter we will show how dispersion and dissipation arise in the finite difference algorithms that we have discussed so far, and show how to derive the numerical dispersion relation for any algorithm.

In Chapter 4 we presented the Yee algorithm for the discretization of Maxwell's equations on a regular, orthogonal grid. In Chapter 6 we discussed the numerical accuracy of the discretization, which depends on the choices of Δx, Δy, Δz, and Δt, as well as the medium of interest (∈, µ, and σ) and the frequencies to be simulated. However, the accuracy of the method is further restricted by how we choose to drive the simulation. The choice of a source is thus integral to accurate simulation in the FDTD method.

Sources can be either internal or external to the simulation. Internal sources are those that are driven at one or more grid cells within the simulation space. An example of an internal source might be a hard antenna, whose parameters we wish to investigate in a particular medium. External sources are those which are assumed to be outside of the simulation space. For example, we may wish to investigate the scattering pattern of an object when excited by a plane wave, but the source of that plane wave may be far from the scattering object. The problem, of course, is that in a numerical simulation, all sources must be applied to real grid cells, which means they are by definition “internal.” However, as we will see beginning in Section 7.2, methods have been derived with great success to make sources appear external.

In Chapter 3 we looked very briefly at the stability of numerical algorithms. By trial and error, we found the condition in Section 3.4.1 that, for a phase velocity of 1 ft ns−1, we required Δt ≤ Δx, or else the solution would rapidly go unstable as time progressed. This condition is known as the Courant-Friedrichs-Lewy or CFL stability condition. In this chapter we will introduce a commonly used method for assessing the stability of finite difference methods. We will also investigate the source of instability and more formally derive the CFL condition.

The most commonly used procedure for assessing the stability of a finite difference scheme is the so-called von Neumann method, initially developed (like many other finite difference schemes) for fluid dynamics related applications. This stability analysis, based on spatial Fourier modes, was first proposed and used by J. von Neumann during World War II at Los Alamos National Laboratory [1, 2].

The von Neumann method is applied by first writing the initial spatial distribution of the physical property of interest (e.g., the voltage or the electric field) as a complex Fourier series. Then, we seek to obtain the exact solution of the finite difference equation (FDE) for a general spatial Fourier component of this complex Fourier series representation. If the exact solution of the FDE for the general Fourier component is bounded (either under all conditions or subject to certain conditions on Δx, Δt), then the FDE is said to be stable.

Our purpose in this text is to provide an introduction to Numerical Electromagnetics, sometimes referred to as Computational Electromagnetics (CEM), a subject much too broad to cover in a single volume. Fifteen years ago, we might have found it difficult to choose between the different techniques to emphasize in our relatively brief coverage. However, due to a number of developments in the late 1980s and the 1990s, partial differential equation (PDE)-based methods, and in particular the so-called Finite Difference Time Domain (FDTD) method, have emerged as the methods with arguably the broadest range of applicability. This is especially true for electromagnetic problems involving complex and dispersive media, photonics applications, and modeling of high-speed circuits and devices. In addition, FDTD modeling of practical problems can now be undertaken with computer resources readily available to individual users. Finally, and quite importantly for our purposes, FDTD methods are relatively straightforward and intuitively follow from a physical understanding of Maxwell's equations, making this topic particularly suitable for both undergraduate and first-year graduate students in the context of a mezzanine-level course. Students with limited or no prior modeling experience will find that the FDTD method is the simplest and most insightful method from which to start their modeling education, and they can write practical and useful simulations in a matter of minutes. With an understanding of the FDTD method under their belts, students can move on to understanding more complicated methods with relative ease.

The FDTD method is widely applicable to many problems in science and engineering, and is gaining in popularity in recent years due in no small part to the fast increases in computer memory and speed. However, there are many situations when the frequency domain method is useful.

The finite difference frequency domain (FDFD) method utilizes the frequency domain representation of Maxwell's equations. In practice, however, rather than Maxwell's equations, the frequency domain wave equation is more often used, as it is a more compact form, requires only a single field, and does not require interleaving. In this chapter we will focus on the latter method, and give a brief overview of the Maxwell's equations formulation toward the end of the chapter. Because the method uses the frequency domain equations, the results yield only a single-frequency, steady-state solution. This is often the most desirable solution to many problems, as many engineering problems involve quasi–steady state fields; and, if a single-frequency solution is sought, the FDFD is much faster than the transient broadband analysis afforded by FDTD.

FDFD via the wave equation

The simplest way to set up an FDFD simulation is via the wave equation. This method is often most desirable as well; since it does not involve coupled equations, it does not require interleaving or a Yee cell, and all fields (E, H, and J) are co-located on the grid.

After a brief exposure to different finite difference algorithms and methods, we now focus our attention on the so-called FDTD algorithm, or alternatively the Yee algorithm [1], for time domain solutions of Maxwell's equations. In this algorithm, the continuous derivatives in space and time are approximated by second-order accurate, two-point centered difference forms; a staggered spatial mesh is used for interleaved placement of the electric and magnetic fields; and leapfrog integration in time is used to update the fields. This yields an algorithm very similar to the interleaved leapfrog described in Section 3.5.

The cell locations are defined so that the grid lines pass through the electric field components and coincide with their vector directions, as shown in Figure 4.1. As a practical note, the choice here of the electric field rather than the magnetic field is somewhat arbitrary. However, in practice, boundary conditions imposed on the electric field are more commonly encountered than those for the magnetic field, so that placing the mesh boundaries so that they pass through the electric field vectors is more advantageous. We will have more to say about the locations of field components later in this chapter. Note that the Yee cell depicted in [2, Fig. 3.1] has the cell boundaries to be aligned with the magnetic field components, rather than with the electric field components as in Figure 4.1.

The methods described in the previous chapter for analytically treating outgoing waves and preventing reflections can be very effective under certain circumstances. However, a number of drawbacks are evident. The reflection coefficient that results from these analytical ABCs is a function of incident angle, and can be very high for grazing angles, as shown in Figure 8.4. Furthermore, as can be seen from the Higdon operators, for higher accuracy the one-way wave equation that must be applied at the boundary becomes increasingly complex, and requires the storage of fields at previous time steps (i.e., n − 1). One of the great advantages of the FDTD method is that it does not require the storage of any fields more than one time step back, so one would prefer a boundary method that also utilizes this advantageous property.

Those methods in the previous chapter should most accurately be referred to as “radiation” boundary conditions. This type of boundary emulates a one-way wave equation at the boundary, as a method for circumventing the need for field values outside the boundary, which would be required in the normal update equation. Strictly speaking, however, they are not “absorbing” boundary conditions.

The family of methods that will be discussed in this chapter are absorbing boundaries. These methods involve modifying the medium of the simulation in a thin layer around the boundary, as shown in Figure 9.1, so that this layer becomes an artifically “absorbing” or lossy medium.

Until now, we have discussed the standard FDTD algorithm, based on the interleaved second-order centered difference evaluation of the time and space derivatives, as originally introduced by Yee in 1966 [1]. This method is an explicit method, in which the finite difference approximations of the partial derivatives in the PDE are evaluated at the time level (n + 1/2), so that the solution at the next time level (n + 1) can be expressed explicitly in terms of the known values of the quantities at time index n.

Explicit finite difference methods have many advantages. For the hyperbolic wave equation, an important advantage is the fact that numerical solutions obtained with the explicit method exhibit finite numerical propagation speed, closely matching the physical propagation speed inherent in the actual PDE. However, an important disadvantage of explicit methods is the fact that they are only conditionally stable. As a result, the allowable time step (vis-à-vis the CFL condition) is generally quite small, making certain types of problems prohibitive simply due to the enormous computation times required. In such cases, which include some important electromagnetic applications as discussed in [2], it is highly desirable to avoid the time step limitation.

Implicit finite difference methods provide the means with which this time step limitation can be avoided. With an implicit method, the finite difference approximations of the partial derivatives in the PDE are evaluated at the as-yet unknown time level (n + 1).

A study of Numerical Electromagnetics must rely on a firm base of knowledge in the foundations of electromagnetics as stated in Maxwell's equations. Accordingly, we undertake in this chapter a review of Maxwell's equations and associated boundary conditions.

All classical electromagnetic phenomena are governed by a compact and elegant set of fundamental rules known as Maxwell's equations. This set of four coupled partial differential equations was put forth as the complete classical theory of electromagnetics in a series of brilliant papers written by James Clerk Maxwell between 1856 and 1865, culminating in his classic paper [2]. In this work, Maxwell provided a mathematical framework for Faraday's primarily experimental results, clearly elucidated the different behavior of conductors and insulators under the influence of fields, imagined and introduced the concept of displacement current [3, Sec. 7.4], and inferred the electromagnetic nature of light. A most fundamental prediction of this theoretical framework is the existence of electromagnetic waves, a conclusion to which Maxwell arrived in the absence of experimental evidence that such waves can exist and propagate through empty space. His bold hypotheses were to be confirmed 23 years later (in 1887) in the experiments of Heinrich Hertz [4].

When most of classical physics was fundamentally revised as a result of Einstein's introduction [6] of the special theory of relativity, Maxwell's equations remained intact. To this day, they stand as the most general mathematical statements of fundamental natural laws which govern all of classical electrodynamics.