Numerical Methods in Electrical Engineering

clappergappawpawUrban and Civil

Nov 16, 2013 (3 years and 8 months ago)

67 views

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.