Finite-Difference Time-Domain Method

1
Finite-Difference Time-Domain Method

• Dennis Sullivan, Ph.D.
• Professor of Electrical Engineering
• University of Idaho
• Moscow, ID USA
• 83844-1023

2
Outline

• FDTD for free space in 1D
• FDTD for biological tissues in 1D
• Calculating the SAR
• FDTD formulation in 3D
• Boundary conditions (the PML)
• Simulation of a dipole antenna
• Modeling biological tissues

3
Outline (continued)
Advanced Topics

• Interpolation across boundaries
• Frequency dependent materials

4
Electromagnetic radiation is governed by the
Maxwells equations
5
In one dimension in free space they become
6
To put these equations in a computer, take the
finite-difference approximations of the partial
derivatives in time and space
7
This is a time-domain method. Each new value of
the electric field E or the magnetic field H is
determined by the previous values
8
The k represents the location in an array in a
computer while n represents time
9
We make a change of variables so E and H have
the same order of magnitude
Once the cell size is chosen, the time steps must
be chose small enough for stability
10
This results in the following two equations of
code in the C program language
exk exk 0.5( hyk-1 - hyk) hyk
hyk 0.5( exk - exk1)
11
nn1
Calculate En1/2
Calculate Hn1
Each time step represents an increment in the
total time T n Dt.
12
The following is a one-dimensional simulation of
an EM pulse propagation in free space. T
represents the number of time steps.
13
(No Transcript)
14
(No Transcript)
15
(No Transcript)
16
(No Transcript)
17
(No Transcript)
18
(No Transcript)
19
(No Transcript)
20
(No Transcript)
21
(No Transcript)
22
(No Transcript)
23
(No Transcript)
24
Media like those found in human tissue
are specified by 1. Relative dielectric
constant 2. Conductivity
25
These are included in the FDTD formulation
Note that the last term is written as the
average over two cells
26
This leads to the following C computer code
Specify the parameters in the cells
eaf dtsigma/(2epszepsilon) cak (1. -
eaf)/(1. eaf) cbk 0.5/(epsilon(1.
eaf)).
Computer code in the main loop
exk cakexk cbk ( hyk-1 - hyk
) hyk hyk 0.5( exk - exk1 )
27
The following simulation shows an EM pulse
propagating in free space and then striking a
material with e 5, s 0.05. (Approximately the
values for human fat or bone.)
28
The cells between 100 and 200 have been assigned
the properties
29
(No Transcript)
30
(No Transcript)
31
(No Transcript)
32
(No Transcript)
33
(No Transcript)
34
(No Transcript)
35
(No Transcript)
36
(No Transcript)
37
(No Transcript)
38
(No Transcript)
39
(No Transcript)
40
(No Transcript)
41
(No Transcript)
42
Most EM sources, like those used in
hyperthermia, produce sinusoidal radiation
43
(No Transcript)
44
(No Transcript)
45
(No Transcript)
46
(No Transcript)
47
(No Transcript)
48
(No Transcript)
49
(No Transcript)
50
(No Transcript)
51
(No Transcript)
52
(No Transcript)
53
(No Transcript)
54
(No Transcript)
55
(No Transcript)
56
(No Transcript)
57
(No Transcript)
58
(No Transcript)
59
(No Transcript)
60
(No Transcript)
61
Once steady state has been reached, the rate of
absorption of energy in the tissue is determined
by the specific absorption rate (SAR)
62
However, detemining SAR in this manner is
difficult for two reasons 1. Sinusoidal
sources are difficult 2. Infromation can only
be obtained for one frequency at a
time.
63
The best method to determine the magnitude of the
E field for a sinusoidal source at a certain
frequency Use a pulse source in the simulation
and then take the discrete Fourier transform at
each cell, at each cell at each frequency of
interest.
64
Remember, this must be done at each cell where
the SAR is to be known. In a 3D simulation, this
is typically many thousands of cells. It would
be impossible to store the time-domain data and
then take the Fourier transforms.
65
However, I can calculate a running Fourier
transform during the simulation, and it only
requires two more computer words per cell, per
frequency
At the end of the simulation, the amplitude
is calculated from
66
The following simulation shows an EM wave
interaction with a section of dielectric
material. The results are calculated for two
different frequencies, 50 MHz and 500 MHz. Note
that a pulse can be used to determine what the
results will be for sinusoidal sources.
67
(No Transcript)
68
(No Transcript)
69
(No Transcript)
70
(No Transcript)
71
(No Transcript)
72
(No Transcript)
73
(No Transcript)
74
(No Transcript)
75
(No Transcript)
76
(No Transcript)
77
(No Transcript)
78
(No Transcript)
79
(No Transcript)
80
The last slide shows the very different patterns
that can occur at different frequencies having
different wavelengths.
81
The Maxwells equations can also be written
The medium characteristics are in the middle
equation
82
nn1
Calculate Dn1/2
Calculate En1/2 from Dn1/2
Calculate Hn1
This adds another step to each increment
83

• There are two main reasons for using this
  formulation
• It is easier to formulate frequency-dependent
  media (We will discuss this under advanced
  topics.)
• It is easier to formulate the perfectly matched
  layer (PML) at the boundaries.

84
Three-dimensional Simulation
The Yee Cell
The E and H fields are interwoven
85
In the FDTD formulation, the Maxwell
equations become six interwoven field calculations
as well as three equation to get the E field from
the flux densities in the x, y, and z directions
86
Up until now, we have not discussed the boundary
conditions at the edges of the problem space.
These are necessary to keep unwanted reflections
from coming back.
87
Probably the best solution is the
perfectly matched layer (PML) which
absorbs out-going waves.
88
The reflection of an outgoing wave is
determined by the reflection coefficient
89
There are two conditions to form a PML

• The impedance going from the background
• medium to the PML must be constant

2. The direction perpendicular to the boundary,
the x direction for example, must be the inverse
of the other directions
90
To implement the PML, we will assume that there
are fictitious values of e and m that we can
attach to the Maxwells equations.
91
We assume e and m are complex. It is the
imaginary part that leads to absorption
92
The following selection of parameters satisfies
these requirements
93
An outgoing wave sees a constant impedance as its
going into the PML. The conductivity causes it
to be absorbed ounce its in the PML.
The values of s are gradually increased as they
go into the PML
94
The following two-dimensional simulation shows
the radiation from a point source. The radiation
is absorbed by the PML.
95
(No Transcript)
96
(No Transcript)
97
(No Transcript)
98
(No Transcript)
99
(No Transcript)
100
(No Transcript)
101
(No Transcript)
102
(No Transcript)
103
(No Transcript)
104
(No Transcript)
105
Remember that the fictitious values of e and m
that are used to implement the PML have nothing
to do with the real values of e that specify
the body being simulated. Therefore, they can
overlap.
106
Applicators can be simulated in FDTD by
specifying the material and the source of energy.
Here is a simple dipole antenna.
107
The portion that is to be metal can be specified
by just holding the E field to zero.
The excitation comes from just specifying the E
field in the gap.
The FDTD method will determine the H field, which
is an indication of the current in the dipole
108
It also calculates the E field that would radiate
out from the antenna.
109
The following simulation is from a
three-dimensional simulation of a dipole
radiating in free space. The first set of
slides shows the H fields next to the metal of
the dipole arms, which are an indication of the
current.
110
(No Transcript)
111
(No Transcript)
112
(No Transcript)
113
(No Transcript)
114
(No Transcript)
115
(No Transcript)
116
(No Transcript)
117
(No Transcript)
118
(No Transcript)
119
(No Transcript)
120
(No Transcript)
121
(No Transcript)
122
(No Transcript)
123
(No Transcript)
124
The next set of slides shows the E field in the
plane of the gap as it radiates away from the
antenna
125
(No Transcript)
126
(No Transcript)
127
(No Transcript)
128
(No Transcript)
129
(No Transcript)
130
(No Transcript)
131
(No Transcript)
132
(No Transcript)
133
(No Transcript)
134
(No Transcript)
135
(No Transcript)
136
(No Transcript)
137
(No Transcript)
138
(No Transcript)
139
(No Transcript)
140
(No Transcript)
141
(No Transcript)
142
(No Transcript)
143
(No Transcript)
144
(No Transcript)
145
Each cell in an FDTD simulation can be specified
as a different tissue type.
146
FDTD does not require an elaborate mesh to
specify the boundaries. Each cell is composed of
a material.
147
If a simple in or out approach is used, a
stair casing effect results.
148
This can be improved by decreasing the cell size,
but that will require more computer resources.
149
Another possibility is to average across the
cells. This gives a better representation
without increasing computer resources.
150
The applicators and the body to be radiated can
be included in the simulation.
151
(No Transcript)
152
(No Transcript)
153
(No Transcript)
154
(No Transcript)
155
Advanced Topics
1. Interpolation to improve accuracy
2. Frequency dependent FDTD formulation
156
The basic FDTD method assumes that the E field is
perpendicular from the plane containing the H
fields.
Ez
Hy
Hx
Hx
Hy
157
In the vicinity of dielectric boundaries, the
actual E field could be substantially different.
Ez
Ez(actual)
Hx
Hx
Hy
158
A High-Resolution Interpolation at Arbitrary
Interfaces for the FDTD Method
J. Nadobny, D. Sullivan, P. Wust, M. Seebass, P.
Dueflhard, and R. Felix
IEEE Transactions on Microwave Theory and
Techniques, Vol. 46, Novmeber 1998.
159
(No Transcript)
160
At the end of the simulation, the FDTD value of
each field, is corrected by the second
term on the right.
161
A program simulates the radiation of a layered
sphere with a plane wave. The results can be
compared with the analytic calculate of the E
fields using a Bessel function expansion method.
162
Accuracy of Ez field simulation on the 45 degree
axis.
163
Accuracy of Ey field simulation on the 45 degree
axis.
164
Frequency dependent methods
The medium characteristics are in the middle
equation
165
Biological tissues can have properties that vary
at different frequencies. The following table
shows the values for human muscle.
Frequency (MHz) er s
166
The dielectric constant and conductivity can
be written as a complex dielectric constant
However, most tissues are frequency dependent and
have one or more additional terms
This could not be incorporated in the previous
method
167
The frequency dependent term must be taken to the
time-domain where it becomes a convolution
This can be formulated in FDTD by the following
168
(No Transcript)
169
The End