See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/3010070

Time-domain finite-element methods

Article in IEEE Transactions on Antennas and Propagation · April 1997

DOI: 10.1109/8.558658 · Source: IEEE Xplore

CITATIONS READS
431 2,015

3 authors, including:

Jin-Fa Lee Andreas C Cangellaris

The Ohio State University University of Illinois, Urbana-Champaign
316 PUBLICATIONS 9,469 CITATIONS 127 PUBLICATIONS 1,748 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Jin-Fa Lee on 19 February 2013.

The user has requested enhancement of the downloaded file.

430 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

Time-Domain Finite-Element Methods

Jin-Fa Lee, Member, IEEE, Robert Lee, Member, IEEE, and Andreas Cangellaris, Member, IEEE

Invited Review Paper

Abstract— In this paper, various time-domain finite-element at using 3-D mesh generator packages. Also, it is not always
methods for the simulation of transient electromagnetic wave phe- possible to get a high-quality mesh for a given geometry.
nomena are discussed. Detailed descriptions of test/trial spaces, Another difficulty is the relative complexity of TDFEM’s
explicit and implicit formulations, nodal and edge/facet element
basis functions are given, along with the numerical stability prop- compared to the standard FDTD method. The time needed
erties of the different methods. The advantages and disadvantages to learn TDFEM’s and write a usable 3-D code is an order of
of mass lumping are examined. Finally, the various formula- magnitude greater than that for the standard FDTD method.
tions are compared on the basis of their numerical dispersion From the above discussion, one is led to the conclusion
performance. that the standard FDTD method is the method of choice when
Index Terms— Electromagnetic transient analysis, finite- modeling geometries of low complexity, while TDFEM’s are
element methods. most appropriate when complicated geometries need to be
modeled for which application of the FDTD method would
I. INTRODUCTION require a very dense grid to resolve accurately the variation
in geometric features. Since TDFEM’s have not received as

N UMEROUS time-domain finite-element methods (TD-

FEM’s) have been proposed for electromagnetic model-
ing in the past ten years; however, the popularity of TDFEM’s
much attention as the FDTD method, they are lacking in both
maturity and variety of applications. Therefore, a comparison
of TDFEM’s with the FDTD method on the basis of current
is very small compared to that for the finite-difference time- capablities would be misleading. Instead, the goal of this paper
domain (FDTD) method. is to present TDFEM’s in a general framework which unifies
There are several reasons why the FDTD method is so pop- much of the previous work done in this area. The merits of the
ular. The most important ones are its programming simplicity, various TDFEM’s are discussed, and some future directions
its minimum bookkeeping complexity, and the simplicity of are proposed.
its numerical integration algorithm. Because of the simplicity
of the method, a novice can develop a three-dimensional (3-D)
FDTD code in a relatively short amount of time.
TDFEM offers some important advantages over the standard II. FAEDO–GALERKIN FORMULATIONS
FDTD method. The most obvious one has to do with its In this section, we discuss the Faedo–Galerkin formulation
use of unstructured grids. Since the grid is unstructured, it as the general procedure for the development of semi-discrete
offers superior versatility in modeling complex geometries. approximations to Maxwell’s time-dependent system. Formu-
Furthermore, the Faedo–Galerkin procedure, used for the lations based on both the hyperbolic system of the two curl
development of the weak statement, provides for a very natural equations and the vector wave equation are discussed. The
way for handling field and flux continuity conditions at mate- derived formulations will serve as a unifying framework for
rial interfaces, thus further enhancing modeling accuracy. In the discussion of specific TDFEM’s in later sections. For the
addition, the use of Galerkin formulation provides a unifying purpose of generality, time discretization is not discussed in
approach for the exploitation of a variety of choices of trial this section. Instead, it is considered in conjunction with the
functions, some of them more appropriate for specific types specific TDFEM schemes in Sections III and IV.
of geometries than others.
Despite the aforementioned advantages there are some ma-
jor reasons why TDFEM’s have not been used widely. The A. Formulation Based Upon the Two First-Order Equations
most important is unstructured 3-D mesh generation. It re- Let , and be the position-dependent, relative permit-
quires a significant amount of effort to develop a proficiency tivity and permeability, respectively, in a volume , and let
Manuscript received May 8, 1996; revised October 2, 1996.
denote an applied current density inside this volume. The
J.-F. Lee is with EMCAD Laboratory, Department of Electrical and initial value problem (IVP) of interest is described by
Computer Engineering, Worcester Polytechnical Institute, Worcester, MA
01609 USA.
R. Lee is with the ElectroScience Laboratory, Electrical Engineering
Department, The Ohio State University, Columbus, OH 43210 USA.
A. Cangellaris is with the Electromagnetics Laboratory, Department of
Electrical and Computer Engineering, University of Arizona, Tucson, AZ
85721 USA. (1)
Publisher Item Identifier S 0018-926X(97)02292-8.

0018–926X/97$10.00  1997 IEEE

LEE et al.: TIME-DOMAIN FINITE-ELEMENT METHODS 431

inside , subject to boundary conditions of the form B. Formulation Based on the Second-Order
Vector Wave Equation
on To derive TDFEM’s in terms of the field, we start with
on the following initial value problem (IVP) description:

where , are perfect electric and magnetic surfaces, in

respectively, inside , and is the unit normal on these
surfaces. For the case of bounded domains, the tangential on
components of or , or an impedance relationship between
on (6)
and on the bounding surface must be specified for
all . For unbounded regions, Sommerfeld’s radiation
condition must hold. Finally, for a unique solution, the initial on
state for all fields inside at must be known.
Assuming the existence of a unique solution to (1), we can for
formulate two weak forms using the Galerkin or weighted
where is the truncation boundary. The formulation in
residual method. The first approach tests (1) by vector func-
terms of the magnetic field can be derived in a similar
tions , viz. square integrable vector functions in
(fields with finite energy). The weak form obtained is thus fashion. We have also adopted, for simplicity, the first-order
absorbing boundary condition (ABC) on to truncate the
infinite domain into a finite region.
The Galerkin form of the IVP (6) can be stated as

find such that

(2)

. By using this weak form, we seek the fields (7)

and
, where , and , where and are the trial and
test spaces, respectively. In general, in the weighted residual
method, the trial and test spaces are not necessarily the same.
One special case is the collocation method, where the test
on space is a finite number of delta functions located at some
chosen collocation points . However, for the collocation
on (3) method to work, the trial function must be smooth enough
at the collocation point such that the integrand in the
bilinear form is bounded. The challenges in the development
In the second approach [1], we multiply the first of (1) by
of a collocation method for unstructured grids are mainly
and integrate over . Similarly, multiplying the
the proper choice of the basis functions (explicitly and/or
second of (1) by , integrating over and
implicitly) and the collocation points. In particular, in the
integrating the curl term by parts, we obtain a weak form as
collocation FEM’s, the trial function must satisfy

(8)

(4) Equation (8) requires to be a smoother function than what

is usually needed for tangential continuity across element
boundaries. The continuity requirements for the trial function
in the collocation FEM’s are
In this case, the admissible trial spaces are
and . Notice that a
(9)
natural boundary condition comes with this second weak form
from the variational principle, namely (10)

on (5) Equations (9) and (10) are also required in the standard weak
Galerkin formulation. However, in the collocation FEM’s, (10)
. has to be a strong boundary condition (essential boundary
432 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

condition) instead of the weak boundary condition obtained are the vector basis functions that span the finite-
through the natural boundary condition in the weak Galerkin dimensional subspace . Note that we have used bold face
formulation. We also comment here that the collocation TD- letters to indicate matrices.
FEM approach, when properly formulated, offers great poten-
tial benefits such as superior accuracy (superconvergence) [2] III. TDFEM’S USING NODAL ELEMENTS
if optimal collocation points are chosen and explicit time step-
ping schemes for higher order (more than linear) interpolation. A. Two First-Order Equations with Point-Matching
However, this is a subject which receives virtually no attention
Among the variety of approaches used for the extension
currently in TDFEM for solving transient wave equation. For of Yee’s FDTD method to nonorthogonal grids, the method
readers who are interested in pursuing this topic further, we of point-matched finite elements was probably the first one
recommend [2] as a starting point. to propose the use of finite-element interpolation practices
A more common approach to the derivation of TDFEM’s is for the development of an explicit time-domain integration
to balance the smoothness requirements between test and trial scheme for Maxwell’s curl equations [3], [4]. The following
functions through the use of Green’s theorem. In particular, a is a concise description of the application of the method to
weak form of the IVP (6) can be stated as the discretization of Maxwell’s hyperbolic system, presented
in the light of the Faedo–Galerkin process.
find such that The point-matched finite-element method maintains the
staggering of the electric and magnetic fields in space;
however, in contrast to the FDTD method, all three
components of the electric field are placed at the same
node. The same is true for the magnetic field. Thus, two
complementary grids are used inside for the development
of the numerical approximation. The two grids are constructed
in such a manner that every electric field node is contained
within the volume of an element in the magnetic field grid, and
every magnetic field node is contained within the volume of an
element in the electric field grid. Using nodal finite elements in
(11) conjunction with these two grids leads to the finite-dimensional
subspaces that are used for the finite-element approximation
, and . of the fields (see Section II-A). More specifically, the fields
A vector function is admissible in this bilinear form if are approximated as follows:
and only if

(12) (15)

The finite-element process is to seek the best solution

of (11) in a finite-dimensional subspace where , are the number of electric and magnetic field
with a finite number of degrees of freedom. Consequently, the nodes in the two grids, and , are the scalar basis
functions used for effecting the finite-element approximations
application of the Faedo–Galerkin process results in a system
in the electric and magnetic grids, respectively. Clearly, the
of ordinary differential equations (ODE’s)
degrees of freedom of the approximation are the temporal
variations , at the nodes of the electric and
(13) magnetic grids, respectively.
To complete the Faedo–Galerkin process, the test functions
where is the coefficient vector of and need to be defined. Different choices of test functions are
possible and lead to either implicit or explicit integration
schemes. The most direct way to the development of an
explicit scheme is through point-matching or collocation. The
collocation points are taken to be the nodes in the two
complementary grids. In other words, delta functions ,
associated with the electric field nodes and
delta functions , associated with
the magnetic field nodes are used as test functions. Substituting
(14) (15) in (2) and testing the first of (2) with the delta functions
associated with the electric field nodes and the second of (2)
LEE et al.: TIME-DOMAIN FINITE-ELEMENT METHODS 433

with the delta functions associated with the magnetic field field nodes inside . The permittivity matrix with elements
nodes yields , , and the permeability ma-
trix with elements , .
The notation is used to denote the 3 3 unit matrix.
Using the aforementioned matrices and vectors, the state
equations (16) assume the compact matrix form

(20)

This step completes the development of the semi-discrete

approximation of Maxwell’s curl equations in the context
of point-matched finite elements. Use of a central difference
(16) scheme for the approximation of the time derivatives, with
and discretized in time in equal intervals and, in such
The above system of state equations may be cast in a a way that the temporal nodes of are shifted one-half time
compact matrix form as follows. First, we introduce the vectors step with respect to those of , results in the simple leapfrog
numerical integration algorithm for the state equations. This
algorithm is conditionally stable. The stability condition is
derived using the procedure of Section V and is
(17)
(21)
where it is
where denotes the maximum eigenvalue of the
matrix defined by .
It is important to point out that if standard finite-element
interpolation functions are used for the approximation of the
spatial variation of the fields, the matrices and are very
and the superscript denotes matrix transposition. Clearly,
sparse due to the compact support of the basis functions,
the dimension of and is , while the dimension of
and the numerical integration algorithm is very efficient.
is . Furthermore, we observe that the cross-product
Even though the use of two complementary grids appears
operations in (16) may be written in matrix form as follows:
cumbersome at first sight, it is actually quite straightforward
if quadrilateral elements in two dimensions and hexahedron
elements in three dimensions are used. The major drawback
of this method comes from the fact that all three components
(18)
for each of the fields are placed at the same node. This
leads to complications in enforcing the appropriate boundary
conditions at surfaces where the electromagnetic properties of
the media exhibit abrupt changes, in the sense that special
for , and treatment is needed for updating the nodes on such surfaces
[3], [4].
To alleviate this difficulty, the use of vector finite elements
has been proposed for the spatial interpolation of the fields [5].
(19) Hybrid finite-element interpolation schemes, implementing
different types of basis functions for the approximation of
the different field components and/or associated fluxes have
been proposed also [6]. As a matter of fact, this is an area of
for . The superscripts in the 3 continuing research that, one might argue, was motivated by
3 array and in the 3 3 array are used the pioneering work of Bossavit [7]. This work is discussed
to indicate that the elements of these arrays are calculated in Section IV.
at points , , and , ,
respectively. These definitions lead to the introduction of two B. Second-Order Equation with Integral Lumping
supermatrices, supermatrix with elements , In this section, we provide a brief review of the Lynch
, , and supermatrix with and Paulsen [8] formulation for a time-domain nodal finite-
elements , , . element method applied to the second-order equation with
Finally, we introduce two diagonal supermatrices to describe integral lumping. Because of the integral lumping, an explicit
the electromagnetic properties at the electric and magnetic time integration scheme can be obtained. The second-order
434 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

equation under consideration is somewhat different than (6) There are several good features about this formulation.
because of Lynch and Paulsen’s concerns about generating First, it is explicit, so no matrix solution is required. Second,
spurious solutions. From the analysis in [9], it is argued that a traditional frequency domain FEM code based on nodal
no spurious modes are produced if the following equation is elements can be easily adapted to this time-domain method.
used: However, there are several things that one must be careful
about in using this method. Because it is a node-based method,
in material discontinuities require special treatment (see [8] for
(22) details). Furthermore, the presence of perfectly conducting
edges and corners can produce large errors in the solution. The
with the same boundary conditions as in (6). The final weak errors are usually not localized near the edge or corner, but are
form of (22) can be written as present throughout the computation domain. Finally, there is
some concern about the accuracy of the method. The numerical
dispersion error for this FEM formulation is compared with
several others in a later section of this paper.

IV. TDFEM’S USING EDGE/FACET ELEMENTS

A. An Explicit WETD Method

It was long advocated by Bossavit [7] that in solving
Maxwell’s equations, Whitney 1-forms (edge elements) should
be used to approximate fields and , whereas Whitney 2-
(23) forms (facet elements) should be used to model the fluxes
and . Since then, many variants (with various interesting
where is approximated in terms of its nodal values by mutations) using this basic notion have been proposed and im-
plemented. Among them, we mention Mahadevan and Mittra
(24) [10], Wong et al.. [11], Lee and Sacks [12], Chan et al. [13],
and Feliziani and Maradei [5]. For Whitney 1-forms, the basis
The traditional central difference in time is used to ap- functions are well known by now. For example, for
proximate the second-order time derivative. Substituting this it is
approximation into (23), we obtain the following form for the
left-hand side of (23) (28)

(25) where is the Lagrange interpolation polynomial at vertex

[14]. The circulation of is one along and zero
where is the electric field at node and time step . along all other edges. Any square integrable vector field, , can
The resulting finite-element equation is implicit, requiring the now be interpolated by a combination of Whitney 1-forms as
solution of a matrix equation.
Lynch and Paulsen apply the concept of integral lumping (29)
(see the discussion on integral lumping later in this paper)
to obtain an explicit time-domain formulation. Thus, (25)
where is the circulation of along . Similarly, the
becomes
vector field for Whitney 2-forms associated with a particular
(26) can be written as

(30)
We can substitute (26) into (23) to obtain
It can be shown that the normal component of is
continuous across element boundaries, and its flux is one
through and zero through all other facets. Such
fluxes relate to the degrees of freedom of a tetrahedral element
by

(31)

where is the flux of vector through .

(27) We shall now present a very interesting dual-field TDFEM
formulation which is an explicit scheme using edge and facet
LEE et al.: TIME-DOMAIN FINITE-ELEMENT METHODS 435

elements. The scheme was first proposed by Bossavit in [7]. They are
It starts with

(38)
(32)
Using collocation, (38) becomes
and

(39)
where

(33)
(40)
Substituting (32) and (33) into Maxwell’s equations, we have
and and are understood to be the numbers
corresponding to the unknowns on and ,
respectively.
Approximation of the time derivatives in (37) can be done in
(34) several ways. Here we choose the central difference scheme
and stagger and in time, evaluating them one
half of a time step apart. Thus
To make (34) operational, we apply collocation method. For
example, on we have

(41)
Equation (41) is now fully discretized and can be used to
update the electromagnetic field distribution in a leapfrog
fashion.

B. General Implicit Two-Step Recurrence Algorithms

To formulate implicit TDFEM’s that involve the solution
of a matrix equation, we start with the semi-discrete equation
(35) derived previously

(42)
From the properties of Whitney 1-forms and 2-forms, (35)
reduce to the simple expressions Also, and are assumed to be given
either from the initial excitation or from previous iteration.
Even though the implicit schemes involve solving a matrix
equation for every time step, it is possible to derive implicit
TDFEM’s that are unconditionally stable. The unconditional
(36) stability is very desirable when solving problems with very
small features that lead to large variations in element size
These can be expressed in matrix form as across the problem domain. Without loss of generality, we
shall concentrate on deriving general two-step algorithms to
obtain using the weighted residual principle.
The approach begins by approximating the function ,
(37) by a second-order polynomial expansion

where indicates a column vector and is the circulation

matrix whose entries are either zero or one [clearly, only add or
subtract operations are involved in performing (37)]. Finally, (43)
to relate fields to fluxes we employ the constitutive relations.
436 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

From (43), the first- and second-time derivatives are approx- A. Stability Condition for the Explicit TDFEM
imated as For the purpose of deriving the stability condition, let us
rewrite (41), assuming in the following form:

(44)

Substituting (43) and (44) into the semi-discrete equation, (42)

results in
(48)

This can be expressed in a more compact matrix notation as

(49)

(45) where
Weighting the residual (45) by a test function , and
defining
and
(50)

A necessary condition for (50) to be stable is [15]

(51)

where is the th eigenvalue of the amplification ma-

trix . Consequently, the stability condition for the explicit
TDFEM of (48) is

(52)
(46)

where is the largest eigenvalue (in terms of mag-

nitude) of matrix . In practical computation, this largest
yields eigenvalue can be estimated efficiently either by the power
method or the Arnoldi method [16]. The stability condition
for the point-matched TDFEM is derived in a similar fashion
and is given in (21).

(47) B. Stability Condition for the Implicit TDFEM’s

Once has been obtained from the above equation, a new We shall follow the general approach outlined in
time step can be started. Implicit algorithms of various kinds Zienkiewicz and Taylor [17] for the second-order wave
can be obtained from (47) by employing different and equation. Assuming a modal decomposition of the system
values. Here, we recommend a particular pair of (13) with no excitation, it suffices to examine the scalar
and that results in an unconditionally stable implicit equation
algorithm with truncation error .
(53)
V. STABILITY CONDITIONS
where is the amplitude for mode in the system. We
In the linear theory of finite methods for the numerical
remark also that and , since is positive
solution of ordinary and partial differential equations, the
definite and and are semi-definite matrices. The two-step
success of numerical schemes is summed up in the well-
recurrence algorithm for takes the form
known Lax equivalence theorem [2]: “consistency and stability
imply convergence.” For general TDFEM’s, the consistency is
self-evident and assured by their formulation. Subsequently,
the question of stability is of paramount importance. In this
section, we present two approaches for deriving the stability
condition for general TDFEM’s. (54)
LEE et al.: TIME-DOMAIN FINITE-ELEMENT METHODS 437

The general solution can be written as this mismatch can be characterized in terms of an error in the
wave number where the numerical wave number is not
a linear function of frequency. Thus, an artificial numerical
dispersion is present in the solution.
(55)
For time-harmonic waves, the effects of numerical disper-
sion is to produce a phase error in the solution. The phase
where is the growth factor for mode . Substituting (55) error is cumulative, so as a wave propagates, the phase error
into (54) provides the second-order characteristic polynomial increases. Consequently, the error increases with the distance
for a wave propagates in the computation domain. The phase error
is especially significant for electrically large problems and
high resonators. A discussion of this phase error for the
FDTD method [18] is given in [19]. Its extension to TDFEM
(56) is straightforward if the mesh is assumed to be regular. In
this section, we compare the phase error for several TDFEM’s
through a two-dimensional dispersion analysis of square el-
For stability, it is necessary that the modulus of the growth
ements. The analysis does not truly model the phase error
factor be bounded by one for all modes, viz. .
for arbitrary unstructured meshes; however, it provides insight
To derive the stability limits, we shall use the so-called
into the relative performances of the various methods. In most
transformation. We introduce the mapping
finite elements, there is always an accumulation of phase error
(57) since the sign of the error is always the same for all angles
of incidence. Thus, although a nonorthogonal or unstructured
where as well as are in general complex numbers. It mesh does not produce the same phase error as a mesh
is easy to show that the aforementioned stability requirement composed of square elements, its overall characteristics does
is equivalent to demanding the real part of to be negative. not change unless the elements is very badly deformed. The
Equation (56) is now rewritten in terms of as one exception for frequency-domain methods is the tetrahedral
edge elements. For this element, the sign of the phase error
changes depending on the incidence angle. The phase error
(58) may cancel especially for unstructured grids [20]. In the time
domain, the phase cancellation can be obtained for methods
Furthermore, through the Routh–Hurwicz condition [17] it can which have implicit time integration, as we will see later in
be shown that in order for the real part of to be negative, this section.
the following conditions must hold To perform the dispersion analysis, let us consider a TE
plane wave propagating in the – plane through an infinite
mesh composed of square elements with sides of length . The
incidence angle of the plane wave is defined such that 0
(59) represent a propagating wave, and 90 represents a
propagating wave. The numerical solution is
The inequalities (59), as well as the fact that and
, lead to the following conditions on and :

and (61)
(60)

for (47) to be unconditionally stable. where we use the discretization , , and .

The fields in (61) can be substituted into the finite-element
VI. NUMERICAL DISPERSION ANALYSIS equations to obtain a transcendental equation for . Examples
of such an equation in the frequency domain are given in [21]
The accurate modeling of electromagnetic phenomena with
and [22]. For the purposes of comparison, we also consider
FEM is strongly dependent on both the shape and electrical
the dispersion analysis for the FDTD method.
size of the elements. Even one poorly shaped element (ele-
In writing out the dispersion relation for the various meth-
ments with large interior angles) can cause large solution errors
ods, we use certain abbreviations. They are ,
throughout the mesh. The element size is also very important;
, and . The numerical wave number
the elements must be small enough to model rapid field
of the FDTD method is given by
variations due to such things as geometrical discontinuities and
highly conductive media. There is an additional source of error
due to the incorrect modeling of the sinusoidal behavior of the
propagating wave. The piecewise polynomial approximation of (62)
FEM does not exactly match a sine function. The error due to
438 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

For the point-matched time-domain node-based FEM [23], the

dispersion relation is

(63)

For the time-domain node-based FEM with integral lumping

[8], we have

(64)

Finally, the dispersion relation for the implicit edge-based

time-domain method in (47) for the lossless case ( ) is
Fig. 1. Plots of phase error per wavelength as a function of the angle of
considered [12]. The numerical wave number, for , is incidence for the FDTD method and three FEM methods (point-matched
the solution to nodal FEM, nodal FEM with integral lumping, and edge-based FEM) with
=
h =20, p = 0:5, and = 23 .

Consequently, the choice results in a point-matched

(65)
FEM dispersion error much larger than what can be achieved
on a uniform grid for values of close to one. However, it
where , and is the objective of TDFEM’s to allow for unstructured grids
and, for general applications with considerable variability in
grid size, the numerical integration for large portions of the
grid will be done at Courant numbers smaller than one. Thus,
(66) the choice appears to be a reasonable one for the
purposes of demonstrating the dispersion error that should
The variable in (66) can be chosen to obtain a range of be expected from the application of point-matched FEM’s in
different time difference schemes. For the scheme is the modeling of realistic structures. Finally, we mention that
unconditionally stable. integral lumping negatively impacts the phase error since other
The phase error for each wavelength that the wave propa- researchers have observed this difference when comparing
gates is given by formulations with and without integral lumping [24].
Another important feature of numerical dispersion is its
anisotropy with respect to the incidence angle. A measure of
Phase Error (degrees) (67) anisotropy is the difference in the maximum and minimum
phase errors over all incidence angles. Let be this difference.
In Fig. 1, the phase error per wavelength ( ) is plotted as a From Fig. 1, we see that FDTD, node-based FEM with integral
function of angle of incidence from 0 to 90 for , lumping, and edge-based FEM have approximately the same
, and . Due to symmetry, the phase error anisotropy ( 0.75 ), while the anisotropy associated with
is periodic with period 90 , so the range that is plotted fully the point-matched node-based FEM is greater ( 6.7 ). The
describes the error for all angles of incidence. Fig. 1 reveals anisotropy of the phase error is important because techniques
some interesting characteristics of various finite methods. The are currently being developed to minimize the phase error for
edge-based FEM shows the least phase error and, in fact, is finite method [25]. These techniques can shift the numerical
the only one with a negative phase error. The point-matched wave number, so that the average phase error with respect to
FEM has the worst error, with the integral lumped nodal is zero; however, the anisotropy is still present.
FEM being the second worst. The poor performance of the In frequency-domain FEM methods, square bilinear first-
two node-based FEM methods compared to the edge-based order elements usually produce the least phase error at 45
method might appear surprising in view of the fact that the and the most at 0 and 90 , since the field is approximated
node-based FEM approximates the field to a complete first by a quadratic function at 45 and a linear function at 0 and
order, while the edge-based FEM only represents the field to an 90 [21]. From Fig. 1, we see that the effect of point matching
incomplete first order. However, it should be pointed out that and integral lumping is to exactly reverse the points of least
the stability condition for the point-matched FEM on a uniform and greatest phase error.
grid is given by , independent of the number of spatial To see how the phase error compares for a different choice
dimensions and, as discussed in [4], the algorithm should be of time step, let us choose with and .
run with as close as possible to one for minimum dispersion. The results are shown in Fig. 2. The phase error increases for
LEE et al.: TIME-DOMAIN FINITE-ELEMENT METHODS 439

Fig. 2. Plots of phase error per wavelength as a function of incidence angle Fig. 4. Plot of phase error per wavelength as a function of incidence angle
for FDTD and three FEM methods (point-matched nodal FEM, nodal FEM for the edge based FEM with h = =20, p = 0:5. The variable is varied
with integral lumping, and edge-based FEM) with h = =20, p = 0:2, and from 0.2–0.8 in increments of 0.2.
= 23 .

shift the level of the phase-error curve without changing the

anisotropy. From these curves, we can deduce that the best
choice of to minimize the maximum phase error is between
0.4 and 0.6.
Increasing the order of approximation within an element
decreases both the phase error and the anisotropy. In [21], it
has been shown that increasing the element order is more ef-
ficient than decreasing the element size for frequency-domain
FEM, and this behavior is expected to hold in the time domain
as well. This behavior provides yet another argument against
mass lumping since, as discussed in the next section, mass
lumping cannot be used in conjunction with higher order
approximations.

VII. THE MASS LUMPING PROCESS:

TO LUMP OR NOT TO LUMP
Formulations of TDFEM’s based on the second-order wave
Fig. 3. Plots of phase error per wavelength as a function of incidence angle equation using Galerkin’s method usually result in implicit
for FDTD and three FEM methods (point-matched nodal FEM, nodal FEM time-marching schemes. Motivated by the successes of finite-
with integral lumping, and edge-based FEM) with h = =10, p = 0:2, and
= 23 . difference methods on regular grids, various approaches have
been proposed to formulate explicit TDFEM’s which do not in-
volve inversion of matrix equation. As mentioned earlier, one
all methods, but the amount of anisotropy remains basically possibility is to employ collocation methods, an approach that
the same as compared to the case. Thus, the effect of has not received substantial attention within the computational
changing the time step is to shift the level of the phase error electromagnetics community yet. Another popular approach is
curves but not its shape. Finally, let us consider a coarser to lump the mass matrix into a diagonal matrix, thus reducing
mesh. The phase error is plotted in Fig. 3 for the case where the algorithm to an explicit TDFEM. In this section, we
, , and . The effect of the coarser mesh outline the lumping process for tetrahedral edge elements and
is to increase both the phase error and the anisotropy. Note comment on its implications in practical implementations. In
that for this grid spacing, the edge-based FEM has slightly general, these comments are also applicable to other types of
smaller anisotropy ( 2.8 ) as compared to FDTD and the finite elements.
node-based FEM with integral lumping ( 3.2 ). The mass lumping is typically achieved by using reduced
Up until this point, the results have been generated with integration. In particular, for a single tetrahedral edge element
for the edge-based FEM. In Fig. 4, the phase error for (see Fig. 5), we require the following mass integration to hold:
the edge-based FEM is plotted for values of from 0.2–0.8
in increments of 0.2. The other parameters are and (68)
. Similar to changing , the effect of changing is to
440 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

We conclude this section with a brief list the most common

problems with mass lumping.
1) It only works for low-order (at most linear) finite ele-
ments.
2) It results in worse numerical dispersion error than the
conventional Galerkin’s method.
3) It sometimes generates a diagonal matrix with zero
and/or negative entries.
In light of these problems as well as the potential benefits
from the use of collocation methods, we do not recommend the
use of mass lumping for the derivation of explicit TDFEM’s.

Fig. 5. A good shaped tetrahedral element. VIII. CONCLUDING REMARKS

In this paper, TDFEM’s were presented for the second-
order vector wave equation and Maxwell’s hyperbolic system.
where is the th diagonal entry in the final diagonal
Both node-based and edge-based elements were considered.
matrix. Thus, essentially, lumping is achieved by the proper
Because of the importance of stability for time-domain meth-
determination of the diagonal entries such that (68) will
ods, the underlying theory for the development of the stability
be exact for any constant field . To proceed
condition for a general TDFEM formulation was presented in
further, we first note that the unknown coefficient for edge
some detail. In addition, the numerical dispersion (or phase)
i is
error performance of various TDFEM’s was compared with
(69) that of the standard FDTD method. Finally, the advantages
and disadvantages of mass lumping were considered.
where it is assumed that edge is the edge from Progress in the advancement of TDFEM’s has been very
to . Substituting (69) into (68), and requiring that slow in comparison to the standard FDTD method. However,
it holds true for any arbitrary choice of , we obtain the there has been an increased interest in the last few years for
following six simultaneous equations for : the development of transient electromagnetic simulators with
geometric modeling versatility superior to the one offered by
the FDTD method. This demand has stimulated a significant
activity toward the development of various methodologies
for enhancing the standard FDTD method to facilitate the
modeling of nonrectangular geometries without the typical
staircase approximation. While many of the proposed numer-
ical approximations appear to be, at first sight, more general
implementations of the standard rectangular FDTD method,
their development is also possible through the Faedo–Galerkin
formulation discussed in this paper. Thus, many of these
methods may be viewed as special examples from the general
class of TDFEM’s. The success of these methods is strongly
dependent on the availability of reliable and robust mesh
generators. As mesh generation techniques mature and reliable
and user-friendly mesh generation packages start becoming
available, research work in development of TDFEM’s and their
(70)
application to transient electromagnetic modeling is expected
to increase significantly.
where is the volume of the tetrahedron. As mentioned in the text, there are several issues in TD-
To demonstrate one of the major undesirable characteristics FEM’s that need to be researched further, and the potential
of mass lumping, let us investigate the element shown in is there for the development of very sophisticated and ver-
Fig. 5. After some algebraic manipulations, it can be shown satile transient electromagnetic simulators. Among them we
that the diagonal entries are and identified the use of collocation methods and the possibility of
. Note that in this case the lumping process using different types of trial functions in different regions of
results in a singular mass matrix. In practical applications, the geometry of interest, appropriately selected so that they are
where many tetrahedral elements are involved (some could be consistent with the anticipated smoothness of the fields in the
very bad shaped), zero and/or negative diagonal entries from various regions. The use of more accurate spatial interpolation
mass lumping are not uncommon. This imposes a very severe schemes for the reduction of dispersion error have started
problem: no time step can be found to make the numerical making their appearance in the computational electromagnetics
algorithm stable. community both in the context of spectral approximations
LEE et al.: TIME-DOMAIN FINITE-ELEMENT METHODS 441

[26], [27] and in the context of wavelet expansions [28]. The [16] W. E. Arnoldi, “The principle of minimized iterations in the solution
use of Galerkin’s method in the development of TDFEM’s of the matrix eigenvalue problem,” Quart. Appl. Mathem., vol. 9, pp.
17–29, Jan. 1951.
offers the ideal framework for the implementation of such [17] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, 4th
diverse types of spatial approximations over different regions ed. New York: McGraw-Hill, 1988.
[18] K. S. Yee, “Numerical solution of initial boundary value problems in-
of the geometry under study in a consistent manner, maintain- volving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas
ing the appropriate field continuity conditions across region Propagat., vol. AP-14, no. 3, pp. 302–307, May 1966.
boundaries. [19] A. C. Cangellaris and R. Lee, “On the accuracy of numerical wave
simulations based on finite methods,” J. Electromagn. Waves Applicat.,
In addition to the aforementioned “hybridization” of the vol. 6, no. 12, pp. 1635–1653, 1992.
trial functions, there are some more topics associated with [20] J. Y. Wu and R. Lee, “The advantages of triangular and tetrahedral edge
TDFEM’s that have not been addressed in this paper and elements for electromagnetic modeling with the finite element method,”
IEEE Trans. Antennas Propagat., to be published.
are the subject of on-going research. One of them is the [21] R. Lee and A. C. Cangellaris, “A study of discretization error in the
development of radiation boundary conditions for TDFEM’s. finite element approximation of wave solutions,” IEEE Trans. Antennas
Research into accurate ABC’s is very sparse. To our knowl- Propagat., vol. 40, pp. 542–549, May 1992.
[22] G. S. Warren and W. R. Scott, “An investigation of numerical dispersion
edge, only the first-order ABC’s have been implemented in in the vector finite element method using quadrilateral elements,” IEEE
conjunction with TDFEM’s. We believe that advancement in Trans. Antennas Propagat., vol. 42, pp. 1502–1508, Nov. 1994.
[23] A. C. Cangellaris, C. C. Lin, and K. K. Mei, “Point-matched time
this area will benefit from the extensive work on ABC’s done domain finite element methods for electromagnetic radiation and scat-
for the FDTD method. More specifically, we expect truncation tering,” IEEE Trans. Antennas Propagat., vol. AP-35, pp. 1160–1173,
boundary conditions based on the perfectly matched layer Oct. 1987.
[24] P. Monk, A. K. Parrott, and P. J. Wesson, “Analysis of finite element
(PML) [29], [30] to find significant application in TDFEM’s time domain methods in electromagnetic scattering,” in 10th Annu. Rev.
also. They have been shown to work well for both the FDTD Progress Appl. Computat. Electromagn., Mar. 1994, vol. 1, pp. 11–18.
and frequency domain FEM, and the recent work in the [25] J. O. Jevtic and R. Lee, “An improved finite difference Helmholtz
equation,” in Joint AP-S/URSI Radio Sci. Meet., Baltimore, MD, July
reformulation of the PML theory in terms of time-dependent 1996, p. 187.
sources [31] is expected to facilitate their incorporation into [26] D. Gottlieb and B. Yang, “Comparisons of staggered and nonstag-
TDFEM’s. gered schemes for Maxwell’s equations,” in Proc. 12th Annu. Rev.
Progress Appl. Computat. Electromagn., Monterey, CA, Mar. 1996, pp.
1122–1131.
[27] A. C. Cangellaris and D. Hart, “Spectral finite methods for the simulation
REFERENCES of electromagnetic interactions with electrically long structures,” in
Proc. 11th Annu. Rev. Progress in Applied Computational Electromagn.,
[1] P. Monk, “A mixed method for approximating Maxwell’s equations,” Monterey, CA, Mar. 1995, pp. 1280–1287.
SIAM J. Numer. Anal., vol. 28, pp. 1610–1634, Dec. 1991. [28] M. Krumpholz and L. P. B. Katehi, “MRTD: New time domain schemes
[2] G. F. Carey and J. T. Oden, Finite Elements III: Computational Aspects. based on multiresolution analysis,” IEEE Trans. Microwave Theory
Englewood Cliffs, NJ: Prentice-Hall, 1984. Tech., to be published.
[3] A. C. Cangellaris, C. C. Lin, and K. K. Mei, “Point-matched time- [29] J. P. Berenger, “A perfectly matched layer for the absorption of
domain finite element methods for electromagnetic radiation and scat- electromagnetic waves,” J. Comp. Phys., vol. 114, pp. 185–200, 1994.
tering,” IEEE Trans. Antennas Propagat., vol. AP-35, pp. 1160–1173, [30] Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly
Oct. 1987. matched anisotropic absorber for use as an absorbing boundary condi-
[4] A. C. Cangellaris and K. K. Mei, “The method of conforming boundary tion,” IEEE Trans. Antennas Propagat., vol. 43, pp. 1460–1463, Dec.
elements for transient electromagnetics,” Progress Electromagn. 1995.
Res.—Finite Element Finite Diff. Methods Electromagn. Scattering, [31] L. Zhao and A. C. Cangellaris, “A general approach for the development
M. A. Morgan, Ed., 1989, vol. PIER-2, pp. 249–285. of unsplit-field time-domain implementation of perfectly matched layers
[5] M. Feliziani and F. Maradei, “Hybrid finite element solution of time for FD-TD grid truncation,” IEEE Microwave Guided Wave Lett., to be
dependent Maxwell’s curl equations,” IEEE Trans. Magn., vol. 31, pp. published.
1330–1335, May 1995.
[6] M. Feliziani and F. Maradei, “An explicit/implicit solution scheme to
analyze fast transients by finite elements,” in Proc. 7th Biennial IEEE
Conf. Electromagn. Field Comp., Okayama, Japan, Mar. 1996, p. 210.
[7] A. Bossavit and I. Mayergoyz, “Edge elements for scattering problems,”
IEEE Trans. Magn., vol. 25, pp. 2816–2821, July 1989.
[8] D. R. Lynch and K. D. Paulsen, “Time-domain integration of the
Maxwell equations on finite elements,” IEEE Trans. Antennas Propagat.,
vol. 38, pp. 1933–1942, Dec. 1990.
[9] K. D. Paulsen and D. R. Lynch, “Elimination of vector parasites in finite
element Maxwell solutions,” IEEE Trans. Microwave Theory Tech., vol.
39, pp. 395–404, Mar. 1991.
[10] K. Mahadevan and R. Mittra, “Radar cross section computation of
inhomogeneous scatterers using edge-based finite element method in
frequency and time domains,” Radio Sci., vol. 28, no. 6, pp. 1181–1193, Jin-Fa Lee (S’85–M’88) was born in Taipei, Taiwan, in 1960. He received
1993. the B.S. degree from National Taiwan University, in 1982, and the M.S. and
[11] M. F. Wong, O. Picon, and V. F. Hanna, “A finite element method based Ph.D. degrees from Carnegie-Mellon University, Pittsburgh, PA, in 1986 and
on Whitney forms to solve Maxwell equations in the time domain,” 1989, respectively, all in electrical engineering.
IEEE Trans. Magn., vol. 31, pp. 1618–1621, May 1995. From 1988 to 1990, he was with ANSOFT Corporation, Pittsburgh, PA,
[12] J. F. Lee and Z. S. Sacks, “Whitney element time domain methods,” where he developed several CAD/CAE finite-element programs for modeling
IEEE Trans. Magn., vol. 31, pp. 1325–1329, May 1995. three-dimensional microwave and millimeter-wave circuits. From 1990 to
[13] C. H. Chan, J. T. Elson, and H. Sangani, “An explicit finite difference 1991, he was a Postdoctoral Fellow at the University of Illinois at Urbana-
time domain method using Whitney elements,” in IEEE Antennas Champaign. Currently, he is an Associate Professor at the Department
Propagat. Symp. Dig., Seattle, WA, July 1994, vol. 3, pp. 1768–1771. of Electrical and Computer Engineering, Worcester Polytechnic Institute,
[14] G. Strang and G. J. Fix, An Analysis of the Finite Element Method. Worcester, MA. His current research interests include analyses of numerical
Englewood Cliffs, NJ: Prentice-Hall, 1973. methods, coupling of active and passive components in high-speed elec-
[15] R. D. Richtmyer and K. W. Morton, Difference Methods for Initial-Value tronic circuits, three-dimensional mesh generation, and nonlinear optic fiber
Problems. New York: Wiley, 1967. modelings.
442 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 3, MARCH 1997

Robert Lee (S’82–M’90) received the B.S.E.E. degree from Lehigh Univer- Andreas Cangellaris (M’86), received the Dipl. degree in electrical engi-
sity, Bethlehem, PA, in 1983, and the M.S.E.E. and Ph.D. degrees from the neering from the Aristotle University of Thessaloniki, Greece, in 1981, and
University of Arizona, Tucson, in 1988 and 1990, respectively. the M.S. and Ph.D. degrees in electrical engineering from the University of
From 1983 to 1984, he worked for Microwave Semiconductor Corporation, California, Berkeley, in 1983 and 1985, respectively.
Somerset, NJ, as a microwave engineer. From 1984 to 1986 he was a Member From 1985 to 1987, he was with the Electronics Department of Gen-
of the Technical Staff at Hughes Aircraft Company, Tucson, AZ. From eral Motors Research Laboratories, Warren, MI. In 1987 he joined the
1986 to 1990 he was a Research Assistant at the University of Arizona. Department of Electrical and Computer Engineering at the University of
In addition, during the summers of 1987 through 1989, he worked at Sandia Arizona, Tucson, where he is now an Associate Professor. His research
National Laboratories, Albuquerque, NM. Since 1990, he has been at The interests include computational electromagnetics, microwave engineering, and
Ohio State University, where he is currently an Associate Professor. His major algorithm development for high-speed digital circuit analysis. In the area
research interests are in the analysis and development of finite methods for of computational electromagnetics his research has emphasized differential
electromagnetics. equation-based techniques, both in the frequency and time domain. He has
Dr. Lee is a member of the International Union of Radio Science (URSI) published more than 50 journal papers in the aforementioned areas of research.
and was a recipient of the URSI Young Scientist Award in 1996. Dr. Cangellaris is a member of the International Union of Radio Science
(URSI). In 1993 he received the International Union of Radio Science (URSI)
Young Scientist Award. He is an active member of the following IEEE
Societies: Antennas and Propagation, Microwave Theory and Techniques, and
Components, Packaging, and Manufacturing Technology.

View publication stats