Modeling sources in the FDTD formulation and their use in quantifying source and boundary condition errors

The modeling of voltage and current sources as either added or replaced sources in FDTD simulations is described and their differences discussed in terms of a transmission line analogy. An infinitesimal current element (ICE) is used to illustrate the validation of added source modeling and to study the errors involved with modeling an infinitesimal element within the finite-sized FDTD grid. This model is also used to illustrate the behavior of radiation boundary conditions as their near-field position with respect to the source is varied. We characterize the errors due to modeling and boundary conditions and give guidelines for obtaining acceptable accuracy in simulations.

810 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 43, NO.4, APRIL 1995 Modeling Sources In the FDTD Formulation and Their Use In Quantifying Source and Boundary Condition Errors Dale N. Buechler, Member, IEEE, Daniel H. Roper, Member, IEEE, Carl H. Durney, Fellow, IEEE, and Douglas A. Christensen, Member, IEEE Abstract-The modeling of voltage and current sources as ei­ther added or replaced sources in FDTD simulations is described and their differences discussed in terms of a transmission line analogy. An infinitesimal current element (ICE) is used to il­lustrate the validation of added source modeling and to study the errors involved with modeling an infinitesimal element within the finite-sized FDTD grid. This model is also used to illustrate the behavior of radiation boundary conditions as their near-field position with respect to the source is varied. We characterize the errors due to modeling and boundary conditions and give guidelines for obtaining acceptable accuracy in simulations. I. INTRODUCTION THE FINITE-DIFFERENCE time-domain (FDTD) method for solving Maxwell's equations [I] has been widely utilized in the analysis of scattering phenomena [2], [3], radiation patterns from antennas [4]-[7] and biomedical applications such as hyperthermia [8]-[11]. For all FDTD applications, proper modeling of sources is essential. However, other than for plane-wave sources [2], [3], FDTD source modeling is not well documented, particularly for current sources. Voltage sources are typically modeled in FDTD formula­tions by either of two methods: 1) replacing the calculated electric field E on a Yee-cell edge by the source E at every time step ("replaced source"), or 2) adding the source E to the FDTD calculated E ("added source"). The replaced source appears to be more commonly used, for example in exciting coax [4] and waveguide structures [6], but knowing which kind of source to use is important because their effects on the system are quite different. For example, replaced sources may cause reflections of waves propagating back to the source location, Manuscript received October 8, 1992; revised July 14, 1994. This work was supported by a grant from the NIH (POI CA 29578) and by a grant of computer time from the Utah Supercomputing Institute, which is funded by the State of Utah and the IBM Corporation. D. N. Buechler is with the Department of Electrical Engineering and the Division of Radiation Oncology, University of Utah, Salt Lake City, UT 84112 USA. D. H. Roper is with Hughes Aircraft Co., Los Angeles, CA 90009 USA. C. H. Durney is with the Department of Electrical Engineering, and the Department of Bioengineering, University of Utah, Salt Lake City, UT 84112 USA. D. A. Christensen is with the Department of Electrical Engineering, Department of Bioengineering, University of Utah, Salt Lake City, UT 84112 USA. IEEE Log Number 9408562. while added sources can be transparent to these incoming waves. The purpose of this paper is to describe both voltage and current added and replaced sources in the FDTD formula­tion, to illustrate transparent source modeling using an added infinitesimal current element (which has a known analytical solution), to determine the extent of the region near this source where errors arise due to the finite grid size, and to illustrate the accuracy of radiation boundary conditions as a function of distance from this nonplane-wave source. The infinitesimal current element (ICE) source is a cur­rent sinusoidally varying in time along a directed line of infinitesimal length. Besides having a known solution for validating code, the ICE is a useful source in its own right. It can be employed, for example, either as a "building block" member of a group of several elements that are weighted with individual amplitudes and phases for modeling a general distributed current source [12], or as a single element by itself for such simulations as antenna feeds [6], radiating dipoles or monopoles [6], [12], [13], or optical emission from fluorescent molecules. The ICE is equivalent to an infinitesimally short dipole of oscillating charge, or Hertzian electric dipole [13], which has been widely used in numerical methods other than FDTD for calculating radiation patterns from such structures as microstrip patch antennas [14], [15] and dipole antennas in complex environments [16]. Errors near the source arise, however, when an infinitesimal current element is modeled within the finite-sized FDTD grid. In later sections, we characterize these errors and give guidelines for obtaining acceptable accuracy in simulations. Absorbing boundaries must be carefully located to keep the model as small as possible to maximize efficiency while maintaining acceptable accuracy. Mur [17] reported a two­dimensional (2-D) study of the performance of radiation boundary conditions (RBC) with nonplanar incident fields from an added isotropic source. We extend that study to 3- D using quantitative comparisons to the known analytical solutions of the ICE and determine the conditions under which the reflections from the boundaries are within acceptable limits. II. ADDED AND REPLACED FDTD SOURCES As explained in the previous section, FDTD sources can be added or replaced sources. The electric field excitation 00 18-9480/95$04.00 © 1995 IEEE BUECHLER et al.: MODELING SOURCES IN THE FDTD FORMULATION FDTD Lumped Transmission Line Analogy Source Cell (a) (b) Fig. I (a) The replacement of the FDTD electric field E by the source field (Es) is analogous to placing a voltage source Vs across a lossless transmission line. Similarly, replacing a current density J by the source current density (J s) is analogous to placing a current source 18 in series with the capacitor in a model of a lossless transmission line. Incoming waves will be reflected by such sources. (b) The addition of the calculated electric field E,. to the source field (Eo<) is analogous to placing a voltage source V s in series with the capacitor in a model of a lossless transmission line. Similarly, adding the calculated current density J c to the source current density (J s) is analogous to placing a current source C across a lossless transmission line. These sources will appear transparent to incoming waves. of a lossless lumped-element transmission line illustrates the different nature of these two kinds of sources (Fig. I). The replacement of the FDTD electric field by the source field (Es) is similar to placing a voltage source (Vs = Esb, where 8 is the cell size) across a capacitance in the transmission line (Fig. I (a)). The FDTD replaced current source is similar to a current source (Is) in series with a capacitance of the transmission line (Fig. lea)). With ideal sources (no internal resistance), these voltage and current sources when deactivated will appear as short and open circuits, respectively, to any incoming waves, thus causing reflections. On the other hand, adding the calculated FDTD electric field (Ec) to the source field (Es) is similar to placing a voltage source (Vs) in series with a capacitance of the transmission line, as shown in Fig. I(b). The FDTD added current source is similar to a current source (Is) in parallel with a capacitance of the the lossless transmission line. These sources, which result in a series combination of an ideal voltage source for the lossless transmission line (Vc + Vs) or the parallel combination of an ideal current source (Ie + Is), will appear transparent (when deactivated) to incoming waves. Modeling FDTD voltage sources is straightforward because the electric field E appears explicitly in the standard FDTD equations and how to add to or replace E with a source E is obvious. Modeling FDTD current sources is perhaps not so obvious because current is usually not explicitly included in the FDTD equations. Where should the current source be placed in the Yee cell, and how should the current source be included in the FDTD equations? We found that the best way 811 8x (a) I Idl liz )- - - 8x ~ Vee Cell (b) Fig. 2. (a) The FDTD representation of the integral form of Maxwell's 'I;" x H equation is shown. The current enclosed by the loop of H fields is equal to the current density J s multiplied by the area f,.l'~y. (b) The FDTD average current density J s is used to represent the current element Idl averaged over the FDTD volume f>.l'tJ y/l ~. to model an added current source in the FDTD formulation is to locate it on the edge of a Yee cell (Fig. 2(a)) and to add it to the current density J in Maxwell's V' x H equation. Starting with the integral form of the V' x H equation, we set J = (Y E + J s, where J s is the z-directed source current density averaged over the entire FDTD source cell. The familiar FDTD equation with one additional source term, C28J s , is then obtained by integrating around the path shown in Fig. 2(a) E~+l(i,j, k) = C1 E:'Ci,j, k) + C2 [H;+O.5(i,j - 1, k) - H~t+05(i,j, k) + H~'+O.5(i,j, k) where - H;+O.5(i - l.j, k) + bJ s], (1) C1 = [~t - ~] / [~t + ~] C2 = [c8[~t + ;c]]-l 8 = 8x = by = bz. Other FDTD expressions are unaffected by the current source. III. THE INFINITESIMAL CURRENT SOURCE MODEL The ICE is a good test case for validating the current-source model because of its readily available analytic solution. First the quantity Idl used in the ICE analytical solution must be related to the source current density J s' The expression for 812 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 43, NO.4, APRIL 1995 the vector potential A generated by an infinitesimal current element is [13] A - JL J J dV' _ JL Idl - 471" --y - 471" -;:-. (2) where the integration is over the volume containing the source, R is the distance from the point of the field A to each source point, and r is the distance from the field point to the center of the source volume, where the current element is located (Fig. 2(b». To model this current source using the FDTD method, we employ a z-directed current density J s centered along the edge of one Yee cell, as shown in Fig. 2(b). In the FDTD formulation, where quantities are assumed uniform within a cell dimension, (3) where t5x, t5y, and t5z are the cell dimensions, and R :::::: r for R » J(t5x)2 + (t5y)2 + (t5z)2. Using (3) in (2), the FDTD source current density J s is related to the infinitesimal current element Idl by Js=~. t5xt5yt5z (4) It can be expected that this representation will be limited in accuracy near the source since the approximation in (3), a consequence of the finite-sized nature of FDTD cells, is less valid as the source is approached. The distance at which the approximation leads to significant errors is investigated in the next section. IV. ERRORS NEAR THE SOURCE AND BOUNDARIES To study errors related to this source model, and also to the proximity of the absorbing boundaries, the ICE source was placed at the center of the free-space region shown in Fig. 3 and employed as an added source in (1). The region was subdivided into cubic Yee cells of. the same size as the cell containing the source. Absorbing boundary conditions (2nd­order Mur [17] with special comer boundary conditions [3]) were used on the outside boundaries of the region. Points for comparison to analytical results were selected along three paths, each beginning at the source and terminating at the face, edge, and comer boundaries respectively. For each point of comparison, the magnitudes of the FDTD electric and magnetic field components were plotted with respect to time to ensure that the steady state had been reached. Because FDTD calculations give transient solutions, the tum-on characteristic of the source strongly affects the FDTD results and sometimes produces dc offsets and/or transients that delay the reaching of the steady-state response. We tested different time functions for the source waveform and found that using r( t) sin wt, where r(t) = 0.5[1 - cos(wtj2a)] is a raised cosine envelope for t between 0 and aT [18], produced the minimum amount of transients and no dc offset. Therefore we used a raised cosine envelope waveform for all simulations. An IBM 3090-6000 computer was used for all FDTD calculations and an HP 9000-850 computer was employed for data display and analysis. A typical run time to reach steady z k:x -----. Face v: Source typical cell Fig. 3. The ICE source is placed at the center of a free-space region subdivided into Yee cells, each having volume /)3. Points for comparison to the analytical results are chosen along three paths, each beginning at the source cell and terminating at the face, the edge, and the comer. Cell Size: 0.0083)" 0.01SH 0.0333)" SO 0.10)" t 40 N W .5 20 g w .~. ~ -20 "c i :::I 0 "ii -40 III II: CD u -60 {f 0.0 0.1 0.2 0.3 0.4 0.5 Distance from Source (Wavelengths) Fig. 4. The relative error in Ez along the face path vs. the distance from the source in wavelengths is shown. Plots are given for four different cell sizes. The overall size of the region is held constant; thus the number of cells in the region varies accordingly. For clarity of presentation, not every cell point is plotted. state for our vectorized code with a 60 x 60 x 60 model was 286 cpu seconds. The analytical results [13] were compared to those obtained with the FDTD simulation as a function of distance from the source normalized to wavelength. Percent relative error was defined as 100 [(Eanalytical - Esimulation) jEanalyticad, where the electric field magnitudes are steady-state peak-to­peak values. The relative errors in Ez and Ex (equal to Ey by symmetry) are plotted in Figs. 4-6 as a function of the distance from the source along the three paths. Fig. 4 shows the relative error in Ez as a function of distance along the face path for four different cell sizes but for a constant overall model size of approximately one wavelength on each side. Fig. 5 shows relative errors in Ex and Ez as a function of position along the face and comer paths for a fixed cell size of 0.0 083'\. In Fig. 6, relative errors in Ex are shown along the face path for four different boundary positions (i.e., for four different model sizes); cubic models with (20)3, (30)3, (60)3, and (120)3 cells of the same cell size (0.0167'\) were utilized for this comparison. BUECHLER et al.: MODELING SOURCES IN THE FDTD FORMULATION ~ !:! Iii Q) > ~ Gi II: 40 20 0 ·20 ·40 ·60 ·80 0.0 ---- EX' Face Path ~ Ez• Face Path --0-- EXI Corner Path -<r-- Ez• Corner Path uaQs) 0.2 0.4 0.6 t E 0 0 0.8 Distance from Source (Wavelengths) Fig. 5. The relative errors in Ez and Ex are shown along the face and corner paths vs. the distance from the source in wavelengths for a cell size of 0.0 083-\. <fi >w< ..5. !:! Iii Q) > i Gi II: 60 30 0 ·30 ·60 ·90 ·120 -----<>- EX' Face Path (20)3model ------ EX' Face Path (30)3model -~ EX' Face Path (60)3model --- EX' Face Path (120)3 model .----,-r-----.------------~M o N ~ 0.0 0.2 0.4 0.6 0.8 1.0 Distance from Source (Wavelengths) Fig. 6. The relative error in Ex is shown along the face path vs. the distance from the source in wavelengths for four locations of the boundaries. The plots are for a (20):3 cell model. a (30):3 cell model, a (GO)3 cell model, and a (120)3 cell model, each with a cell size of 0.0167-\. The face boundary location for each case is indicated. V. DISCUSSION OF RESULTS For the simulations shown in Fig. 4, the boundaries were kept at a fixed distance (0.5>') from the source as the cell size was changed. This was done to determine the ICE source model accuracy near the source as a function of cell size, and to study reflections from the RBC's. Note that errors in Ez near the face boundary are small regardless of cell size. Note, however, that cell size significantly affects the accuracy of the results near the source. As the cell size decreases, the region of appreciable error around the source shrinks, but even for small cell size, accuracy immediately adjacent to the source is still poor. This is understandable since the FDTD method calculates only the spatially averaged field values for each Yee cell, and the rapid variation of the fields near the infinitesimal current element cannot be adequately modeled with finite cell dimensions. As the cell size is made smaller, this region of high error shifts closer to the source. For convenience of 813 comparison, we defined a ten-percent relative error point. This point can be seen to occur at about six cells from the source for cell sizes of 0.0 167>' and 0.0 083>' (slightly closer for b 2 0.0 333>.). We have found that as cell size is decreased even more (not shown), this ten-percent error point will shift closer to the source in terms of wavelengths but will never get closer than approximately six cells. In Fig. 5, the error near the source is larger for the face path than for the comer path; this is understandable since the analytical result for an ICE source shows more rapid spatial variation in the fields along the face path than along the comer path near the source. At the face boundary (normal to the x­direction), the Ex field component has a larger relative error than the Ez component. This is not surprising since the Mur boundary conditions, which were developed to minimize the reflections of normally incident waves, minimize the error of Ez at the face boundary. At the comer, the order is reversed, with Ez having the larger relative error. The edge path errors (not shown) for Ex and Ez were each found to be a few percent below the respective extremes of Ex at the face and Ez at the comer. Fig. 6 shows an increase in the relative error along the face path as the boundaries are placed closer to the source in terms of wavelengths. It demonstrates that the plane-wave behavior assumed in the RBC's is increasingly invalid for boundaries located in the near field closer than about 0.5>. to the source. High boundary errors for the closer boundary locations (i.e., 0.167>. and 0.250>') were also seen along the edge and comer paths (results not shown). We found that placing the boundaries further than 0.5>. from the source reduces these errors at the boundaries. Fig. 6 also shows that the results immediately adjacent to the source are not changed significantly by the closer proximity of the boundary. The lack of interaction between the boundary location and the source is indicative of a transparent (added) source. Otherwise, it is expected that the reflections caused by the boundaries would have had a more noticable effect on the results immediately adjacent to the source. VI. SUMMARY AND CONCLUSION We found that the ICE model is valid in FDTD simulations at distances of about six cells or more from the source, independent of cell size (assuming, of course, the conventional limit of b < AI 10). Closer in than six cells, however, the rapidly varying fields immediately adjacent to the source were difficult to accurately model. This is apparently due to the fact that the FDTD method inherently models a current density which is uniformly distributed throughout a source cell volume, while the analytical expressions are for a precisely located infinitesimal element, leading among other things to a difference between the magnitudes of Rand r in (3). Reduc­tion of the cell size shifts the error curve physically closer to the source, but our studies indicate that the point of ten-percent error is not closer than approximately six cells from the source for small cell sizes. Within this limitation, the ICE source can be used with confidence in FDTD applications. For general current sources where the current is distributed over finite 814 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES, VOL. 43, NO.4, APRIL 1995 dimensions, the above limitation may be eased, but that part of the error which results from the FDTD method's inability to model rapidly varying fields over distances comparable to the cell size will remain. The errors near the boundary substantially increased when the boundaries were placed closer than about one-half wave­length from the infinitesimal current element source. This increased error is consistent with the findings in the 2-D study by Mur [12] using an isotropic source. Placement of the boundaries further than this one-half wavelength point results in a significant reduction of the errors at the boundaries. For some models, however, it may be inconvenient or impossible to place the boundaries one-half wavelength away due to computer memory limitations. This applies to large-area or low-frequency simulations at high resolution, e.g., treatment­planning models for regional heating hyperthermia devices [9], [19]. It is therefore important that work continues on the development of improved boundary conditions and new algorithms which incorporate variable cell sizes [20] and shapes [21], [22] to allow more efficient use of computer resources. REFERENCES [I] K. S. Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. Antenn. Propagat., vol. AP-14, pp. 302-307, 1966. [2] A. Taflove and M. E. Brodwin, "Numerical solution of steady-state elec­tromagnetic scattering problems using the time-dependent Maxwell's equations," IEEE Trans. Microwave Theory Tech., vol. MTT-23 , pp. 623-630, 1975. [3] A. Taflove and K. R. Umashankar, "The finite-difference method for numerical modeling of electromagnetic wave interactions with arbitrary structures," M. A. Morgan, Ed., in PIER 2 Finite Element and Finite Difference Methods in Electromagnetic Scattering. Amsterdam, The Netherlands: Elsevier, ch. 8, 1990. [4] R. Luebbers, L. Chen, T. Uno, and S. Adachi, "FDTD calculation of radiation patterns, impedance, and gain for a monopole antenna on a conducting box," IEEE Trans. Antenn. Propagat., vol. AP-40, pp. 1577-1583, 1992. [5] A. Reineix and B. Jecko, "Analysis of microstrip patch antennas using finite difference time domain method," IEEE Trans. Antenn. Propagat., vol. AP-37, pp. 1361-1369, 1989. [6] P. A. Tirkas and C. A. Balanis, "Finite-difference time-domain method for antenna radiation," IEEE Trans. Antenn. Propagat., vol. AP-40, pp. 334-340, 1992. [7] J. G. Maloney, G. S. Smith, and W. R. Scott, Jr., "Accurate computation of the radiation from simple antennas using the finite-difference time­domain method," IEEE Trans. Antenn. Propagat., vol. AP-38, pp. 1059-1067, 1990. [8] J. A. Shaw, C. H. Durney, and D. A. Christensen, "Computer-aided design of two-dimensional electric-type hyperthermia applicators using the finite-difference time-domain method," IEEE Trans. Biomed. Eng., vol. BME-38, pp. 861-870, 1991. [9] D. Sullivan, "Three-dimensional computer simulation in deep regional hyperthermia using the finite-difference time-domain method," IEEE Trans. Microwave Theory Tech., vol. MTT-38, pp. 204-211, 1990. [10] M. J. Piket-May, A. Taflove, W. C. Lin, D. S. Katz, V. Satiaseelan, and B. B. Mittal, "Initial Results for automated computational modeling of patient-specific electromagnetic hyperthermia," IEEE Trans. Biomed. Eng., vol. BME-39, pp. 226-237, 1992. [II] J. Y. Chen and O. P. Gandhi, "Numerical simulation of annular-phased arrays of dipoles for hyperthermia of deep-seated tumors," IEEE Trans. Biomed. Eng., vol. BME-39, pp. 209-216, 1992. [12] S. Ramo and J. R. Whinnery, Fields and Waves in Modern Radio, 2nd ed. New York: Wiley, 1953, p. 496. [13] E. C. Jordan and K. G. Balmain, Electromagnetic Waves and Radiating Systems. Englewood Cliffs, NJ: Prentice-Hall, 1968. [14] D. R. Jackson and N. G. Alexopoulos, "Simple approximate formulas for input resistance, bandwidth, and efficiency of a resonant rectangular patch," IEEE Trans. Antenn. Propagat., vol. AP-39, pp. 407-410, 1991. [15] E. Tonye, J. B. Temdenou, and J. Ntonme, "Computer simulation of multimicrostrip patches antenna problem," lEE Proc., pI. H, vol. 139, pp. 59-64. [16] I. Y. Hsia and N. G. Alexopoulos, "Radiation characteristics of hertzian dipole antennas in a nonreciprocal superstrate-substrate structure," IEEE Trans. Antenna. Propagat., vol. AP-40, pp. 782-790, 1992. [17] G. Mur, "Absorbing boundary conditions for the finite-difference ap­proximation of the time-domain electromagnetic field equations," IEEE Trans. Electromagn. Compat., vol. EMC-23, pp. 377-382, 1981. [18] D. H. Roper and J. M. Baird, "Analysis of overmoded waveguides using the finite-difference time domain method," in IEEE MTT-S Int. Microwave Symp. Dig.,1992, pp. 401-404. [19] D. M. Sullivan, D. Buechler, and F. A. Gibbs, "Comparison of measured and simulated data in an annular phased array using an inhomogeneous phantom," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 600-604, 1992. [20] S. S. Zivanovic, K. S. Yee, and K. K. Mei, "A subgridding method for the time-domain finite-difference method to solve Maxwell's equations," IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 471-479, 1991. [21] M. Fusco, "FDTD algorithm in curvilinear coordinates," IEEE Trans. Antenn. Propagat., vol. AP-38, pp. 76-89, 1990. [22] J.-F. Lee, R. Palandech, and R. Mittra, "Modeling three-dimensional discontinuities in waveguides using nonorthogonal FDTD algorithm," IEEE Trans. Microwave Theory Tech., vol. MTT-40, pp. 346-352, 1992. Dale N. Buechler (S'83-M'86) was born in In­glewood, CA, on May II, 1962. He received the B.S. and M.S. degrees in electrical engineering from the University of Arizona in 1984 and 1986 respectively, including completion of the Clinical Engineering Option. While working on his Master's Degree he was involved in research and clinical work in hyper­thermia and continued this work as a research associate in the Division of Radiation Oncology at the University of Wisconsin-Madison (1987-1988) and at the University of Utah Division of Radiation Oncology (1988-1991). He continued to work there on research in Clinical Radiation Therapy until 1993. He is currently working on the Ph.D. degree in electrical engineering at the University of Utah, where his research interests include biomedical engineering and numerical techniques in electromagnetics. Daniel H. Roper (S'91-M'92-S'92-M'93) was born in Sierra Madre, CA, on Dec. 7, 1966. He received the B.S. and M.S. degrees in electrical engineering from the University of Utah in 1990 and 1992, respectively. From 1990 to 1993 he worked as a research assistant in the Department of Electrical Engineering on numerical analysis of microwave components and on mode measurement systems for high-power circular waveguides. Since 1993 he has been with the Corporate Engineering Rotation Program of Hughes Aircraft Co., Los Angeles, CA. His current research interests are in active array antennas and microwave integrated circuits. Carl H. Durney (S'60-M'64-SM'80-F'92), photograph and biography not available at the time of publicaiton. Douglas A. Christensen (M'69), photograph and biography not available at the time of publicaiton.