数学建模英文论文(l.2a)
(l.2b)
One-Dimensional Simulation with the FDTD Method
This chapter is a step-by-step introduction to the FDTD method. It begins with the simplest
possible problem, the simulation of a pulse propagating in free space in one dimension. This
example is used to illustrate the FDTDformulation, Subsequent sections lead to the formulation
for more complicated media
1.1 ONE-DIMENSIONAL FREE SPACE FORMULATION
The time-dependent Maxwell's curl equations in free space are
aE 1 - = -V x H (l.la) at 80
en 1
- = --V x E. (l.lb) at /-to
E and H are vectors in three dimensions, so in general, Eq. (l.la) and (l.lb) represent three
equations each. We will start with a simple one-dimensional case using only Ex and By, so
Eq. (l.la) and (l.lb) become
aEx __~ aHy
at - 80 az
en, 1 ss,
=--- at /Lo OZ
These are the equations of a plane wave with the electric field oriented in the x direction, the
magnetic field oriented in the y direction, and traveling in the z direction.
Taking the central difference approximations for both the temporal and spatial derivatives
gives
E;+1/2(k) - E;-1/2(k) = _~ H;(k + 1/2) - H;(k - 1/2)
~t ~ ~x
H;+l(k + 1/2) - H;(k + 1/2) 1 E;+1/2(k+ 1) - E;+1/2(k) =-----------
~t /Lo ~x
(l.3a)
(l.3b)
1
2 Chapter 1 • One-Dimensional Simulation with the FDTD Method
In these two equations, time is specified by the superscripts, i.e., "n" actually means a time
t = Ilt · n. Remember, we have to discretize everything for formulation into the computer.
The term "n + 1" means one time step later. The terms in parentheses represent distance,
i.e., "k" actually means the distancez = Ilx · k. (It might seem more sensible to use Ilz as
the incremental step, since in this case we are going in the z direction. However, Ilx is so
commonlyusedfor a spatialincrement that I will use Ilx.) The formulation of Eq. (1.3a)and
(1.3b) assumes that the E and H fields are interleaved in both space and time. H uses the
arguments k + 1/2 and k - 1/2 to indicatethat the H field values are assumedto be located
between the E field values. This is illustrated in Fig. 1.1. Similarly, the n + 1/2 or n - 1/2
superscript indicatesthat it occurs slightlyafter or beforen, respectively.
E n-ll2
x k-2 k-l k k+l k+2 L
\ I Hj k - 1 112 Ik - 112
~ ~ Ik+ 112 k+2112 Ik + 11/2
\ I
E n+ll2 ~ ~ x I k-2 k-l k k+ 1 I k+2 L
Figure 1.1 Interleaving of the E andH fieldsin spaceandtimein theFOTDfonnulation.
Tocalculate Hy(k + 1/2), for instance, the neighboring values of Ex at k and
k +1areneeded. Similarly, to calculateEx(k + 1), the valueof Hy atk +1/2
and k + 1 1/2 are needed.
Eq. (1.3a)and (1.3b)can be rearranged in an iterative algorithm:
En+l/2(k) =En- I/2(k) - ~[Hn(k + 1/2) - Hn(k - 1/2)] (1.4a)
x x 80 • Ilx y y
Hn+l(k + 1/2) = Hft(k + 1/2) - At [E ft+l/2 (k +1) - En+I
/
2 (k)] . (1.4b)
Y Y J.Lo · Ilx x x
Noticethat the calculations are interleaved in both spaceand time. In Eq. (1.4a),for example,
the newvalueof Ex is calculatedfrom the previous valueof Ex and the most recent valuesof
By. Thisis thefundamental paradigmofthefinite-difference time-domain (FDTD)method[1].
Eqs. (1.4a) and (1.4b) are very similar, but becauseeo and #Lo differ by severalorders
of magnitude, Ex and By will differby severalorders of magnitude. This is circumvented by
makingthe following changeof variables [2]:
E= r!i.E . (1.5) ,,~
Substituting this into Eqs. (1.4a)and (1.4b)gives
£n+l/2(k) =En- 1j 2(k )__1_ dt [Hn(k + 1/2) - Hft(k - 1/2)] (1.6a)
x x JeO#Lo Sx y Y
H n+1(k + 1/2) = Hn(k + 1/2) - _1_ dt [En+I /2(k + 1) - Eft+I / 2(k )] (1.6b)
y y JeoJLo dx x x
Section 1.1 • One-Dimensional Free Space Formulation 3
Once the cell size ~x is chosen, then the time step ~t is determined by
~x
~t=-2
· Co
(1.7)
where Co is the speed of light in free space. (The reason for this will be explained later.)
Therefore,
1 ~t ~x/2 . Co 1
---=Co· =
Jeo~o ~x ~x 2
Rewriting Eqs. (1.6a) and (1.6b) in C computer code gives the following:
(1.8)
ex[k] = ex[k] + 0.5*( hy[k-l] - hy[k] )
hy[k] = hy[k] + 0.5*( ex[k] - ex Ik- t l ).
(1.9a)
(1.9b)
Note that the n or n + 1/2 or n - 1/2 in the superscripts is gone. Time is implicit in the
FDTD method. In Eq. (1.9a), the ex on the right side of the equal sign is the previous value
at n - 1/2, and the ex on the left side is the new value, n + 1/2, which is being calculated.
Position, however, is explicit. The only difference is that k + 1/2 and k - 1/2 are rounded off
to k and k - 1 in order to specify a position in an array in the program.
The program fdld_l.1.c at the end of the chapter is a simple one-dimensional FDTD
program. It generates a Gaussian pulse in the center of the problem space, and the pulse
propagates away in both directions as seen in Fig. 1.2. The Ex field is positive in both
directions, but the By field is negative in the negative direction. The following things are
worth noting about the program:
1. The Ex and By values are calculated by separate loops, and they employ the interleaving
described above.
2. After the Ex values are calculated, the source is calculated. This is done by simply
specifying a value of Ex at the point k = kc, and overriding what was previously
calculated. This is referred to as a "hard source," because a specific value is imposed
on the FDTD grid.
T= 100
~~ 0 ..---------~
20 40 60 80 100 120 140 160 180 200
-1 '----__~__..lo____~__.....4.-__~ __01___~____L.__~
o
80 100 120 140 160 180 200
FDTDcells
20 40 60
-1 ~-_ _'_____....l__...__~__--'--____L_____L____...L______'______L_____J
o
::t:~ 0t-----------.
Figure 1.2 FDTD simulation of a pulse in free space after 100 time steps. The pulse
originated in the center and travels outward.
4 Chapter I • One-DimensionalSimulationwith the FDTDMethod
PROBLEM SET 1.1
1. Get the program fdId_I.t.e running. What happens when the pulse hits the end of the array?
Why?
2. Modify the program so it has two sources,one at kc - 20 and one at kc + 20 (Noticethat
ke is the center of the problemspace). Whathappenswhenthe pulses meet? Explain this from
basic EM theory.
3. Instead of Ex as the source, use By at k = kc as the source. What difference does it make?
Try a two-point magnetic source at kc - 1 and kc such that hy [kc-l] = -hy [kc].
What does this look like? What does it correspond to physically?
1.2 STABILITY ANDTHE FDTDMETHOD
Let us return to the discussion of how we determine the time step. An electromagnetic wave
propagating in free space cannot go faster than the speed of light. To propagate a distance
of one cell requires a minimum time of l1t = Sx Lc«. When we get to two-dimensional
simulation, we have to allow for the propagation in the diagonal direction, which brings the
time requirement to 6.t = 6.x/(,J'ico). Obviously, three-dimensional simulation requires
6.t = 6.x/(,J3co). This is summarized by the well-known "Courant Condition" [3,4]:
~x
!!it -< -In.- (1.10) co'
where n is the dimension of the simulation. Unless otherwise specified, throughout this book
we will determine l1t by
~x
6.t = --.
2· Co
This is not necessarily the best formula; however, we will use it for simplicity.
PROBLEM SET 1.2
(1.11)
1. In fdId_l.l.c, go to the two governingequations, Eq. (l.9a) and (1.9b), and change the factor
0.5 to 1.0. What happens? Change it to 1.1. Now what happens? Change it to .25 and see
what happens.
1.3 THE ABSORBING BOUNDARY CONDITION IN ONE
DIMENSION
Absorbing boundary conditions are necessary to keep outgoing E and H fields from being
reflected back into the problem space. Normally, in calculating the E field, we need to know
the surrounding H values; this is a fundamental assumption of the FDTD method. At the edge
of the problem space we will not have the value to one side. However, we have an advantage
because we know there are no sources outside the problem space. Therefore, the fields at the
edge must be propagating outward. We will use these two facts to estimate the value at the
end by using the value next to it [5].
Suppose we are looking for a boundary condition at the end where k = O. If a wave is
going toward a boundary in free space, it is traveling at Co, the speed of light. So in one time
step of the FDTD algorithm, it travels
. ~x ~x
distance =Co . l1t =Co . - = -.
Co 2
Section 1.4 • Propagation in a Dielectric Medium 5
This equation basically explains that it takes two time steps for a wave front to cross one cell.
So a common sense approach tells us that an acceptable boundary condition might be
(1.12)
It is relatively easy to implement this. Simply store a value of Ex(l) for two time steps, and
then put it in Ex (0). Boundary conditions such as these have been implemented at both ends
of the Ex array in the program fdld_I.Z.c. (In the printout labeled fdld_I.Z.c at the end of the
chapter, the entire program hasn't been reproduced, only those parts which are different from
fdld_I.I.c. Furthermore, key points are highlighted in boldface.) Figure 1.3 shows the results
of a simulation using fdld_I.Z.c. A pulse that originates in the center propagates outward and
is absorbed without reflecting anything back into the problem space.
~~ 0.5
01------:-------:---
20 40 60 80 100 120 140 160 180
~~ 0.5
0
20 40 60 80 100 120 140 160 180
1 ~
I I I I I I I
~~ 0.5 \ T=250 ~ ) 0
20 40 60 80 100 120 140 160 180
FDTD cells
Figure 1.3 Simulation of an FDTD program with absorbing boundary conditions. Notice
that the pulse is absorbed at the edges without reflecting anything back.
PROBLEM SET 1.3
1. The program fdld_1.2.c has absorbing boundary conditions at both ends. Get this program
running and test it to ensure the boundary conditions are completely absorbing the pulse.
1.4 PROPAGATION IN A DIELECTRIC MEDIUM
In order to simulate a medium with a dielectric constant other than one, which corresponds to
free space, we have to add the relative dielectric constant e, to Maxwell's equations:
aE 1
-=-VxH ot e.e«
an 1
-=--VxE. at /-Lo
(1.13a)
(I.13b)
6 Chapter 1 • One-Dimensional Simulation with the FDTD Method
(1.14a)
(1.14b)
=---
Wewill stay with our one-dimensional example and make the change of variables in Eq. (1.5),
oEx(t) 1.oHy(t)
-o-t- =ErJEoJ-Lo 8z
oHy(t) 1 8Ex (t ) --=---.--,
ot JeoJ-Lo 8z
and then go to the finite difference approximations
E;+1/2(k) - i;-1/2(k) H;(k + 1/2) - H;(k + 1/2)
6.t erJeo~o ~x
H;+l(k + 1/2) - H;(k + 1/2) 1 E;+1/2(k + 1) - E;+1/2(k)
= -----------
t1t J-Lo Sx
From the previous section
1 6.t 1
---=-,
Jeo~o ~x 2
so Eq. (1.14) becomes
E;+1/2(k) = i;+1/2(k) + 1:2 [H;(k + 1/2) - H;(k + 1/2)]
H;+l(k + 1/2) = H;(k + 1/2) - ~[i;+l/2(k + 1) - i;+l/2(k)].
From these we can get the computer equations
(1.15a)
(1.15b)
ex [k]
hy[k]
ex[k] + cb[k] * (hy[k-l] - hy[k] )
hy [k] + O. 5 * (ex [k] - ex [k+1] ),
(1.16a)
(1.16b)
where
cb[k] = .S/epsilon (1.17)
over those values of k which specify the dielectric material.
The program fdId_I.3.c simulates the interaction of a pulse traveling in free space until
it strikes a dielectric medium. The medium is specified by the parameter cb in Eq. (1.17).
Figure 1.4 shows the result of a simulation with a dielectric medium having a relative dielectric
constant of 4. Note that a portion of the pulse propagates into the medium and a portion is
reflected, in keeping with basic EM theory [6].
PROBLEM SET1.4
1. The program fdId_I.3.c simulates a problem partly containing free space and partly dielectric
material. Run this program and duplicate the results of Fig. 1.4.
2. Look at the relative amplitudes ofthe reflected and transmitted pulses. Are they correct? Check
them by calculating the reflection and transmission coefficients. (See the appendix at the end
of this chapter.)
3. Still using a dielectric constant of 4, let the transmitted pulse propagate until it hits the far right
wall. What happens? What could you do to correct this?
40 60 80 100 120 140 160 180 200
I I I L_._.I_._.~._.~._._l_. __
T=440 Eps=4
I I I
40 60 80 100 120 140 160 180 200
FDTD cells
0.5
~>o(
0
-0.5
20
1 ~ I
0.5 ~
~><
°V·
-0.5
20
T=320 Eps=4
Figure 1.4 Simulation of a pulse striking a dielectric material with a dielectric constant of
4. The source originates at cell number 5.
1.5 SIMULATING DIFFERENT SOURCES
In the first two programs, a source is assigned a value to Ex; this is referred to as a hardsource.
In fdld_I.3.c, however, a value is added to Ex at a certain point; this is called a soft source.
The reason is that with a hard source, a propagating pulse will see that value and be reflected,
because a hard value of Ex looks like a metal wall to FDTD. With the soft source, a propagating
pulse will just pass through.
Up until now, we have been using a Gaussian pulse as the source. It is very easy to
switch to a sinusoidal source. Just replace the parameter pulse with the following:
pulse
ex [S]
sin[2*pi*fre~in*dt*T]
ex[S] + pulse.
The parameterfreq_in determines the frequency of the wave. This source is used in the program
fdld_I.4.c. Figure 1.5 shows the same dielectric medium problem with a sinusoidal source.
A frequency of 700 MHz is used. Notice that the simulation was stopped before the wave
reached the far right side. Remember that we have an absorbing boundary condition, but it
was only for free space!
Note that in fdld_l.4.c the cell size ddx and the time step dt are specified explicitly.
We do this because we need d t in the calculation of pulse. The cell size ddx is only specified
because it is needed to calculate dt from Eq. (1.7).
8 Chapter 1 • One-Dimensional Simulation with the FDTD Method
'[j= 150 Eps=4
. - . - . - . - . ~----"'~----------------i
-1
20 40 60 80 100 120 140 160 180
._'_'_0_._0_0_'_'_0_._.
T=42 Eps=4
~~ 0
-1
20 40 60 80 100 120 140 160 180
FDTD cells
Figure 1.5 Simulation of a propagating sinusoidal wave of 700MHz striking a medium
with a relative dielectric constant of 4.
PROBLEM SET 1.5
1. Modify your program fdld_1.3.c to simulate the sinusoidal source (see fdld_1.4.c).
2. Keep increasing your incident frequency from 700 MHz upward at intervals of 300 MHz. What
happens?
3. A type of propagating wave function that is of great interest in areas such as optics is the
"wave packet," which is a sinusoidal function in a Gaussian envelope. Modify your program
to simulate a wave packet.
1.6 DETERMINING CELLSIZE
(1.18)
Choosing the cell size to be used in an FDTD formulation is similar to any approximation
procedure: enough sampling points must be taken to ensure that an adequate representationis
made. The number of points per wavelengthis dependent on many factors [3, 4]. However,
a good rule of thumb is 10 points per wavelength. Experience has shown this to be adequate,
with inaccuracies appearing as soon as the sampling drops below this rate.
Naturally, we must use a worst-case scenario. In general, this will involve looking at
the highest frequencieswe are simulating and determiningthe correspondingwavelength. For
instance, suppose we are running simulations at 400 MHz. In free space, EM energy will
propagate at the wavelength
Ao = Co = 3 X 10
8
m/sec =.75 m.
400 MHz 4 x 108 sec- 1
If we were only simulating free space we could choose
~x =Ao/IO = 7.5 em,
(1.19)
However, if we are simulating EM propagation in biological tissues, for instance, we must
look at the wavelengthsin the tissue with the highest dielectricconstant,because this will have
the corresponding shortest wavelength. For instance, muscle has a relative dielectric constant
of about 50 at 400 MHz, so
Am = col.;so = .424 X 10
8 m/see = 10.6 em,
400 MHz 4 x 108 sec"!
and we would probably select a cell size of one centimeter.
Section 1.7 • Propagation in a Lossy Dielectric Medium
PROBLEM SET 1.6
1. Simulate a 3-GHz sine wave impinging on a material with a dielectric constant of e, =20.
1.7 PROPAGATION IN A LOSSYDIELECTRIC MEDIUM
9
So far, we have simulated EM propagation in free space or in simple media that are specified
by the relative dielectric constant e.. However, there are many media that also have a loss
term specified by the conductivity. This loss term results in the attenuation of the propagating
energy.
Once more we will start with the time-dependent Maxwell's curl equations, but we will
write them in a more general form, which will allow us to simulate propagation in media that
have conductivity:
aE
Ba-t =V x H-J
aH 1
-=--VxE. at f-Lo
J is the current density, which can also be written
]=u·E,
(1.20a)
(1.20b)
where a is the conductivity. Putting this into Eq. (1.20a) and dividing through by the dielectric
constant we get
aE 1 a
-=-VxH--E. at BoBr BoBr
We now revert to our simple one-dimensional equation:
aEx(t) = __1_ . oHy(t) _ ~Ex(t),
at BrBo aZ BrBO
and make the change of variable in Eq. (1.5), which gives
8ExCt)
-8-t- = __1_ . _8H_y_Ct_) _ _u-Ex(t)
Br-jBOf-LO 8z BrBO
CI.21a)
CI.21b)
(1.22)
8HyCt) 1 8ExCt)
-8-t- = - -jBof-Lo •az
Next take the finite difference approximations for both the temporal and spatial derivatives
similar to Eq. (1.3a):
E~+1/2(k) _ E~-1/2(k) H;(k + 1/2) - H;(k - 1/2)
=
ti.t Br-jBof-LO ti.x
a E~+1/2(k) + E;-1/2(k)
BrBO 2
Notice that the last term in Eq. (1.2Ia) is approximated as the average across two time steps
in Eq. (1.22). From the previous section
1 ti.t 1
-jBof-LO ti.x = 2'
(1.23a)
(1.23b)
10 Chapter 1 • One-Dimensional Simulation with the FDTD Method
so Eq. (1.22) becomes
i;+1/2(k) [1 + D.t · u] = i;-1/2(k) [I _~t · u]
2crcO 2erco
1/2
- -[H;(k + 1/2) - H;(k - 1/2)]
Cr
or
(
1 D.t · u)
£n+l/2(k) = - ~ £n-l/2(k) _ 1/2 [Hn(k + 1/2) - Hn(k - 1/2)].
x (1 + D.t · u) x e; .(1 + fl.t · u) y y
2crcO 2cxco
From these we can get the computer equations
ex [k] ca [k] *ex [k] + cb [k] * ( hy [k-l] - hy [k] )
hy [k] = hy [k] + O. 5 * ( ex [kl - ex [k+1] ) I
where
eaf = dt*sigma/ (2*epsz*epsilon) (1.24a)
ca[k] = (1. - eaf)/(l. + eaf) (1.24b)
cb [k] =0.5/ (epsilon* (1. + eaf)) . (1.24c)
Theprogramfd1d_I.5.c simulatesa sinusoidalwavehittinga lossymediumthat has a dielectric
constantof4 andaconductivityof 0.04. Thepulseis generatedatthefar left sideandpropagates
to the right. This is illustratedin Fig. 1.6. Notice that the waveform in the mediumis absorbed
before it hits the boundary, so we don't have to worry about absorbing boundary conditions.
-1
Eps=4
Cond = 0.04
20 40 60 80 100 120
FOTD cells
140 160 180
Figure 1.6 Simulation of a propagating sinusoidal wave striking a lossy dielectric material
with a dielectric constant of 4 and a conductivity of 0.04 (S/m). The source is
700 MHz, and originates at cell number 5.
PROBLEM SET1.7
1. Run program fdId_I.S.c to simulate a complex dielectric material. Duplicate the results of
Fig. 1.6.
2. Verify that your calculation ofthe sine wave in the lossy dielectric is correct, i.e., it is the correct
amplitude going into the slab, and then it attenuates at the proper rate. You can best do this by
writing a little program that calculates the parameters given in the appendix of this chapter.
3. How would you write an absorbing boundary condition for a lossy material?
4. Simulate a pulse hitting a metal wall. This is very easy to do, if you remember that metal has a
very high conductivity. For the complex dielectric, just use (J' = 1.e6, or any large number. (It
does not have to be the correct conductivity of the metal, just very large.) What does this do to
the FDTD parameters ca and cb? What result does this have for the field parameters Ex and
By? If you didn't want to specify dielectric parameters, how else could you simulate metal in
an FDTD program?
References
APPENDIX 1.A
11
When a plane wave traveling in medium 1 strikes a medium 2, the fraction that is reflected is
given by the reflection coefficient I', and the fraction that is transmitted into medium 2 is given
by the transmission coefficient r. These are determined by the intrinsic impedances 71 I and 712
of the respective media [6, p. 398]:
r = Ere! = 112 - 111 (1.A.l)
Einc 712 + 711
Etrans 2712
t' = - = (I.A.2) s.; 712 + 711
The impedances are given by
71= f-!i
BOB;
where B; is the complex relative dielectric constant
(I.A.3)
(1.A.5)
(I.A.4)
e*; a =Br + -.- ·
JWBo
For the case where JL =JLo, Eqs. (I.A.I) and (I.A.2) become
1 1
~-~ 1Ef-/Bf
r = _1_ + _1_ = IEf +/Bf
/Bf1Ef
2//Bf 2· IEf t'=----- _1_ + _1_ - IEf +/Bf'
/Bf1Ef
The amplitude of an electric field propagating in the positive z direction in a lossy
dielectric medium is given by
where Eo is the amplitude at z =O. The parameters ex and fJ are determined by the dielectric
constant Br » the conductivity a, and the radian frequency W = 21Tf of the propagating wave,
via the following two formulas [6, p. 420]:
a= (i) re; [ 1+ (_a)2 _1] 1/2 (Npjm)
CoV2: WBOBr
fJ = (i) re; [ 1+ (_a)2 + 1] 1/2 (rad/m),
CoV2: WBOBr
REFERENCES
(I.A.6)
(I.A.7)
[1] K. S. Vee, Numerical solution of initial boundary value problems involving Maxwell's equations in
isotropic media, IEEETrans. Antennasand Propagat., vol. 17, 1966, pp. 585-589.
[2] A. Taflove and M. Brodwin, Numerical solution of steady state electromagnetic scattering problems
using the time-dependent Maxwell's equations, IEEETrans. Microwave TheoryTech., vol. 23, 1975,
pp. 623-730.
12 Chapter 1 • One-Dimensional Simulation with the FDTD Method
[3] A. Taflove, Computation Electrodynamics: The Finite-Difference 'lime-Domain Method. Boston,
MA: Artech House, 1995.
[4] K. S. Kunz and R. J. Luebbers, The Finite Difference 'lime Domain Method/or Electromagnetics.
Boca Raton, FL; CRC Press, 1993.
[5] G. Mur, Absorbing boundary conditions for the finite-difference approximation of the time domain
electromagnetic field equations, IEEE Trans. Electromagn. Compat., vol. 23,1981, pp. 377-384.
[6] D. K. Cheng, Field and Wave Electromagnetics, Menlo Park, CA: Addison-Wesley, 1992.
C Programs 13