IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 48, NO. 2, FEBRUARY 2000 231 Numerical Stability of NonorthogonalFDTD Methods Stephen D. Gedney, Senior Member, IEEE, and J. Alan Roden, Member, IEEE Abstract—In this paper, a sufficient test for the numerical stability of generalized grid finite-difference time-domain (FDTD) schemes is presented. It is shown that the projection operators of such schemes must be symmetric positive definite. Without this property, such schemes can exhibit late-time instabilities. The origin and the characteristics of these late-time instabilities are also uncovered. Based on this study, nonorthogonal grid FDTD schemes (NFDTD) and the generalized Yee (GY) methods are proposed that are numerically stable in the late time for quadri- lateral prism elements, allowing these methods to be extended to problems requiring very long-time simulations. The study of numerical stability that is presented is very general and can be applied to most solutions of Maxwell’s equations based on explicit time-domain schemes. Index Terms—FDTD methods, numerical stability. I. INTRODUCTION T HE finite-difference time-domain (FDTD) method has been highly successful for the analysis of a plethora of electromagnetic interaction problems [1]. Yet one principal limitation of the classical FDTD method is the restriction to orthogonal grids. A number of techniques have been proposed to develop FDTD methods based on conformal meshing such as the contour path FDTD (CPFDTD) method [2], the nonorthogonal FDTD (NFDTD) method [3], and the discrete integral equation (DSI) [4] and generalized Yee (GY) methods [5]. The advantage of these techniques is that through the introduction of more generalized discretizations, error due to boundary discretization can be alleviated. Unfortunately, these and similar methods can sometimes suffer from late-time instabilities [6]. It was demonstrated by Craddock et al. that the source of in- stability of the CPFDTD method was due to the nonreciprocal nature of the original algorithm [7]. By analyzing the FDTD method as a passive circuit, an alternative stable and accurate solution was proposed. Such an extension is not directly appli- cable to the NFDTD, DSI, and GY methods because of the com- plexity of the projection operations required to project the fields normal to the grid faces onto the dual edges passing through the faces. In this paper, it is demonstrated that the source of late-time instabilities in the NFDTD and DSI/GY methods is due to the basic definitions for the projection operators. Consequently, Manuscript received March 25, 1998; revised November 15, 1999. This work was supported in part by the National Science Foundation under Grants ECS-9624628 and CDA-9502645 and the Army Research Office under Grant DAAH04-94-G-0243. S. D. Gedney is with the Department of Electrical Engineering, University of Kentucky, Lexington, KY 40506-0046 USA. J. A. Roden is with IBM Personal Systems Group, IBM Corporation, Re- search Triangle Park, NC 27709 USA. Publisher Item Identifier S 0018-926X(00)01650-1. these methods can lead to formulations that are unstable in the late time independent of the time step. Alternative methods based on the NFDTD and the DSI/GY algorithms are proposed that are accurate and stable when employing quadrilateral prism elements (i.e., elements that are orthogonal in the vertical di- rection and irregular quadrilaterals in the horizontal direction). II. GENERAL FORMULATION The FDTD methods described herein assume a dual-stag- gered grid with an edge-based discretization of Maxwell’s equa- tions. Specifically, the electric field intensity is projected onto the edges of a primary grid and the magnetic field intensity is projected onto the edges of the secondary grid. Implementing Maxwell’s curl equations in their integral form, the flux den- sities normal to the primary and secondary grid faces are nat- urally updated given the circulation of the dual field about the faces. Before performing the field update, the normal flux den- sity vector must be projected onto the dual edge passing through the face. This is accomplished in the NFDTD algorithm [3] via the local Jacobian tensor and local interpolation and a local in- terpolation scheme in the DSI and GY methods [5], [8]. The de- tails of these algorithms are not repeated here. Rather, the reader is referred to the original articles as well as [9] for a detailed summary of the algorithms. Both the NFDTD and the DSI/GY methods result in explicit update schemes that can be expressed in general form as a cou- pled set of first-order difference equations (1) (2) where and are vectors of the discrete vector flux densities, the superscripts refer to discrete time, and represent the dis- crete contour integrals of the electric and magnetic fields about primary and secondary cell faces, respectfully, is a diagonal matrix with entries representing the inverse of the relative per- mittivity, and and are the projection matrices. Note that for simplicity the domain is assumed to be lossless and nonmag- netic. However, the following analyses are not limited by these simplifications. III. STABILITY ANALYSIS The coupled difference equations in (1) and (2) are explicit in nature and are conditionally stable. To derive a sufficient sta- bility condition, the coupled difference equations are reposed as a single first-order difference equation. To this end, first substi- tute (1) into (2). Then, introducing the vector (3) 0018–926X/00$10.00 © 2000 IEEE