Finite-Difference Time-Domain Methods

Publisher Summary

This chapter discusses finite-difference time-domain (FDTD) methods. During the 1970s and 1980s, a number of researchers realized the limitations of frequency-domain integral-equation solutions of Maxwell's equations. This led to early explorations of a novel alternative approach: direct time-domain solutions of Maxwell's differential (curl) equations on spatial grids or lattices. FDTD and related space-grid time-domain techniques are direct solution methods for Maxwell's curl equations. These methods employ no potentials. Rather, they are based upon volumetric sampling of the unknown electric and magnetic fields within and surrounding the structure of interest, and over a period of time. For simulations where the modeled region must extend to infinity, absorbing boundary conditions (ABCs) are employed at the outer lattice truncation planes, which ideally permit all outgoing wave analogs to exit the region with negligible reflection. Phenomena such as induction of surface currents, scattering and multiple scattering, aperture penetration, and cavity excitation are modeled time-step by time-step by the action of the numerical analog to the curl equations.

Introduction

Prior to about 1990, the modeling of electromagnetic engineering systems was primarily implemented using solution techniques for the sinusoidal steady-state Maxwell's equations. Before about 1960, the principal approaches in this area involved closed-form and infinite-series analytical solutions, with numerical results from these analyses obtained using mechanical calculators. After 1960, the increasing availability of programmable electronic digital computers permitted such frequency-domain approaches to rise markedly in sophistication. Researchers were able to take advantage of the capabilities afforded by powerful new high-level programming languages such as Fortran, rapid random-access storage of large arrays of numbers, and computational speeds orders of magnitude faster than possible with mechanical calculators. In this period, the principal computational approaches for Maxwell's equations included the high-frequency asymptotic methods of Keller, 1962, Kouyoumjian and Pathak, 1974 and the integral-equation techniques of Harrington [1968].

However, these frequency-domain techniques have difficulties and trades-off. For example, while asymptotic analyses are well suited for modeling the scattering properties of electrically large complex shapes, such analyses have difficulty treating nonmetallic material composition and volumetric complexity of a structure. While integral equation methods can deal with material and structural complexity, their need to construct and solve systems of linear equations limits the electrical size of possible models, especially those requiring detailed treatment of geometric details within a volume, as opposed to just the surface shape.

While significant progress has been made in solving the ultra-large systems of equations generated by frequency-domain integral equations (see, for example, Song and Chew [1998]), the capabilities of even the latest such technologies are exhausted by many volumetrically complex structures of engineering interest. This also holds for frequency-domain finite-element techniques, which generate sparse rather than dense matrices. Further, the very difficult incorporation of material and device nonlinearities into frequency-domain solutions of Maxwell's equations poses a significant problem as engineers seek to design active electromagnetic/electronic and electromagnetic/quantum-optical systems such as high-speed digital circuits, microwave and millimeter-wave amplifiers, and lasers.

During the 1970s and 1980s, a number of researchers realized the limitations of frequency-domain integral-equation solutions of Maxwell's equations. This led to early explorations of a novel alternative approach: direct time-domain solutions of Maxwell's differential (curl) equations on spatial grids or lattices. The finite-difference time-domain (FDTD) method, introduced by Yee [1966], was the first technique in this class, and has remained the subject of continuous development (see Taflove and Hagness [2000]).

There are seven primary reasons for the expansion of interest in FDTD and related computational solution approaches for Maxwell's equations:

• FDTD uses no linear algebra. Being a fully explicit computation, FDTD avoids the difficulties with linear algebra that limit the size of frequency-domain integral-equation and finite-element electromagnetics models to generally fewer than 10⁶ field unknowns. FDTD models with as many as 10⁹ field unknowns have been run. There is no intrinsic upper bound to this number.

• FDTD is accurate and robust. The sources of error in FDTD calculations are well understood and can be bounded to permit accurate models for a very large variety of electromagnetic wave interaction problems.

• FDTD treats impulsive behavior naturally. Being a time-domain technique, FDTD directly calculates the impulse response of an electromagnetic system. Therefore, a single FDTD simulation can provide either ultrawideband temporal waveforms or the sinusoidal steady-state response at any frequency within the excitation spectrum.

• FDTD treats nonlinear behavior naturally. Being a time-domain technique, FDTD directly calculates the nonlinear response of an electromagnetic system.

• FDTD is a systematic approach. With FDTD, specifying a new structure to be modeled is reduced to a problem of mesh generation rather than the potentially complex reformulation of an integral equation. For example, FDTD requires no calculation of structure-dependent Green's functions.

• Computer memory capacities are increasing rapidly. While this trend positively influences all numerical techniques, it is of particular advantage to FDTD methods which are founded on discretizing space over a volume, and therefore inherently require a large random access memory.

• Computer visualization capabilities are increasing rapidly. While this trend positively influences all numerical techniques, it is of particular advantage to FDTD methods which generate time-marched arrays of field quantities suitable for use in color videos to illustrate the field dynamics.

An indication of the expanding level of interest in FDTD Maxwell's equations solvers is the hundreds of papers currently published in this area worldwide each year, as opposed to fewer than ten as recently as 1985 (see Shlager and Schneider [1998]). This expansion continues as engineers and scientists in non-traditional electromagnetics-related areas such as digital systems and integrated optics become aware of the power of such direct solution techniques for Maxwell's equations.

FDTD and related space-grid time-domain techniques are direct solution methods for Maxwell's curl equations. These methods employ no potentials. Rather, they are based upon volumetric sampling of the unknown electric and magnetic fields within and surrounding the structure of interest, and over a period of time. The sampling in space is at sub-wavelength resolution set by the user to properly sample the highest near-field spatial frequencies thought to be important in the physics of the problem. Typically, 10–20 samples per wavelength are needed. The sampling in time is selected to ensure numerical stability of the algorithm.

Overall, FDTD and related techniques are marching-in-time procedures that simulate the continuous actual electromagnetic waves in a finite spatial region by sampled-data numerical analogs propagating in a computer data space. Time-stepping continues as the numerical wave analogs propagate in the space lattice to causally connect the physics of the modeled region. For simulations where the modeled region must extend to infinity, absorbing boundary conditions (ABCs) are employed at the outer lattice truncation planes which ideally permit all outgoing wave analogs to exit the region with negligible reflection. Phenomena such as induction of surface currents, scattering and multiple scattering, aperture penetration, and cavity excitation are modeled time-step by time-step by the action of the numerical analog to the curl equations. Self-consistency of these modeled phenomena is generally assured if their spatial and temporal variations are well resolved by the space and time sampling process. In fact, the goal is to provide a self-consistent model of the mutual coupling of all of the electrically small volume cells constituting the structure and its near field, even if the structure spans tens of wavelengths in three dimensions and there are hundreds of millions of space cells.

Time-stepping is continued until the desired late-time pulse response is observed at the field points of interest. For linear wave interaction problems, the sinusoidal response at these field points can be obtained over a wide band of frequencies by discrete Fourier transformation of the computed field-versus-time waveforms at these points. Prolonged “ringing” of the computed field waveforms due to a high Q-factor or large electrical size of the structure being modeled requires a combination of extending the computational window in time and extrapolation of the windowed data before Fourier transformation.

Current FDTD and related space-grid time-domain algorithms are fully explicit solvers employing highly vectorizable and parallel schemes for time-marching the six components of the electric and magnetic field vectors at each of the space cells. The explicit nature of the solvers is usually maintained by employing a leapfrog time-stepping scheme. Current methods differ primarily in how the space lattice is set up. In fact, gridding methods can be categorized according to the degree of structure or regularity in the mesh cells:

• Almost completely structured. In this case, the space lattice is organized so that its unit cells are congruent wherever possible. The most basic example of such a mesh is the pioneering work of Yee [1966], who employed a uniform Cartesian grid having rectangular cells. Staircasing was used to approximate the surface of structural features not parallel to the grid coordinate axes. Later work showed that it is possible to modify the size and shape of the space cells located immediately adjacent to a structural feature to conformally fit its surface (see, for example, Jurgens, Taflove, Umashankar and Moore, 1992, Dey and Mittra, 1997). This is accurate and computationally efficient for large structures because the number of modified cells is proportional to the surface area of the structure. Thus, the number of modified cells becomes progressively smaller relative to the number of regular cells filling the structure volume as its size increases. As a result, the computer resources needed to implement a fully conformal model approximate those required for a staircased model. However, a key disadvantage of this technique is that special mesh-generation software must be constructed.

• Surface-fitted. In this case, the space lattice is globally distorted to fit the shape of the structure of interest. The lattice can be divided into multiple zones to accommodate a set of distinct surface features (see, for example, Shankar, Mohammadian and Hall [1990]). The major advantage of this approach is that well-developed mesh-generation software of this type is available. The major disadvantage is that, relative to the Yee algorithm, there is substantial added computer burden due to:

  • memory allocations for the position and stretching factors of each cell;

  • extra computer operations to implement Maxwell's equations at each cell and to enforce field continuity at the interfaces of adjacent cells.

  Another disadvantage is the possible presence of numerical dissipation in the time-stepping algorithm used for such meshes. This can limit the range of electrical size of the structure being modeled due to numerical wave-attenuation artifacts.

• Completely unstructured. In this case, the space containing the structure of interest is completely filled with a collection of lattice cells of varying sizes and shapes, but conforming to the structure surface (see, for example, Madsen and Ziolkowski [1990]). As for the case of surface-fitted lattices, mesh-generation software is available and capable of modeling complicated three-dimensional shapes possibly having volumetric inhomogeneities. A key disadvantage of this approach is its potential for numerical inaccuracy and instability due to the unwanted generation of highly skewed space cells at random points within the lattice. A second disadvantage is the difficulty in mapping the unstructured mesh computations onto the architecture of either parallel vector computers or massively parallel machines. The structure-specific irregularity of the mesh mandates a robust pre-processing algorithm that optimally assigns specific mesh cells to specific processors.

At present, the best choice of computational algorithm and mesh remains unclear. For the next several years, we expect continued progress in this area as various groups develop their favored approaches and perform validations.

For computational modeling of electromagnetic wave interaction structures using FDTD and related space-grid time-domain techniques, it is useful to consider the concept of predictive dynamic range. Let the power density of the primary (incident) wave in the space grid be P₀ W/m². Further, let the minimum observable power density of a secondary (scattered) wave be P_S W/m², where “minimum observable” means that the accuracy of the field computation degrades due to numerical artifacts to poorer than n dB (some desired figure of merit) at lower levels than P_S. Then, we can define the predictive dynamic range as 10log₁₀(P₀/P_S) dB.

This definition is well suited for FDTD and other space-grid time-domain codes for two reasons:

• It squares nicely with the concept of a “quiet zone” in an experimental anechoic chamber, which is intuitive to most electromagnetics engineers;

• It succinctly quantifies the fact that the desired numerical wave analogs propagating in the lattice exist in an additive noise environment due to nonphysical propagating wave analogs caused by the imperfect ABCs.

In addition to additive noise, the desired physical wave analogs undergo gradual progressive deterioration while propagating due to accumulating numerical dispersion artifacts, including phase velocity anisotropies and inhomogeneities within the mesh.

In the 1980s, researchers accumulated solid evidence for a predictive dynamic range on the order of 40–50 dB for FDTD codes. This value is reasonable if one considers the additive noise due to imperfect ABCs to be the primary limiting factor, since the analytical ABCs of this era (see, for example, Mur [1981]) provided outer-boundary reflection coefficients in the range of about 0.3–3% (−30 to −50 dB).

The 1990s saw the emergence of powerful, entirely new classes of ABCs including the perfectly matched layer (PML) of Berenger [1994]; the uniaxial anisotropic PML (UPML) of Sacks, Kingsland, Lee and Lee, 1995, Gedney, 1996; and the complementary operator methods (COM) of Ramahi, 1997, Ramahi, 1998. These ABCs were shown to have effective outer-boundary reflection coefficients of better than −80 dB for impinging pulsed electromagnetic waves having ultrawideband spectra. Solid capabilities were demonstrated to terminate free-space lattices, multimoding and dispersive waveguiding structures, and lossy and dispersive materials.

However, for electrically large problems, the overall dynamic range may not reach the maximum permitted by these new ABCs because of inaccuracies due to accumulating numerical-dispersion artifacts generated by the basic grid-based solution of the curl equations. Fortunately, by the end of the 1990s, this problem was being attacked by a new generation of low-dispersion algorithms. Examples include the wavelet-based multi-resolution time-domain (MRTD) technique introduced by Krumpholz and Katehi [1996] and the pseudo-spectral time-domain (PSTD) technique introduced by Liu, Q.H., 1996, Liu, Q.H., 1997. As a result of these advances, there is emerging the possibility of FDTD and related space-grid time-domain methods demonstrating predictive dynamic ranges of 80 dB or more in the first decade of the 21st century.

Using FDTD and related methods, we can model electromagnetic wave interaction problems requiring the solution of considerably more than 10⁸ field-vector unknowns. At this level of complexity, it is possible to develop detailed, three-dimensional models of complete engineering systems, including the following:

• Entire aircraft and missiles illuminated by radar at 1 GHz and above;

• Entire multilayer circuit boards and multichip modules for digital signal propagation, crosstalk, and radiation;

• Entire microwave and millimeter-wave amplifiers, including the active and passive circuit components and packaging;

• Entire integrated-optical structures, including lasers, waveguides, couplers, and resonators.

A key goal for such large models is to achieve algorithm/computer-architecture scaling such that for N field unknowns to be solved on M processors, we approach an order(N/M) scaling of the required computational resources.

We now consider the factors involved in determining the computational burden for the class of FDTD and related space-grid time-domain solvers.

• Number of volumetric grid cells, N. The six vector electromagnetic field components located at each lattice cell must be updated at every time step. This yields by itself an order(N) scaling.

• Number of time steps, n_max. A self-consistent solution in the time domain mandates that the numerical wave analogs propagate over time scales sufficient to causally connect each portion of the structure of interest. Therefore, n_max must increase as the maximum electrical size of the structure. In three dimensions, it can be argued that n_max is a fractional power function of N such as N^1/3. Further, n_max must be adequate to step through “ring-up” and “ring-down” times of energy storage features such as cavities. These features vary from problem to problem and cannot be ascribed a dependence relative to N.

• Cumulative propagation errors. Additional computational burdens may arise due to the need for either progressive mesh refinement or progressively higher-accuracy algorithms to bound cumulative positional or phase errors for propagating numerical modes in progressively enlarged meshes. Any need for progressive mesh refinement would feed back to factor 1.

For most free-space problems, factors 2 and 3 are weaker functions of the size of the modeled structure than factor 1. This is because geometrical features at increasing electrical distances from each other become decoupled due to radiative losses by the electromagnetic waves propagating between these features. Further, it can be shown that replacing second-order accurate algorithms by higher-order versions sufficiently reduces numerical dispersion error to avoid the need for progressive mesh refinement for object sizes up to the order of 100 wavelengths. Overall, a computational burden of order(N·n_max)=order(N^4/3) is estimated for very large FDTD and related models.