Copyright: David Thiel 2009
Numerical Methods in Electrical
Engineering
Professor David Thiel
Centre for Wireless Monitoring and
Applications
Griffith University, Brisbane Australia
Copyright: David Thiel 2009
Purpose
•
This is a follow

on workshop from the
Optimization workshop.
•
The aim is to introduce two methods of
electromagnetic modelling:
–
The Method of Moments (MoM): frequency domain,
open boundaries.
–
The Finite

Difference Time

Domain (FDTD): time
domain, closed boundaries.
•
To outline the validation processes essential to
obtaining reliable solutions.
Copyright: David Thiel 2009
My interest?
•
I wished to solve some fundamental problems in
electromagnetic geophysics where subsurface
targets were small and the earth was both highly
conductive and semi

infinite.
•
I have worked with many numerical EM
techniques including
–
MoM
–
FDTD
–
Impedance method
–
FEM
–
FDM
–
Spice (Circuit simulator)
Copyright: David Thiel 2009
Topics
1.
Overview of Numerical Modelling for
electromagnetics.
2.
The Method of Moments.
3.
The Finite

Difference Time

domain
method.
4.
Review.
Copyright: David Thiel 2009
Topic 1
Overview of the Modelling Process
Copyright: David Thiel 2009
Why numerical modelling?
•
To reduce the number of prototypes
constructed.
•
To use as the forward solver in an
optimization routine to gain the best
possible design.
Copyright: David Thiel 2009
What are the alternatives?
•
Design and build lots of physical models
–
Time consuming
–
Expensive
–
Measurement difficulties (cables, calibration,
mechanical precision, interference, mutual
coupling, field confinement, etc)
•
Analytical solutions
–
Limited to simple models because of the
computational complexity
Copyright: David Thiel 2009
But for Verification ….
•
You MUST make a physical model!
•
You MUST undertake some simple theoretical
modelling for use in your numerical modelling
analysis!
•
If you don’t, you have no verification of your
model and the results are not acceptable for
publication!
•
This is true EVEN if you use “proven”
commercial software.
Copyright: David Thiel 2009
Numerical Modelling Procedure
•
Select/write your numerical modelling
code.
•
Find a simple analytical model and divide
into segments/pixels/voxels.
•
These MUST have every dimension less
than
l
/10 where
l
is the wavelength in the
material being modelled. (Note that this is
NOT the free space wavelength).
Copyright: David Thiel 2009
Let’s do the calculation of
l
•
The complex propagation coefficient
g
is given by
the equation
where
s
is the conductivity,
e
is the absolute permittivity,
m
is the absolute magnetic permeability,
w
is the angular radiation frequency, and
j
is the
imaginary number,
a
is the attenuation constant
b
is the phase constant
)
(
we
s
wm
b
a
g
j
j
j
Copyright: David Thiel 2009
The wavelength in the material
•
At this frequency we have
•
The wavelength depths on the electromagnetic
properties of the material.
•
Remember NO segment/pixel/voxel can be
larger than in any dimension
l
/10
.
•
This can make life very difficult when you have
electrically large conductive materials.
b
l
2
Copyright: David Thiel 2009
Solution
•
MoM and FDTD both use perfectly
conducting materials.
•
To introduce finite conductivity, we can
use lumped impedance elements.
Copyright: David Thiel 2009
Modelling process
1.
Enter the model
2.
Pre

processing: Divide the model into
segments/pixels/voxels. You need to know the
coordinates of very line and corner.
3.
Apply the solver
4.
Post

processing: Extract the important
parameters (gain, front to back ratio, beam
width, bandwidth, input impedance, radiation
pattern, near

field strength etc)
Copyright: David Thiel 2009
Numerical Modelling Procedure
•
Select/write your numerical modelling
code.
•
Find a simple analytical model and divide
into segments/pixels/voxels.
•
Solve the simple model analytically and
numerically.
•
Compare your solutions.
•
Alter the segments/pixels/voxel size and
recalculate the numerical model.
Copyright: David Thiel 2009
Standard numerical modelling test
•
Recalculate use different sized
segments/pixels/voxels.
•
Do not use a factor of 2, 5 or 10.
•
Why?
–
To ensure that your modelling routine is robust
–
To assess what accuracy is needed to get good
results
–
To minimize the number of segments/pixels/voxels in
the model to reduce the computational time and the
memory requirements.
Copyright: David Thiel 2009
Numerical Modelling Procedure
•
Enter your full model and check results.
•
Alter the segments/pixels/voxel size and
recalculate the numerical model.
•
Only when you have robust answers can
you conclude that the model is correct.
•
It is often a good idea to re

enter the
model to ensure conductor connectivity,
no air gaps, etc.
Copyright: David Thiel 2009
Reducing the model size
•
Look for one or more lines of symmetry
and apply a boundary (perfect conductor
or perfect magnetic conductor).
•
One line of symmetry can halve the
computational effort.
•
For circular or cylindrical symmetry, you
can change a 3D problem to a 2D or 1D
problem.
Copyright: David Thiel 2009
Topic 2
The Method of Moments
Copyright: David Thiel 2009
Background
•
Numerical technique used to solve integral
equations.
•
The model common code is the Numerical
Electromagnetics Code (NEC) developed at
Lawrence Livermore Laboratories to check the
effect of EMP on antennas.
•
Originally created in the 1980’s, NEC2 is freely
available.
•
Now commercially available in more user

friendly form.
•
The method is always in 3D.
Copyright: David Thiel 2009
MoM Basics (1)
•
The model is constructed of infinitely thin,
perfectly conducting wires.
•
To construct a plane, you need a 2D
intersecting grid of wires.
•
To construct a volume, you need a 3D
intersecting grid of wires.
•
To add loss, each segment must contain a
lumped impedance resistance.
Copyright: David Thiel 2009
MoM Basics (2)
•
To account for finite radius, a perturbation
technique is used in the code
–
BUT, the
current on the wire is always axially
symmetric.
•
If this approximation is not good enough,
you need to construct a cylindrical grid of
wires to represent the conductor of finite
diameter.
Copyright: David Thiel 2009
MoM Basics (3)
•
You must add at least one source (voltage or
current) to one of the conducting segments to
provide energy to the system.
•
NEC2 provides an option for calculations in free
space or calculations in the vicinity of a ground
plane (infinite extent but finite or infinite
conductivity).
•
The ground plane calculation is done using
image theory and the Sommerfeld surface wave
(if selected).
Copyright: David Thiel 2009
MoM Basics (4)
•
The MoM computational engine places all
conductive segments in a matrix and calculates
the induced current in every segment using
mutual coupling equations.
•
The current on a single straight wire is fitted by a
basis function which smooths the current
transition between segments
•
The near electric field, the near magnetic field
and the electromagnetic far field (radiated field)
is calculated using the standard Hertzian dipole
formulation.
Copyright: David Thiel 2009
Pitfalls: Inaccurate results if
•
Wire segments are too long then incorrect
results (phase errors).
•
Wire segments are too short then incorrect
results (truncation errors).
•
Wire segments are too close and not
connected (axial symmetry is not a valid
approximation).
Copyright: David Thiel 2009
Pitfalls: Inaccurate results if
•
The feed point(s) is not on a conducting
segment.
•
The wires do not exactly meet. An air gap
is a major problem.
Copyright: David Thiel 2009
NEC code
•
The code was designed for computer card
entry.
•
Each line contains one piece of
information.
•
Most codes have a user

friendly interface
and the codes are created semi

automatically.
Copyright: David Thiel 2009
NEC Coordinate System
x
y
z
q
f
The
xoy
plane (
z
= 0) is where the ground plane is located if used.
Copyright: David Thiel 2009
Coding Examples
•
Comment Cards
–
CM
comments

can only be at the start of
the program
–
CE
end comments

only one required.
Copyright: David Thiel 2009
Coding Examples
•
Structure Geometry Cards
–
GA
wire arc specification
–
GE
end geometry

0 means no ground plane
–
GF
use numerical Green's function
–
GM
shift and duplicate structure
–
GR
generate cylindrical structure (symmetry)
–
GS
scale structure dimensions

0, 0, 1
–
GW
specify wire

number, number of segments in
wire, x1,y1,z1, x2, y2, z2, radius
–
GX
reflect structure
–
SP
specify surface patch
–
GH
generate helix
Copyright: David Thiel 2009
Coding Cards
•
Program control cards
–
FR
frequency specification

0, # of steps, 0, 0, start
f
, D
f
–
GN
ground parameter specification

0 for free space, 1 for
perfect gnd, 2 finite gnd
–
LD
structure impedance loading

0 series
RLC
or 1
parallel
RLC
, wire #, start seg #, end seg #,
R
,
L
,
C
–
EX
structure excitation

0, wire #, segment #, 00, volts
real, volts imag
–
NT
two

port network specification
–
TL
transmission line specification
–
EN
end of data flag
–
GD
additional ground parameter specification
–
NE
near electric field request
–
NH
near magnetic field request
–
RP
radiation pattern

0, # of
q
steps, # of
f
steps, 1000,
q
start,
f
start, D
q
, D
f
Copyright: David Thiel 2009
Coding Example
NEC Coding example: Dipole antenna in free space
CM Simple dipole antenna in Free Space
CM Optimized for resonance at 300 MHz
CE
GW 1, 9, 0,

.2418, 0, 0, .2418, 0, .0001
GS 0, 0, 1
GE 0
EX 0, 1, 5, 0, 1, 0
FR 0, 1, 0, 0, 300, 1
RP 0, 181, 1, 1000,

90, 0, 1, 1
RP 0, 1, 360, 1000, 90, 0, 1, 1
EN
Copyright: David Thiel 2009
NEC Activity
Copyright: David Thiel 2009
MMANA
•
http://www.smeter.net/antennas/mmana

tutorial.php
•
This is based on Mininec
–
a preNEC
program which is very efficient
computationally.
Copyright: David Thiel 2009
References
•
Iskander, M.F., 1992.
Electromagnetic fields and
waves
, Prentice Hall. (Section 4.9)
•
Harrington, R.F., 1993.
Field computation by
moment methods
, IEEE Press.
•
Guru, B.S., and Hiziroglu, H.R., 1998.
Electromagnetic field theory and fundamentals.
PWS Publishing Company. pp. 530

533.
•
NEC manual. (on the web)
•
Sadiku, N.O., 1992.
Numerical Techniques in
Electromagnetics
, CRC Press. (Section 3.8,
Chapter 5)
Copyright: David Thiel 2009
Topic 3
Finite

Difference Time

Domain
Copyright: David Thiel 2009
Background
•
Method was developed by Yee in 1966 to
perform a direct solution of Maxwell’s
equations.
•
Using a Finite

difference scheme for
Maxwell’s equations, the method
calculates the magnetic and then the
electric field components through space
containing the objects(s) after each time
step.
Copyright: David Thiel 2009
Background
•
By directly calculating electric and
magnetic fields through space, the
technique inherently includes all EM
phenomena including:
–
Surface waves
–
Far field values (Radiated waves)
–
Near field values
–
Mutual coupling
Copyright: David Thiel 2009
Background
•
By recording the time dependent field
strength at various places in the solution
space, it is possible to calculate
parameters such as:
–
Input impedance
–
Radiation pattern
–
Reflection coefficients
–
Effective dielectric constant
Copyright: David Thiel 2009
Background
•
The system of equations can be solved in
1D (a line), 2D (a surface) and 3D.
Copyright: David Thiel 2009
Applications
•
All electromagnetic systems including
–
EM Couplers
–
Antennas
–
UWB systems
–
Radar crossection
–
Discontinuities in transmission lines
(microstrip lines, waveguides, etc)
–
electromagnetic geophysics
–
etc
Copyright: David Thiel 2009
Basic formulation
•
The differential form of Maxwell’s
equations are solved using the central
difference formulation of the partial
differential operator.
•
If
u
is a field (
H
or
E
), then
x
x
x
u
x
x
u
dx
x
du
2
)
(
)
(
)
(
0
0
0
2
0
0
0
2
0
2
)
(
)
(
2
)
(
)
(
x
x
x
u
x
u
x
x
u
dx
x
u
d
Copyright: David Thiel 2009
3D PDE equations
x
z
y
x
H
y
E
z
E
t
H
'
1
m
y
x
z
y
H
y
E
z
E
t
H
'
1
m
z
y
x
z
H
y
E
z
E
t
H
'
1
m
x
y
z
x
E
z
H
y
H
t
E
s
e
1
y
z
x
y
E
x
H
z
H
t
E
s
e
1
z
x
y
z
E
y
H
x
H
t
E
s
e
1
Copyright: David Thiel 2009
2D equations
z
x
y
z
y
z
y
x
z
x
E
y
H
x
H
t
E
H
x
E
t
H
H
y
E
t
H
s
e
m
m
1
'
1
'
1
z
y
x
z
y
z
y
x
z
x
H
x
E
y
E
t
H
E
x
H
t
E
E
y
H
t
E
'
1
1
1
m
s
e
s
e
TM case
TE Case
Copyright: David Thiel 2009
1D equations
z
y
z
y
z
y
E
x
H
t
E
H
x
E
t
H
s
e
m
1
'
1
z
y
z
y
z
y
H
x
E
t
H
E
x
H
t
E
'
1
1
m
s
e
TM case
TE case
Copyright: David Thiel 2009
1D Difference Equations
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
H
x
H
H
t
E
E
H
x
E
E
t
H
H
s
e
m
1
1
1
1
1
'
1
Copyright: David Thiel 2009
Time step equations
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
H
t
E
E
x
t
H
E
t
H
H
x
t
E
m
m
e
s
e
'
1
1
1
1
1
1
Copyright: David Thiel 2009
•
The
E
and
H
fields are not calculated at exactly
the same point so the time at which the
calculations are performed must reflect this.
•
Chose
E
and
H
to be offset by
D
x
/2
and so the
time spacing between the
E
and
H
calculations
must be
D
t
/2
. We must therefore write the
equation as:
1
2
/
1
2
/
1
2
/
1
1
1
2
/
1
2
/
1
2
/
1
2
/
1
2
/
1
1
'
1
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
E
t
H
H
x
t
E
H
t
E
E
x
t
H
e
s
e
m
m
Copyright: David Thiel 2009
Leap

frog stepping
•
E
is calculated in the half spaces (
i
+1/2
) and at
integer time (
n
).
•
H
is calculated in the integer spaces (
i
) and at
half times (
n+1/2
).
1
2
/
1
2
/
1
2
/
1
1
1
2
/
1
2
/
1
2
/
1
2
/
1
2
/
1
1
'
1
n
i
n
i
n
i
n
i
n
i
n
i
n
i
n
i
E
t
H
H
x
t
E
H
t
E
E
x
t
H
e
s
e
m
m
Copyright: David Thiel 2009
Leap

frog stepping
–
update
equations
x
t
n
+1
n
+1/2
n
n

1/2
n

1
i

1
i

1/2
i
i
+1/2
i
+1
E
field
H
field
Copyright: David Thiel 2009
Source Field
•
This is a time domain method so need a
time domain source.
•
Simple source is a Gaussian pulse
•
The pulse is centred at time
n
0
and has a
1/e
decay constant of
n
decay
time steps. The
source field is given by:
2
0
/
)
(
0
decay
n
n
n
e
E
E
Copyright: David Thiel 2009
Gaussian pulse
0
10
20
30
40
50
60
70
80
90
100
0.2
0
0.2
0.4
0.6
0.8
1
1.2
Time step
E field
Gaussian pulse source
n
0
= 4n
decay
n
decay
= 10
Copyright: David Thiel 2009
Source Field
•
This can be applied at any time and any
E
field place.
•
An
H
field source can be used at an
H
field
place.
•
This is called a “hard source” as the field is
fixed by the equation.
•
Usually need to use a “soft source” which
includes energy from other nodes.
Copyright: David Thiel 2009
Soft Source Field
•
E
i
n
is calculated from the update
equations.
n
i
n
n
n
n
si
E
e
E
E
decay
2
/
1
/
)
(
0
2
/
1
2
0
Copyright: David Thiel 2009
Boundary Conditions
•
The FDTD method is a closed method.
•
It is essential to treat every boundary in
some way.
•
At PEC (perfect electric conductor) we
have:
•
In 1D we have
E
= 0
at the boundary.
incident
reflected
E
E
tan
tan
Copyright: David Thiel 2009
Boundary Conditions
•
The FDTD method is a closed method.
•
It is essential to treat every boundary in
some way.
•
At PMC (perfect magnetic conductor) we
have:
•
In 1D we have
H
= 0
at the boundary.
incident
reflected
H
H
tan
tan
Copyright: David Thiel 2009
Absorbing Boundary (ABC)
•
This is quite difficult to manage when the
angle of incidence is oblique.
•
Mur’s ABC in 2D
•
In 1D we have
E
i+1/2
n+1
=
E
i

1/2
n
•
Berenger used a multilayered PML
n
j
i
n
j
i
n
j
i
n
j
i
E
E
x
t
c
x
t
c
E
E
2
/
1
,
2
/
1
1
2
/
1
,
2
/
1
2
/
1
,
2
/
1
1
2
/
1
,
2
/
1
Copyright: David Thiel 2009
Courant Condition
•
Time step increments must be sufficiently
small to ensure a stable process.
•
In 2D
•
In 1D
2
2
2
1
1
1
1
z
y
x
c
t
2
c
s
t
3
c
s
t
Copyright: David Thiel 2009
Workshop
•
Run the program FDTD1D.m
•
Read the program and identify the
components.
•
Change the two boundaries
•
Change the material properties
Copyright: David Thiel 2009
Transform to Frequency domain
•
Record
E
(
t
)
at one point in the solution
space for the duration of the pulse (There
must be energy loss in the system either
through an ABC, PML or conductive
material).
•
E
(
w
)
the Fourier transform of
E
(
t
).
Copyright: David Thiel 2009
Student Evaluation of Course
•
Survey!
Copyright: David Thiel 2009
References
•
Taflove, A., 1995.
Computational Electrodynamics: the Finite

difference Time

domain
Method,
Norwood, MA: Artech House.
•
Kunz, K.S., and R.J. Luebbers, 1993.
The Finite Difference Time Domain Method for
Electromagnetics,
Boca Raton, FL: CRC Press.
•
Stutzman, W.L., & Thiele, G.A., 1998.
Antenna theory and design
, 2nd ed., Wiley
(Chapter 11)
•
Berenger, J.P., 1994.
A perfectly matched layer for the absorption of electromagnetic
waves
. J. Computational Physics, vol 114, pp. 185

200.
•
Isaacson, E., & Keller, H.B., 1967.
Analysis of numerical methods
. Wiley, New York.
•
Mur, G., 1981. Absorbing boundary conditions for the finite

difference approximation
of the time

domain electromagnetic field equations
. IEEE Trans. Electromagnetic
Compatability
, vol. 23, pp. 377

382.
•
Sadiku, M.N.O., 1992.
Numerical techniques in electromagnetics
, CRC Press, Boca
Raton. Section 3.8.
•
Thiel, D.V., & Mittra, R., 1997. Surface impedance modelling using the finite

difference time

domain method.
IEEE Trans. Geoscience & Remote Sensing,
vol 35,
pp. 1350

1356.
•
Yee, K.S., 1966. Numerical solution of initial boundary value problems involving
Maxwell's equations in isotropic media.
IEEE Trans
AP

14, pp 302

307.
Comments 0
Log in to post a comment