# Numerical Methods in Electrical Engineering

Urban and Civil

Nov 16, 2013 (4 years and 5 months ago)

90 views

Numerical Methods in Electrical
Engineering

Professor David Thiel

Centre for Wireless Monitoring and
Applications

Griffith University, Brisbane Australia

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.

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)

Topics

1.
Overview of Numerical Modelling for
electromagnetics.

2.
The Method of Moments.

3.
The Finite
-
Difference Time
-
domain
method.

4.
Review.

Topic 1

Overview of the Modelling Process

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.

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

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.

Numerical Modelling Procedure

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).

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

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

Solution

MoM and FDTD both use perfectly
conducting materials.

To introduce finite conductivity, we can
use lumped impedance elements.

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
pattern, near
-
field strength etc)

Numerical Modelling Procedure

code.

Find a simple analytical model and divide
into segments/pixels/voxels.

Solve the simple model analytically and
numerically.

Alter the segments/pixels/voxel size and
recalculate the numerical model.

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.

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.

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.

Topic 2

The Method of Moments

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.

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.

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.

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).

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.

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).

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.

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.

NEC Coordinate System

x

y

z

q

f

The
xoy

plane (
z

= 0) is where the ground plane is located if used.

Coding Examples

Comment Cards

CM

-

can only be at the start of
the program

CE

-

only one required.

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

Coding Cards

Program control cards

FR

frequency specification
-

0, # of steps, 0, 0, start
f
, D
f

GN

ground parameter specification
-

perfect gnd, 2 finite gnd

LD

-

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

NE

near electric field request

NH

near magnetic field request

RP

-

0, # of
q

steps, # of
f

steps, 1000,
q

start,
f

start, D
q
, D
f

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

NEC Activity

MMANA

http://www.smeter.net/antennas/mmana
-
tutorial.php

This is based on Mininec

a preNEC
program which is very efficient
computationally.

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)

Numerical Techniques in
Electromagnetics
, CRC Press. (Section 3.8,
Chapter 5)

Topic 3

Finite
-
Difference Time
-
Domain

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.

Background

By directly calculating electric and
magnetic fields through space, the
technique inherently includes all EM
phenomena including:

Surface waves

Near field values

Mutual coupling

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

Reflection coefficients

Effective dielectric constant

Background

The system of equations can be solved in
1D (a line), 2D (a surface) and 3D.

Applications

All electromagnetic systems including

EM Couplers

Antennas

UWB systems

Discontinuities in transmission lines
(microstrip lines, waveguides, etc)

electromagnetic geophysics

etc

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

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

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

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

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

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

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

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

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

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

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

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.

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

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

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

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

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

Workshop

Run the program FDTD1D.m

Read the program and identify the
components.

Change the two boundaries

Change the material properties

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
).

Student Evaluation of Course

Survey!

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.

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.