A new method for solving Maxwell's equations in the time domain for arbitrary values of permittivity, conductivity, and permeability is presented. Spatial derivatives are found by a Fourier transform method and time integration is performed using a second order, semi-implicit procedure. Electric and magnetic fields are collocated on the same grid points, rather than on interleaved points, as in the Finite Difference Time Domain (FDTD) method. Numerical results are presented for the propagation of a 2-D Transverse Electromagnetic (TEM) mode out of a parallel plate waveguide and into a dielectric and conducting medium

Shebalin, John V.

English

NASA Technical Reports Server

NASA Technical Memorandum 104162
sl 7U _J'
A Fourier Collocation Time Domain Method
for Numerically Solving Maxwell's Equations
(NASA-TM-I04162) A FOURIER COLLOCATION
DOMAIN METHOD FOR NUMERICALLY SOLVING
MAXWELL'S EQUATIONS (NASA) 19 p CSCL
TIME N92-I1742
12A
Unclas
G3/64 0051754
John V. Shebalin
October 1991
[¢/LqA
National Aeronautics and
Space Administration
Langley Research Center
Hampton, Virginia 23665-5225
https://ntrs.nasa.gov/search.jsp?R=19920002524 2020-03-17T15:15:37+00:00Z

Abstract
A new method for solving Maxwelrs equations in the time domain for arbitrary values of
permittivity, conductivity, and permeability is presented. Spatial derivatives are found by a
Fourier transform method and time integration is performed using a second-order, semi-
implicit procedure. Electric and magnetic fields are collocated on the same grid points,
rather than on interleaved points, as in the finite difference time domain (FDTD) method.
Numerical results for the propagation of a 2-D TEM mode out of a parallel plate wave guide
and into a dielectric and conducting medium is presented.
\
\
\
\
\
\
\
\
Introduction
Solving Maxwell's equations by finite difference procedures in the time domain has been
steadily gaining popularity since the pioneering work of Yee [1] and Taflove [2]; a good
example (with more references) of this f'mite difference time domain (FDTD) method is
given in the recent work of Britt [3]. Here, we present the beginnings of a new method, a
Fourier collocation time domain (FCTD) approach. We will concentrate here on two-
dimensional (2-D) simulations, as these provide a good compomise between reality and
resolution; 1-D and 3-D simulations are straightforward extensions of the results presented
here.
First, we will define our coordinate system and model problem, and second, set up the 2-D
system of Maxwelrs equations which we wish to solve. Third, we will explicitly discuss
the FCTD method of solution. Fourth, we will present numerical results, and lastly, give a
summary and conclusion.
2
Problem Definition
The physical problem we are attempting to solve by numerical simulation consists of a
macrowave horn emitting a single frequency, continuous wave signal through a dielectric
slab and into a plasma. The motivation for this problem is the planned use of a microwave
reflectometer to determine critical electron number densities in the stagnation region of
hypersonic reentry vehicle.
The 2-D model problem space is shown in Figure 1; the horizontal coordinate is z, the
vertical coordinate is x, and the (third dimensional) y axis is perpendicular (and outward) to
the plane of the figure. The initial conditions will consist of a sinusoidal TEM mode
(truncated to one wavelength) in the parallel plate part of the horn. The plasma to the right
of the dielectric slab has a gaussian profile in both the z and x directions. The physical
3dimensions of the (square) 2-D space are 10 cm on a side and the TEM mode has a
wavelength of 2.5 cm.
,w t
2-D Maxwell's Equations
For now, the magnetic permeability will be assumed to have only its free space value;
this case, the 2-D Maxwelrs equations are (SI format):
in
TM eqs.:
E _tEx = - _ _zBy - jx
3tBy = 3xEz- 0zEx
8 DtEz = g_ 0x By- jz
TE eqs.:
0tBx = 0zEy
 0tF_,= l(0zBx- xBz)-jy
_tBz = -_xEy (1)
where the current j (and conductivity o) are:
j (t) = f',. o (t-x) E ('0 dx ; MO(t) = _ _(_) e-i_t do3
(2)
Some initial conditions which can be used for the parallel plate wave guide of Figure 1 are:
TM modes: TE modes:
TEM: Ex = By, Ez = 0
TMI: k_='/gk, k= m
2 C
Ex = COs(kl z) cos(_ x)
By= Ex
_c
Ez = _3 sin(k, z) sin(_ x)
TE_: k_=_-k, k=m
2 C
Bx = _ sin(k1 z) sin(_k x)
_2_.¢.
F-,y= /g. Bx
Bz = _ cos(k,z)cos(_kx)
(3)
In the frequency domain, the complex 'permittivity' is defined in terms of the dielectric
permittivity ¢ and conductivity o as:
i o(c0)
(4)
4
where ofco) is the Fourier transform of o(t) (and vice versa):
t
_ (t) = 2-_nI_ _ (co) eic° t dco
(5)
but e _ the Fourier transform ofe(co) (i.e., e is always positive while Re(efco_) need
not be). For a plane wave of angular frequency co propagating in a plasma, a useful
expression for o(c0) is [4]:
_(o_) = i
co (6)
In the time domain, assuming a single frequency, we use for e whatever the coresponding
value is for that frequency; for o(t) in (2), we will use a constant value equal to I ¢c0)lat
the frequency of interest, co = me:
o(t) = Oo (t), Go="-v
020 (7)
Thus, in (1), j= ooE. In our numerical example, we will choose
COp= 2n x 25 GHz (maximum) coo = 2n x 12 GHz
so that at maximum electron density, Cro= 3 mho/m.
The FCTD Method
The FCTD Method consists of a semi-implicit time integration scheme and a Fourier
transform technique ("spectral method") for determining spatial derivatives. Here, spatial
grids are NxN, where N is a power of 2; in the numerical example we will use N=128.
Since j= (_oE, the first TM equation in (1) can be written as, for example:
5
e 0rEx = - _ 0zBy - ffoEx (8)
When numerically integrating this, we can treat the last term on the right-hand-side
implicitly, i.e.,
8E_x l:EE_x_At[l.tj(OzB_)+t_oE_x q _ E_x l=E_x-At _J(_zB_)
e+At _o (9)
This first-order time integration is absolutely stable with respect to the conductivity; in fact,
any non-negative value of _ is allowed on the model grid.
Following this example, a second-order time integration mehtod for the TM equations in
(1) is:
Predictor
8_+ Atao
2 e°
2 l_o
Corrector
e_x- At c2(0zB_ -la)1-
8_+ At °°
Eo
B_ 1 = B_ + At (_)xEz_+la- _)zE_xm)
+Atc:(OxBWWI=
Er + At_
(10)
The 1/2timestepvaluesof thefieldsareonly usedasintermediatevaluesandarenotsaved
from onetime stepto thenext.
6
The spatial derivatives in (10) are evaluated by a Fourier transform method; first, we
transform to k-space, f'md the Fourier coefficients of the derivatives of the fields, and then
transform back to x-space:
E(k) = !r_ E(x)e -ik'x ==¢, VE(x) =lN_ ikE(k)e ik'x
x k (11)
In (11) (and (12) below) we form a dyadic which can be specialized to a dot or cross
product; also, x=2rt(j,k)/N (O<j,k <N-I), and k=(m,n), (N/2- l<m,n<N/2); j, k, m,
and n are integers.
The Fourier transform method is equivalent to an N th order f'mite difference method, if we
view a generalized difference method as a convolution:
VE(x) = _ D(x-x') E(x') D(x-x') -- _ -k sin [k-(x-x')]
x' N k (12)
However, instead of taking O(N 2) operations for x-derivatives per point on an NxN grid,
as the corresponding (N th order) f'mite difference operation would, the number of
operations (per point) is proportional to NlogN. The ratio of accuarcy to time spent
differentiating is thus much higher than the corresponding f'mite difference method. (The
generic accuracy of spectral methods is well known [5], and will not be explicitly
considered here.)
An FDTD method [1,2,3], has only second order spatial differencing, which makes it
faster, but less accurate than the FCTD method described here. Another difference is that
in theFCTD method,bothE andB aredefinedatthesamephysicalgrid point,while in the
FDTD method,E andB arenotdef'medatthesamephysicalgrid point,but rather
'leapfrog'overoneanotherandlie on interleavedgrid points.
7
Thus,in theFCTD method,E andB arecollocated at all grid points. Collocation,
however, also has another meaning in FCTD. In numerical analysis, a collocation method
is a member of the methods of weighted residuals, and is also referred to as a
pseudospectraI method [6]. A collocation method assumes a solution in terms of known
expansion functions and requires that the solution be exact (with respect to the function
expansion) at the grid points; here, this means that we determine our spatial derivatives
exactly (in terms of trigonometric functions). Since we are solving time-dependent
equations, we have a Fourier-collocation time-domain, or FCTD, method. (If the
equations were linear, which they are not because of the spatial variation of the material
constants c, It , and _, we would have a Fourier-Galerkin time-domain method [6].)
There are two kinds of boundaries in the model space. First, there are interior boundaries,
i.e., boundaries between various media within the model space. Second, there is the
exterior boundary of the model space itself. In the numerical example to be presented here,
the interior boundaries are linear, and at either O, 45, or 90 degrees to the horizontal axis.
Material properties (such as _) are allowed to change discontinuously across the interior 0
or 90 degree boundaries, while on the 45 degree boundaries, material properties are set to
a value halfway between those on either side of the interface (although this is not an
essential part of the method). The exterior boundary, because of the use of a Fourier
transform in evaluating derivatives, is periodic; this poses few problems, however, as long
as the bulk of the model electromagnetic energy is contained naturally within the exterior
boundary (as occurs in the numerical example to be presented). A problem with an open
exterior boundary will require an algorithm for absorbing energy, or expansion functions
other than trigonometric (or both), and is not addressed here.
Numerical Results
As a preliminary test case, the 2-D geometry shown in Figure 1 was set up. There are three
essentially separate regions: a parallel plate microwave "horn" on the left of the figure, a
dielectric slab of a relative permittivity of four in the middle of the figure, and a "plasma" to
the fight of the dielectric slab. The size of the numerical grid corresponding to Figure 1 is
128x128; the physical size of the grid is 10 cm on a side. The dielectric slab begins on
horizontal grid point 52 and and ends on grid point 84; this corresponds to about 2-1/2 cm.
The lower comers of the horn are at (32,48) and (48,32), while the upper corners are at
(32,80) and (48,96). The material of the horn is assumed to be solid and have a
conductivity of 10,000 mhos/m at the transmission frequency of 12 GHz (although any
value can be chosen). The plasma has a conductivity of
[ 2 2o(z,x) = o o exp -(Z-Zo) ¢2Sz- (X-Xo)Z/2s
-1
Oo=3 (_m) , Zo=8.44 cm, Xo=5.0Ocm, Sz=7.81 cm, Sx = 1.56 cm(13)
8
and is conf'med to the fight-hand-side of the dielectric slab.
The results are shown in Figures 2 through 7. The simulation time step size was At - 10 14
sec, the total number of time steps was 48,000, and the cpu time/time step was 0.164
seconds/At (on a Cray-2). The total energy at the beginning of the simulation has a relative
value of unity, and decreases as time wears on. The electromagnetic energy density, which
is displayed in the Figures, is defined as
TZ-M = 2(eE2 + I.tH 2) (14)
TE-Mcan also be integrated over a selected region in physical space. This is done for the
three disjoint regions of the model space: to the left of the dielectric slab, within the
dielectricslab,andto therightof thedielectricslab;
of theFigures.
thisregional energy is given on each
g
It should be noted that there is a reflected wave in Figures 3 and 4 which barely shows up
since its energy, which is small relative to the transmitted wave's energy, falls mostly
below the lowest contour value in those Figures. The various reflections show up much
better in Figures 5-7, as the primary wave has been greatly attenuated in the plasma to the
right of the dielectric slab, and the relative energy of these reflections increases.
Conclusion
The results, as shown in the Figures, appear to be an accurate representation of the
propagation of an electromagnetic wave through dielectric and conducting media (a 1-D test
case with a linear grid of 2048 points, in which a plane wave was normally incident on a
wide dielectric slab of relative permittivity 4, reproduced analytic predictions for reflection
and transmission amplitudes to one part in 104). Material constants may be arbitrarily
specified on a computational grid, so that practically any physical situation can be
simulated. The FCTD method is of intrinsically higher accuracy than the FDTD method, as
the FCTD method has (essentially) N th order spatial differencing, while the FDTD method
has only 2 nd order spatial differencing. While a 2-D example was presented here, it is a
straightforward procedure to use the FCTD method for either 1-D or 3-D problems. (We
should also note, at this point, thatfinite element techniques utilizing unstructured grids
offer an alternative which may prove to be more efficacious, in general, than either FDTD
or spectral methods.)
Although a Fourier transform method was used to evaluate derivatives, other function
expansions (which may be more suitable for modelling certain boundary conditions) can be
utilized (e.g., Chebyshev polynomials); the FCTD method is thus only one example of a
spectral method [7]. In the future, we plan to investigate the use of alternative function
expansions,aswell asto investigatetechniquesfor approximatingabsorbingboundary
conditionsandtechniquesfor modellingthepropagationof signalswith broadfrequency
content.Applicationswhichcanbeaddressedincludethesimulationof microwave
reflectometersinteractingwith theplasmageneratedin thebowshockof ahypersonic
reentryvehicle.
10
Acknowledgements
Thenumericaldata presented here came from computer runs performed on the NASA
Ames Cray-2. The efficient transfer and visualization of numerical data was accomplished
using the National Center for Supercomputing Applications (NCSA) Telnet and Imagetool
software, respectively.
References
[1] K. S. Yee, "Numerical Solution of initial boundary value problems involving Maxwelrs equations in
isotropic media", IEEE Trans. Antennas & Propagation, AP-14, 302-307, 1966.
[2] A. Taflove, "Application of the finite-difference time-domain method to sinusoidal steady-slate
electr0magnetic-_fietration problems", IEEE Trans. Electromagn. Compat_, EM_-22, 19 i-202, 1980.
[3] C.L. Britt, "Solution of electromagnetic scattering problems using time domain techniques", IEEE
Trans. Antennas & Propagation, AP-37, 1181-1192, 1989.
[4] J.D. Jackson, Classical Electrodynamics, 2nd Ed., p 287, Wiley, New York, 1975.
[5] D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods." Theory and Application,
pp 135-8, SIAM, Philadelphia, 1977.
[6] C. A. J. Fletcher, Computational Galerkin Methods, p 27 & p 205, Springer-Verlag, New York, 1984.
[7] E. A. Coutsias, F. R. Hansen, T. Huld, G. Knorr, & J. P. Lynov, "Spectral Methods in Numerical
Plasma Simulation", Physica Scripta, 40, 270-279, 1989.
11
o = 10 4 (f2 m) 1
X
Z
£r=4
Regional energy = 1.000 = 0.000
Total energy = 1.000 at 0.00 nsec
= 0.000
Figure 1. 2-D model space with initial energy distribution, along with pcrmittivity and conductivity.
12
Regional energy = 0.696 = 0.301 0.000
Total energy = 0.997 at 0.08 nsec
Figure 2. Electromagnetic energy distribution at 0.08 nsec.
13
<::D
Regional energy = 0.128 = 0.858
Total energy = 0.987 at 0.16 nsec
= 0.001
Figure 3. Eleclromagnetie energy distribution at 0.16 nsec.
14
<=>
<_
Regional energy = 0.094 = 0.785
Total energy = 0.953 at 0.24 nsec
= 0.074
Figure 4. Electromagnetic energy distribution at 0.24 nsec.
\
15
Regional energy = 0.082 = 0.171
Total energy = 0.306 at 0.32 nsec
= 0.053
Figure 5. Elecu'omagnetic energy distribution at 0.32 nsec.
16
Regional energy= 0.031 = 0.167
Total energy = 0.206 at 0.40 nsec
= 0.008
Figure 6. Electromagnetic energy distribution at 0.40 nsee.
17
c::::::z::>
Regional energy = 0.086 = 0.093
Total energy = 0.183 at 0.48 nsec
- 0.004
Figure 7. Elecuomagnetic energy distribution at 0.48 nsec.
Form Approved
REPORT DOCUMENTATION PAGE OMBNo.0704-0188
Pub c reporting burden for this toltec%ton of information is estimated to average 1 hour per respOnse, including the time for reviewing instructions, searching existing data sources.
gather ng and maintaining the data needed, and completing and reviewing the collecl=on of _nformat=On, Send comments regarding this burden estimate or any other aSl_eCt of this
collection of =nformation, including suggestions for reducing this burden, lo Washington Headquarters Services, Directorate for information operations and Reports, 1215 Jefferson
Dav=s Highway. Suite 1204 Arlington, VA 22202-4302. and to the Office of Management and Budget. Paperwork Reduction Project (0704-0 _88), Washington, DC 2050].
1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED
October 1991 Technical Memorandum
4. TITLE AND SUBTITLE 5. FUNDING NUMBERS
A Fourier Collocation Time Domain Method for
Numerically Solving Maxwell's Equations WU 505-64-70-01
6. AUTHOR(S)
John V. Shebalin
7. PERFORMING,ORGANIZATION NAME(S) AND ADE)RESS(ES)
NASA Langley Research Center
Hampton, VA 23665-5225
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)
National Aeronautics and Space Administration
Washington, DC 20546-0001
8. PERFORMING ORGANIZATION
REPORT NUMBER
10. SPONSORING / MONITORING
AGENCY REPORT NUMBER
NASA TM-I04162
11. SUPPLEMENTARY NOTES
12a. DISTRIBUTION/AVAILABILITY STATEMENT
Unclassified - Unlimited
Subject Category 64
12b. DISTRIBUTION CODE
13. ABSTRACT (Max_um2OOwords)
A new method for solving Maxwell's equations in the time domain for arbitrary values
of permittivity, conductivity, and permeability is presented. Spatial derivatives
are found by a Fourier transform method and time integration is performed using a
second-order, semi-implicit procedure. Electric and magnetic fields are collocated
on the same grid points, rather than on interleaved points, as in the Finite
Difference Time Domain (FDTD) method. Numerical results for the propagation of a
two-di_ensional Transverse Electromagnetic (TEM) mode out of a parallel plate wave
guide and into a dielectric and conducting medium is presented.
14. SUBJECT TERMS
computational methods, electromagnetics, Fourier analysis
17. SECURITY CLA'$SIFICATION 18.
OF REPORT
Unclassified
NSN 7540-01-280-5500
SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION
OF THIS PAGE OF ABSTRACT
Unclassified
15. NUMBER OF PAGES
18
16. PRICE CODE
A03
20. LIMITATION OF ABSTRACT
Standard Form 298 (Rev 2-89)
Pfescrtbed by ANSI Std Z3g-18
29B.I02


A Fourier collocation time domain method for numerically solving Maxwell's equations

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server