5:.

lfl"./"1' “

‘ . “rainy
‘ "3-5.: V;

v 5 5 - "V;
Zv'zgki'fiztfig Ev "

 

 

..
V. 613-.

’1'.‘ "I—

.. .,. ....o.vvs,r -=<‘:<t'r'-"-"

r
v.
‘v

I,”

[ |
~‘m v..
¥} 1)

.J, ..
. . " - rv‘v ‘
"333:" J 51". " ‘
H , . m .' ,LA r ‘
23% . V» savvy“: 1W???“
- 3%)?) M509»; ’ ‘ 33¢“
. ‘10 55556 and“ .

p- .)aA'I ' I "x"?

v .
.h ‘ -, “J '7 Y"“’
I.- . , t‘II‘IH. ‘J ' .01; L'l'ww"va~C€r\L.",“:r?\-d"g
V‘ I.a ‘nh'm'xax‘w ‘\-\’\‘:J. ’ v, '
.u n u' ". ; 0‘6“.
. a.»
x

. \ .
“it“ 3*: <¢

"A

- 1.
1% z 5.
“w u; 3:? ~

=5?
5 mm
(1’ 1

IL
,_ A A.
"‘:5‘¢:2ii,:3~:‘u n
‘ (‘3’;
15‘.

5%.!
TI
V v
5!" ‘;£;::..‘:%
,v. my; '
a a" \h.

g A
’l .‘ .‘
1' m
7 ,‘uX‘j 3%)
. . lg' .HH
4-: " I}
.

.
ma

1

F .1'
*1?“

!~

7.

p «
Sig-3%

‘1‘ ~ 4‘4
‘31: Hi‘sva‘.

I

up .,y <
3“"77‘“ .

u

:4 $$

.vw u 1"“:
_ 33:“ \ 3 Vi
\_v \

1|? R‘.

6

N“ fivr‘
‘9 . "u z

.1533). ‘
' . jwy 5‘

i :33"

w, .‘ 7 ‘, ~ {m
5%? we»

4.

m '-'r 18' a,"
. ‘ “n 1:»?
. wfix‘ft‘b‘fgzfiééim
1.}, Xng‘“
I ' "far". I“ _ '

m .3 ..

'._ >3! J'- , , _ .
‘$%§"ag Ti“
,5 1 v. -' ‘
.‘~ I: ‘8‘"; ..I l {J I." .‘
oug- $3K§ 511,: -. ‘
I? §$§
95%}:

I"

. , w‘
\. - 22“"
135313371”, fi‘é‘y":

‘ V
. #"é “m
{95: 22‘": e'

 

 

 

Ci WHG73WOO

lllllilillliiilll‘iil‘iiliiliiliitilliiiiliill‘iil 1/

3 1293 00612 5763

 

 

 

 

 

 

 

This is to certify that the

dissertation entitled

 

lIBRARY

Michigan State
University

 

Numerical Solutions of Electromagnetic
Problems by Integral Equation Methods and
Finite-Difference Time-Domain Method

presented by

Xiaoyi Min

has been accepted towards fulfillment

of the requirements for

Ph.D degree in Electrical Engineering

 

 

/é/%

 

Major professor

Date 8/! /;7

MS U i: an Affirmative Action/Equal Opportunity Institution

0-12771

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

  

DATE DUE DATE DUE DATE DUE

 

 

 

HM

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(if:

MSU Is An Aflirmdive ActiorVEqual Opportunity Institution

 

 

 

NUMERICAL SOLUTIONS OF ELECTROMAGNETIC
PROBLEMS BY INTEGRAL EQUATION METHODS AND

FINITE-DIFFERENCE TIME-DOMAIN METHOD

By

Xiaoyi Min

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering

1989

ABSTRACT

NUMERICAL SOLUTIONS OF ELECTROMAGNETIC
PROBLEMS BY INTEGRAL EQUATION METHODS AND
FINITE-DIFFERENCE TIME-DOMAIN METHOD

By
Xiaoyi Min

This thesis first presents the study of the interaction of electromagnetic waves
with three-dimensional heterogeneous, dielectric, magnetic, and lossy bodies by surface
integral equation modeling. Based on the equivalence principle, a set of coupled sur-
face integral equations is formulated and then solved numerically by the method of
moments. Triangular elements are used to model the interfaces of the heterogeneous
body, and vector basis functions are defined to expand the unknown current in the for-
mulation. The validity of this formualtion is verified by applying it to concentric
spheres for which an exact solution exists. The potential applications of this formual-

tion to a partially coated sphere and a homogeneous human body are discussed.

Next, this thesis also introduces an efficient new set of integral equations for
treating the scattering problem of a perfectly conducting body coated with a thin mag-
netically lossy layer. These electric field integral equations and magnetic field integral
equations are numerically solved by the method of moments (MOM). To validate the
derived integral equations, an alternative method to solve the scattering problem of an
infinite circular cylinder coated with a thin magnetic lossy layer has also been

developed, based on the eigenmode expansion. Results for the radar cross section and

current densities via the MoM and the eigenmode expansion method are compared.
The agreement is excellent. The finite difference time domain method is subsequently
implemented to solve a metallic object coated with a magnetic thin layer and numeri-

cal results are compared with that by the MoM.

Finally, this thesis presents an application of the finite-difference time-domain
approach to the problem of electromagnetic receiving and scattering by a cavity-backed
antenna situated on an infinite conducting plane. This application involves
modifications of Yee’s model, which applies the difference approximations of field
derivatives to differential operators in the Maxwell’s curl equations, and applies the
radiation boundary condition on a truncated boundary surface. The modifications are
based on the integral forms of the Maxwell equations and image theory. The effects
of an infinitely thin impedance sheet on receiving and scattering characteristics of the

antenna are investigated.

ACKNOWLEDGMENTS

The author wishes to express her heartfelt appreciation to her academic advisor,
Dr. Kun-Mu Chen, for his sincere guidance and encouragement throughout the course
of this work. She also wishes to appreciate Dr. D. P. Nyquist and Dr. E. J. Rothwell

for their generous support and constructive suggestions during the period of this study.

A special note of thanks is due Dr. Byron C. Drachman for his special assistance

and concern to the author.

In addition, the author thanks her parents, Mr. and Mrs. Longqiu Min, and her
parents-in-law, Mr. and Mrs. Jiantang Sun for their continued love, sacrifice and

encouragement in her research years.

The research reported in this thesis was supported by the National Science Foun-

dation under Grant No. ECE 8509213.

iv

Table of Contents

LIST OF FIGURES .................................................................................................

CHAPTER 1 INTRODUCTION ...........................................................................

CHAPTER 2 INTERACTION OF ELECTROMAGNETIC FIELDS

WITH THREE-DIMENSIONAL, HETEROGENEOUS,
DIELECTRIC, MAGNETIC AND LOSSY
BODIES ............................................................................................

2.1 Introduction ..............................................................................................

2.2 Derivation of Coupled Surface Integral Equations for
a Heterogeneous Body ...........................................................................

2.2.1 The Preliminary Theorems ........................................................

2.2.2 A Heterogeneous Body Without Any Perfect
Conductors Inside ........................................................................

2.2.3 A Heterogeneous Body Enclosing a Perfect Conductor ..........
2.2.4 Matellic Body with Partial Coating ..........................................
2.3 Numerical Algorithms .............................................................................
2.3.1 Basis Function .............................................................................
2.3.2 Properties of Vector Basis Function ...........................................
2.3.3 Testing Procedure and Matrix Equation .....................................
2.4 Numerical Results ....................................................................................
2.5 Comparision with Other Methods ...........................................................
2.6 Some Comments ......................................................................................
2.7 Conclusion ................................................................................................

CHAPTER 3 FUNDAMENTALS OF FINITE DIFFERENCE TIME

DOMAIN METHOD IN SOLUTIONS OF
ELECTROMAGNETIC SCATTERING AND
ANTENNA PROBLEMS ...............................................................

3.1 Introduction .............................................................................................

3.2 BASIC FD-TD ALGORITHMS ............................................................
3.2.1 Maxwell’s Curl Equations .........................................................
3.2.2 Discretization of the scalar Maxwell equations

viii

10
18
21
28
28
3O
32
43
71
71
72

74
74
77
77

by using Yee’s model ................................................................. 79

3.3 Radiation Boundary Conditions ............................................................. 85
3.3.1. Derivation of Radiation Boundary Condition by
Wave Equation ........................................................................... 88
3.3.1.a. The First method ............................................................ 88
3.3.l.b. The Second Method ....................................................... 94
3.3.2 Finite-difference Approximation of Radiation
Boundary Conditions ................................................................... 95
3.3.2.a. Three Dimensional case ................................................. 95
3.3.2.b. Two Dimensional case .................................................. 101
3.3.3 Radiation Boundary Condition for the 2D and
3D Comer Points ....................................................................... 105
3.4 The Stability Analysis of the FD-TD Method Schemes ....................... 108
3.4.1 The Methods of Stability Analysis ............................................ 110
3.4.2 Summary ...................................................................................... 118
3.5 Total Field Region and Scattered Field Region .................................... 118
3.6 Integral Interpretation of the FD-TD Algorithms ................................. 120
3.7 Numerical Irnplememtation and Validation of FD-TD Modeling
.......................................................................................................................... 123
CHAPTER 4 APPLICATION OF FINITE DIFFERENCE TIME
DOMAIN METHOD AND MOMENT METHOD TO
TWO DIMENSIONAL ELECTROMAGNETIC
SCATTERING PROBLEMS ......................................................... 132
4.1 Introduction ............................................................................................. 132
4.2 Integral Equations for Perfectly Conducting Cylinders
Coated with Thin Magnetic Materials ................................................... 133
4.2.1 Derivation of Integral Equations for three-dimensional
structures .................................................................................... 133
4.2.2 Integral Equations for Two-Dimensional Geometries .............. 141
4.2.3 Some Properties of New Surface Integral Equations ............... 148
4.3 Eigen-mode Expansion Solution to an Infinitely Long
Circular Cylinder .................................................................................... 149
4.4 Finite Difference Time Domain Method ............................................... 154
4.5 Numerical Results ................................................................................... 161

4.5.1 Perfectly Conducting Circular Cylinder Coated with

Magnetic Thin Layer ................................................................. 161

4.5.1.a Completely Coated Cylinder ......................................... 161
4.5.1.b Partially Coated Cylinder .............................................. 162
4.5.2 Perfectly Conducting Square Cylinder Coated with
Magnetic Thin Layer ................................................................. 163
4.6 Extension to Three Dimensional Case ................................................... 164

CHAPTER 5 EFFECTS OF AN IMPEDENCE SHEET ON THE

CHARATERISTIC PROPERTIES OF CAVITY
BACKED ANTENNA BY FINITE DIFFERENCE

TIME DOMAIN METHOD .......................................................... 182

5.1 Introduction ............................................................................................. 182

5.2 Basic Finite Difference Time Domain Scheme ..................................... 184

5.3 Integral Interpolation of Maxwell’s Curl Equation ............................... 184

5.4 Modification of Radiation Boundary Condition .................................... 196

5.5 Fields in the Scattered Field Region ..................................................... 202

5.6 Backseattered Fields ............................................................................... 202

5.7 Numerical Results ................................................................................... 206

5.7.1 Open Cavity situated in the Ground Plane ............................... 206

5.7.2 Cavity Backed Antenna ............................................................. 207

CHAPTER 6 SUMMARY ..................................................................................... 228
APPENDIX A ........................................................................................................... 230
APPENDIX B ............................................................................................................ 232
APPENDIX C ........................................................................................................... 239
BIBLIOGRAPHY ..................................................................................................... 243

vii

2.1
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
2.10

2.11

2.12
2.13

2.14

2.15

2.16

2.17

List of Figures

Geometry of a heterogeneous body with a plane wave excitation ................ 8
Geometry for illustration of the equivalence principles ................................. 9
A heterogeneous body without any perfect conductors ................................. 12
A heterogeneous body with a perfect conductor inside ................................. 19
A perfect conductor partially coated with lossy material .............................. 23
Sources maintaining fields in region 1 ........................................................... 24
Sources maintaining fields in region 2 ........................................................... 25
Local coordinates associated with an edge ..................................................... 29
Geometry for normal component of basis function at common edge ........... 31
Coordinate for calculating centroids and moments of basis

vectors .............................................................................................................. 33
Local coordinates and edges for source triangle 1" with

observation point in triangle 7” ....................................................................... 40
Definitions of areas used in defining area coordinates .................................. 42

O-component of equivalent magnetic current on the outer surface

SI of a concentric sphere with: (koaz = 0.0595; 8,2 = 16.0;

pa = 4.0; tan(57) = 0.39 ) ( koal = 0.13: En = 9.0: uri = 1.0;

tan(8,) = 0.0) ..................................................................................................... 45

¢-component of equivalent magnetic current on the outer surface

51 of a concentric sphere with: ( koaz = 0.0595; 8,; = 16.0;

u,2 = 4.0; tan(62) = 0.39 ) ( koal = 0.13; 8,1 = 9.0; u,1 = 1.0;

tan(81) = 0.0 ) .................................................................................................... 46

e-component of equivalent electric current on the outer surface

S, of a concenuic sphere with: ( koaz = 0.0595: 8,2 = 16.0;

1,1,2 = 4.0. tan(87) = 0.39) ( koal = 0.13; 6,1 = 9.0; llri = 1.0;

[311(51) = 0.0 ) .................................................................................................... 47

¢~component of equivalent electric current on the outer surface

S, of a concentric sphere with: ( k002 = 0.0595; 8,; = 16.0;

[in = 4.0; tan(87_) = 0.39 ) ( koa, = 0.13; 8,1 = 9.0: “71 = 1.0;

tan(5,) = 0.0 ) .................................................................................................... 48

G-component of equivalent magnetic current on the inner surface

52 of a concentric 'sphere with: ( koaz = 0.0595; 8,2 = 16.0;

Ila = 4.0; tan(87) = 0.39 ) ( [coal = 0.13; Eri = 9.0; llrr = 1.0;

tan(51) = 0.0 ) .................................................................................................... 49

viii

2.18

2.19

2.20

2.21

2.22

2.23

2.24

2.25

2.26

2.27

2.28

2.29

¢-component of equivalent magnetic current on the inner surface

S; of a concentric sphere with: ( koa; = 0.0595; 6,; = 16.0;

Ila = 4.0; tan(6;) = 0.39 ) ( koa, = 0.13; 8,, = 9.0; 1171 = 1.0;

tan(8,) = 0.0 ) ..................................................................................................

e-component of equivalent electric current on the inner surface

S; of a concentric sphere with: ( koa; = 0.0595; 8,; = 16.0;

ufl = 4.0; tan(8;) = 0.39 ) ( koa, = 0.13; 8,, = 9.0; 11,, = 1.0;

tan(8,) = 0.0 ) ....................................................................................................

¢~component of equivalent electric current on the inner surface

S; of a concentric sphere with: ( koa; = 0.0595; 6,; = 16.0;

u,; = 4.0; tan(8;) = 0.39 ) ( koa, = 0.13; 6,, = 9.0; ”71 = 1.0;

tan(6,) = 0.0 ) ....................................................................................................

e-component of equivalent magnetic current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.
( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 1.0; tan(5) = 0.0 ) ...............................

«ti-component of equivalent magnetic current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.
( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 1.0; tan(8) = 0.0 ) ...............................

G-component of equivalent electric current on the outer surface

S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 1.0; tan(8) = 0.0 ) ...............................
«ti-component of equivalent elecuic current on the outer surface

S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 1.0; tan(5) = 0.0 ) ...............................
G-component of electric current on the inner surface

S; of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:: koa, = 0.51:; e, = 1.0; u, = 1.0; tan(5) = 0.0 ) ...............................

tit-component of electric current on the inner surface
S; of a perfectly conducting sphere coated with a lossy layer.
( koa; = 0.41:; koa, = 0.51:; e, = 1.0; it, = 1.0; tan(5) = 0.0 ) ...............................

O-component of equivalent magnetic current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.
( koa; = 0.41:; koa, = 0.51:; e, = 4.0; u, = 1.0; tan(8) = 0.0 ) ...............................

o-component of equivalent magnetic current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.
( koa; = 0.41:; koa, = 0.51:; e, = 4.0: u, = 1.0; tan(8) = 0.0 ) ...............................

B-component of equivalent electric current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.

ix

50

51

52

53

54

55

56

57

58

59

60

2.30

2.31

2.32

2.33

2.34

2.35

2.36

2.37

2.38

3.1
3.2
3.3
3.4
3.5
3.6
3.7
3.8

( koa; = 0.41:; koa, = 0.51:; e, = 4.0; u, = 1.0; tan(8) = 0.0 ) ..............................

¢~component of equivalent electric current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 4.0; u, = 1.0: tan(8) = 0.0 ) ..............................

6-component of electric current on the inner surface
S; of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; c, = 4.0; u, = 1.0; tan(8) = 0.0 ) ..............................

¢-component of electric current on the inner surface
S; of a perfectly conducting sphere coated with a lossy layer.

( [(002 = 0.43; kOal = 0.53; 8, = 4.0; Ll, = 1.0; {311(5) = 0.0 ) ..............................

O-component of equivalent magnetic current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 4.0; tan(5) = 0.0 ) ..............................

¢-component of equivalent magnetic current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:: koa, = 0.51:: e, = 1.0; u, = 4.0; tan(8) = 0.0 ) ..............................

O-component of equivalent electric current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 4.0; tan(8) = 0.0 ) ..............................

¢-component of equivalent electric current on the outer surface
S, of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; ll, = 4.0; tan(5) = 0.0 ) ..............................

O-component of electric current on the inner surface
S; of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; ll, = 4.0; tan(8) = 0.0 ) ..............................

o-component of electric current on the inner surface
S; of a perfectly conducting sphere coated with a lossy layer.

( koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 4.0; tan(5) = 0.0 ) ..............................

Lattice unit cell ---- Yee’s model .....................................................................
Two-dimensional lattice unit cell .....................................................................
Lattice unit cell near outrnost surface ( y=0 plane ) ........................................
Graphic interpretation of Eq. 3.3.9b .................................................................
Points used in the Mur’s second order approximation ....................................
Two-dimensional Mur’s second order approximation .....................................
Total field and scattered field regions ..............................................................
Examples of spatially orthogonal contours in free space. (a) Ampere’s

. 61

62

63

64

65

66

67

68

69

70

81
84
87
91
96
102
119

3.9

3.10

3.11

3.12
3.13

3.14

3.15

4.1

4.2
4.3

4.4

4.5

4.6

4.7

4.8

4.9

4.10

law for E,; (b) Farady’s law for H2 . .............................................................

Current distribution on an infinite long, perfectly conducting
rectangular cylinder without coating in the case of the TM
excitation ( koa = 21: ) .......................................................................................

Comparison of results on the induced current on a perfectly conduct-
ing strip using integral equations and the FD-TD method ............................

Induced current on a strip of thin film with variable conductivity:
n = 1/(orZo)= 2(x/a)2 ........................................................................................

A perfectly conducting wing-shaped cylinder ...............................................

Amplitude distribution of the z-component electric field in the
total field region with koa = 20 ........................................................................

Amplitude distribution of the y-component magnetic field in the
total field region with koa = 20 ........................................................................

Amplitude distribution of the x-component magnetic field in the
total field region with koa = 20 ........................................................................

The illustration of equivalent currents and equivalent charges .......................
Total equivalent charges on a thin layer ..........................................................

Separation of the total fields into incident fields and scattered
fields maintained by equivalent sources ..........................................................

Application of the surface E-field integral equation on the surface
of a rectangular cylinder ...................................................................................

The cross section of a circular cylinder coated with a magnetically
lossy layer .........................................................................................................

Truncation of an infinite two-dimensional space into a finite
region .................................................................................................................

A basic cell of the Yee’s model adjacent to the coating with TM
excitation ...........................................................................................................

A basic cell of the Yee’s model adjacent to the coating of a cylinder
with TE excitation ............................................................................................

Radar cross sections of an infinitely long conducting circular cylinder
with or without a magnetic coating in the case of TE excitation
( koa = 51:. u,t/a = 0.01 - 10.03 ) ........................................................................

Amplitude distribution of the e-component current on the surface

of an infinitely long conducting circular cylinder with or without a
magnetic coating in the case of TE excitation

( Icoa = 51:, tut/a = 0.01 - j0.03 ) ........................................................................

xi

122

125

126

127
128

129

130

131

135
135

137

146

150

156

158

160

165

166

4.11

4.12

4.13

4.14

4.15

4.16

4.17

4.18

4.19

4.20

4.21

Radar cross sections of an infinitely long conducting circular cylinder
with or without a magnetic coating in the case of TM excitation
( Icon = 21:, u,t/a = 0.01 —10.03 ) .......................................................................

Amplitude distribution of the z-component current on the surface

of an infinitely long conducting circular cylinder with or without

a magnetic coating in the case of TM excitation

( koa = 21:, u,t/a = 0.01 -10.03 ) .......................................................................

Radar cross sections of an infinitely long conducting circular cylinder
with or without a magnetic coating in the case of TM excitation

( koa = 51:, Ma = 0.01 -10.03 ) .......................................................................

Amplitude distribution of the z-component current on the surface
of an infinitely long conducting circular cylinder with or without a
magnetic coating in the case of TM excitation

( Icon = 51:, urt/a = 0.01 -— 10.03 ) .......................................................................

Radar cross sections of an infinitely long conducting circular cylinder
partially coated with a magnetically lossy thin layer in the case
of TM excitation ( koa = 21:, u,t/a = 0.01 - j0.03 ) ..........................................

Amplitude distribution of the z-component current on the surface
of an infinitely long conduCting circular cylinder partially coated
with a magneticallt lossy thin layer in the case of TM excitation
( koa = 21:, u,t/a = 0.01 - 10.03 ) .......................................................................

Radar cross sections of an infinitely long conducting circular cylinder
partially coated with a magnetically lossy thin layer in the case of
TE excitation ( koa = 21:, tut/a = 0.01 - j0.03 ) ................................................

Amplitude distribution of the O-component current on the surface

of an infinitely long conducting circular cylinder partially coated

with a magneticallt lossy thin layer in the case of

TE excitation ( koa = 21:, ug/a = 0.01 — j0.03 ) ................................................

Amplitude distribution of the z-component current on the surface
of an infinitely long conducting rectangular cylinder in the case
of TM excitation ( koa = 21:, u,t/a = 0.0 ) ........................................................

Amplitude distribution of the z-component current on the surface

of an infinitely long conducting rectangular cylinder coated with

a magnetically lossy thin layer in the case of TM excitation

( koa = 21:, u,tla = 0.01 - 10.03 ) .......................................................................

Amplitude distribution of the z-component current on the surface

of an infinitely long conducting rectangular cylinder partially

coated with a magnetically lossy thin layer in the case of

TM excitation ( koa = 21:, tut/a = 0.01 - j0.03 ) ..............................................

xii

167

168

169

170

171

172

173

174

175

176

177

4.22

4.23

4.24

4.25

5.1
5.2

5.3
5.4.a
5.4.b
5.4.0
5.5
5.6
5.7
5.8
5.9
5.10

5.11

5.12

Radar cross sections of an infinitely long conducting rectangular

cylinder partially coated with a magnetically lossy thin layer in the

case of TM excitation ( koa = 21:, pig/a = 0.01 — j0.03 ) ..................................
Amplitude distribution of the z-component current on the surface

of an infinitely long conducting rectangular cylinder partially

coated with a magnetically lossy thin layer in the case of TM

excitation ( koa = 21:, u,t/a = 0.01 —10.03 ) ......................................................

Radar cross sections of an infinitely long conducting rectangular
cylinder partially coated with a magnetically lossy thin layer in the
case of TE excitation ( koa = 21:, LL,t/a = 0.01 - j0.03 ) ...................................

Amplitude distribution of the 0-component current on the surface
of an infinitely long conducting rectangular cylinder partially coated
with a magnetically lossy thin layer in the case of TE excitation

( [too = 21:, u,t/a = 0.01 —10.03 ) .......................................................................

Cavity-backed antenna with an impedance sheet on the aperture .................

Radiation boundary condition and boundary separating total field
region and scattered field region .....................................................................

Yee’s model near the impedance sheet ...........................................................
Integral path for E, ..........................................................................................
Integral path for H, ..........................................................................................
Intergal paths ....................................................................................................
Points used in the Mur’s second order approximation ...................................
Mur’s difference scheme near the gound plane .............................................
Sources and their images .................................................................................
Decomposition of total fields in the scattered field region ............................
Apeture fields and their images .......................................................................

(a) x-component of scattered electric field on the empty aperture

of an opened cavity in the ground screen; (b) x-component of

scattered electric field on the aperture covered with an impedance

sheet. ...............................................................................................................

(a) y—component of scattered electric field on the empty aperture

of an opened cavity in the ground screen; (b) y-component of

scattered electric field on the aperture covered with an impedance

sheet. ..............................................................................................................
Radar cross section of an opened cavity ( koa = 21:, d/l = 1 )

when a plane wave is normally incident on it ( 9 = 0.0 ): (a) an

xiii

178

179

180

181

183

185
187
189
189
193
197
199
200
203
204

5.13

5.14

5.15

5.16

5.17

5.18

5.19

5.20

5.21

5.22

5.23

opened cavity with its aperture empty; (b) an opened cavity with its
aperture covered by a thin film ( or = 0.01 ). ............................................... 212

Comparison of radar cross sections in the X-Z plane of an opened
cavity ( koa = 21:, (II): = 1 ) when a plane wave is normally
incident on it ( 0 = 0.0 ). ................................................................................ 213

z-component of electric field on the empty aperture of an opened
cavity in the ground screen: (a) Ez at half cell above the
aperture; (b) E, at half cell below the aperture. ............................................ 214

Total electric field distribution at the plane 2 = - d + Sdz which

is 5 cells away from the bottom of the cavity ( koa = 21:, d/l = 1 )

when a plane wave is normally incident on it ( 6 = 0.0 ) and

no film is covered on the aperture. ............................................................... 215

Total magnetic field distribution at the plane 2 = — d + Sdz which

is 5 cells away from the bottom of the cavity ( koa = 21:, d/l = 1 )

when a plane wave is normally incident on it ( 6 = 0.0 ) and

no film is covered on the aperture. ............................................................... 216

(a) x-component of total electric field in the X-Z plane inside the
cavity with its aperture empty; (b) x-component of total electric
field in the Y-Z plane inside the cavity with its aperture empty. ................ 217

(a) y-component of total magnetic field in the X-Z plane inside the
cavity with its aperture empty; (b) y-component of total magnetic
field in the Y-Z plane inside the cavity with its aperture empty. ................ 218

(a) z-component of total magnetic field in the X-Z plane inside the
cavity with empty aperture; (b) z-component of total magnetic field
in the Y—Z plane inside the cavity with its aperture empty. ........................ 219

(a) x-component of total electric field in the X-Z plane inside the

cavity with its apertutre covered by an impedance film; (b) 1:-

component of total electric field in the Y-Z plane inside the cavity

with its aperture covered by an impedance film. .......................................... 220

(a) y-component of total magnetic field in the X-Z plane inside the

cavity with its aperture covered by an impedance film; (b) y-

component of total magnetic field in the Y-Z plane inside the cavity

with its aperture covered by an impedance film. .......................................... 221

(a) x—component of scattered electric field on the empty aperture
of a cavity-backed antenna; (b) x—component of electric field on
the impedance-film covered aperture of a cavity-backed antenna. .............. 222

(a) y-component of scattered electric field on the empty aperture
of a cavity-backed antenna; (b) y-component of electric field on
the impedance-film covered aperture of a cavity-backed antenna. .............. 223

xiv

5.24 Comparison of radar cross sections of a cavity-backed antenna
(koa = 21:, d/l = 1, location of the antenna MA = 17/21, Z", = 50 0)
when a plane wave is normally incident on it ( 6 = 0.0 ):
(a) Radar cross section in the X-Z plane; (b) Radar cross section
in the Y-Z plane. ............................................................................................. 224

5.25 (a) x-component of total electric field in the X-Z plane inside the
cavity produced by a cavity-backed antenna with its aperture empty;
(b) x-component of total electric field in the X—Z plane inside the
cavity produced by a cavity-backed antenna with its aperture covered
by an impedance film. .................................................................................... 225

5.26 (a) y-component of total magnetic field in the X-Z plane inside the
cavity produced by a cavity-backed antenna with its aperture empty;
(b) y-component of total magnetic field in the X-Z plane inside
the cavity produced by a cavity-backed antenna with its aperture
covered by an impedance film. ...................................................................... 226

5.27 Comparison of the total current density distribution on the
cavity-backed antenna : ( centrally fed impedance Z", = 50 Q;
the antenna is (/1. = 17/21 from the bottom of the cavity;
receiving power in the case of no film covered P = — 54.4085 db;
receiving power in the case of an impedance covered aperture
= — 65.6098 db ) ............................................................................................ 227

XV

CHAPTER I

INTRODUCTION

In recent years, due to the increasing complexity of electromagnetic problems and
advances in computer capability, much more attention of the research community is
focused on seeking numerical solutions of various engineering problems in the field of
electromagnetics. A number of numerical techniques have been introduced and
applied to a variety of specific problems. However, there remains the demand for

finding more efficient numerical methods to meet growing needs.

There are several numerical approaches available to study electromagnetic scatter-
ing by specific three-dimensional objects [1-7]. These approaches can be categorized
into frequency-domain techniques, time-domain techniques and hybrid techniques.
Frequency-domain techniques, including the method of moments (MoM), uni-moment
method, finite-element method and iterative methods, treat electromagnetic scattering
as a boundary value problem by deriving and then solving an integral equation for the
electric field and/or magnetic field with specified boundary conditions. Time-domain
techniques, including the finite-difference time—domain method (FD-TD), potential
integral approach, and the Singularity Expansion Method, strive to examine the tran-
sient behavior of fields before the onset of a steady state and to model problems where
non-sinusoidal excitations lead to scattering responses that may have no sinusoidal
steady state. The hybrid techniques, including the method of moments/high frequency
method, the FD-TD/method of moments and the FD-TD/high frequency method, are

rooted in the combination of techniques in the other two categories.

In the past decade, most problems involving electromagnetic scattering have been
handled by frequency-domain methods, especially by the method of moments. The

basic steps of the moment method in solving integral equation formulations involve

converting the integral equation into a matrix equation and then solving the matrix

equation by standard routines.

The method of moments is based on establishing an integral equation for a given
geometry and composition. To solve a heterogeneous or inhomogeneous problem, a
volume integral equation is normally used and discretized by ten'ahedral modeling with
a vector basis [8]. The surface integral equation approach is suited to analyzing homo-
geneous objects[12] or heterogeneous objects. The coupled surface integral equations
are set up in terms of equivalent electric and magnetic currents on the interfaces of a
heterogeneous body. In Chapter 2, this method is extended to analyze an arbitrarily
shaped heterogenuous body. Based on the equivalence principle, a set of coupled sur-
face integral equations is derived in section 2.2 and then solved numerically by the
method of moments in section 2.3. Triangular elements are used to model the inter-
faces of the heterogeneous body, and vector basis functions are introduced to represent

the unknown surface currents in the integral equations.

Because of its characteristics in matrix operation, the moment method becomes
very inefficient and impractical when objects become electrically large. The finite
difference time domain method provides a good alternative to the traditional moment
method. The FD-TD method is a direct solution of Maxwell’s time-dependent curl
equations. The main algorithm of the FD-TD method includes Yee’s model, which
applies the simple second-order central-difference approximations of both spatial and
temporal derivatives of the electric and magnetic fields directly to the differential
operators of Maxwell’s curl equations. It also applies the radiation boundary condition
on the outer truncated boundary surfaces to simulate the outside extension when
infinite space must be truncated. The system of equations developed by Yee [?] to
update the field components is fully explicit such that the required computer storage

and running time is proportional to the electrical size of the volume modeled. This

sets a remarkable difference from traditional moment-method which requires the

, inversion of a matrix.

The FD-TD method has been used to solve problems which include two and three
dimensional electromagnetic wave scattering, electromagnetic wave penetration and
coupling for both two and three dimensions, inverse scattering reconstructions in one
and two dimensional cases, and microsuip and microwave circuit models. However,
little effort has been made in the application of the FD-TD method to metallic objects
coated with thin material and to transmitting and receiving antennas. Chapter 3 intro-
duces the fundamental algorithm of the FD-TD method. It examines the derivation of
Yee’s model, the radiation boundary condition, and stability analysis, to provide a
basis for extending the FD-TD method to important applications in Chapter 4 and

Chapter 5.

It is known that a thin layer of electrically lossy material on a perfectly conduct-
ing body cannot efficiently reduce its radar cross section. The tangential component of
electric field is very small near the surface of a perfect conductor and consequently the
induced current and dissipated power in the coating layer are very small. On the other
hand, if a thin magnetically lossy layer is used to coat the body, its radar cross section
can be significantly reduced. The tangential component of magnetic field is very large
on the surface of a conducting body, resulting in a large equivalent magnetic current

and a high dissipated power in the coating layer.

In Chapter 4, a new set of coupled integral equations is derived for treating the
scattering problem of a perfectly conducting body coated with a thin magnetically
lossy layer. These electric field integral equations and magnetic field integral equa-
tions are numerically solved by the method of moments. To validate the derived
integral equations, an alternative method to solve the scattering problem of an infinite

circular cylinder coated with a thin magnetic lossy layer has also been developed,

based on the eigenmode expansion. Results. for the radar cross section and current
densities via the MoM and the eigenmode expansion method are compared. A special
application of the FD-TD method to the thinly coated cylinder is also presented. A
new algorithm for the FD-TD is developed based on the integral representation of

Maxwell’s equations.

In some applications, it is desirable to hide an airplane from the detection of radar
systems. Since an antenna on the airplane is an efficient scatterer, it is necessary to
cover the antenna with a lossy layer to reduce its radar cross section. In Chapter 5,
the effects of an impedance sheet covering a cavity backed antenna, on the scattering
and receiving characteristics, is studied. The integral equation technique is difficult to
apply to this problem due to its complex geometry and the infinite conducing surface
involved. Therefore, the finite difference time domain method is applied. Based on
the fundamental theory described in Chapter III, Yee’s model and the radiation boun-
dary condition are modified for treating an infinitely thin impedance sheet and the

infinite structure.

CHAPTER II

INTERACTION OF ELECTROMAGNETIC FIELDS WITH
THREE-DIMENSIONAL, HETEROGENEOUS, DIELECTRIC,
MAGNETIC AND LOSSY BODIES

2.1 Introduction

The interaction of electromagnetic waves with three-dimensional heterogeneous(or
piecewise inhomogeneous), dielectric, magnetic, and lossy bodies has been extensively
studied recently because of its relevance in the modification of the radar cross section
of a metallic body by magnetic coating. Other motivations for the study are due to the
need for quantifying the EM power deposition in human bodies and the application of
EM fields in medical diagnostics. In the literature, a number of methods have been
developed and applied to three-dimensional scattering problems [1-7]. Some of these
methods have been demonstrated only for homogeneous bodies, while others are res-
tricted to bodies of revolution. Even though the tetrahedral modeling, based on the
volume integral equation [13], is applicable to arbitrarily shaped heterogeneous

bodies, a more efficient and accurate numerical method is desirable.

More recently, the hybrid finite element method (HFEM ) has received increasing
attention and has been applied to two-dimensional electromagnetic scattering problems
[16,18]. In this technique, the finite element method has been used for the near field
region and the boundary-element method has been used for the exterior far field
region. Objects to be modeled can simultaneously contain conductors, lossy dielectrics
and lossy magnetic materials with complicated high aspect ratio geometries. This
method can theoretically be applied to three—dimensional heterogeneous bodies, but so
far only one such effort has been reported [16]. In that paper, the boundary conditions

across the interfaces of a heterogenuous body must be accounted for during the assem-
bly of the matrix equations. Thus, the advantage of the HFEM for seeking a sys-
tematic solution to the complex geometries is lost, and this method becomes more

difficult to apply.

With the rapid advance of supercomputers, the finite difference method in the
time-domain (FD-TD) has become quite popular. The FD—TD method is a direct solu-
tion of Maxwell’s time-dependent curl equations with boundary conditions and initial
values. For an open-boundary problem, the infinite space is truncated into a finite
space by introducing an artificial boundary condition or the radiation boundary condi-
tion. Space and time discretizations are selected to bound errors in the sampling pro-
cess, and to insure numerical stability of the algorithm. The finite difference equations
at each time step are fully explicit, so that the computer storage and running time are
proportional to the electrical size of the volume to be modeled. This makes it promis-
ing when a very large body is involved, such as a man model with 50,000 volume
cells [19,20]. However, for this method it is not easy to model complicated curved
surfaces, and numerical artifacts such as instability and nonphysical wave reflections

have to be discussed [17].

The surface integral equation approach is very well suited to analyzing homo-
geneous objects [12] or heterogeneous objects. The usual procedure in this method is
to set up coupled surface integral equations in terms of equivalent electric and mag-
netic currents on the interfaces of a heterogeneous body. Instead of considering the
whole domain, only boundaries have to be involved. We extend this method to
analyze an arbitrarilly-shaped heterogeneous body by choosing an efficient and simple

numerical scheme.

In this chapter, a method of solving this problem efficiently has been developed.

Based on the equivalence principle, a set of coupled surface integral equations is

formulated in section 2.2, and then solved numerically by the method of moments in
section 2.3. Triangular elements are used to model the interfaces of the heterogeneous
body, and vector basis functions are defined within each triangular element to insure
that the normal components of equivalent currents are continuous across triangular
edges. The validity of this method will be established by applying it to a concentric
sphere, with or without a perfectly conducting sphere inside which an exact solution
exists. This method is applicable to a partially coated sphere and a homogeneous
human model. This method will be compared with other existing methods in section

2.6.

2.2 Derivation of Coupled Surface Integral Equations for a Heterogeneous Body

The geometry of a heterogeneous body is shown in Figure 2.1. In this section,
we will only consider two different cases: 1) A heterogeneous body without any per-
fect conductor inside and 2) A heterogeneous body which includes a perfect conductor.
By methods similar to the formulations of these two problems, a set of coupled surface

integral equations for a more complicated geometry can be easily derived.

2.2.1 The Preliminary Theorems

Before the coupled, surface integral equations are derived, we first review the

general solutions of Maxwell’s equations in terms of sources and surface fields.

Consider the geometry shown in Figure 2.2. An infinite region V is bounded by
an infinite spherical surface S... and within the region, there are volume EM sources,
J,p,J,,,,p,,, and a closed surface S, which may enclose some other EM sources. rt is the
unit normal vector pointed outward from V. We can find E and H fields at an obser-
vation point r(x,y,z) maintained by all EM sources ( including the sources inside S, )

and express them in terms of J.p.J,,..p,,, and the surface fields ( E and H ) on S, as

Ei

 

F'
1g
. 2
.1
G
e
om
CU'Y
of
a
h
ete
ro
gen
eo
us b
od
y
w
ith
a
pl
an
e
w
av
e
ex
cit
ati
o
n

     
   

“if: .L ........

:ifiz-z':-'.1;?:- t-Jr-r-fig"
W"Q"‘¢¢".?W~""

Sso

Figure 2.2 Geometry for illustration of the equivalence principles

10
follows:

E(r) = 41.] {—jqud) — meV’cb + P-V'cb] dV'
1! 8
7% j [—ju)u(rb<H)<l> + (r’tXE)xV’<I> + (Ii-E) V’<I>] dS’ (2.2.1a)
S,+S_
_ _T_ _- , r ,
H(r)— my )meJmcr>+JxV<r>+ updefl dV'

__7_'_ j [jwe(fiXE)<I>+(ft><H)xV’<I>+(fi-H)V’<I>] as" (2.2.1b)

47‘ s,+s_

The factor T comes from a singular integral surrounding the observation point :

l r inside V
T = 2 r on S,
0 otherwize

/ . . . . .
with (D = e R and R = Ir — r’l, where r’ rs the posruon vector locating the source pornt

 

and r represents an arbitrarily—located observation point. I: = 0)er is the wavenumber.

The complex permittivity e is defined as e = 80(8, — j&) where c, is the dielectric

constant, c is the conductivity of the body and u is its permeability. The integral

over S... vanishes by the radiation condition, and

fi-E = J—V’-(m<H) (2.2.2a)
(08

1111 = -,—1—V’-(nb<E) (2.2.2b)
1th

The above relations between the normal components and tangential components

of E and H fields are valid on a smooth surface and they are derived in appendix A.

2.2.2 A Heterogeneous Body Without Any Perfect Conductors Inside

We first consider the formulation of coupled surface integral equations for a

heterogeneous body which does not enclose a perfect conductor. A heterogeneous

11

body as shown in Figure 2.3 consists of two different regions, V, and V;. The surface
S; encloses the inner region V; and 11; is the unit vector pointed outward from V;. The
surface S, encloses the whole body and fl, is the unit outward normal vector on 5,. An
electric field E‘ and a magnetic field H‘, defined to be the fields due to the EM sources
in the absence of the scatterer, are incident upon the body. The time dependence factor
of equmt) is assumed and suppressed throughout the analysis. Based on the
equivalence principle as discussd in section 2.1, the fields near the interfaces in the

different regions can be determined in terms of equivalent surface currents.
In Region 1 V0 (Eo,llo) :

The medium in the region external to S, is assumed to be free space with permea-
bility lie and permittivity 80. According to the theory discussed in section 2.2.1, the
total electric or total magnetic field just outside the surface S, can be expressed in
terms of the incident E‘ or the incident Hi field and tangential components of E and H

fields on the surface 5,. Using Eq.(2.2.1a), we can write
E(r) = {Id [—j(ou<I> — meV’CD + five] dV'
7%] [-jtou(r‘b<H)<I> + (fixE)xV’<I> + (an) V'cb] dS’ (2.2.3a)
1
Where the volume integral in the expression is equal to E‘. From Eq.(2.2.1a),
H(r) = .4—T1t-l [-jcoeJ,,,<I> + JxV’cp + TIL-pmV’CD] dV’
.72! [itoe(ri><E)<I> + (riXH)xV’<I> + (fi-H)V’<I>] dS’ (2.2.3b)
1

where the volume integral in the expression is equal to H‘. When the observation
point r approaches the surface S, from region 1, T equals to 2. The fields on S, in
region 1 are denoted by (E, H1). 11 is directed as defined in section 2.2.1. Due to
different definitions of fi,r’i,,r’i;, there are sign changes of fi when applying the above

formula to our problems. Finally the expressions will have following forms:

12

( Region 1) V0

 

Fig. 2.3 A heterogeneous body without any perfect conductors

13

E1(r)= 23% 711-:- i [—j(ouo(fi,xH1)<Do + (a,xE‘)xv'<t>o + BLEEV’-(fi,xH1) Who] (15’ (2.2.4a)
I
1110) = 211‘ +-Zl;sj[jmeo(rt,xE1)<po + (ri,le)xV’<l>o + fiv-(mel) Who] (18’ (2.2.4b)
.l
l

e-jkoR
where (1),, =

 

is the free space Green’s function, and k0 = (ox/tron, is the free space
wavenumber.
In Region 2 V, (e,,o,,u,) :

In the region between surface S, and S;, with permeability u, and complex per-
mittivity 8,, the E and H fields just inside the surface S, or just outside the surface S;

can be expressed by tangential components of E and H fields on the surfaces S, and S;.

E(r) = 41; [-ju)|J.J<I> - meV’CD + P-V'cb] dV'
1: 8
—% j [—jwu(rp<H)<I> + (IVE)XV’<I> + (II'E) V’cb] dS' (2.2.5a)
S, + S;

The volume integral in the above expression is equal to zero because of no source

within region 2. Similarly,

911m = 711? j [—ju)eJ,,,<I> + JxV’(l> + ipmV’CD] dV'

—-1— j [ime(ri><E)<D+(rixH)xV’<D+(ri-H)V’<I>] dS’ (2.2.5b)

47‘ s, + 5;
where the volume integral is zero also. The factor T is equal to 2 when r is on S, in
region 2 or on S; in region 2. Let ( E2, H2 ) denote the electric and magnetic fields on
S, in the region 2 side, and ( E3, H3 ) denote the fields on S; in the region 2 side. Not-

ing the sign changes in 11, and we can express fields in the forms:

When r is on S, ,
1320-) = 371:! [—jwu,(fi,xH2)<I>, + (ri,xE2)xV’<l>, + iv-(mxnz) V’tb,] dS’
, ' 1

+ 31;! t-jwulmzxnecbl + (fisz31xV'<b, + "w‘Le-V'ifizxm) W11 dS’ (2.2-6a)
2 1

14

H2(r) = :1. Uwe,(ri,xE2)<I>, + (ri,xH2)xV’<b, + -—.-1—V’-(r’i,xE2) V'¢,] dS’
27‘ 1 10411
+ —l- UC081(fi2XE3)‘D1 + (52XH3)XV'¢1 + fiL-V’(ri;XE3) V'q’ll 615’ (2-2~6b)
21: 2 103111
When r is on S;, the expressions for E,H essentially remain the same in form.
3 __ -1 . A 2 A 2 I . I A 2 I I
E (r) _ 3;![-j(0ll,(n,XH )6, + (n,xE )xV (I), + —J—w8 V '(n,xH ) V (1),] as
l
l
+ 31; J [-ju)u,(ri;xH3)<I>, + (fi;><E3)xV’<I>, + Bé—V’~(fi;xH3) V’CD,] as (2.2.7a)
l
2
Wm = :3- Uwe,(fi,xE2)<I>, + (fi,xH2)xV’<I>, + -.—1—V’-(fi,xE2) V’<D,] dS’
21: I Jmlit

+ i J U081(fi2XE3)‘D1 + (42XH3)XV"D1 + '.—1"V"(’32XE3) V’q’il d5, (2.2-7b)
21: 2 100%

 

where (13,: e is the Green’s function for region 2, and k, = 0) u,e, is the
complex wavenumber in this region.

In Region 3 (e;,o;,u;) :

Similarly in region 3 where the permeability and permittivity are u; and s;
respectively, the region is enclosed by S; with an outward unit normal ii; . The E and
H fields just inside the surface S; can be determined in terms of the tangential com-

ponents of E and H on the surface S; as follows.
E(r) = {Fl [—j(ouJ<I> - meV’CD + Jim” W
—4—:J[-jmu(m<H)d> + (fixE)xV’<I> + (ii-E) V’m] dS’ (2.2.8a)
2
H(r) = 1‘] [—ja)e.l,,,<l> + JxV’<I> + -1-pmV'<I>] dV’
41: u
77;! Unre(rp<E)<I> + (moxv'cb + (fi-H)V’<l>] (13' (2.2.8b)
2

The factor T equals to 2 when the fields on S; in region 3 are evaluated. The
volume integral terms are equal to zero because of source free condition. Let ( E‘,H‘ )

denote the electric and magnetic fields on S; in region 3. Then we can express these

15

fields in the forms:
1240-) = —311-c-J[—jwu;(fi;xH4)<b; + (fi2xE4)xV’<I>; + —w-ie—V'-(a2xn4) V’CD;] as (2.2.9a)
2
2

V’-(r’t‘;xE4) V’<D;] dS’ (2.2.9b)

 

0

11%) = —-2};![ime;(fi;xE4)<D; + (fi;xH4)xV’<I>; +
2
2
e—jkzk
WhCI'C (D2 =

 

is the Green’s function for region 3, and k; = (mime; is complex

wavenumber in this region.

Up to this point, we have the total induced fields (E,H) at the interfaces of the
three different regions. Before going any further, we should introduce the boundary
conditions which must be satisfied at the surfaces S, and S2. The boundary conditions
require the tangential components of the electric and magnetic field be continuous

across S, and 8;.
Boundary Condition :

When the observation point r is on S ,, we may write

r~:'(r)l,,n = E2(r)lm (2.2.10a)

H1(r)lm = H2(r)lm (2.2.10b)
For r on S; :

E3(r)ltan = E4(r)ltan (2.2.11a)

H3(r)ltan = H4(r)ltan (2.2.11b)

Equivalently, we can write

r‘i,xE'(r) = fi,xE2(r) (2.2.1221)
r‘i,xH1(r)= fi,xH2(r) (2.2.1211)
fi2xE3(r) = fi;xE‘(r) (2.2.12c)
ii;xH3(r) = ii;xH4(r) (2.2.12d)

By forcing the electric fields to satisfy the boundary conditions on the two inter-

faces, we have a set of so-called electric field integral equation(EFIE) as follows:

16

For r on S,:

E‘(r)l,,,, = 152ml,an
47:90),” = JUGXfi,XHl)(u0<Do + fli‘pi) ’ (file‘)><V’(<I>o + (D1) -

- ch (1)
.LV’-(ri,le) V'(—9- + -—‘)1 (15’ + jt—jwu,(a;xn3)<t>, +
to 80 81 ;

(fisz3)><V’<I>1 + -m£-V’-(fi;xH3) mp], dS’ltan (2.2.13a)
I

For r on S; :

E3(r)ltan = E4(r)ltan

. A 3 A 3 I _L I A 3 I (D1 (D2 I
JUMHZXH )(u1¢1 '1' u2¢z) - (H2XE )XV (Cpl '1' (D2) — 0) V '(NzXH ) V (? + '8—)] d5 [tan
2 l 2

= Juwufifilxfl‘w, —(fi1><E‘)><V’<1>, — —w£—V’-(ri,xH1) V’<I>,] dS'llan (2.2.13b)
1
1
Similarly by forcing the magnetic fields to satisfy the boundary condition, a set of
so-called magnetic field integral equation(MFIE) is formulated.
For r on S,:
1110011,”, = 112mlum
4KHi(r)lm=J[—jw(filXEl)(€0¢0+€1¢1)—(fi1XH1)XV’((D0 '1" (DO—Jig) V"(fi1XEl)
l
(D (D
V’(Elf-+71)]dS’+J[/'u)e,(fi;xE3)<D,+(fi;xH3)xV’<D,+-fil—V’-(ri;xE3)V’<I>,]dS’ltan (2.2.13c)
r 2 1
For r on S; :

H3(r)|tan = H4(r)ltan

(I) (I)
JUMfizXEBXEIqH + 82¢z) + (fizXH3)XV’(¢1 “'1' (D2) + 7:;V"(fi2XE3) VT]; + fHdS’ltan
2 l 2

= IUcoe,(fi,xE1)<D,i (fi,xH1)xV’<I>, + j—oil—V'tmn‘) V’¢,] dS’ltan (2.2.13d)
l
1
To simplify the notations, let’s define the equivalent surface current as follows:

-h‘,xE‘(r)

= M1(r) (2.2.14a)
1lo

17

_. 3
____”2"E 0’) = M2“) (2.2.14b)
1lo
ri,xH1(r)= K‘(r) (2.2.140)
r’i;xH3(r) = K2(r) (2.2.14d)

After matching boundary conditions and scaling the physical dimensions by k0,

we have a set of coupled surface integral equations:

r on S,:
41: .- __ 1 1 , . , 1 , (Pr ,
TE(r)ltan ._ JUK (<r>O + u,,<b,) + M xV (<1)0 + ¢,) —,V -K V ((1),, + E.” as
0 , 71
+ JP—jfllezqh - MZXV’q), + ELVAKZVQDI] dS'ltan (2.2.153)
2 I]

r on S2:

 

2 I . I 2 I ¢1 ¢2 ’
J [K (u,,<r>, + u,;<r>;) + szV ((D, + o?) — ,V -K V (:- + 8 )1 as Itan
2 rl r2

Jungle, + M‘xV’cb, — ELV’K‘ V'¢,] dS’ltan (2.2.15b)
, 71

r on S,:

. <r>
41:H‘(r)lm = JDM‘GDO + 2.14%) - K‘xV’(<I>o + <D1) + -}V’~M‘ V’(<I>o + -1-)] dS’

1 rl

 

 

+ j[-je,,M2<r>, + szV’cb, - Tll—V'M2 V’<D,] dS’ltan (2.2.15c)
2 J r1
I' 01'! S2:
2 I 1 I 2 I (D1 (Dz I

JUM (e,,<l>, + e,;<l>;) — szV (CD, + o9 + —_V -M V (— + )1 as Itan

2 .l “'71 ur2
=errtM1¢1 - K‘XV’¢, + 1,: WM1 V’<D,] dS’ltan (2.2.15a)

1 I]

ll: 8." . o . . .
where u,,- = -— and 6,,- = — = e,,-’ - j—. Due to the normalization wrth res ect
#0 80 (060 p
am

to k0, the Green’s functions have the forms (1),: e

 

and [6,, =

5
k0 .

18

In the above four coupled surface integral equations, (K1,M‘,K2,M2) are four unk-
nowns which can be determined by solving by these four coupled surface integral

equations.

2.2.3 A Heterogeneous Body Enclosing a Perfect Conductor

As shown in Figure 2.4, a perfect conductor covered by a dielectric, magnetic,

and lossy medium is located in free space. The medium has the complex permittivity
of e = Eda-1.398;) , and permeability u. A similar derivation to the previous case,
with the exception of vanishing tangential components of elecnic field on the surface
of the perfect conductor, will be presented.

In Region 1 (50410) :

The expressions for E‘,H1 remain the same as Eq.(2.2.4).

 

E1(r)=—-2E‘%J [—jwu0(zi,xH‘)<bo+(ri,xE1)xV’<Do+TD-l; V’-(fi,xH‘)V’<DO]dS’ (2.2.16a)
l
H1(r)=2Hi+?1£J[imeo(r‘i,xE1)<I>o+(ri,xH‘)xV’<DO+ , u V'-(a,xE‘)V'<r>O]as (2.2.l6b)
, 0

1“ Region 2 (81,01,111) 2

Due to the zero tangential component of E field on the perfect conductor surface
S;, the integrands involving ri;xE3 over S; disappear. For the fields (E2,H2) on the sur-

face S, in the region 2 side, we have the following simplified expressions:
E2(r) = €11?! [—jtou,(ii,xH2)<b, + (fi,xE2)xV’<I>, + —wj£—V"(fi,xH2) V'cp,] dS’
, l
+ -1—$[[—ja)u,(ii;xH3)<I>, + —LV'-(a;xH3) V'¢,] dS’ (2.2.17a)
21‘ 2 (1)81
H2(r) = :1 Ume,(r‘i,><E2)<b, + (ii,xH2)xV’<I>, + -,—1—V'-(a,xE2) V’<D,] dS’
21: ‘ - lelr

+ 31; j, + (fi3xH3)xV’<I>,] aS’ (2.2.17b)
2

19

   
  

( Region 1)

 
  
   
 

( Region 2 )

53/“ ----------

(81 411,01)

Fig. 2.4 A heterogeneous body with a perfect conductor inside

On the surface of the perfect conductor 5;, only the equation for E3 field in region

2 side is required.
E3“) = ijl-jwuimixflzm +(fi1XE2)XV’(D, + T1L)J‘;3—V’(r1,xnz) V'cp,] as’
l
I
+ ifl-jw‘mzxmmt + J-V’tfiszh V’<t>,] as’ (2.2.18)
21'! 2 (08,

Next we will apply the boundary condition to set up the coupled integral equa-
tions. On the surface S,, tangential components of E and H should be continuous,
while on the surface of the perfect conductor 5;, the tangential component of E3 is

equal to zero.

r on S,:
1310):“, = E2(r)lm (2.2.l9a)
H‘(r)lum = H2(r)l,an (22.1%)
r on S; :
E3(r)|tan = 0 (2.2.l9c)

Alternatively, we may write:

ri,xE1(r) = r’i,xE2(r) (2.2.2021)
r’i,xH‘(r) = ri,xH2(r) (2.2.10b)
Ii;xE3(r) = o (2.2.20c)

with the equivalent surface currents defined as:

—fi,xE1(r)

= M1(r) (2.2.21a)

llo
r’i,xH‘(r) = K‘(r) (2.2.21b)
I‘r‘;xH3(r) = K2(r) (2.2.21c)

After matching boundary conditions and scaling by k0 ,we have another set of

coupled surface integral equations as follows:

For r on S,:

21

_ (I)
glacial, = juxlteo + un<1>1>+ Mle'tcbo + «m -1V'-K‘ V'<<t>o + e—‘n d5’
0

1 rl

+ J {-ju,1K2<r>, + ELV’KZV’dh] dS’ltan (2.2.22a)
2 71
For r on S;:

<r>
J UKZu,,<I>, — jV'-K2V'—‘] dS’ltan
2

 

£71
= JujinK‘d), + Mle'cp, — ELV’Je V’d),] dS’Itan (2.2.2213)
1 r1
for r on S,:
. (I)
mutt-Man = j UM1(¢0 + 5,14,)- K‘xV'tcbo + chi) + iVM‘ V’(<I>o + ‘ )1 dS’
1 rl
l
+ J[K7'XV’<D,] dS’ltan (2.2.22c)
2
h “i and e, e,’ ’ ‘ o,- D e to the normali ation ith respect to k0
wereur1=-— .-=—=e,.--J—-. u z w .
lie 60 (050

-jk.eR k-
the green’s functions have the forms (1),: £7?— and k,,- = 7;; , the same forms as that
o

of section 2.2.1.
On the surface 5;, only K2 is unknown because M2 equals zero. In this case, we

have only three unknowns, (K‘,M1,K2) . These three equations are sufficient to solve

for the three unknowns.

2.2.4 Metallic Body with Partial Coating

A more general geometry than that of section 2.2.3 is a perfect conductor par-
tially coated with a layer of dielectric, magnetic, and lossy material as shown in Fig.
2.5. K1 is the electric current on the uncoated area S, of the perfect conductor in
region 1, and K3 denotes the electric current on the coated region S; of conducting

body in region 2. K2 and M2 denotes the equivalent electric and magnetic currents on

22

the interface S; between free space and coated region in region 1. 1? denotes the unit
outward vector of S, , S; and 5; respectively as in Fig. 2.5. As S, shrinks to zero, we
have a uniform coated conducting body which is the case of section 2.2.3. The study
of the EM fields interacting with partially coated metallic body has its significance in
the modification of the radar cross section of a radar target, and this topic has received

considerable attentions recently.

We now consider the derivation of a set of coupled surface integral equations for
the geometry of Fig. 2.5. Following the procedure previously, we will write the fields
in the different regions and then use the boundary conditions to set up the coupled sur-
face integral equations.

Define K} = MH, on S, in region 1, K? = fixH; on S; in region 1 and
M1; = —r’r><F; on S; in region 1 as shown in Fig.2.6. We can then write the fields in

region 1 (ascetic) :
i T . 1 . I l I I
(r) 4“ j I was, 0 #0380 1 01615
l
— {dummxioo + M7,'><V’<I>, — —w-£LV'-K% V'ct>,] aS' (2.2.2321)
2 1

1 VIM? v1¢01 dSI
H0

 

H1(r) = r H‘ + fifl-jmeoMffbo + Kij'cbo +
2
+ 771; jmngeo, aS' (2.2.23b>
1

Define K; = 113<H; on S; in region 2, M3 = -rb<E; on S; in region 2 and K3 = mm,

on S; in region 2 as Fig. 2.7. We can formulate the fields in region 2 (e,,o,,u,) as :
E2(r) = «£- [-jcou,t(§<r>, — ngV'cb, + —-LV'-K§ V'<r>,] (15’
47C 2 (081
— l—j[-jtou,K3<D, + —-LV’-K3 V’<I>,] dS’ (2.2.24a)
41C 3 (08,
H2(r) = 1- ymmgo, - K§xV’<l>, — -,—1—V’~M§ V'¢,1 dS’
47‘ , 103%

+ 3%,, + 1(ngde as (22.24!»
3

23»

 

( Region 2

S3

( Region 1 )

Ei

Hi

Fig. 2.5 A perfect conductor partially coated with lossy material

24

Bi
1
, 3:1,
H1
K2/
M21
32 A
“\
1
K
(€0,1-10,00) 1

Fig. 2.6 Sources maintaining fields in region 1

=1)

25

 

Fig. 2.7 Sources maintaining fields in region 2

26

where K}, K? and M? indicate equivalent currents on Sl and S; in region 1, while

K2, M§ and K3 indicate equivalent currents on S; and S3 in region 2.

Since K? = Kg and M2 = M3 , we have only four unknowns, K2, M2, K} and K3
For convenience, the subscripts of these four unknowns are omitted and they are
referred to as K2, M2, K1 and K3 . To determine these unknowns, four boundary con-

ditions are required on the three interfaces, S1, 5;, S; .

For 1' on S;, tangential components of E and H are continuous across S; , i.e. E‘Ium
= Ezlm and H‘Ium = Hzlum or alternatively M1: M2 and K1: K2 . Applying these

boundary conditions, we have:

41tE‘(r)ltan =

- <1> (D
JquZQJOcpo + “1(1)” + M2XV’((DO + (D1) - 6V"K2 V’(—£O£ + 24)] d5,
1
2
+ JUquK‘CDO + J-V’K‘V'cbo] (15’
1 (D60
+ Il-jwu1K3¢1 + -w'%-V’-K3V’(I>1]dS’Itan (2.2.25a)
1
3
' - 2 I _L I 2 I (DO (D1 I
4nH'(r)Im=5[uooM (80¢0+81<D1)—K2><V (<p0+¢,)- m V -M V (3+?)] dS
2 l
+ J[—K‘xV’<DO] (15’
l
+ J[K3XV’<I>1] dS’ltan (22.2513)
3
For r on S; , the tangential component of E is zero. Thus, we have :
0 = J— J’ [—jo3p.1K2<D1 — M’xV'cb1 + -L-V"K2 V'cbl] dS’
41E 2 (1)81
— —1—J[-jtou1K3¢l + —j—V'-K3 V'cbl] dS’ltan (2.2.26)
47‘ 3 (1)81

For r on S, , again‘the tangential component of E is zero. From Eq.(2.2.23a) , it

leads to:

27.

41:E‘(r)lum = JUtouoK2<Do + szV’CDO — Tot-V’Kz V’cbo] dS’

2

+ JUwuoK‘¢o + M‘xV'cpo — Riv-K1 V’<D0] (15’ (2.2.27)
1

We can scale these equations by the free space propagation constant [to and sum-

marize them as follows:

 

For r on S; :
filflflltfln = JUKZM’O + Hrl‘b1) + M2xV’(<I>o + (D1) - IV"K2 V’W’o + (D: )1 (15'
2 r
+ JUKOCDI +jv”K1V’¢o] ‘15,
l
+ J [-ju,1K3(I>1 + EL‘V’K3V'CD1] dS’Itan (2.2.28a)
3 r

 

. (I)
maxrym = Jun/12m, + e,1<1>,) — szv'(<r>o + (bl) + -}-V’-M2 V’(<I>0 + ‘ )1 (15’
2 71
+ J[—K1xV'¢O] dS’
l
+ J[K3xV’<D1] dS’ltan (2.2.28b)
3

For r on 5;:

0 = i J [-jtt,1K2<I>1 — szv'cpl + -LV’-K2 V’CD,] (15’
2

8r:
- :1; lat—Mm», + évam mpg dS’ltan (2.2.29)
For r on $1 :
:‘T’SEiaymn = J [,chpo + szv'cpo — ,V'-K2 V’<1>o] (15’
2
+ JUK‘d’o -1V’-K‘ V’¢ol dS’Itan (2.2.30)
1

The above equations can uniquely determine the four unknowns K1, K2, K3 and
M2 . As S1 or S; shrinks to zero, two special cases are obtained as the cases of a per-

fect conductor with uniform coating or a perfect conductor without coating. Special

28

attention should be paid to the interface S; when a numerical scheme is applied, and

this will be discussed later.

2.3 Numerical Algorithms

In the previous section, we have developed three sets of coupled surface integral
equations. In this section, we will solve these coupled surface integral equations by
the method of moments. Triangular elements are chosen to model arbitrarily—shaped
surfaces. The advantages of triangular patch surface modeling have been elaboratly
discussed in references[9,10]. Vector basis functions for the equivalent electric and
magnetic currents are defined in each of the triangular surface patches, and a Galerkin
method is implemented to solve for the unknown surface equivalent electric and mag-

netic current distributions.

2.3.1 Basis Function

Vector basis functions [8,10] are defined in this section both for electric current
and magnetic current distributions. Detailed derivations of the basis functions are dis-
cussed in reference [10]. Figure 2.8 shows two triangles, T; and T; with the It“ com-
mon edge. The electric and magnetic currents flow along radial direction, (5;, in trian-
gle 7;, and similarly flows along radial direction, 9;, in triangle 7;. Referring to Fig-
ure 2.9, if I; is the base lenth of common edge, then height lenths of triangles T; and
T; are respectively given by 2.4;”; and 2A;/l,,, where A: represents the areas of Tfi.
Any point in triangles 7; can be conveniently defined either with respect to the global
origin, 0, or with respect to vertices 0;. The plus or minus designation of the triangles
is determined by the choice of a positive current reference direction [10] for the n’“
edge, which is assumed to be from T; to T;. We define a vector basis function associ-

ated with the n’” edge as

29

11th edge

 

Figure 2.8 Local coordinates associated with an edge

30

f

 

 

 

1,, +
2A: p" r in T;
< l
f;(r)= " ; r in 1; (2.3.1)
2A;p
; 0 otherwise

The vector basis function stated above is used to represent surface electric
current, K, and surface magnetic current, M, on the triangulated surfaces of a given
heterogeneous scatterer. Further mathematical properties of vector basis functions
have been discussed in reference [10], and it gives some of the elegant properties in

detail. In the following section, the main properties are summarized.

2.3.2 Properties of Vector Basis Function, f,,(r)

(a) Along boundary edges, the currents flow basically parallel to the edges.
Hence, they have no normal components to boundary, and no line charges exist along

the boundary.

(b) The component of current normal to the It” edge is constant and continuous
across the edge as can be seen in Fig. 2.9, which shows that the normal components of
pi along n‘” is just the heights of triangles Tfi with edge n as the base and the height
expressed as (ZAf)/l,,. Using this term to normalize the vector basis, its normal com-
ponent to the It” edge, is unity. This result, together with a), implies that all edges of

T; and T; are free of line charges.

(c) The surface divergence of current basis function , which is, in fact, propor-

tional to the corresponding surface charge density, is given by,

31

 

 

 

 

 

 

 

Fig. 2.9 Geometry for normal component of basis function at common
edge

32

 

V In
E r in T;
i 1,,
(V.-f..(r)) = Z: r in 7;; (2.3.2)
; 0 otherwise

The surface charge density is obviously constant in each triangle, total charge associ—
ated with the triangle pair T; and T; is zero, and the basis functions for the charge evi-

dently have the form of pulse doublets.

(d) The surface integral of basis function over adjacent triangles represents

moment given by

1,.
l l W = ""195? + 95? (2.3.3)
7: + 7: 2
= ln(r:+-r;)
AS shown in Fig 2-10. p? is the vector between the centroid of T; and 0; and pg“
is the vector between the centroid of T; and 0;. rf,+ and 1‘? are the distance vector to

centroids of triangles T; and T; from the global reference point, 0.

2.3.3 Testing Procedure and Matrix Equation

Either (K‘,M1,K2,M2) in section 2.2.2 or (K‘M‘.K2) in section 2.2.3, which need
to be solved, are expanded in terms of vector basis functions defined in the previous
section. If interfaces (S,.S;) are discretized into triangular patches with numbers of
edges (N1,N;) respectively, then the equivalent currents can be represented by

N,- .
K‘(r’) = )3 1:,r,(r') (2.3.4)

n=1

and

Ni
M‘(r’) = Z Mf,f,,(r’) (2.3.5)

n=l

‘v'
..

'V‘

33

r (2)

Sub domain
T rrangle \\ (1 Subdomain

Tn+ :P/ Triangle
a". 41 T -
Centr01d of P - \ Centroid of

 

Il

 

 

Fig. 2.10 Coordinate for calculating centroids and moment of basis
vectors

34

where 1': 1,2 for both K and M in section 2.2.2, and 1': 1,2 for K, i: 1 for M in sec-
tion 2.2.3. If, and M; are unknown coefficients to be determined. Since the normal
component of f; at the n‘” common edge connecting T; and T; is unity, each coefficient
of 1:, and M, can be interpreted as the normal components of the electric and magnetic
current densities flowing through the It” edge. For a given triangular face, there exists

three edges corresponding to three vector basis functions.

Substituting Eq.(2.3.4) and Eq.(2.3.5) into Eq.(2.2.15) and choosing the weighting
function the same as the expansion function, Galerkin method, a matrix equation is
formulated. Testing is enforced based on the following symmetric product to reduce

the operator type integral equations to the corresponding functional type equations,

<f,g> = fl; f'g ds (2.3.6)

For the convenience to write the matrix equation into compact form, we define

the electric and magnetic vector and scalar potentials as following:

A1: £K <1>,-as (2.3.7a)
F,- = 1M <1),- as (2.3.7b)
v, = 3[VP-K V’d), as (2.3.7c)
U,- = £V’-M V’<D,- as (2.3.7d)
and
P,- = in V’ ch,- as (2.3.7e)
Q.- = £Mx V’ o, as (2.3.70

where i may vary from 0 to 2 in the case of section 2.2 or from O to 1 in the case of

section 2.3.

The matrix equation for the coupled surface integral equations given in section

2.2.2 is:

 

L

35

 

2‘1: 2'35 C21]: Cg; 11 E
Zr. 52: C2). C21. ’2 = [01
”it. Di:- Y” ”5., 3‘ [’6’]
021 022 Vii 22 .2N>0N 2

(2.3.8)

where ’1’ indicates the outer surface or S1 and ’2’ for the inner surface or S; , ’E’

stands for EFIE and ’M’ for MFIE. Z,Y,C and D are submauices with size N10,; by

N10,; The system matrix size is 2N by 2N with N = N1 + N;.
tors with lengths of N1 or N;, and [0] is N; - zero vector .

which will be solved for.

expressions :

(ZEIEmn

+

212nm:—

(Z22 )Inn=

 

 

           

 

1 eh
m + pm _
=fl 1m;{+ 2 ”riAimn +7 EuriAimn
'=O (=0
1 1 V
.023 8,
pS.‘ _
. _. A
+ 2 url Inn
1
__(Vi'mn— :nlmn)}
€71
. of: pf; _
Jlm{ 2 ' urlAlmn+-_2'-' urlAlmn
pt: 2 of; 2 -
1 _. _. ~A-
-18'IM{W—2 i=21urlA 2 i=21url mm

1, M, V and H are vec-
I, M are unknown vectors

The various matrix elements are given by the following

(2.3%)

(2.3.9b)

(2.3.9c)

(2.3.9d)

(2.3.9e)

36

0+ C"
(YMIZ M)mn=“—jh%}filn eflFTm +'E_“ eflr

2 2
+ —-1—(UT,,.. _ 1m) } (2.3.90
”d
(jg?5mn=:jk%?%%n EHFTWI+-%§m EHFEmr
+ ELM” - HM} (2399
M pm 2 p: 2
(333)mn=h_jlm -E—. 2:€%F+ 4'-§_. EiefiFEm
1‘] i=1
2
+ )2 u —(U7m,.— Um} (2.3.9b)
i=1 n
l
(Cfiw)mn={ 2Fm+ Erma}
i=0 i=0

= _ (0%“)... (2.3.9i)

(CEMM)

all
r—J‘—\

:3

E

+

J

3
L_~__J

= _ (0mm (2.3.9j)

(67% M)mn==‘{fqmn+'PEmi}

 

— - (03’1“)... (2.3.910

EM =_ 2 P’-’ 2

(C32)mn 2: mm'F 2: mm

i=1 i=1

= _ (02425»... (2.3.91)
_ 41:1 _p__$.’ _ - p_f: ,-.
Em— Tlo [— 2 E‘+ 2 -E ] (2.3.10a)
u- c—

H,,, = —4n 1,, [ ”Tm-n” + ”TWE- 1 (2.3.10b)

For convenience, we write the matrix elements by testing at the centroids of the

triangles, which is a special case of the general Galerkin method discussed.

37

Next, we consider the set of coupled surface integral equations in section 2.2.3
when a perfect conductor is inside the body. Similarly for the Eq.(2.2.22), the matrix

equation can be formed as follows:

 

 

2ft“ 2%" Cfl‘ 1, E
2%? 2%? oil” 12 = [21 (2.3.11)
Dt’F DEE Y’tl“ M1

Only a small modification of matrix elements in previous case will be needed to

obtain the matrix elements for this problem.

 

 

       

. pf: ‘ of; 1 _
(leffzrmzjlmzi{ A+ . ZuriAimn

 

 

 

= 2 i=0
1
+ 2 (V5,... — OV'W)} (2.3.12a)
= oer.-
pc+ pc-
(212 mn = -j1m TD ° “rlAi-mn + ; ° urlAlmn
+ _me _ VIM} (2.3.12b)
Err
. of: of; _
(251E = Jlm{—2' urlA-i-mn + T urlAlmn
+ —1-(Vlm _ W} (2.3.12c)
Err
c+ 1 c— 1
(25525)... = -fl..{ 92’" ~ 2 mat... + p; - 2 was...
i= 0 i: 0
1
+ 2 —1-(V*.-,,.,. - ma} (2.3.12d)
i=0

 

+ ‘2‘. —lj(U?‘..... — Um} (2.3.12e)

'=0 i=0

1 1
<Cil’...={ £Pt....+ EFT-m}

= _ (0’11")m (2.3. 12f)

38

 

= _ (0mm (2.3.12g)
4“ 1m pic:- 1+ pic;- i-
0
H _4 of: .. p; ..
m = 751ml 7‘“ + TH ] (2.3.1313)

Now we move to partial coated perfect conductor as described in Eq.(2.2.28-

2.2.30). The matrix structure can be written in a form as below:

 

 

 

 

2%”: [0] 25328 Cg" I3 [0‘]
[0] Zsubl 155 2%; C1152! ’1 = E; (2 3 14)
2% 25:? 25525 C52“ ’2 E‘; ' '

l 03435 034,5 03423 r35" M2 L H‘;

where the lower subscript indicates the numbering of the surfaces and upper subscript
for distinguishing EFIE and MFIB, for example, 2%: denotes source points and field
points are both on S; by EFIE. Special attention should be paid for the fields on S;
when a numerical scheme is applied. In our case, the boundary conditions for perfect
conducting surface got to be applied to the triangle edges on S; which have connec-

tions with perfect conductor surface as illustrated in Fig.2.5.

The matrix elements have similar expressions as previous two cases, and will not

exhibited here to avoid repetition.

The vector and the scalar potential integrals take the following forms:

A3“:- I j r,(r)<r>,(rf:,r3as' (2.3.15a)
mm
=de

V5,“: 1 j [V'-r,(r)]<t>.(r$r3ds’ (2.3.15b)

am
= mm

39

 

1,;-

Piimn = “2%;ij I If;(r’)xV’<I>,(r,‘f,t,r')dS’ (2.3.15c)

"' Ti. (7:77.)
= Qi
and
e-jkfii
Writ) =
= Ir}; — r’l
e-jkfi-

 

Mar r')— - (r... — r ’)(1 +jk.R*)(R,

The above integrals are in a convenient form for numerical evaluation. However,
each face has three edges associated with it. Identical integrals would be recomputed
nine times if the elements were computed sequentially by edges. To avoid the costly
and inefficient recomputation of integrals, we instead compute the various matrix ele—
ments by considering faces. This cuts down by approximately ninefold computer the
time required to generate matrix elements. Morever, we note that elements of Z and Y
only differ from each other by the appropriate coefficients which can be conveniently
incorporated while filling matrix elements. Similarly for the elements of C and D,
they are similar except by a multiplying constant such that elements of D can be

directly obtained from the elements of C and vice versa.

We note also that the expressions of the submauix elements contain terms
belonging to different regions. These expressions of integrals are identical except for
the characteristic parameters which appear in the propagation constants and in the mul-
tiplication constants. For numerical efficience, we use the same routine to generate the
integrals in the different regions simultaneously which principally, make up various

matrix elements.

In accordance with the above discussion, we now consider the evaluation of the

various integrals. As illustrated in Fig.2.ll, assume that an observation point in face p

4O

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Fig. 2.11 Local coordinates and edges for source triangle '1'q with
observation point in triangle Tp.

41

and a source point residing in face q. Let’s use a " local indexing scheme " for con-
venience. We number the edges of face q by 1, 2, 3 with edge length l,,l;, and I3, and
opposite vertices at r1,r;,r3, respectively. Face q will be denoted simply as triangle 7"
with area A" , and face p as T" with area A". There are three integrals which essen-

tially make up various matrix elements:

 

 

_kRP
J’rj- A3as= j r, j f,(r’) ‘31 as (2.3.16a)
17 7? 1" RP

1' li _kRP
I .qu= I (—’—) J(—) 3’ dS' (2.3.16b)
1P 1? AP 1" Aq RP
I"... "',_ I pi:[Hf.<r“)><V’<D.(r”.r’)dS’] (2.3.16c)
m 7’ 7"
and
RP = lr” — r’l

These three integrals are most conveniently evaluated by transforming to a local
system of area coordinates [14]. As shown in Fig. 2.12, the vector p, divide T" into
three regions of areas A1, A;, andA; which are constrained to satisfy A1+A1+A3 =A" .

The normalized area coordinates are defined as

_ Ar _ Az _ A3
a- Evn-th—Fv (2.3.17)

and

§+n+§= 1, (2.3.18)

(§, 11, C) vary between zero and unit in T", and the triangle comers r1, r; and r3,
corresponding to the area coordinate (£319 takes the values ( 1, 0, 0 ), ( 0, 1, 0 ), and
( 0, 0, l ) , respectively. The transformation from Cartesian to area coordinates may

be written in vector form as

r’ = fir] + nr; + §r3 (2.3.19)

and surface integrals over TV can be transformed as follows:

42

 

Fig. 2.12 Definitions of areas used in defining area coordinates

43

ll—n
[grads = 2A?! ]; gain + Tlrz + <1 - t - n)r3)d§dn. (2.3.20)
7"
Numerical evaluation of the integrals in Eq.(2.3.16) may be accomplished by
using numerical quadrature techniques specially developed for triangular domains [15]
together with the procedures discussed in reference [10]. Appendix B carries out these

three integrals (2.3.16) numerically in details.

2.4 Numerical Results

In this section, numerical results are presented for the equivalent electric and
magnetic current distributions on some selected scatterers under plane wave illumina-
tion. The geometries considered are heterogeneous concentric spheres and a perfect
conducting sphere with coating. To demonstrate applicability of the above formulation
and to validate computer algorithms, numerical results for the concentric spheres are

presented and compared with the exact solution of Mie series.

As the first example, an electrically-small two-layered sphere is investigated. The

electric size of the inner sphere is [coal = 0.0595 where [to = 3}? is free space propaga-

tion constant and a1 is the radius of the inner sphere; For the outer surface, koa; = 0.13
where a; is the radius of the outer sphere. The inner sphere has the relative dielectric
constant of 8,1 = 16 , lossy tangent tan(8) = 0.39 and relative permeability u, = 4. In the
region between S, and S;, the relative dielectric constant is 8,; = 9. The inner and
outer surfaces are modeled by 92 triangular patches respectively. There are in total
192 triangles consisting of 288 edges which result in a matrix size of 576 by 576. A
plane wave is axially incident along the positive z direction. In Fig.2.13 to Fig. 2.20,
equivalent electric and magnetic currents are shown along two circumferential priciple

arcs that means that the 6 component of K and the 4) component of M are plotted

44

along the arc in q) = 0 plane, and the 4) component of K and the 9 component of M are
plotted along the arc in q) = 12‘— plane. In Fig. 2.13 to 2.16, (K5. Kl, M3,, M3,) are the

four components of equivalent currents on the outer surface S1. In Fig. 2.17 to 2.20,
(K2, K3, M2, M3) are the four components of equivalent currents on the inner surface
S;. The numerical solutions of all these eight components have very good agreement

with the Mie series solutions.

Consider next a perfectly conducting sphere coated with a layer of dielectric,
magnetic, and lossy material as the second example. The elecuic size of a perfectly
conducting sphere is koa = 0.41: and the thickness of the layer is 0.11: . First, we con-
sider a perfectly conducting sphere coated with a layer of air. Figs. 2.20 to 2.24 show
the four equivalent current distributions (K3,, K}; M3,, M3) on the outer surface 5,, and
two current distributions (K3,. K5,) on the perfectly conducting sphere are shown in Fig.
2.25 &2.26. This special case is used to check the coupled surface integral equations
and the results are the same as that from electric field integral equation for a perfectly
conducting sphere. All the results also agree very well with exact solutions. Next, a
perfectly conducting sphere coated with a dielectric layer a, = 4 is considered. Again,
all the six components on the outsuface and on the perfectly conducting sphere are
plotted in Fig.2.27 to 2.32, and the results agree very well with exact solution. In the
last case, the numerical results for a perfectly conducting sphere coated with a layer of
magnetic material [1, = 4 are shown in Fig.2.33 to 2.38. Again all the six components

have good correspondence to the exact solutions.

45

Equivalent Magnetic Surface Currents

 

 

 

 

 

 

08 A A numerical sol.

. — exact solution
0.6-

IJmO I I IE0 I

0.4 -
0.2— t '

0 l I I

0 1 2 3
0 (radian )

Figure 2.13 G-component of equivalent magnetic current on the outer surface S, of a
concentric sphere with:

( koa; = 0.0595; 6,; = 16.0; u,; = 4.0; tan(8;) = 0.39 )
(Icon, = 0.13; e,, = 9.0;m,=1.0. tan(8,) = 0.0)

Equivalent Magnetic Surface Currents

 

 

 

 

 

A A numerical sol.
0-3 " exact solution
0.6 —
|Jm¢ I / IE0 I
0.4 —
0.2 -—
0 [ l I
O 1 2 3

9 ( radian )

Figure 2.14 ¢ocomponent of equivalent magnetic current on the outer surface S, of a
concentric sphere with:

( koa; = 0.0595; 8,; = 16.0; 11,; = 4.0; tan(5;) = 0.39 )
(koa, = 0.13; 8,]: 9.0; [1,1=1.0; tan(5,) = 0.0)

47

Equivalent Electrical Surface Currents

 

 

 

 

 

 

A A A A
0.8 —
0.6 -
If e l / l H o I
0.4 —4
0-2— A A numerical sol.
exact solution
0 T l r
0 l 2 3
9 ( radian )

Figure 2.15 O-component of equivalent electric current on the outer surface S, of a
concentric sphere with:

( [(002 =-’ 0.0595; 872 = 16.0; “72 = 4.0; {311(6)} = 0.39 )
(koa, = 0.13; 2,, = 9.0; u” = 1.0; tan(8,) = 0.0)

48

Equivalent Electrical Surface Currents

 

 

 

 

 

 

A
0.8 —
0.6 —
If,» I / IH o I
0.4 —(
0.2 - A A numerical sol.
exact solution
0 l l l
0 l 2 3

6 ( radian )

Figure 2.16 (tr-component of equivalent electric current on the outer surface S, of a
concenuic sphere with:

( koa; = 0.0595; 2,; = 16.0; M = 4.0; tan(5;) = 0.39 )
(koa, = 0.13; 8,, = 9.0. #31 = 1.0; tan(6,) = 0.0)

 

49

Equivalent Magnetic Surface Currents

 

 

 

 

 

 

A A numerical sol.
0'8“ exact solution
0.6 q
IJmO I I IE0 I

0.4 —
0.2 — a a. ‘3 a :. a

O I I l

0 l 2 3

9 ( radian )

Figure 2.17 O-component of equivalent magnetic current on the inner surface S; of a
concentric sphere with: '

( koa; = 0.0595; 8,; = 16.0; u,; = 4.0; tan(8;) = 0.39 )

(koa, = 0.13; 8,, = 9.0: u,,=1.0, tan(8,) = 0.0)

50

Equivalent Magnetic Surface Currents

 

 

 

 

 

A A numerical sol.
0'8 ‘ exact solution
0.6 -
IJm¢ I I IE0 I
0.4 —
0.2 «—
O I I I
0 1 2

0 ( radian )

Figure 2.18 (tr-component of equivalent magnetic current on the inner surface S; of a
concentric sphere with:

( koa; = 0.0595; 8,; = 16.0; m = 4.0. ran(t‘>;) = 0.39 )
(koa, = 0.13; E,1= 9.0; 11,, = 1.0; tan(5,) = 0.0)

51

Equivalent Electrical Surface Currents

 

 

 

 

 

 

 

0.8—
0.6-—
A A
IJel/lHol . . . M
0.4a
0'2“ A A numerical sol.
exact solution
0 I I I
O 1 2 3

0 ( radian )

Figure 2.19 O-component of equivalent electric current on the inner surface S; of a
concentric sphere with:

( koa; = 0.0595; e,; = 16.0; a,; = 4.0, tan(5;) = 0.39 )
(koa, = 0.13; 8., = 9.0. I», = 1.0. tan(51) = 0.0)

 

527

Equivalent Electrical Surface Currents

 

 

 

 

 

A A numerical sol.
0.8 -, exact solution
0.6 —
A
IJ 4) I / l H o | A

0.4 - A
0.2 -I

0 I I I

0 1 ‘ 2 3

0 ( radian )

Figure 2.20 (to-component of equivalent electric current on the inner surface S; of a
concentric sphere with:

( koa; = 0.0595; 8,; = 16.0; 13,; = 4.0; tan(5;) = 0.39 )
(koa! = 0.13; 871 = 9.0;].171: 1.0; m(81) "'-' 0.0)

53

Equivalent Magnetic Surface Currents

 

 

 

 

0.8 a
0.6 —
A
I M 10 I / I E "w I A
0.4 -—
A
02 d — exact solution
AA numerical sol.
0 I I I
0 1 2 3

0 (radian )

Figure 2.21 O-component of equivalent magnetic current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(Icon; = 0.41:; koa, = 0.51:; 8, = 1.0; u, = 1.0; tan(& = 0.0 )

54

. Equivalent Magnetic Surface Currents

 

  
 

 

 

 

0.8 —
0.6 -—
A
IM it) | / I E "w I
0.4 —
0.2 _ —— exact solution
A A
AA numerical sol.
O I ‘P I
0 1 2 3
0 ( radian )

Figure 2.22 o—component of equivalent magnetic current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(koa; = 0.4x; koa, = 0.51:; 8, = 1.0; u, = 1.0; tan(6) = 0.0)

55

- Equivalent Electrical Surface Currents

 

 

 

 

2 ._
A
1.5 _.
inc
IK 10 | / I” l A
l _
A
0.5 — .
— exact solution
AA numerical sol.
0 I l l

9 ( radian )

Figure 2.23 O-component of equivalent electric current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.
(koa; = 0.41:; koa, = 0.51:; e, = 1.0; u, = 1.0; tan(8) = 0.0)

56

Equivalent Electrical Surface Currents

 

 

 

 

2 _.
A
1.5 -— A
' A
|K1¢|/IH‘"°|
1 _
0.5 - .
-— exact solution
AA numerical sol.
0 I I I
0 l 2 3
9 ( radian )
‘ Figure 2.24 Mmponent of equivalent electric current on the outer surface S, of a

perfectly conducting sphere coated with a lossy layer.
(Icon; = 0.41:; koa, = 0.51:; 8, =1.0; ,1, =1.0; tan(8) = 0.0)

57

Equivalent Electrical Surface Currents

 

 

 

 

A A
2 _
1.5 —
A
I K 29 | / I H "w I
1 - A
0.5 -— .
—— exact solution
AA numerical sol.
o I I I
0 l 2 3

0 ( radian )

Figure 2.25 e-component of electric current on the inner surface S; of a perfectly
conducting sphere coated with a lossy layer.

(koa; = 0.4x; Icon, = 0.51:; 8, =1.0; u, = 1.0; tan(8) = 0.0)

58

Equivalent Electrical Surface Currents

 

IKw/IH’MI

 

— exact solution

AA numerical sol.

 

I I
0 l 2

9 ( radian )

94—1

Figure 2.26 o-component of electric current on the inner surface S; of a perfectly
conducting sphere coated with a lossy layer.

(koa; = 0.41:; too, = 0.51:; 8, =1.0; u, =1.0; tan(5) = 0.0)

-n-» ___ ____.r

 

 

2...”... “..-__-.-..—..».-

59

Equivalent Magnetic Surface Currents

 

  
 

 

0.8 -—
0.6 -
I M 19 I / | E "‘c I
0.4 —
02 q —— exact solution
AA numerical sol.
0 I I I

 

 

O 1 2 3
k 9 ( radian )

Figure 2.27 O-component of equivalent magnetic current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(koa; = 0.41:; koa, = 0.51:; 8, = 4.0; u, =10; tan(5) = 0.0)

Equivalent Magnetic Surface Currents

 

   

 

 

 

0.8 --
0.6 —
inc
IM 14) I / IE I
0.4 -—
0.2 _ — exact solution
AA numerical sol.
0 ' I I I
O 1 2 3

0 ( radian )

Figure 2.28. ,- («component of equivalent magnetic. current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(koa; = 0.41:; koa, = 0.5x; e, = 4.0, ,1, =10, tan(8) = 0.0)

61

Equivalent Electrical Surface Currents

 

IKreI/IHI"°I

 

 

 

 

05 - .
— exact solution
AA numerical sol.
O I I I
0 1 2 3

6 ( radian )

Figure 2.29 O-component of equivalent electric current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(koa; = 0.41:; too, = 0.51:; 8, = 4.0; u, =1.0; tan(5) = 0.0)

62

Equivalent Electrical Surface Currents

 

IK,.I/IH"'“‘I

 

 

 

 

0.5 — .
— exact solutron
AA numerical sol.
0 I I l
0 1 2 3

9 ( radian )

Figure 2.30 @component of equivalent electric cru'rent on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(koa; = 0.41:; koa, = 0.51:; 8, = 4.0; u, = 1.0, tan(8) = 0.0)

63

Equivalent Electrical Surface Currents

 

 

IK2¢|/|H‘"°I

0.5 - .
—— exact solution

AA numerical sol.

 

 

 

0 l l
0 1 2

9 ( radian )

m—

Figure 2.31 e-component of electric current on the inner surface 52 of a perfectly
conducting sphere coated with a lossy layer.

( [(002 = 0.41:; koa, = 0.51:; e, = 4.0; u, = 1.0; tan(5) = 0.0 )

. Equivalent Electrical Surface Currents

 

 

 

 

 

3
2-
IKzel/lHincl
1—
— exact solution
AA numerical sol.
0 I I l
0 l 2 3
9 ( radian )

Figure 2.32 ¢-component of electric current on the inner surface 82 of a perfectly
conducting sphere coated with a lossy layer.

( 15002 -.- 0.41:; [(001 = 0.5m 8, = 4.0. I» = 1-03 tan(8) = 0'0 )

65

Equivalent Magnetic Surface Currents

 

 

 

 

1.5
A A
A
l _
A
lMieI/IE‘MI
0.5 — A
A
— exact solution
AA numerical sol.
0 T I l
0 l 2 3

9 ( radian )

Figure 2.33 e-component of equivalent magnetic current on the outer surface S, of a
perfectly conducting sphere coated with a lossy layer.

(koaz = 0.41:; koal h= 0.5x; e, = 1.0; u, = 4.0; tan(5) = 0.0 )

66

Equivalent Magnetic Surface Currents

1.5

 

|M1¢|/|Ei"c|

 

A — exact solution

AA numerical sol.

 

 

 

l l I
O 1 2 3

9 ( radian )

Figure 2.34 o-component of equivalent magnetic current on the outer surface S1 of a
perfectly conducting sphere coated with a lossy layer. ’

(koaz = 0.43; too, = 051:; e, = 1.0; u, = 4.0; tan(8) = 0.0 )

67

Equivalent Electrical Surface Currents

 

 

 

 

 

3
2 —-4
IKloI/IHWI
1 _
A — exact solution
A
AA numerical sol.
O l I I
0 l 2 3

9 (radian )

Figure 2.35 e-component of equivalent electric current on the outer surface SI of a
perfectly conducting sphere coated with a lossy layer.

( koaz = 0.41:; [coal = 0.51:; e, = 1.0, u, = 4.0, tan(5) = 0.0 )

68

Equivalent Electrical Surface Currents

 

   
 

 

 

 

3
2-
|K1¢|/|Hi"c|
1-
exact solution
A AA numerical sol.
0 I I I
0 l 2 3

9 (radian )

Figure 2.36 «tr-component of equivalent electric current on the outer surface S] of a
perfectly conducting sphere coated with a lossy layer.

( koaz = 0.41:; [coal = 0.51:; e, = 1.0, u, = 4.0; tan(8) = 0.0 )

69

Equivalent Electrical Surface Currents

 

 

 

 

3
A
A
2 —+
leeI/IH'""|
l ‘—l
A
A .
— exact solution
AA numerical sol.
O I I I
0 1 2 3

9 (radian )

Figure 2.37 e-cornponent of electric current on the inner surface S 2 of a perfectly
conducting sphere coated with a lossy layer.

(too; = 0.4x; Icon, = 0.51:; e, = 1.0; u, = 4.0; tan(8) = 0.0 )

70

Equivalent Electrical Surface Currents

 

 

 

 

 

3
2 _
leal/IH‘“ I

1 '1
— exact solution
AA nurr‘ierical sol.

O I I I

0 l 2 3

9 (radian )

Figure 2.38 (Ir-component of electric current on the inner surface S; of a perfectly
conducting sphere coated with a lossy layer.

( koaz = 0.41:; too, = 0.51:; e, = 1.0; u, = 4.0; tan(5) = 0.0 )

71.

2.5 Comparison with Other Methods

One of the existing methods to quantify the electromagnetic scattering by an arbi-
trarily shaped heterogeneous lossy body is to use the tetrahedral modeling method
based on volume integral equations[13]. However, this method is not numerically
efficient when the scatterer has a simple structure. The main advantage of the pro-
posed method in this chapter is that instead of discretizing the whole domain as in the
volume integral method, only the interfaces need to be discretized, thus, leading to
much less unknowns. For example, let’s take a concentric sphere with parameters of
£1: 16,k0a1= 0.0595,€2 = 9,k0a2 = 0.13 , the same as the first example in section 2.4.,
and compare the numerical results between these two method. For the tetrahedral
modeling method, the sphere was modeled by 512 tetrahedral cells leading to 1088
unknowns which led to the results in reference [13]. While for the present triangular
element method, only 576 unknowns associated with 192 triangular cells are used to
obtain the results which agree very well with the exact solution. Furthermore, the
accuracy of the latter method is superior to that of the former method. On the other
hand, these two method provide the same modeling flexibility because the advantages
of triangular elements for surface modeling are analogous to that of tetrahedral ele-
ments for volume modeling [8-9]. However, it is obvious that surface modeling is

much simpler.

The proposed method may be used in many situations where the evaluation of the

scattering by a composite object or the analysis of a coated structure is involved.

2.6 Some Comments

In section 2.2, a set of coupled surface integral equations for a heterogeneous

body either with or without a perfect conductor inside has been developed. However,

72.

the coupled surface integral equations are not suited for a heterogeneous body coated
with a very thin layer on it. It can be seen that as the thickness of the layer
approaches to zero, the equivalent electric and magnetic currents on the two interfaces
of the layer become the same which leads to the resulted matrix become singular. The
failure for this extreme case is due to the formulation of the integral equations. An

alternative approach for this extreme case will be proposed and discussed in chapter 4.

Another thing we observed from the numerical results is that the convergence rate
become pretty slow as the electric size of a body is increased. This is mainly resulted
from the characteristics of moment method which needs to solve a matrix equation.
For example, considering a perfectly conducting sphere coated with a layer of lossy
material, we have six unknowns currents in the resulted coupled surface integral equa-
tions. Supposing each unknown current is expanded into vector basis functions and
represented by N unknown coefficients, then we have 6N total unknowns which leads
to 6N by 6N complex matrix. To get a good resolution of an electrically large body,

the requirement of the storage easily exceeds the memory size of computer system.

For a more complex body, for example, a multilayer sphere or a heterogeneous
body with more than one different media, a similar derivation of a set of coupled sur-
face integral equation can be deduced. However, more equations or more terms in the
set of integral equations may be resulted and the computer codes for solving the prob-
lem has to be changed. This prevents a systematic way to handle an arbitrary hetero-

geneous body.

2.7 Conclusion

To analyze electromagnetic scattering by arbitrarily shaped three dimensional
heterogeneous objects, a set of coupled surface integral equations has been developed

based on the equivalent principle. Using the method of moments, th coupled surface

73-

integral equations are solved by an efficient and simple numerical algorithms. Concen-
tric spheres are used as a test case to validate the theory and computer algorithms.
Numerical results indicate that the proposed method is more efficient than the existing
tetrahedral modeling method based on volume integral equations for simple composite

SU’UCtllI'CS.

CHAPTER III

FUNDAMENTALS OF FINITE DIFFERENCE TIME DOMAIN
METHOD IN SOLUTIONS OF ELECTROMAGNETIC
SCATTERING AND ANTENNA PROBLEMS

3.1 Introduction

The demand for solving electromagnetic problems involving electrically large
bodies and complex structures is increasing in engineering designs. The problems may
involve large composite bodies, cavities, anisotropic media and the bodies may have
dimensions of a few wavelengths. To solve these large complex problems, the finite
difference time domain method provides a good candidate over traditional moment
method. The main idea of the finite difference time domain method(FD-TD) is quite
straight forward. It is a direct solution of Maxwell’s time-dependent curl equations.

The main steps in this method will be discussed here.

The Yee’s model [22] has been used almost universally. It applies simple,
second—order central-difference approximations for both spatial and temporal deriva-
tives of the electric and magnetic fields directly to the differential operators of the
Maxwell’s curl equations. Electric and magnetic field components are interleaved in
space to permit a natural expression of Farady and Ampere’s laws. Space and time
discretizations are selected to bound errors and insure numerical stability of the algo
ritlrm. In addition, the system of equations developed by Yee[22] to update the field
components is fully explicit such that the required computer storage and running time
is proportional to the electrical size of the volume modeled. This sets a remarkable

difference from traditional moment method which needs inversion of matrix.

74

75_

When the fields must be computed is unbounded, as the case of scattering prob-
lems, the infinite space must be truncated into a finite region since it is impossible to
store an unlimited amount of data into computers. A special technique, the radiation
boundary condition, is proposed on the outer truncated boundary surfaces to simulate
the outside extension. For example, the second order radiation boundary condition
introduced by Mur[28] allows all outgoing scattered wave analogs ideally propagate
through the lattice truncation planes with negligible reflection to exit the sample

region.

Overall, the FD-TD method is a marching-in-time procedure which simulates the
continuous actual waves by sampled-data numerical analogs propagating in a data
space stored in a computer. Electromagnetic phenomena such as induction of surface
currents, scattering and multiple scattering, penetration through apertures, and cavity
excitation are modeled time-step by time-step by the action of the curl equations ana-
log.

Recently, the FD-TD method has received more and more attentions because it
has advantages over the moment method in solving problems involving electrically
large bodies and complex structures. Most attentions can be addressed in two aspects:
a) Improvement of the technique theoretically by investigating better radiation boun-
dary condition. The most recent work proposed by Feng [27] is a new method to
improve the radiation boundary condition by introducing a correcting factor. It seems
that this method works more efficiently than others by raising the order of radiation
boundary condition. However, it has been demonstrated only for two dimensional
cases. b) Completing the algorithm for wider applications. By employing the integral
interpretation of Yee’s model, a new method [29,30] was developed for a simple but
efficient modeling of thin-slot coupling, thin-wire coupling, and smomhly curved sur-

faces , thus, overcomes the main drawbacks of the FD-TD method.

76.

The FD-TD method has been widely applied in electromagnetics area. It has
been used to solve problems which include two and three dimensional electromagnetic
wave scattering, electromagnetic wave penetration and coupling for both two and three
dimensions, very complex three-dimensional structures like human bodies and anisotro-
pic objects, inverse scattering reconstructions in one and two dimensional cases, and
rrricrostrip and microwave circuit models. However, little effort has been made in the
application of the FD-TD method to metallic objects with thin material coating and to

transmitting and receiving characteristics of antenna problems.

When solving the electromagnetic scattering problems of metallic objects with
material coating by using moment method or the finite element method, these methods
become very inefficient when the objects become electrically large or sharp curvature
geometries are involved. For some antenna problems, such as cavity-backed antenna
in a infinite ground plane, integral equation or modes matching techniques also become

unpractical.

Main effort of this chapter is devoted towards providing a basis of expanding the
FD-TD method to two important electromagnetic topics by using most advanced algo-
rithms: two dimensional partial coating of metallic objects, and transmitting and

receiving antennas.

In section 3.2, we start from Maxwell’s curl equations to develop two dimen-
sional basic FD-TD algorithm details. It includes Maxwell’s equations, the Yee’s
algorithm. In section 3.3, radiation boundary conditions for two and three dimensional
cases are studied intensively. Different approaches of deriving a radiation boundary
condition are discussed and both first and second order radiation boundary conditions
of two and three dimensional cases are derived. The radition boundary condition for
corner points is also developed. In section 3.4, a systematic way of analyzing the sta-

bility is developed and this method is illustrated through a few examples of different

77

schemes. In section 3.5, dividing the truncated region into scattered field region and
total field region is discussed. In section 3.6, the basic idear of implementing the
integral intepretation of Maxwell’s curl equations to conform the integration paths is
introduced and the application of the idear will be expained in the chapter 4 and

chapter 5 when a thin impedance sheet is taken into account.

As we mentioned above, the theories discussed in this chapter will provide a basis
for the future applications. As a part of chapter 4, basic FD-TD algorithms will be
modified to handle the two dimensional metallic objects coated with thin magnetic
material. In chapter 5, the FD-TD method is employed and a few modifications are
made to study the effects of an impedance sheet on the receiving and scattering

characteristics of a cavity backed antenna.
3.2 BASIC F D-TD ALGORITHMS

3.2.1 Maxwell’s Curl Equations

The FD-TD method is a direct implementation of the time-dependent Maxwell

equadons:

3‘) _ _ .323 __ _
a: _VxH J at _ VxE J,,, (3.2.1)

Consider a source-free region with constituent electrical parameters which are

independent of time, Maxwell’s curl equations can be rewritten into a form:

8—13 = .1..va _ 2E (3.2.2a)
t 8 E

9!. = .ivXE _ an (3.2.2b)
at u u

Assuming e, o, p’ and u are isotropic, where e is the electrical permittivity; o is
the electrical conductivity; p’ is an equivalent magnetic resistivity and u is the mag-

netic permeability; E is the electric field and H is the magnetic field .

78.

If writing the above equations into scalar forms in the rectangular coordinate

systems(x, y, 2), we have a set of scalar equations for three dimensional cases:

3H, 1 3E), 8Ez

 

 

 

 

 

 

 

a: WETE' 9
an 1 as as
__l __. _ __z._ x _ '
a: u( 8x 32 9”?)
8H, 1 8E, 8E
_ = _( __l _ ' )
a: u ay ax ’
as, _ 1(8H, an, E
a: _ a By a: O ")
BEy _ 1 8H, 81-12 E)
at " e 32 ax ° y
as, 1 an, 311,
= — —- - 015,)

_aTeax ay

(3.2.33)
(3.2.3b)
(3.2.30)
(3.2.3d)
(3.2.36)

(3.2.30

Now, consider two dimensional EM scattering problems. If we assume that nei-

ther the incident wave excitation nor the modeled geometry has any variation

in thez

direction, all the derivatives with respect to z are equal to zero ( i = 0 ) such that

82

Maxwell equations are simplified. Due to the linearity of Maxwell equations, any

polarization of incident wave excitation can be decomposed into a linear composition

of a transverse magnetic (TM) mode and a transverse electric (TE) mode. Maxwell

equations finally are simplified to two sets of scalar equations according to TM mode

and TE mode to describe two-dimensional wave interaction with objects. The relevant

equations for these two modes are as follows;
TM mode with H,=0,E,=0, and Ey=0z

3H,_ 1(8E,+ ’H
a: ‘ u ay 9 9

 

TEmode withE,=0,Hx=0, and H,=0:

(3.2.4a)
(3.2.4b)

(3.2.40)

79

 

EiEJ,_1(BHz E (325)
at — 6 8y 0 ,) . . a
BE), -1 3H,
"—3: _ '7: (__8x + 05,) (3.2.5b)
8H, _ 1 8Ex BEy ,
‘a‘: - wire—x " P ”2) (32-50)

3.2.2 Discretization of the scalar Maxwell equations by using Yee’s model

The method of solving Maxwell equations directly is equivalent to the mathemati-
cal problem of solving a set of linear partial differential equations plus boundary con-
ditions with initial values. The derivatives for both spatial variation and temporal vari-
ation will be approximated by using finite difference. The best numerical model avail-
able so far is the one proposed by Yee [22]. The reason is that from the mathematical
point of view , it achieves the second order accuracy in the space and time increments
respectively and also it is a natural geometrical interpolation of Maxwell equations.

These facts will become very clear as we go through this chapter.

In 1966, Yee introduced a set of finite-difference equations for the systems of
Eq.(3.2.3). Following Yee’s notation, we denote a space point in a rectangular lattice
as

(131'. k) = (iAx,jAy, kAZ)
F'(i, j, k) = F(iAx, jAy, kAz, nAt)
where Ax, Ay, and A2 are respectively, the lattice space increments in the x , y, and z

coordinate directions; At is the time increment; and i, j, k, and n are integers.

To approximate the derivatives, central difference formula is used here.

. 1 . . 1 .
amt, .’ k) = F"(t+3.1. k) - F“(z——.J. k)

ax Ax + 0(Ax2) (3.2.6a)

 

. 1 l
815*(1'. j. k) = F 20“.}. k) - F 201, ’0 +0013) (3 26b)
3: At . .

80

Applying above central difference equations to the scalar forms of Maxwell equa-
tion (3.2.3), one may observe that in order to achieve the second order accuracy in
spatial derivative, the components of E and H about a unit cell of the lattice are natur-
ally positioned as shown in Fig. 3.1 , and also that in order to achieve the second
order accuracy in temporal derivatives, E and H are evaluated at alternate half time
steps. Using Eq(3.2.6) to approximate the derivatives of Eq.(3.2.3), we obtain a set of

difference equations for three dimensional cases:

With the definitions of:

22011.02 1'. k) - 9’01 j. 10th
220mm} 1'. k) + 9’0". 1'. k)<:oAt
zzollrU» is 16)
ZZou,(i. j. k) + 9’03 1'. [OCOA‘

 

CA(i, j, k) =

C303 1'. k) =

 

We have:

n+1..11 ..11n-‘..1 1
H, 70,117,“? = CAO, 1+3,k+-7:)H, 70, 1+3, k+—2-) (3.2.7a)
1 1

contain, j+-2-.k+-§) 1 1
_ [E,‘(t j+l, HE) — E20, j. 1:47)]

. . 1 l
lira: J+—’ kl?)AyZO

2
CoAiCBa. j+_; sk+'%') 1 1
+ 1 1 [530,143, k+1) — 530'. j + -2-. k)l
Hr“. 1+3, “30/3120

 

 

l l
”+2'.1.1_.1.l"'z'.1.1
H). (ii-“'2', j, k+—2-)— CA(l+-2-, j, k+-E')Hy (z+—2-, j, [(+3) (3H27b)

coAtCB(i+—;-. j. k+-%-) 1 1
+ [E',‘(i+1, j, k+-2—) — E',‘(i. j. k+3)l

”14%, j, 15%)szo

 

coAtCB(i+-21-. j. k+-;-) 1 1
_ 1 1 [5307, j, k+1) — E;'(i+—2-. j. k)]
11,(i+-2-, j, k+—2-)AzZO

 

81

 

 

 

 

 

 

 

 

 

 

 

 

Ey
1421/
Z
A E Ex
: A El
By /
E Hy
EZA II
EX
(1.1,k) Ey
: y

Fig. 3.1 [Lattice Unit Cell ----Yee’s Model

82.

+ 1 - 1
H: 7 (14%, j+—;-, k): 010%, 14%, k)H: 7 (14%, 14%, k) (3.2.7c)

coAtCB(i+%. 14%. k) 1 1
+ [E;(i+3, j+1, k) — E;(i+—2-. j. k)]

. 1 . l
IMO-['3 . 1+3. k)A)’Zo

 

coAtCB(i+%. 14%. k) 1 1
_ [5;(i+1, j+3, k) — E’;(i, j+-2-. k))

11.04%. 14%. Iowa

 

1 25(z+%, j, k)—Ato(i+%, j, k) 1
gym-7,1, k) = 1 1 520+? 1'. k) (3.2.7d)
2804—3, j, k)+AtO(i+3, j, k)

 

Ml n+ l
+ 2A! [Hz I(l'+l, 121.;1- k) - Hz I(i+—;'.j‘%: k)]

Aytze<i+-;-. j. lemma-is j. k))

 

M 1 8+ 1
_ 2“ [Hy I(i+-L, j, k+-;-) — Hy 7 (14%. j. k—%)l

Az(2€(i+-%- . j. k)+AtG(i+-;-. j. k)) 2

 

zen, j+-1-, k)—Ato(i, j+-1-, k)

l 2 2 l

 

 

 

 

 

 

5;“(1, 1+3, k) = 1 E;(i, 1+3, k) — (3.2.7c)
250', 121;, k)+Ato(i, 1+3, k)
l l
1 2’” 1 [HT7<i+-;—,j+%.k)—HT’<i--;‘Z-.n%. 101
l l
1 2’1“ 1 1H2” 30.14%. k+%)—H:+ ’(z‘. 14%. k—%>]
AZ(2€(i, j+_2-’ k)+At6(i, 1+3: k))
.. 1 . . 1
280,], k+3-)—Ato(t, j, k7) 1
5;“(i.j.k+-)= 1 1 £?(i.j.k+3) (3.2.70
22(i,j, k+3HAto(i. j. lei-3)
l 1
+ 12’” 1 [Hf7(i+—;-, j. k+—;)-H;+r (1%, j. k+—;_)]
Ax(2£(i, j, H; Whoa. 1. k7»
2At ”in 1 Mi.» 1
- H t -a +— -' v _9k+—
1[x(11+2 k 2) Hz (ll-2 2)l

Aytzeo'. 1'. mg >+Aro<zz 1. k7»

83

In a similar way, two dimensional models for the TM or TE mode can be
described in the Fig. 3.2 a or Fig. 3.2 b. Following Fig.3.2 , Eq(3.2.4) and Eq(3.2.5),

finite difference equations for both TM and TE cases result in:

Define:

CAU- J) = 2110.1) + P'(i. DA!
C8011) = 741031) - p’(i. DAt
For the TM case:

E?+l(i.D= 28(1J2—At0(lofl E:(l,j)

 

 

 

 

 

 

 

28(i, j)+Ato(i, j)
2At Mi. 1 "+4. 1
H e —, —H ~ —,
+ M28(i.1)+At0(i.D)[ ’ (”2 D , 0 2 1)]
2A: "+412 . l M}. . 1
- H ~ , — —H, , — ..
Ay(28(i./)+Ato(i.D)[ x (I 1+2) (1 #2)] (3288)
. . 1
n+1 CBO’f'E) "-1 1 2131
”1: 2.0.143): —_'1—(Hx 2-(i.j+—) - [5:03 j+1) -E?(i» DI (32.8w
CA(i, it?) CA(i, j+3-my
n+1 1 Cali-id) "-1 1 2A:
H, Inez-.11: ——1——H, 7<i+-2—.D+ 115’;<i+1.n—£:'(i.m (3.2.8c)
CA(i+—2-. 1) 010+? pm
For the TE case:
. 1 . 1
1 2804-3, j)—Ato(z+3, j) 1
Ez“<i+-2-. n= 1 1 520+? J) + (3.2.921)
220+? mama? j)
1 "+1
1 2A‘ 1 [HT7(i+%, j+2i — Hz I (14%, j—%)]
Ay(2£(i+-2-, puma-5, 1)) )
. . 1 . . 1
1 28(1.J+'2-)fA10(l.J+-2') 1
E;'*‘(i.j+3)= 1 1 E2(i.j+3)- (3.2.9b)
26(1'. ft-i-HAIUU. #3)

84

 

 

 

 

 

 

El 11 HY Ez
55 Ex

1 ”y
El 13Z

( a ) For TM Mode

 

 

 

 

 

 

E
Hz A Y Hz
"Ex ‘Ex
11
BY
H2 Hz

(b) For TE Mode

Fig. 3.2 Two Dimensional Lattice Unit Cells

85-

2A:
M220, j+%HAtO(i. 14%))

»+1 n+1
[.1 . l 7.1 .1
[20212) 20212)]

 

. 1 . 1
CB(‘+3”+2) -%. 1 . 1
1 1 Hz (Pr—2',j+3) (3H29C)

 

n+2- . '
Hz (PI—’17):

+ 2’” [E:(z+- .n— E:<i-—. m

0 1 I 1 2
CA +—, +— A

 

fl . . 1 n . . 1
— 12’” 1 [E,<z.1+3)-E,(z.j-3)1
CA(t+-2-, 143m

 

The increments should satisfy At 5 for two dimensional case

 

l
1 +_1_]2

W? M

l
i
+— -1— +—1—]2

“[3 M A22

for above equations. This will be discussed in details in section 3.4.

and At 5

 

for three dimensional cases to insure the stability

3.3 Radiation Boundary Conditions

Considering the open " problems where the domain of the computed field is
ideally unbounded, that is the whole space, it is obvious that the infinite space needs to
be truncated into a finite space to fit in the finite storage of a computer. The computa-
tion zone must be large enough to enclose the structure of interest, and an artificial
boundary condition must be posed on the outrnost truncated surfaces to simulate the
outside extension of the infinite space. Based on the physical meaning, this boundary
condition should be transparent to outgoing wave from the computation zone theoreti-

cally and numerical implementation of the boundary condition should limit the

86-

reflections of outward propagating waves to some acceptable level. This boundary
condition has been called either the radiation boundary condition or the artificial boun-
dary condition.

From a mathematical point of view, boundary conditions are necessary for solving
a set of partial differential equations. For the " open " problem, Sommerfeld radiation
boundary condition at infinity plays a role as the boundary condition. For a finite
computation zone, the radiation boundary condition at outer lattice truncation interfaces
is needed instead of Sommerfeld condition. We can see that the radiation boundary
condition must be used through the detailed finite difference schemes. From the Yee’s
model and the numerical algorithm of Maxwell curl equations as described in section
3.2.2, the radiation boundary condition can not be directly obtained. As illustrated in
Fig. 3.2, a central difference scheme requires knowledge of the field one-half space
cell to each side of an observation point. We can not use the same scheme on the out-
most lattice plane as interior point because no information for the field one-half space

cell outside the outmost lattice plane is available which is illustrated in Fig.3.3.

In 1966 when Yee first employed the FD-TD method to solve an EM problem, he
used perfect conductor condition, tangential components of E field being equal to zero,
on the outmost surface which definitely caused a lot of reflection. Perhaps that was
the reason that the FD-TD method didn’t become popular because of limited applica-
tions. Later on, a number of better methods were developed to compute the fields at
boundary nodes: Taylar et a1. [31] used simple space extrapolation ; Taflove &
Brodwin [23] used averaged process in an attempt to account for all possible angles of
propagation of outgoing waves; Merewether [32] and Kunz and Lee [33] used far field
approxmation at large distances from the center of the scatteres to obtain an absorbing
boundary condition. All these methods mentioned above are equivalent to the first

order approximation of the radiation boundary condition. They have disadvantages of

87

0 Plane

ace at Y

“““““““““““““

rf

“Outmost Su

““““““““““““““

 

 

 

 

 

 

 

.“-“---“‘-‘--“‘- ‘
O

 

 

.11’:--.----..----
' H

 

 

 

ll
.j.1<)

z

 

E

---‘--““---“--“---

f

2

.................

“““““““““““““““““““““““““““““““““““““““““““““““““““““““““

Lattice Unit Cells Near Outmost Surface ( Y=0 Plane )

88-

causing considerable reflections when the fields near the boundary of mesh do not pro-
pagate in a specific direction ( either the direction normal to the boundary or radial

direction from the center of the scatter ).

Only after Mur [28] in 1981 employed the Engquist and Majda’s first and second
order absorbing boundary conditions [24], which are first one and two terms of Pade
expansion of the exact absorbing boundary condition for electromagnetic field prob-
lems, it was made possible for the FD-TD method to be employed to much wider
applications of solving EM problems. So far the second order Engquist and Majda’s
absorbing boundary condition is the best among the methods mentioned above. This
section will develop the theory and numerical implementation of the radiation boun-

dary condition through different view points based on Engquist and Majda’s theory.

3.3.1. Derivation of Radiation Boundary Condition by Wave Equation

Through Eq.(3.2.3) of section 3.2 and Yee’s model as shown in Fig.3.l, radiation
boundary conditions for the Maxwell’s equations on the outmost mesh are only
required for the three components of electric fields if the outmost surface is placed as
shown in Fig.3.3. The reason is that to evaluate the magnetic fields in the computa-
tional zone, the electric fields needed in Eq.(3.2.3) are also inside computational region
which are available and that to calculate electric fields on the truncated lattice mesh,
we need the magnetic fields outside of the computational region for which we have no
information about them. For the same reason, if the outmost surface is placed half cell
away from that in Fig.3.3, the radiation boundary condition is only needed for the

magnetic field components.

a. The First method

Maxwell equations permit wave propagation in all directions, and after eliminat-

89

ing H field from Maxwell’s equations, we obtain

(83 + a; + a} - c528?)E = 0 (3.3.1)
That means each component of electric fields independently satisfies the three-
dimensional wave equation. If we denote the each component of electric fields as W,

then

(83 + a; + 83 — cazafw = 0 (3.3.2)

As stated in the previous discussion, we need radiation boundary condition at out-
most lattice. This means that a partial differential equation which permits wave propa-
gation only in certain directions, which is called a "one-way wave equation", is

needed.

Consider a function W( x, y, z, t) and the Fourier transform over the space
domain of it as w(k,,, k,. k,, t). Taking the Fourier transform of Eq.(3.3.2), which
yields:

-(k?. + lei + Ié)w(k.. k, k.. t) — cazfiwacp k, k,. z) = 0 (3.3.3)
Solving Eq. (3.3.3) w(k,, Icy, k,, t) , we have:

wag. Icy. k. t) = Mk. k, of” + 3a.. 1., In)?“ (3.3.4)
where lc2 = k? + k3 + k} , A and B are coefficients with parameters (lg, Icy, kz) . We
can see that the general spectral solution of wave equation (3.3.2) are two waves trav-
eling along opposite directions. This can be seen more clearly if we write W(x, y, z.t)
in terms of its spectra:
4,...

W(x, y, 2,0 = i3- [ woe, k). k,. 06—173? dk, dk, dk, (3.3.5)

Denoting 1?: (k,, Icy, k,) and 7’: (x, y,z) , Substituting Eq.(3.3.4) into Eq. (3.3.5),

we have

W(x, y, z,:) = 271:,— ] [A(k,. k, k,)eic°k"jz?+ B(k,, ky, k,)e"c°"'"P?] dk, dk, dk, (3.3.6)

90-

From Eq.(3.3.6), an arbitrary W(x, y, z, t) is a superposition of two waves pro-
pagating in opposite direction , I? and 47 , over the whole spectral domain. Now let’s

consider only the wave in one direction along 1?: Eq.(3.3.6) becomes:

we, y, z,t) = 31—3. 1A(k,, Icy, [egg/CW"? dk, dky dk, (3.3.7)
1: —oo

Denoting A(k,, Icy, k,)/C°“_jk'? by wa(k,, ky, k,, t) , if represents a wave propagating

 

along 7?, let’s consider 9%? and the derivative of wa along the direction 71’, 3512—0 :
83:0 =jc0kwa (3.3.8a)
awa = dwaix + awagz + Bwa 82 (3.3.8b)

8n 8x 8n By an 82 ‘3';

8x av 82

n + '98,: + [93,2]
533-]

an

Comparing these two equations, we observe that by multiplying proper

= —jwa[kJr

= —-jwal?[fix- + y-a-y- +
an an

coefficients, one equation can be related to another:

Mfibffl + jg}; .1.. 2222—] + C0 8W0 = 0 (3.3.93)

 

at [c an an (in

If we denote C(—)1_:_ = is, + 73, + is, , then 53 + s? + sf = C62 . Equation(3.3.7)

becomes:
Bwa 8x av 82 awa _
a: [“3" + Stan + slan] + 3' ‘ 0 (3391’)

This can be interpreted physically as shown in Fig.3.4. For any wa which
satisfies above partial differential equation, waves will exit ideally the plane with its

normal ii propagating along direction 1?. If I? is parallel to n , Eq.(3.3.9) is reduced to :

_awa alive-
an +60 8: _0 (3.3.10)

This is called the first order approximation of the radiation boundary condition

which describes a wave propagating along ii . In other words, by using Eq.(3.3.10), it

 

91

 

 

Fig. 3.4 Graphic Interpretation Of Eq. ( 3.3.9b )

92-

causes no reflection for waves propagating along r1 .

If , in particular but without loss of generality, taking r1 as —-JE direction and
assuming the mesh is located in the region O<x, we give the boundary condition for

the plane x=0. Eq.(3.3.9) is changed to a form:

(__., - -—a—)wal.._.o = 0 (3.3.11)
8x

By using relation 3,2, «1» xi + sf = C62 and noting s,<0 , we have an equation which

is consistent with Eq.(8) of Mur’s[28]:

1

3.3.
a:

For given values of s,, s, , the solution of above equation will determine a wa on

(58; - (c52 — s3 - s3) )waI-o = 0 (3.3.12)
the outer surface which is consistent with an outgoing wave. Since the angles of
incidence of the wave approaching the surface is unknown, an approximation of Eq
(3.3.12) must be made. By approximating (c52 — s; - sf) with Taylor expansion, it will

result in different order approximation of the radiation boundary condition.

By taking first term of Taylor expansion:

1
<1 — (can)2 - (c.2162 = 1 + 0<<cos,>2 + (00592) (3.3.13)

We obtain as a first order approximation:

.3- _ 41 _
(ax Co a: )Walpo — 0 (3.3.14)

It is also a special case of Eq.(3.3.11) when Tr’is -x’.

If we keep the first two terms of Taylor expansion, then:

1
(1 - (cosy)2 — (cos?) 2 = 1 - 0.51<cos,)2 + «0:921 + 0(((Covy)2 + (coax?) (3.3.15)
After substituting Eq(3.3.15) into Eq.(3.3.12), the second order approximation is

yielded:

93—

(c5133, — c5233 + —;—(a§ + 83))Wl-0 = 0 (3.3.16)

The second order approximation of radiation boundary condition for waves
towards other planes ( x = 0, h, y = 0, h. z = 0, h ) can be derived through a similar pro-

cess. They are summarized as follows:

(c5133, — c523} + %(33 + 513))w1,=0 = o (3.3.17a)
(c513; + c523} — &(ag + 83))Wl-h = o (3.3.17b)
(c5183, — c5283 + $433 + abut/1,3, = 0 (3.3.17c)
(c5133, + c523} — gag + 33))WI,..,. = o (3.3.17d)
(c5133, — c5233 + é—(ai + 33))w1g0 = o (3.3.17e)
(c5133, + .5233 — $33 + 83))Wl,=,, = 0 (3.3.170

For a two dimensional problem, we can simplify the Eq.(3.3.17) by taking
83W = 0 or it can be further simplified in such a way : For TM mode , we only have

E,H,J-I, and Eq.(3.2.4a) becomes :

1103,11, = «3,5, (3.3.18)

substituting ( 3.3.18 ) into (3.3.16) with w = E, and noting % = 0, after integrat-

ing with respect to t with E, = 0 for t < O we obtain:

(8.5. - ca‘afi. - (Colic/Dayflxhw = 0 (3.3.19)
For the TE mode, we only have H,, E,, Ey . In this case, we can go through a

similar derivation as the TM case or directly use the equation for the TM mode by
changing E, to H, and E, to -H, based on the complementary properties of Maxwell

curl equations. We then obtain:

. (a,H, — calaH, + (capo/mam, = 0 (3.3.20)

94.

b. The Second Method
Another approach to derive the radiation boundary condition is to use operator

factoring . Rewritting wave Eq.(3.3.2) as:

(3% + 33 + a} - cazafnv = 0 (3.3.21)

Define partial differential operator L as follows:

1

L = 0,2, + 03 + D} —- 70,2 (3.3.22)
Co
Where
82 82 82 82
D§=—;DZ=—;Df=—;oz=— 3.3.23
8x2 ’ ay2 az2 ‘ 3:2 ( )
Now Eq (3.3.21) can be written as:
LIV = O (1.324)
Observe that L can be factored into L+ and L":
.1. .1.
L = U1: = ((03 + 03 + of) 2+:1-D,)((Df + D; + of) Z—Ciop (3.3.25)
0 0
where L“ and L‘ are :
l
L+ = ((133 + D; + of) 2+2142) (3.3.263)
0
i
L’ = ((03 + D; + of) 2—-c-1-D,) (3.3.26b)
0

Applying L‘ to the wave function W will result in an analytical radiation boun-
dary condition which will absorb a plane wave propagating towards the boundary , for
example, towards x=O plane from the right at any angles. On the other hand, operator
L” will perform on the same function as L“ does except for the wave propagating

towards opposite direction, for example, towards x=O plane from the left .

If, as shown in Fig.3.4, we consider the waves out of x=0 plane from the right

half space where the meshes are located, the approximations are made by expanding L‘

95.

into Taylor series:

D 2 D _1_ 2 2
D,(1+(—5’-) +(Bi)2)2 =D,(1+o.5(%-) +0.5(B£)) (3.3.27)

1 I I

D D
For 31<1 and 132-<1 , Substituting it back into L’ = 0:

I X

—

Dl D’2+D’2310W—01+050’205 ‘2
[:(+(Dx) (DIM-C0114; '(Dx)+'(3:))
_ _1_Dt]w = 0 (3.3.28)

Co

Multiplying both sides of the equation by D, , we then have:

((033 + 0.503 + 0.503)—ch,D,))w = 0 (3.3.29)
0

Replacing DZ by using wave equation (3.3.21), it leads to:

(@5203- — 0.50; — 0.505—71ngw = 0 (3.3.30)
0
After rearranging above equation, we have the same second order equation as
Eq.(12) of Mur’s [28] :

(Clef, — c528? + 0.5(33 + 83))WI,._O = 0 (3.3.31)
0

3.3.2 Finite-difference Approximation of Radiation Boundary Conditions

In this subsection, we will do the numerical implementation for both second and
first approximations of the radiation boundary condition and we will discuss this for

two and three dimensional cases separately.

a. Three Dimensional case

First, the numerical scheme of the second order approximation introduced by Mur
will be illustrated here. -According to Mur’s formula, this scheme can be explained

graphically as shown in Fig.3.5 for the three-dimensional grid case at the x=O grid

96

W(O , j, k+1 )

 

 

 

 

W(1, j, k+1
W“), j'lo k) W(Ot j, k ) W(O, j+1, k)
W(O.5, j, k)
. W l, ' l, k
W(1, J‘l, k) W(1, j, k) ( 1+ )
W(Ov j) k‘l )
W0. j. k'1 )
X

Fig. 3.5 Points Used in the Mur’s Second Order Approximation

Y

97.

boundary. Referring to Fig. 3.5, W"(i, j, k) represents the component of E or H field
which is located at the x=0 grid plane and is tangential to it. For example, W“(O, j, k)
can be determined by the previous values of points surrounding it as illustrated in Fig

3.5. In the following, the Mur’s scheme will be derived step by step:

Rewritting Eq.(3.3.16):

(c518; — c5283 + %(33 + 83))Wlpo = 0 (3.3.32)
In the above equation, aiw, 82W, 83W and 83W must be approximated. Mur
employed a finite difference scheme to write Eq.(3.3.32) as:

x 1 1 ..
D‘aDoW’iM - EDwf—(Wa‘jjt + My) + ZDXDXOVBJL '*' Wm.)
2

+ 31-0103) 035. + with) = 0 (33-33)
where the subscript of D indicates different scheme to approximate derivatives, 0
means central difference, + means forward difference, and - means backward
difference. The superscript of D indicates which variable the differential operator is
operating on. For example, 0; means using central difference formula to approximate

the temporal derivative. We will look at the terms of the Eq.(3.3.33) one by one.

The first term aiw is approximated by using central difference for temporal

derivative at r=nAt and by central difference for spatial derivative of x at x = i .

 

 

 

 

2
n+1 n-l
W2:(_2_, j, k) = 2A! (3.3.348)
[ W'"‘(1.j. k) - w"*‘(0.j. 1.)] _ [ “No.1. k) — W*‘(o.j. Io]
___ Ax Ax
2A:

Second term 33w can be approximated by averaging 8,2,W at x=1 and x=0 and then

using formula for second derivative to t :

82w" . 32w"
0, ,k +
( 1 ) a?

 

[
1 . _ a?
W:(-2-o.,9 k)- 2

 

(1. j. k)]
(3.3.34b)

98-

= W"*"(0. j, k) - 2W"(0, L k) + wr-‘(q L k)
A12

 

+ w"+1(1,j, k) — 2W"(1,l, k) + w"-1(1. j. k) ,2
A?

Similarly for the third term 83yW , first take the average of agw at x=0 and x=1

 

and then approximate agw at x=0 and x=1:

2 2
3 Wm. k) + aTgfl—(u. k)]

2

 

[
Wag. j. k) = 8’2

 

(3.3.34c)

 

= W"(O, j+1, k) — 2W"(O, j, k) + W"(0.j—1, k) +
Ay2

W"(l,j+1, k) - 2W"(1,j, k) + W"(1,j—l, k) /2

Ay’

 

The same procedure for 832W .
82w" . 82W" .
—— 0, ,k + 1, ,k
I a 2 ( J ) 822 ( j )l

l . _ z
W'z'.(—2-. ). k) — 2 (3.3.34d)

_ W“(0.j. k+1) - 2W"(O.j. k) + w"(0, j, k—l) +

_ A22

W"(1,j, k+1) — 2W"(1,j, k) + W"(1,j, k—l) /2
A22

Substituting Eq.(3.3.34 ) into Eq.(3.3.32) and solve for W"+'(O, j, k) , we obtain the

 

 

 

 

following time stepping difference scheme at the x=0 grid plane:

coAt-Ax

l . 1 .
CoAH'Ax W” (1’ 1' k) + W"‘ (0.1. k)] (3.3.35a)

w"+‘(o, j, k) = —W"‘1(1, j, k) +

 

ZAx
coAt+Ax

(COAth
+
2Ay2(c0At+Ax)

+

 

[W"(1. j. k) + W"(0. j. k)]

 

[W"(O, j+1, k) - 2W"(O, j, k) + W"(0, j—l, k)

+ W"(l,j+1, k) — 2W"(1,j, k) + W"(1,j—1, k)]

99

(coAt)2Ax
2A22(coAt+Ax)

 

[W"(O, j, k+1) - 2W”(0, j, k) + W“(0, j, k—l) +

W"(1,j, k+l) - 2W"(1, j, k) + W"(1,j, k—1)]

For a cubic cell, Ax = Ay = A2 = 5 . The above formula can be reduced to that is
given in Mur’s paper.

The second order finite difference equations for waves propagating along the
other planes can be derived in a similar way by using the same numerical scheme as

described in Eq.(3.3.33) and Eq.(3.3.17). Those equations are summarized as follows:

For the wave traveling in the direction of increasing x towards the x=1 plane (

only for 5),, E2 ) :

 

 

+1 . 1 Com Ax n+1 n—l
W" (1. j. k): -W"_ (0.]. 10+ [W (0 J. k) + W (1 j. k)] + (3.3.35b)
C0A+ H‘AX
2Ax . .
60AM” [W"(1. j. k) + W"(O. j. k)]
(COAt)2Ax

 

W"0,'l,k-2W"O,',k+W"0,'—1,k
+2Ay2(c0At+Ax)[ (1+ ) (J) (1 )

+ W"(1,j+1, k) — 2W"(1,j, k) + W"(l,j—l, k)]

(comm
2Azz(coAt+Ax)

 

[W"(O, j, k+1) — 2W"(0, j, k) + W"(0, j, k—l) +

W"(1,j, k+1) — 2W"(1,j, k) + W"(1,j, k—l)]

For the wave traveling in the direction of decreasing y towards the y=0 plane (

for 15,, E2 ):
Wn+l - _ n—l - CIGA- A)’ 1 "-1
(l, O, k) — -W (t, 1, k)+ coAt-l-A -——[W’|+ (i, l, k) + W (i, 0, k)] + (3.3.350)
_2_AL_
coAt+Ay [W"(i,l, k) + W"(i, 0, k)]
(CoA‘)2A)’

 

W"°1,0,k —2W" ’,0,k W" '-1,0,k
2sz(coAt+Ay)[ (1+ ) (I )+ (I )

100

+ W"(i+1, 1, k) - 2W"(i, 1, k) + W"(i-1, 1, k)]

(com)sz

[W"(i, 0, [(+1) - 2W"(i, 0, k) + W"(i, 0. k-l) +
2Azz(coAt+Ay)

 

W"(i, 1, k+l) — 2W"(i, 1, k) + W"(i, 1, k—l)]

For the wave traveling in the direction of increasing y towards the y=1 plane (

only for E1, E, ):

coAt—Ay
coAt+Ay
+ ——2—AX—[W"(i, 1, k) + W"(i, 0, k)]
coAt+Ay
(coAt)2Ay
+
2M(COAI+Ay)

 

w"+‘(i, 1, k) = —w"-1(1, o, k) + [w"+‘(i, 0, k) + w"*1(1, 1. k)] (3.3.35d)

 

[W”(i+1, 0, k) - 2W"(i, O, k) + W”(i—l, O, k)

+ W”(i+1, 1, k) - 2W"(i, 1, k) + W"(i—1, l, k)]

(coAt)2Ay
2A22(coAt+Ay)

 

[W"(i, 0, k+1) — 2W"(i, 0, k) + W"(i, 0, k-l) +

W"(i, 1, k+1) - 2W"(i, l, k) + W"(i, l, k—1)]
For the wave traveling in the direction of decreasing 2 towards the 2:0 plane (
only for E,, E, ):

coAt—Az
coAt+Az

 

W”‘(i.j. 0) = —W"“(i.j. 1) + [W"*‘(i, j, 1) + w""(z, j, 0)] (3.3.35e)

2A2
coAt-l-Az

+ (comfm
2Ax2(coAt+Az)

+

 

[W"(i. j. 1) + W"(i. j. 0)]

 

[W"(i+1, j, 0) - 2W”(t', j, 0) + W"(i-1, j, 0)

+ W"(i+l,j, 1) — 2W"(i, j, 1) + W"(i—1,j, 1)]

(coAt)2Az
+
2Ay2(coAt+Az)

 

W"(i, j+1, 1) - 2W"(i, j, 1) + W"(i,j—1, 1)]

101

For the wave traveling in the direction of increasing 2 towards the 2:1 plane (

only for E,, Ey ):

coAt-Az

W"+1 ', ',1 =—W"'l ', 30+
(I1) (I) ) COAHAZ

[W""1(i, j, 0) + W"'1(i, j, 1)] + (3.3.351)

 

2A2
coAt-t-Az
(coAt)2Az
+
2Ax2(coAt+Az)

 

[W03 1'. 1) + W"(I'. j. 0)]

 

[W"(i+1, j, 0) — 2W"(i, j, 0) + W"(i—1, j, 0)

+ W"(i+l,j, 1) - 2W“(i, j, 1) + W"(i—1,j, 1)]

(comm
+
2Ay2(coAt+Az)

 

[W"(i, j+1, 0) — 2W"(i, j, 0) + W"(i, k—l, 0) +

W"(i, j+1, 1) — 2W"(i, j, 1) + W"(i,j—1, 1)]

It seems that with Eq (3.3.35) we can now handle the grid points on the outer
surfaces. However, we observe that at comer points of the outer boundary as illus-
trated in Fig. 3.5, the formulae for the second approximation just derived also use the
points outside computational region. Unfortunately, we can not apply the second
approximation to corner grid points . Instead of the second approximation, a formula
of first order approximation will be derived for the points at grid corners of the out-

most boundary and this will be discussed in section 3.3.3.

b. Two Dimensional case

Derivation of finite difference equation for the two dimensional second order

approximation is straight forward from Eq.(3.3.35) by just setting terms related to

82W =

822 0:

 

As shown in Fig.3.6, for the wave traveling in the direction of decreasing x

towards the x=0 plane:

coAt-Ax
coAt-t-Ax

 

w"+‘(o, j) = -w"“(1, j) + [w"+1(1, ,) + W"“(0. 1)] + (3.3.36a)

102'

 

 

 

 

YA
W(0,j+1) W(1,j+1)
W(0.j) W(1.j)
W(0,j-1) W(1.j-1>

 

Fig. 3.6 Two Dimensional Mur’s Second Order Approx.

103

2Ax
coAt-I-Ax

(coAt)2Ax
2Ay2(coAt+Ax)

 

[W"(1. J) + W"(0. 1)]

 

[W"(0, j+1) - 2W"(0, j) + W"(O, j—l) +

W"(1, j+1) -- 2W"(1, j) + W“(1,j—-1)]

For the wave traveling in the direction of increasing x towards the x=1 plane:

coAt—Ax
coAt+Ax

 

W"+1(1, j) = —W"+‘(0, j) + [w"+1(0, j) + w""(1, 1)] + (3.3.36b)

2Ax
coAt+Ax

(COAt)2Ax
2Ay2(c0At+Ax)

 

[W"(1. J) + W"(0. 1)]

 

[W"(0, j+1) - 2W"(O, j) + W"(0, j—l) +

W"(1, j+1) — 2W"(1,j) + W"(1,j—1)]

For the wave traveling in the direction of decreasing y towards the y=0 plane:

coAt-Ay

n+1 - -1 -
coAt+Ay [W (1, 1) + WI (1, 0)] + (3.3.36c)

W’”1(i,0)= —W”‘1(i, 1) +

 

W"'1,0—2W"',0 W”‘-1,0
2Ax2(coAt+Ay)[ (1+ ) (I )+ (I )+

W"(i+1, 1) — 2W"(i, 1) + W"(i—1, 1)]
For the wave traveling in the direction of increasing y towards the y=1 plane:

coAt—Ay
coAt+Ay
2A2 . n .
coAt+Ay [W"(z, 1) + W (t, 0)]
(CoAI)2A)’
2Ax2(coAt+Ay)

W“l(i, 1) = -W"+‘(z', 0) + [W””(i, 0) + W"’1(i, 1)] + (3.3.36d)

 

[W"(i+1, 0) - 2W"(i, 0) + W"(i—l, 0) +

W"(i+1, 1) — 2W"(i, 1) + W"(i—l, 1)]

104

For the two dimensional case, W represents E2 or Hz for TM mode or TE mode,

respectively.

If we use the further simplified equations(3.2.4) of the second order approxima-
tion, the algorithm is developed for TM mode as follows:

Rewriting Eq(3.2.4)

(8.5. — c5125. - mono/mafia...) = 0
The central difference is used to approximate the derivatives of Eq.(3.3.19) and

8E, BE, ,
the derivatives are averaged. For example, 797 and E- are approxrmated as fol-

lows:

BE. = [521104) - Him. 1) + E20. 1) - Exo. m

 

 

“a: 2.. (3.3-37)
8E 15"“ 1, -E;'1, 15:“ 0, -E;o,

. = 1 z ( 1) ( D+ ( J) ( 1)] (3.3.38)
8: 2A:

In the same way, the other derivatives of Eq(3.3.l9) are approximated such that a

set of difference equations are obtained:

For wave propagating in the direction of decreasing x towards x=0 plane:

n+1 _ COAt _ AX n+1 _ n
E: (0.1) — 530» J) + __th + Arm. (1. J) 52(0. 1)) (33.393)
CZOHOAIAI "+21"

 

n+1 1 1
— H 70, '+— -H. 0,'—-—
2(COAHMAyl . ( J 2) ( J 2)

n+1 1 n+l 1
+ H. z(1.1+3-1- H. 70.17)]

For the wave propagating in the direction of increasing x towards x=1 plane:

1 _ coAt-Ax
5;“(1.1)-E:<o.n+ —c0At+Ax
A céqurAx

— 2(c0At + Ax)Ay

(arm 1) — 8:0. 1» (3.33%)
n+ 1 1 n+ 1 1
(H. Ito. 1+3) - H. {(0.1-3)

M1 1 M1 1
+ H. 70.14;) — H. fur—311

105-

For wave propagating in the direction of decreasing y towards y=O plane:

n+1 . _ - coAt—Ay n+1
E: (to 0) - E20! 1) + CoAt+ A)’ (E:
C2 [A n l
_ 0110A Y [Hy+2'(
2(c0At + Ay)Ax

(i, 1) - E20, 0)) (3.3.39C)

1
. 1 "+7. 1
17,0)-Hy (17,0)

 

M1 +1
2' . 1 " 2' . 1
+Hy ((+3, 1) -Hy (l-E', 1)]

For wave propagating in the direction of increasing y towards y=1 plane:

ET'U, 1) = E'z'(i, O) + wfi’f’la, O) — E'z‘(i, 1)) (3.3.39d)
c + Ay

0A:

_ CéquIAy
2(c0At + Ay)Ax

 

n+1 "+1
2' . 1 2' . 1
[Hy (Ii—5'. 0) - Hy (17. 0)
n+2). 1 13+}. 1
+Hy (1+3, 1) -Hy (1-3, 1)]

Equations (3.3.39), two dimensional second order approximation, are used as radi-
ation boundary conditions in solving two-dimensional problems in this chapter. For
the three dimensional case, this two dimensional second order approximation can not
be applied to the comer points on the outer boundary as shown in Fig. 3.6 so that the

first approximation has to be used. This point will be discussed in section 3.3.3.

3.3.3 Radiation Boundary Condition for the 2D and 3D Corner Points

As discussed before, the second order approximation can not be applied to the
comer points of two and three dimensional outmost boundaries. This can be easily
seen from Figs 3.S&3.6 since the formula also uses points outside the boundary which
are not available. For the corner points, we have to use other schemes which are
equivalent to the first order approximation. This prevents further improvement of the
radiation boundary condition by raising the order of it, thus, searching for a better radi-

ation boundary condition remains an interesting topic.

106

Back to section 3.3.1, the general form of the first order approximation can be

written into a form:

awa _1 Bwa
8n +60 8:

If we assume scattered waves propagating spherically from the origin, then the fi

= 0 (3.3.40)

in the above equation would be taken in the outgoing radial direction for the corner

points. The derivative along the radial direction becomes:

awa awafii+8wa§l+awa§£

 

 

8r = 8x 8r ay 3r 3. ar (33'4”
Where
3J5 = ———"——; (3.3.42a)
3’ V(x2+y2+22)
8v 2
-“ = ; (3.3.42b)
3’ V(f+y2+zz)
82 z
— = ———; (3.3.420)
8r V(x2+y2+zz)
r = V(x2+y2+22) (3.3.42d)
Substituting Eq.(3.3.41) & (3.3.42) back to Eq.(3.3.40), we have:
Bwa x dwa v awa z + c’1 Bwa = 0 (3.3.43)

 

 

 

HWW W76 W 0..
To discretize above equation to difference equation, special attention should be
paid to the partial derivatives of wa to ensure the stability of the schemes which will
be discussed in details in the next section. Concerning the stability, the first order

approximation for the comer points in different regions have the following forms:
In the region (x>O, y>0, z>0):
wa"+‘(i, j, k) = wa"(i, j, k) (3.3.44a)

— coAdr{-i[wa"(i, j, k) — wa"(i—1, j, k)] + i[wa"(i, j, k) - wa"(i, j-1, k)]
Ax . Ay

+ 322mm j, k) - wa"(i. 1'. HM}

1 07.

In the region (x<0, y>O, z>0):

wa’*1(i, j, k) = wa"(i, j, k) - COAt/{— ix—[wmnr j, k) — wa”(i, j, k)] (3.3.44b)
+ L[wa"(i, j, k) — wa”(i, j—1, k)] + —z—[wa"(i, j, k) — wa"(i, j, k—1)]}
Ay Az
In the region (x>0, y<0, z>0):
mafia, j, k) = wa"(i, j, k) — com/{fimfia j, k) — wa"(i—1, j, k)] (3.3.44c)
- -Y-[wa"(i, j+1, k) - wa"(i, j, k)] + l—[wa"(i, j, k) - wa”(i, j, k—1)]}
Ay Az
In the region (x<0, y<0, z>0):
wa"+1(i, j, k): wa"(i, j, k) - coau+ fimflm, j, k) — wa"(i, j, k)] (3.3.44d)
- L[wa”(i+1. j, k) — wa"(i, j, k)] + i-[wa"(i, j, k) — wa”(i, j, k—1)]}
Ay Az
In the region (x>0, y>0, z<0):
wa"+1(i, j, k): wa"(i, j, k) - COAflr{-:—x[wa"(i, j, k) - wa"(i—1, j, k)] (3.3.446)
+ J—[wa"(i, j, k) — wa”(i, j—l, k)] - i[wd'(i, j, k+1) — wa"(i, j, k)]}
Ay Ax
In the region (x<0, y>O, z<0):
wa”+l(i, j, k): wa”(i, j, k) — coAdr{— -/:x—[wa”(z+1, j, k) — wa"(i, j, k)] (3.3.441)
+ j-[WflII j. k) - wa"(i.j-1. k)] - -?-[wa"(I'. 1'. H1) - wa"(I'. j. k)]}
Ay Az
In the region (x>O, y<0, z<0):
wa’”l(i, j, k): wa"(i, j, k) -_ coAt/{éMa’Kz} j, k) — wa"(i—1, j. k)] (3.3.44g)

__z_ . . __ . . __z_ . . _ . .
Ay[wa"(t,j+1,k) wa"(1, j, k)] AZ[wa"(t, j, k+1) wa"(r, j, k)]}

108-

In the region (x<0, y<O, z<0):

wa”+l(t', j, k): wa"(i, j, k) - coAt/{— 7::[wa"(i+1, j, k) — wa"(i, j, k)] (3.3.44h)

__z_ ,,. . _ . . __z_ . . _ . .
Ay [wa (1+1, j, k) wa"(t, j, k)] Az [wa”(t, j, k+1) wa"(t, j, k)]}

We have used the above scheme to handle comer points in our numerical applica-
tions. The results show that this scheme is stable as it is analyzed theoretically in sec-

tion 3.4.
The first approximation for two-dimensional cases can be developed in a similar

way by letting 15";— = 0

There are other ways to handle corner points, one way we want to mention here
is that used by Taflove et.al.[30]. They simulated the outgoing radial waves in a

different way and a stable scheme who also proposed.

3.4 THE STABILITY ANALYSIS OF THE FD-TD METHOD SCHEMES

The finite difference time domain method (FD-TD) has been increasingly used to
solve a variety of open bounded electromagnetic problems. Its stability criterion has
also been established based on the iterative behavior of fields discretized in an interior
space. In fact, the stability of a scheme is attributed to the behavior of fields discre-
tized both in the interior space and on the boundary of the truncated space. This sec-
- tion introduces a systematic way to analyze the stability of a FD-TD scheme when an
artificial boundary is introduced. As examples, a few well known schemes used in
electromagnetic problems are analyzed. It is shown that the conclusion by the theory
is consistent with the empirical rule in numerical practice. But the application of the

introduced method is applicable to general schemes.

109-

In the early years of applying the FD-TD method to electromagnetic problems,
Yee [22] proposed a finite difference scheme and established a stability criterion which
was based on the behavior of fields at the points inside a closed region. Later Taflove
[23] conducted a mathematical derivation and modified this criterion. Nowadays,
engineering designs are involved in many complex models which may include lossy,
inhomogeneous, and anisotropic systems and of which the structure may involve aper-
tures, cavities, and antennas [21]. In most realistic cases, the basic algorithms of the
FD-TD method have to be modified and the study of stability becomes very important

for a successful difference scheme.

As an example, when an open bounded electromagnetic scattering problem is con-
sidered, the infinite space has to be truncated and the truncation is implemented by
enforcing a radiation boundary condition on the outer surface of a truncated finite
space as discussed in section 3.3. Thus the open bounded problem is changed into a
bounded one, which is comprised of Maxwell’s curl equations satisfied at the interior
points and the radiation boundary condition matched on the outer surface. The variant
boundary conditions enforced on the truncated space and different approximations of
derivatives yield various finite difference schemes. The established criterion is no

longer suitable to all numerical schemes.

The stability is a vital factor to determine the applicability of a new modification
of the FD—TD method scheme, when a modification of well known difference schemes
has to be made for a specific problem. Different approximations on the boundary con-
ditions may result in different numerical schemes, thus result in different stability
requirements. It is necessary to introduce a systematic way to study the stability of a

FD-TD scheme before its numerical solution is attempted.

In fact, the stability of a scheme is attributed not only to the property of fields

discretized at the interior points of a truncated space but also to that of fields

110-

discretized at the points on the boundary. . This section introduces a theoretical
approach to analyze the stability of a whole system. One of the popular FD-TD
numerical model is taken as an example which is constructed by applying the Yee’s
model to the interior points, the Mur’s second order radiation boundary condition to
the boundary points except for the comer points, and the first order radiation boundary

condition to the comer points. But the analysis is applicable to any specific schemes.

3.4.1 The Methods of Stability Analysis

Without the loss of generality, the stability is discussed for the model which is
constructed by applying the Yee’s model to the interior points, the Mur’s second order
radiation boundary condition to the boundary points, and the first order radiation boun-
dary condition to the corner points. In the interior region, the set of difference equa-

tions derived by using Yee’s model is stable in iteration if

 

1
A 2 3.4.1
t c[1/Ar2 + My2 + l/Azz]“2 ( )

holds. Eq. (3.4.1) has been mathematically verified by Taflove. However, Eq. (3.4.1)
is not generally true for the stability of the finite difference equations discretized from
a radiation boundary condition. An alternative approach can be introduced here to
reach the same stability restriction for the interior finite difference equations, and to
find the new criterion for the finite difference equations from a radiation boundary con-
dition.
Based on the Yee’s model, a set of explicit finite difference equations was
developed [22]. These equations can be written into a matrix form as
X"+1 = A X" + 8" (3.4.2)
Assume an error E0 is somehow introduced at the initial step. After 11 steps the error

is accumulated up to

111-

E" = A"Eo (3.4.3)
where A is the matrix of the detailed finite difference equations, which are finite

difference approximations on the Maxwell equations and the radiation boundary condi-
tions. The sufficient and necessary condition for the stability of Eq.(3.4.2) can be
stated as [25]

p(A) = maxlkil S 1 (3.4.4)
if A has a full set of eigenvectors with p(A) being the spectral radius of A. If A does

not possess a full set of eigenvectors ( such as A has a repeated eigenvalue ), the

sufficient condition for the stability of Eq.(3.4.2) is then

p(A) = maxlkil < 1 (3.4.5)
If the eigenvalues of the matrix A can be found, the stability criterion of Eq.(3.4.2) is

consequently obtained. The stability study by a matrix method [25] is based on seek-

ing the eigenvalues of the matrix A.

Hereafter the Fourier method (Von Neumann) [25] is introduced for the general
analysis. The Fourier method can be applied to borh the interior finite difference equa-
tions and the finite difference equations on the radiation boundary. It is based on the
separation of time variable and space variable in the Fourier transform domain. The
applicability of the method is justified by considering a function fit, x,y,z), which is a
general function of time and space, it can be represented by its Fourier spectrum as

fir, x,y,z) = (31;? j F(t, e, t), y)e99‘e9°’e‘i*’ dedody (3.4.6)

where cj = {—7, F(t, 0.41.7) is the Fourier spectrum of fit, x,y,z) in spatial domain, 0, o, y
are the space frequencies. If the function flt, x,y,z) is stable as time is stepped, the
functions of F(t, 0,¢,y)ec’éxec’¢yecm should be stable at any time and space frequencies.
Thus, the stability analysis of fit, x,y,z) is resorted to the stability analysis of

F(t, 0,11),y)ec’e"‘ec’¢y cc” at all frequency values.

112

Instead of working on the whole set of the finite difference equations grouped
from both the interior region and the radiation boundary, the Fourier method can be
applied respectively to finite difference equations of interior region and those of the
radiation boundary. Since the criterion of Eq.(3.4.1) has been established and vali-
dated true for the finite difference equations based on Yee’s model in an interior
region, the emphasis here is put on the analysis of the first order radiation boundary
conditions. The Mur’s second order radiation boundary condition can also be analyzed

via the same method.

The first order radiation boundary condition is often used at the comer points on

the radiation boundary [24,26]. It is mathematically represented by:

2W. 4219.-
an H, a. .0 (3.4.7)

where w is a scalar wave function, Co is the light speed in free space, and ii is the nor-

mal direction of the outgoing wave-front. If the scattered wave is approximated as a

spheric wave radiated from the origin, the normal derivative in Eq.(3.4.7) is given by

W Wi+§1§x+éxfia

'37? = ‘3; M ay are a. are (3'43)
and it can be shown that Eq.(3.4.7) leads to
211__x__+fl_2__+§11_z__+caliw_=o (3,4,9)
ax V(x2+y2+22) 32’ V(x2+y2+zz) 3’ V(12+yz+zz) 3‘

A discrete form of Eq.(3.4.9) can be created by substituting finite differences for
the derivatives in Eq.(3.4.9). A proper forward, or backward finite difference is used
to avoid the possible occurrence of the instability. For example, in x > 0 and y > 0 and

z > 0 region the discrete form of Eq.(3.4.9) is written as:

W‘a. j. k) = W”(I. j. k) —_ coAt/R{§Mi. j. k) — “MI-1.1. HI (3410)

.1. . . _ . ._ _z_ . . _ . . _
+ Ay [#0.]. k) w"(z.1 1.k)1+ min/Kw. k) WWI/c 1)]}

113

where the superscript n indicates the n"' step of time, and R = V)? + y2 + 22. As dis-

cussed previously the stability analysis by a Fourier method is alternated to analyze the

. ex . . .
wave functions of F (t, 6, (11,7) ec’ ecmec’fl. Therefore, the separatron of variables of w

can be assumed as
W(i, j, k) = rwa, j, k) (3.4.11)
where T' is the function of time at the 11’” step, and W is the function of space. Substi-

tuting Eq.(3.4.11) into Eq.(3.4.10) yields:

 

Tn+l
= 3.4.12
7" It ( a)

+ (W(i, j. k) — W(i. H k»? + (W021. k) — Wu. 1. lem—’3‘?
where Ax = Ay = A2 = 8 is taken, 11 is a constant to be determined from (12.b), and 1:
equals coAt. It can be deducted from the linearity of Eq.(3.4.12) that if an error e is

introduced at the initial step, the error function obeys the same equation as w, thus

the error function may be represented in the same form of Eq.(3.4.11)

e: j, ,, = r'wa, j, k) (3.4.13)
from (3.4.12.1) the function r' is related to the initial 7° by

T" = T'RTO (3.4.14)
7° can be assigned to unity without the loss of generality. The boundness of the error

function Eq.(3.4.13) is thus turned to the establishment of Inl s 1, which is determined

by Eq.(3.4.12b).

Consider all the possible solutions of n in Eq.(3.4.12b) which are represented in a

Von Neumann’s form

W(i, j, k) = ec’m em ec/q (3.4.15)

where 9, o, y are real space frequencies in the Fourier frequency domain, and CI = 11:1- .

114

Now 11 in Eq.(3.4.12b) is dependent on 9, (1), y and % since different values of n will

be resulted from Eq.(3.4.12b) at different frequencies.

It is seen that the numerical scheme Eq. (3.4.9) is stable for a given % if

111(1), t1), 7, %)l s 1 (3.4.16)
holds for all real 6, o, y .

At the initial step of the iteration, the error function is

e2 j. k = W(i, j, k)?“ = W(i, j, k) (3.4.17)
Using the Von Neumann’s form Eq. (3.4.15) for W(i, j, k) in Eq. (3.4.12b) results in:

n = 1 -- RM [0' + j + k) - (1.396 + jecl" + ke‘”)]

 

 

 

 

r(t' + j + k) r . .
= (1 — + [10059 + jcoso + kcosy]
W12 +j2 + k2) Vuz +j2 + k2)
r . . .
— c [isrne + jsmo + ksrny] (3.4.18)
jW?+fl+H)

where r is defined as coAt/S. The amplitude of 11 should be bounded by:

 

 

|n|2=r(l— r(i+j+k) + r

2
(isine + jsino + ksiny)] S 1 (3.4.19)

2
(icosG + jcoso + kcosy)]

 

 

01'

2
[(1102 + f + k2) — r(i + j + k)) + r(icos9 + jCOS¢ + kcosy)]
+ [r(isin9 + jsin¢ + ksiny)]2 s (i2 + 1'2 + k2) (3.4.20)
should hold for any combination of 9, t), 7. After a few steps of algebraic manipula-

tion, an equivalent form of Eq.(3.4.20) can be obtained:
r(i + j + k)2 + 12 + 1'2 + k2 + 2r(ijcos(9 - (1)) + ikcos(9 - y) + jkcos(¢—y))

+ 2(iCOSO + jcos¢ + kcosy)[\r(i2 + f + k2) - r(i + j + k)]
52 11(12 + 1'2 + k2)(i + j + k) (3.4.21)

 

115-

It can be shown that Eq.(3.4.21) holds if

\/(i2 +1"2 + k2) - r(t' + j + k) .2 0 (3.4.22)
To find the solution of Eq.(3.4.22) for r, it is helpful to introduce the relation

 

(~1— (3.4.23)
m i=1 m
It is deduced:
v? +12 + k2 2 713—0 +j + k) (3.4.24)
The possible range of r for Eq.(3.4.22) to be valid is restricted to
l
g _ .4.
r )5 (3 25)

then the satisfaction of Eq.(3.4.21) is provided by Eq.(3.4.25). Finally it can be stated

_1_

that if r5 ‘5

, the me, o, y, r)| s 1 holds, then the equations created from Eq.(3.4.10)

are stable.

It is informative to see that the criterion (3.4.1) established from the interior finite
difference equations by the Yee’s model is reduced to (3.4.25) when Ax = Ay = A2 is
forced. This criterion also holds for the Mur’s second order approximation on the
radiation boundary condition. At this point, it can be concluded that (3.4.25) is the
stability criterion for the system specified in the introduction section. As a matter of
fact, the conclusion has already been verified by the numerical practice of applying the

FD-TD to three dimensional scattering problems in transient or steady state.

Next two different approximations on the first order radiation boundary condition
for two dimensional comer points are presented to investigate the stability. The first
scheme to be presented is the analysis on the scheme we used in the research. The

first order radiation boundary condition was proposed in the region x<0 and y<0 as

_13w 131+fl

Co E _ -‘j-2:( ax ay ) = 0 (3.4.26)

116-

If the forward finite differences are used to replace the derivatives in Eq.(3.4.26),
Eq.(3.4.26) will lead to a stable scheme. But a different approach may result in an

instable scheme.

First, Eq.(3.4.26) is discretized by forward finite differences into the form

ficallW‘ti. f) - Ma. mm: =
[w"(i+1, j) - W(i, 1)]/Ax + [w"(i, j+1) — W(i, j)]/Ay (3.4.27)
Using Eq.(3.4.11) in Eq.(3.4.27) yields

r nwa, j) = [(r - 2)W(i, j) + W(i+1,/) + W(i, j+l)] (3.4.28)

where W is a two dimensional function, and r is defined as r = NEAx/(coAt). Subse-

quently the W is replaced by its Von Neumann’s form, The equation

0 = [(r - 2 — m) + e619 + e°i°1 (3.4.29)
is obtained. The stability Eq.(3.4.26) requires l n I s 1. It is equivalent to require

[r — 2 + c056 + COS¢]2 + [(sine + sint11)]2

Inl2 = r2

s 1 (3.4.29)

The solution for r is shown to be

r (2 - (case + cos¢)) 2 3 - 2(cosO + cos») + cosOcoso + sinesino (3.4.30)
When cost) = 1 and cost) = 1 , both sides become zeros and equality holds. Other-
wise

, > 3 — 2(cose + coso) + cosOcoscp + sinesinq)
- 2 — (c059 + cost)

It can be proven that the right hand side of Eq.(3.4.31) is bounded by 2, thus r 2 2

(3.4.31)

guarantees the satisfaction of Eq.(3.4.29). Note that r 2 2 implies coA S 5N2 which is

the same criterion for the finite difference equations of the interior region points.

However, if Eq.(3.4.26) is discretized by the finite differences suggested in [24],
an unconditionally instable scheme results. This can be demonstrated as follows.

Assume the finite differences are taken in Eq.(3.4.26) as denoted in [24]

117-

(D‘ + Dy )(uf; + 14’1“) — —‘J2COD'+(u?J+l + 111:.” + u3)- — (3.4.32)
The operator Di denotes forward finite difference on 1: variable, and the similar nota-
tions are used for D1 , Di, . This algorithm seems better because derivatives w.r.t. the
space coordinates are averaged on the two time instants t= nAt and t = (n+1)Ar, and
the derivative respect to space is averaged over three locations (1', j), (1+1, j), and

(i, j+1). But it creates an instable scheme.

The explicit form of Eq.(3.4.32) is written as

 

 

 

W(I11- — -E—'—;2) +2) m 11+ EL— ’3 (Wan 11+ W(I 1+1»
+-E———1+ + 2; ———(w"(i+l, j) + W(i, j+1)) (3.4.33)
Using of Eq.(3.4.11) gives rise to the equations
1
11.: I}; (3.4.34a)
i—lw WI.( 11+ El—+—'l(W(I+1 11+ W(I.1+11)
( + 2)
T1= (3.4.34b)
we: 1. k1- Lamar. 11 + W. 1+1»
(r + 2)
If the W is substituted by its Von Neumann’s form, then 11 is given by
691 11'
_(r - 2) + (1 + r)(e + e ) (3.4.35)
(1+2) _(1 - r)(e’8 + e”)
For the scheme to be stable, 11 should be bounded by
W2 _[(r- 2) + (1 + r)(cose + cos¢)]: + (sine + smb):(1 + 0:3 <1 (3.436)
[(r + 2) - (1 - r)(cos9 + cos¢)] + (sine + sine) (1 - r)
Unfortunately the further expansion of Eq.(3.4.36) yields
—(0039 + coso) + 2(sinosin9 + cosecoso) S 0 (3.4.37)

which is independent of r, and the inequality is apparently faulty at many values of 0
and (p. The criterion of Eq.(3.4.36) is no longer satisfied for any 0 and 11. Thus the

scheme constructed by Eq.(3.4.32) is not a stable one. This conclusion agrees with the

118

numerical practice in which the practical numerical application of Eq.(3.4.32) experi-

enced the instability.

Several other modifications used by Taflove et al.[26] have also been verified by
the proposed method, the analysis results are consistent with their empirical conclu-

sions.

3.4.2 Summary

The Fourier method has been introduced to analyze the stability of a FD-TD
scheme. Instead of investigating the stability of a whole numerical scheme, the intro-
duced method can be applied respectively to the interior region, the radiation boundary
and the corner points. It is useful when a modification on the boundary condition, or
on the finite difference formula for the interior points is required. The stability
analysis becomes extremely important when a FD-TD method is used to solve steady
EM problems since a steady state takes long time to be reached, and the knowledge of

the iterative behavior of the FD-TD scheme helps save numerical computation efforts.

Several popular schemes used in the application of the FD-TD to electromagnetic
problems have been analyzed via the method. The deducted conclusions for stability

are consistent with the numerical practice.

3.5 Total Field Region and Scattered Field Region

Due to the linearity of Maxwell’s equations, the numerical algorithms derived
from the Maxwell’s equations can be applied to incident EM fields, scattered fields or
the total fields. On the outmost surface, radiation boundary condition has to be applied
to the outgoing waves where in most cases they are scattered fields. While on the
interfaces of different mediums, tangential components of total B and H fields must be
continuous across the interfaces. One technique to treat these two cases is to zone the

numerical space lattice into two distinct regions, as shown in Fig.3.7, separated by a

119

Radiation Boundary

 

Scattered Field Region

 

Total Field Region

 

 

 

 

 

 

 

Fig. 3.7 Total Field And Scattered Field Region

120

rectangular surface which serves to connect the fields in two regions[34].

Total field region is the inner region of the FD-TD lattice. It includes all interact-
ing structure of interest. The finite difference system for the Maxwell’s equations

operates on the total field vector components.

The outer region of the FD-TD lattice is denoted as the scattered field region. In
this region, the finite different system for the curl equations operates only on the scat-
tered field vector components. Radiation boundary condition can be applied directly to

the points on the outmost truncated surface.

Dividing the FD-TD space lattice into two regions, it provides some convenience.
In addition to the advantages of applying radiation boundary conditions and boundary
conditions across the interfaces, the incident plane waves can be generated on the con-
nection surface of two regions to insure the consistence of the fields in using the finite

difference systems.

3.6 Integral Interpretation of the FD-TD Algorithms

In most realistic cases, many problems involve thin wires, slots and curved sur-
faces. However, the Yee’s algorithm for the FD-TD method was originally interpreted
as a direct approximation of the pointwise derivatives of Maxwell’s time dependent
curl equations. This resulted in a staircase approximation of the curved surfaces which

might limit the predictive powers of the FD-TD method.

A few work has been done in this area to come up with new FD-TD algorithms
for curved surfaces. Recently, a simple but efficient technique has been used success-
fully based on the integral interpolation of Maxwell’s curl equations. Maxwell’s curl
equations can be written in integral forms which are called Ampere’s law and Farady’s

law.

111.11: 11.12.11. $11.11 (3.6.1a)

121

113-171 = 11,113 — 337111-113 (3.6.1b)

As shown in Fig.3.8, the above integral equation can be implemented in a
geometrical interpretation. These contours intersect in the manner of links in a chain.
With this geometrical interpolation, surface curvature can be conformed by deforming
contour paths. Specifically, if we take the contours as shown in Fig.3.8, we find out
that it yields the identical Yee’s scheme. Here we take the Ampere’s law as an exam-

ple. The equivalence of Yee’s formula and contour integral can be seen as follows:

iJD-dS=¢EH-dl
at 1 1

As shown in Fig.3.8, we apply the Ampere’s law along contour C1 . Assuming
that the field value at a midpoint of one side of the contour equals the average value of

that field component along that side, the right side of the equation becomes:
n+2!» 1 "+21- 1 ”+21' 1
J H-JI = H, (1, j—E, k)Ax + H, (1+3, j k)Ay — H, (1, 1+3, k)Ax (3.6.2)
1

n+ l 1
- Hy 7(1—3, j, k)Ay
If we further assume that E,(i, j k) equals the average value of E, over the surface

S1 , and Ax: Ay = A2 = 5 . The time derivative can be numerically realized using a

central-difference expression, and left side of the equation also can be simplified:

 

.3. . _ 2 E7103], k)-Eg(i,j, k)
at JD ‘15 “ ‘05 [ A, (3.6.3)

After equating (3.6.2) and (3.6.3), the equation becomes:

82 $10.91.! k) - 5:0: j, k)
80 At

 

 

 

n+ l 1 11+ 1 1
H, I(i, j—-2-, k) + H, Tat-3, j k)

n+21- 1 12+} 1
_ Hz (i1j+-2_1 k) — Hy 0-31.1.1 k) 8 (3°6°4)

where the superscript indicate field values at time steps 11, 111%.- and 11+] . Rearranging

122

Ey(i. j+1/2, 1t+1/2)

 

 

 

  
 
 

    

 

 

 

 

 

 

 

    
 

 

    

 

 

( a)
11 i . .
EZ(I.J+1.1<)
Hx(i, j+1/2. 1:)
d1; 7
Hy(i-1/2.j.k) AZ 1) Ez(i.j.k)
Ay (HI/2.1.10
S1 T tis1 Hy

A c1

r7 = \

11x6. 1:112. k)
Ey(i. j+1/2.k-1/2)
1
11 y(i. j+1/2. k+1/2)
( b )
1 ll . .
1 HZ(I.1+1.k)
Exa' j+1/2.k)
d1;
Ey(i-1/2.j.k) Al 1) 1126.1. 1:)
A . .
82 y I d8 Ey(l+1’20 Jo k)
2
Ax +— C2
F] = \

 

Ex(i. j-1/2. k)
HyG. j+1l2. k-l/Z)

 

 

Figure 3.8 Examples of spatially orthogonal contours in free space. (a) Ampere’s
Law for E, ; (b) Farady’s Law for H, .

123

the above equation:

15;“(1. 1. k1: E21131. k) + —A‘—,
208

 

11+1 1 n+1 1
H, I(1, j—E, k) + Hy r(143,1 k)

11+ 1 n+ l
_ H, 7(1, 1%, k) — H, 704%, j, k)]fi (3.6.5)

We get exactly the Yee’s expression for E, , for the free space case, which was

derived directly from the Maxwell’s curl H equation.

Similarly, we can apply the Faraday’s law along contour C2 in Fig.3.8b. which

yields the same expression as Yee’s from the Maxwell’s curl E equation.

By using integral interpretation of the curl equations, the Yee’s model can be
modified to handle the surface curvature, infinite thin sheets. In next Chapter, we are
going to use integral interpretation to develop algorithms for two dimensional scatter-

ing problems.

3.7 Numerical Implementation and Validation of FD-TD Modeling

The main applications of the FD-TD method we concern can be stated as: a) Two
dimensional metallic objects coated with magnetic material; b) Three dimensional cav-
ity backed antenna with an impedance sheet over it. Based on the theory and algo—
rithms discussed in this chapter, a few modifications will be made to solve above prob-
lems. Since a lot of details are involved, these two topics will be investigated

separately in chapter 4 and chapter 5.

In this section, only a few well known examples have been chosen to validate the
computer program and the standard FD-TD algorithms. As the first example, a
infinite long perfectly conducting square cylinder is tested to verify the program and
algorithm. As shown in,Fig.3.9. , each side of the square cylinder equals to koa = 21: .
EM wave is propagating along y direction with TM polarization and both solutions of

FD-TD and MOM are plotted in Fig.3.9. Another simple example is a perfectly

124

conducting strip as shown in Fig.3.10 and a strip of thin film with variable conduc-
tivity n = 1/(orZ0) = 2(1r/a)2 as shown in Fig.3.11, where a is the half length of a strip.

It can be seen that the results of both FD-TD and MoM agree with each other very

well.

The capability of FD-TD in solving the electrically large bodies can be seen in
the following example of a two dimensional electrically large airplane wing. The air-
plane wing is approximated as shown in Fig.3.12 and field distribution is plotted in
Fig.3.13-3.15. The length of the airplane wing is about six wavelength and one end of
it is very sharp. We can see the field singular behaviors around sharp comers and
some small discontinuity of fields resulted from the staircase approximation of the sur-

face of an airplane wing.

125

 

 

 

 

   

 

 

 

3.4—

3 q
2.6—
2.2‘ b 1

11, l/IH‘ I

1.8-
1.44 .

1...
0.6—

FDTD solution
0.2 - -— MOM solution
I - ! I L -
0 60 120 180 240 300 360
9 (in degree)

Figure 3.9 Current distribution on an infinite long, perfectly conducting rectangular
cylinder without coating in the case of the TM excitation ( koa = 21: )

126

 

 

 

 

Currents on the Thin Film
5
9 (D
' I
4 —: AA Integral Equation 1
I 00 mm MCIhOd :
i :
1 i
: 1
3 .1.) .'
i :
1 I
I
I1, 1, ”10' ‘1 :
I I
1
2 _ f. 1
1| :
1'. ‘~ '
.t *4
1-
0 I I I
-2 0 2
x/A

Figure 3.10 Comparison of results on the induced current on a perfectly conducting
strip using integral equations and the FD-TD method

127

 

 

 

 

Currents on the Thin Film
2.5
an Integral Equation
2 ..
oo FD-TD MCIhOd
f \ '
I
1.... ,1 ‘1 .
I K
I \
I] IH l ‘ ‘
2 ll 0 ’1 \
t
i- ,t )I
.1 V
I, \
‘8
x1" 3‘ 1’
05 — ‘. I '
1r - 1! h
0 l I F
_2 0 2
x/k

Figure 3.11 Induced current on a strip of thin film with variable conductivity:
11 = l/(OtZo) = 20:11:)2

128‘

 

 

Figure 3.12 A perfectly conducting wing-shaped cylinder

129'

  
   

     
  

‘I “I .
11313111111 .
"WWI“ ”le \\\ “:11 M11)" M

 
       
   
 

     
 

  

  
  
  

y", \\\) I\

.7 '.‘

5.39543 “O;

‘4'“ 'I :.-.. ii I'II"'\\\ "i \M“

,;,}:.~;.'~,,,,..,,..‘..,,,., IIIII’I 11,111,“. \ 23 529
'II II)" I" "NIH,” ," i." 11,114.41, 1 Elk)“
a +1622 7 (1&1, ‘.~,.‘;: ' "',".'I,'I ’".,‘ZIII ”III” "5" ‘
\\\\\\ ‘\.\ \\\\.'o' "If? h ‘i‘ ‘ "J"
(IQ \\\ " ‘ ". “HE/11’" ,111' .1 765

    
  

     

 

" \‘.\\\\
\\.’\\\\\

  
  

 

‘1‘) t" I \\\\\\"\\\ "9
\\\\\\\ M“ . ‘\\\\\\\, “(\\\

           
    
     

    
        
    

III)
I

          

     

         
 
  

\\ \\\) \\\,\\\\ \\\i‘ 1 ’ .
0 0000J W)» 1)“ EH \Mii‘. ’19". III“ 0 0000
[le WEIWAIIE i111“, .9
— 56176 1\ _1, 765 \s
11 484 M" \"hk‘ \\\:ka
1096* 23 529 -23 :29

Figure 3.13 Amplitude distribution of the z-component electric field in the total field
region with koa = 20

Hg/YU

130

             
  
    
 

I ~'

.‘I‘ \ 1'\,\\\

I“ “I0 \I“ III“ “I

IIII‘N III II“
\\\III "\\\”. ‘1’.

I“

“I“

\I“\\,~’\ \
I‘I‘NIII“ .1
“I “I.“

   

          
    
    
       
    

  

 
  
    
   

     

  
 

    

  

o \IHW
:\¢.'I\‘\\I1MIIIII’\\ I'IIIII“ “II‘ “:1“
\\\;I \\ .0
. .3 “.~“ ”II“ I“ \\\I \IIIIWI“ “3'
, 1'1"“. l;s‘.o'..',..’:0:1;.:..
8 80 38 j \\\‘\\\\.\ ',21’1‘173’513 ".",I"I.1'"ll',':l’ 33 589

   

“\\\“\\\“\\\\\\\\\\w§i

“\\\

.\\\,\ .,
III“
0 0000 J18 3

I \,i
III“

' I
'".-fifi'II'Ii‘Il‘IIII

    

'1:'1;::'\;"11‘1/II'II”/ I 152:.“
III‘i' I'll/III! III

    
  

Figure 3.14 Amplitude distribution of the y-component magnetic field in the total
field region with hp = 20

Hx/YO

131'

    

 
       
   
    
       
  

.«9. I‘
'I1‘ .51,
”'1

   

 
  

      
 
 
 
  

     
 

  
  

      
     

     
  
 
  
     

  

   
   
   
 
       
   
    
     

     
  
   
  
 
 
  
 

 

   
  

 
 
  

jam ~
' {..éfl.‘ '1 o‘ 1"“ "‘
1'1.'1I'10“'0\\ ll“
"(3.1/5:11 M231}! :« 'WI H. I"
3‘37”? W" ;7'l n!”’1“.‘l.‘l¢&l'&‘f'$"“ 41‘4““
,,Il_lll ‘5’ '"...." “H‘ “$3!“ 9‘5“}, 11,1. 1_ $\\ \ IO; \Ii)‘\\v:”'
$7,," 'I'I'I"‘1‘\:: 7&3; ’ "9'13"" s"1\\\"II':‘\ 9,1”,
:5...- III'.’~"'I'1.I1 '1‘.” ‘1'?)“(55' “\\\I ‘ \\II ‘}1‘.I'I" "I!”hul’l'lb’lt
' .1.'.~.\ '.',’,";f1' ,,‘\ . “"., "‘11.“: 7". \\\ «C8 \D1‘ *1. “v. ’,I lull/Ii] I
MW”, ll” 'I'IIN' " 337’” “"""' I'I: 23 235
IIiI'IiiiI ”W H "N {am ' u
"f L. ’1',“ 3 In":
, VIII?"
:m‘u I’.’ w. II,','II ‘11“ M.
”3'? 1\.\\:3,.,’!,‘.”" 0““ ‘ 1 1 618
‘\\""1\\v {:1}. 15%.; I ‘ \ ‘53.’...‘ 1‘ ‘b ‘
’95“: Guy.“ ”'0 ‘33};1'291": 3‘
'II'I '
I,

\l I '\\\:;':':’l"
1

V;.. '33“: I1 "“...ZIW/
\\' I. ‘ .1., ‘
.., ~ «m III

—12 607 \\\: \\
\\} N

 

    

       
 
 

   

   

   
      

“"3 "\\\ II\\/Il1:s;'lm,\;\\, “'l O 0 000
c.”
0.1% \\\ II1\\\‘§“/55‘~~!o :III'IL ' >9
— qe176 ‘ ' \
\a‘tl' 555%“ —11 618 Q
11 L+8L+ "I
4096* 23 529 —23 235

Figurc 3.15 Amphtudc dlstnbuuon of the x-component magnetic field in the total
field region with koa = 20

CHAPTER IV

APPLICATION OF FINITE DIFFERENCE TIME DOMAIN
METHOD AND MOMENT TO TWO DIMENSIONAL
ELECTROMAGNETIC SCATTERING PROBLEMS

4.1 Introduction

It is known that radar cross section of a conducting body can be reduced if it is
coated by an electrically or magnetically lossy layer. In practice, it is desirable to
make the coating layer thin. A thin layer of electrically lossy material on a perfectly
conducting body can not reduce its radar cross section because the tangential com-
ponent of electric field is very small near the surface of a perfect conductor and conse-
quently the induced current and dissipated power in the coating layer are very small.
On the other hand, if a thin magnetically lossy layer is used to coate the body, its
radar cross section can be significantly reduced because the tangential component of
magnetic field is very large on the surface of a conducting body, resulting in a large

equivalent magnetic current and a high dissipated power in the coating layer.

In this chapter, a new set of coupled integral equations is derived for treating the
scattering problem of a perfeCtly conducting body coated with a thin magnetically
lossy layer. These electric field integral equation and the magnetic field integral equa-
tion are deduced and numerically solved by the method of moments (MoM). To vali-
date the derived integral equations, an alternative method to solve the scattering prob-
lem of 'an infinite circular cylinder coated with a thin magnetic lossy layer has also
been developed based on the eigenmode expansion. The results of the radar cross sec-
tion and currents via the MoM and the eigenmode expansion method are compared.

The agreement is excellent. The finite difference time domain method is then implen—

132

133,
mented to solve a metallic object coated with a magnetic thin layer and numerical
results are compared with that by the MoM.

4.2 Integral Equations for Perfectly Conducting Cylinders Coated with Thin

Magnetic Materials
4.2.1 Derivation of Integral Equations for three-dimensional structures

A new set of integral equations is derived based on the equivalent principle,

focusing on the extreme case of a perfect conductor with a thin magnetic coating.

Let’s start from Maxwell’s Equations:

VxH = -82 +Jc (4.1a)
a:
BB
V = — — — ,
XE at J,,. (4 1b)
VD = p: (41C)
V-B = pm (4.1d)
and continuity equations of
V-J + 1p = 0 (4.1e)
e at e
B
V- — = ,

In the frequency domain, the fields are preassumed to be harmonic dependence of

e’m‘. Thus the Maxwell’s equations become:

VxH = 1'er + J, (4.2a)
VxE = - jtouH - J,,. (42b)
V'(8E) = Pe (4.2C)
V-(uH) = 9». (42d)

where E and H are the electric field and the magnetic field, 8 and u are permittivity

and permeability of the media, J, and J,,, are the electric current and magnetic

134

current, and p, and p", are the electric and magnetic charges. The above equations

can also be rewritten as:

VxH = jweoE + jo)(e - £0)E + J, (4.3a)
VxE = -j<ouoH -j<o(u - uom - J... (43b)
V-(EOE) = - V-((e - €o)E) + 9. (43¢)
V-(uoH) = — V-((u - item) + 9... (43d)

When the permittivity and permeability of a medium are different from that of

free space, equivalent sources can be defined:

J?! = J". + jaw - #0)}! (4.4a)
J:" = J. +jw(e - so»: (4.4b)

5" = — V-((e - so)E) + p. (4.4c)
pf: = - V-((u - item) + pm (4.4d)

In the real-source free region where Jm, J,, p,, pm are zeros, the equivalent sources

become:
J3? = jwm - now = 0...}! (4.5a)
Ji" =1w<e - so»: = oE (4.5b)
i" = - V-((e - €0)E) (4.5c)
pi? = - V-((u - llo)H) (4.5a)

The equivalent sources are depicted in Fig. 4.1, and they obey the continuity

equations of

V-Jj": — jwpgq (4.6a)

However, It can be shown that the total equivalent charges per unit surface area are
zero in the thin film in the limiting case. As shown in Fig. 4.2, the equivalent charges

on two surfaces of the coated film are:

p, = fi-(e — eon), (4.7a)
92 = - fi-(e - eo)E2 (4.7b)

135

 

 

/
chc= j(0(8—80)E
\
(“0960) ++ +F eq (“0980)

++ + _ p mz‘V'((“‘“o)H)

++ + —

eq :_ . ..

++ + _ p e V((8 80)E)
(“0:80)

Fig. 4.1 The illustration of equivalent currents and equivalent charges

(“'0’ 60)

Am

92=n:(e-eo )E2

++++

 

 

Fig. 4.2 Total equivalent charges on a thin layer

136.

When the film thickness approaches zero, the total equivalent charge per surface area

diminishes to zero since E2 is the same on two surfaces of the film.

p‘°‘“’ = 91 + P2 = 0 (47¢)
Following the same argument, it can be shown that the total equivalent magnetic

charge per surface area diminishes too.

pi?“ = O (4.7d)
As shown in Fig. 4.3 the total fields inside region 1 can be considered to be

summation of the incident fields and the scattered fields maintained by the equivalent
sources in the thin layer and the electric current on the surface of the perfect conduc-
tor. Based on the equivalence principles, the general integral equations for E, H fields

can be established:

E'O‘“’(r) = Ty — jwuorgcp — $.wi + Peiv'mdv

“I
+ T3 — jcouojfiqd) — Jfng'cp + —e;—V'<r>]dv
+ T][ — jwu0(r’ixH)<D + WE)xv'<r> + (fi-E)V’<D] (15’ (4.8a)
S

H'o‘“’(r) = T y — jmeoJfinCD + ngv'cp + p—ug'v'omv‘

«I
+ rfl - jwengfcp + Jj"xV’<I> + fiV'cdeV'
+ T? jmeo(n><E)<r> + (riXH)xV’d> + (fi-H)V’<D] as (4.8b)

1
hr<b,’=e
wee (rr) 4M

 

, R = Ir — r’I, and T is 2 when r is on the surface of the perfect

conductor and is 1 when r is in the space outside the perfect conductor.

The second volume integrals are contributed by the equivalent sources in the
coated magnetic layer, and the third surface integrals are contributed by the equivalent

sources on the surface of the conductor.

137

Jeqm =jto(Ll—llo)H

Je = n x H
+

Scattered Fields

It

pe = ( n - E ) / 80

( Region 1 )

 

Fig. 4.3 Seperation of the total fields into incident fields and scattered
fields maintained by equivalent sources

138,

The first volume integrals, which are contributed by real sources to maintain the

incident fields, is denoted respectively by TE‘ in the electric field integral equation and

THi in the magnetic field integral equation.

“I
e

V’<I>]dV'
so

 

E‘o’“’(r) = TE‘ + T3 — jwuojfiqd) — Jfng'cp +
+ TE - jtou0(rixH)<l> + (n3<E»<V’¢ + (ri'E)V’<D] dS’
. «1
H'O‘a’(r) = TH‘ + T1} - jmeorgfcb + £ka + fivmav
+ T? jweo(rb<E)<l> + (r’zXH)><V’<D + (ii-H)V’<I>] as

On the surface of a perfectly conducting body

fixE = 0
fi-H = O
and
.13" = (SE = O

(4.9a)

(4.9b)

(4.10a)
(4. 10b)

(4.11)

By assuming the null electric conductivity, 0 = 0 . Note that the tangential component

of electric field is not continuous across the magnetic current sheet due to

fix(E1 - E7) = —Jf,fi’t . Consequently Jfi" = 0 is not deduced from fixE = O on the surface

of the perfect conductor. In fact, ME is not zero in the thin layer when a magnetic

current sheet exists . It is also interesting to note that the nullity of the currents in an

electrically lossy layer is provided by the null tangential component of electric field on

the surface of the perfect conductor and its continuation across the lossy layer. How-

ever, when a magnetic sheet exists electric current Jfi" = CE in the thin layer may con-

tribute to the scattered fields and then radar cross section if c is not equal to zero. For

the time being, only the case of o = 0 is considered. It is also easy to justify that

fi"=0, t—~>O

pf:=0, t—ao

(4. 12a)
(4.12b)

139.

where t is the thickness of the coated layer. Using above conditions, equations 4.9 are

simplified to:

E‘0'°’(r) = 713‘ + Ty - me'mdv (4.13a)

+ T3“ — jmu0(n3<H)<l> + (fi-E)V’<I>] dS’

H‘O‘“(r) = TH‘ + Ty — jweomowv (4.13b)
+ T3 (mme'cm dS’

Using r’i-H=0 on the surface of a perfect conductor , and shrinking the volume

integrals into surface integrals, when t is very thin, leads to:

Emma) = TE“ + T? — (omt)HxV’<l> — jmu0(rb<H)<D + (ii-E)V’<I>] dS’ (4.14a)
H‘°‘°’(r) = TH‘ + T? — jmeoothCD + (WIDXV’CIH dS’ (4.14b)
where (4.5.a) is used and
EW’: E’ + E" (4.15a)
H‘O‘a’ = H’ + H" (4.15b)
If the currents are related to the tangential components of the fields and proper

boundary conditions are applied, both the electric field integral equation and the mag-

netic field integral equation can be established.

The tangential component of total E field is zero on the surface of perfect con-
ductor, that is fixE‘m’ = 0 . However, matching the boundary condition on the surface
of perfect conductor will result in an incorrect integral equation. E‘“"’(r) on the left
side of (4.14a) should be understood as the total field on the outer surface of the mag-
netic thin layer, which will be illustrated in Appendix C. If the observing point is on
the outer surface of the thin layer, boundary condition is n3<E‘°‘°’ = — .17.? t . Matching
the boundary condition on the outer surface of the thin layer and setting T: 2 result

in:

140

330:: n... = Biol... “‘13 - (omz)HxV’<I> — jwuowm‘” + (“Wm (15,}

tan

or

tan

—E‘(r)lmn=—-;- fixothlm+{£[-(omt)HxV’<b—ja)uo(fixH)<b+(fi-E)V’<I>]dS’} (4.16)

By using fi-E= fiV-J , an integral equation in terms of the unknown of rb<H is

obtained.

-E‘(r)lm=—% fixothImo{£[-(omt)HxV’<b—jwuo(rixH)QH—E.% V’-(rb<H)V’<l>]dS’} (4.17)

tan

This integral equation is the electric field integral equation (EFIE) for arbitrarily
shaped three—dimensional bodies coated with thin magnetic layers, and it is the second

kind of Fredhelm integral equation.

Similarly, the magnetic field integral equation can be derived. The tangential
component of the magnetic field, which is continuous across the magnetic sheet when
0 is zero, is related to the electric current on the surface of a perfect conductor as

fixH = J, .
Using the boundary condition:

I‘lXH = J, (4.183)

and

fill = 0 (4.18b)
yields a magnetic field integral equation (MFIE)

&H‘m’afim = nil” + {g - jo)e0(omt)H<I> + (rb<H)xV’¢] dS'} (4-19)

tan

where H is the unknown to be solved. It is observed that MFIE does not have the
term of V’-(I’D<H)V’<I> which involves the divergence of fixH . Thus, the MFIE has an

advantage of numerical simplicity over the EFIE.

141.

Both the MFIE and the EFIE are applicable to the perfectly conducting bodies
coated with a thin magnetic lossy layer. It is expected that magnetic lossy coating can
significantly change the scattered fields and the radar cross section of a conducting
body. This is verified when the developed integral equations are applied to two-

dimensional problems.
4.2.2 Integral Equations for Two-Dimensional Geometries

Equations (4.17) and (4.19) can be simplified when they are applied to two
dimensional cases. Only the EFIE is taken as an example to show the derivation pro-
cedure. In a two-dimensional problem, the tangential component of H is related to
the electric current by rb<H = I and the electric current I(r) is assumed to have a propa-

gation constant B in the z direction as

mm) = I(r) = Kane-113z (4.20a)
E‘(r) = vane-13* (4.20b)
H‘(r) = H’(B)e”B’ (4.20c)
where
fi=fi+fi

Substituting the above two relations into the EFIE and integrating over the 2 variable

yield:
— E‘(fi)e”lem = — gammmpw/B’Im + 4 fl - (amoncme', j e-mz' <D(r,r’) dz’
c —oo
A” -sz’ a r - a -}B2' I
— (omt)H(p’)xz j e 574%”) dz - )wany) j e <D(r,r’) dz
+ $6”, — jag-rcmvgie'lwmrx') dz’

.9.

_L ' _ - . .. 432' ' ,
+ 8000(V, 1B2“) Jqfifie a2,430.1") dz 1 6” hm (4°21)

142.

It is well known that the two-dimensional Green’s function is:

G = $118140? - Barr? - m)

 

 

: Te‘jfia’ " 2) exp—jBJ'fi - 5'2 + (Z, _ Z)2 d2, (4 22)
- War—mu (2 - z’)2 '
and
1° flax — 2) ix exam/fr? — F" + (2' - z)2 W
—~ 32 \hp’ - m2 + (z — {)2
= ,we - z) (exp-1W1? — 6"2 + (z’ - 2f ,1
‘1'? - 3'2 + (z - 2’)2
_ fire — ., wax exg-jfi‘lrfi - m2 + (2' - 29 W
- Vrfi—m2+<z—z)2
: jB I 6'71“"-I - 2) (”(22/13er _ pflz + (Z, - Z)2 )dzl
«w ire—69+ <2- .32
= jBI-erSWIc’ - [32rd - WI) (4.23)

The integrals over 2 variable in the integral equation can be represented by the two-

 

 

 

 

 

 

 

 

dimensional Green’s function:

_ 4nE‘Cm'tan = - ZKGmIWHCm'm + {ll - (Omtmfmxw': + 113230040c2 - [52W - T3")

- jwuonxnamdwz - Barri - m)
+ $ng — 1132)me + Jfiaccmra - rm] ctr} (4.24)

tan

where H0?) is the unknown, but it can be related to the electric current I by H = -rb<I
and n3<H = I since rb<fixH = - H + fi(fi~H) and fi-H = 0. Thus the EFIE is established

in terms of the electric surface current:

- 4nE‘Cfi'fim = + 21t0mtI©|m + {it (amt)rb<ICfi')><(V’r + J'l3z‘)G(‘J(k2 - [3%l - E")

— jwuoltmcoftkl - 52w - rm
+ 25:67: — jBfl-I@(V’: + 1132300108 - [52W — fi' 01 dl’} (425)

tan

143

which is an electric field integral equation for a two-dimensional arbitrarily shaped
body. Compared to the pure perfect conductor case, it is seen that there is an extra
term contributed by the magnetic current. The contribution from the magnetic current
will give rise to a different scattering property from that of a perfect conductor without
coafing.

If TM polarized fields of E = z‘Ez , H = Hg? 4» Hy)? illuminate a circular cylinder,
the integral equation has only the 2 component

I = fixH = z“! (4.26a)

Note that

3.272 = 0 (4.261))

In the cylindrical coordinates.

mama - m) = (é—é— + OgiflGM/cz - B’W - 15") (427a)

 

pae'
with
(5 £00108 - W6 — F7") = :50th - Barri - [3") (k2 — [35335.46 - 3' (4271’)

 

since rp - p'l = R = \[p2 + p,2 -— 2pp’cos(9 — 9'), we have:

3376 — fi'l = 2(p’ - pcos(6 — 9’))/R (4.28)
Equation (4.25) can be decomposed into scalar components in cylindrical coordinates
for a cylinder of arbitrary cross seetion by calculating each term as above in the
integral equation. In particular if a circular cylinder is considered, then p = p’ and

above expression can be simplified:

 

R = p‘ll - cos(9 - OWE (4.2%)

and

 

3:76 — T3" = m1 — cos(9 - 9'» (4.2%)

144

In addition
G'o/(Ic2 — 52w - 17') = —:$-H6’>'N(k2 — Bard — m)
= - -’jF-H?>< (k2 - leffi - m) (430)
Thus for a circular cylinder, the first integrand in (4.25) is reduced to

(amomcmxtvz + jBaGo’tkz — B’m? - 3|)
= .@<m<zs<o)§;0d<k2 — BZW — WI)

 

= — 1,(m2—£H<,2)(«/(k2 — B’W - rm ‘1th — B2) ((1 - cos<e — e'wi (4.31)
J
The simplified electric field integral equation in a cylindrical coordinate is obtained:

— 4nE‘ca) = + 2nomzlca) (4.32)

2::

+ l [ j<omt>H32>’(\/<k2 - Barr? - WW - B’Ml — cos<e - e'»N§

 

— wuoImHBzW/cz - WW - (5'01 p’de’}

It is equivalent to

— 4E‘(p) = + Zomthi')

 

21:
+ {I - Elfin—O'Hiz’TW/Cz - [32W - 3|)sz - BZMI — c08(9 - 9W5

- WWWEEWWW - (7')] r'de'} (4330

Generally speaking, equation (4.25) can be written either in a rectangular coordi-

nate or a cylindrical coordinate. In some cases, its representation in a rectangular

coordinates is easer to use. For example, the integral equation of (4.25) applied to a
rectangular cylinder with a TM polarized illumination can be written as:

E: 0'! [‘0 font 2 {—3. .
— = — '— _H - - I - . ._ .
20 2201 + 4 £{Zo'p- DVI 1( 0‘2 B )Ifi ff )Slgn(x ’9’an IO]

145.

+ [1 — (%)2]H3(V(k2 — B2)? — my}! (4.34)

where x,’=1’ and x,=x ifp' is on the side Iyl=a , andx,’=y’ and xi=y ifp" is on the
side Ix! = b as depicted in Fig. 4.4. sign(x—x,-’) is a function which is +1 if (x—x,’)=>0 or
—1 if (x—x,’)=<0 . Calculation of the diagonal elements of the moment matrix for (4.33)

and (4.34) needs special care for singular integrals.

In the case of TE polarization, due to the divergence of unknown 1 , the EFIE
becomes difficult to solve numerically . An alternative approach is to use MIFE for

treating the TB polarization. The proper MFIE can be summarized as follows:
2nH®lm = 41rH‘(fr’)Ium + {ll — jtr)eo(o,,,r)H(‘(3‘)G(‘J(k2 - [32ml - fil)

+ WHCWXW': +J'132’)G(V(k2 - [32W - 30] d5’} (4.35)

tan
In the case of TE polarization, HQ?) = EHZCH) and EC?) =£E,+ y‘Ey . The MFIE

with TE excitation reduces to an integral equation with only 2 component:

2nz‘H’Cm = 41:2ng + {gt — jweommwmcdu? — we — T3")?

+ Hz®m®<<vz>cd<k2 — we - r701 49'} (4.36)

The second integrand in (4.36) is similar to the first one in (4.25), which can be

conducted in a cylindrical coordinate as:

memxwzmdaz - 62w — 17 I)

= _ H,(p)z*-’J-F-H?)’(\J(k2 - 52m - 3:) W — (32) v'(1 — cos(e - 9'))N§ (4.37)

The MFIE for a circular cylinder with TE excitation then becomes:

 

21%) = 411:0?) +{1l — weotomrwztmfIb’kW/é - Barr

 

- 15" bi + jflzchHizl'Ntkz - Barr? — 15" I) 40:2 - £32) ((1 — cos<e - awn/i] 49’} (4.38)

146

 

 

 

 

 

 

Fig. 4.4 Application of the surface E-field integral equation on the surface
of a rectangular cylinder

147—

and the MFIE for a square cylinder with TB excitation is:

 

C

—T)"I)I + sign(x—x--(i——') 1? JHZCsz-fiv' BOW—301%! (4.39)
which can be solved by the moment method.

When the currents are known, radar cross section can be evaluated by using the

asymptotic expressions of Bessel functions.
Hf?)(2) -> 2/7tze"("”m’m) ( z —> oo ) (4.40)
For the TM case, radar cross seetion of a circular cylinder is represented by

E; 2
0(9) = 21th—.—|
El

2

= 2an éfljmngzkkp) + jomtcos(6 - 8’)HP(kp)]l(8’)d9’l2
C

ka ’JU‘D ' —) 2) 2
= 2 I — — -— + m 9- 8 I 6’ de’ I
Itp 4], Mp e {I no 16,, os< )e M )
—k“—2-Ij[ 1 — -°—"‘-cos(e— 6’)]l(9’)d9 I2 (4.41)
=4llo c 110

and the normalized radar cross section is by

0(9) = ka2
7» 81"]0

 

I? 1 — g"'—tcos(9 - 0’)]l(6')dO’l2 (4.42)
'00

The radar cross section of a rectangular cylinder is given by

 

2 = "02 I i (_Gmta + 1)ef*p'°°‘<°‘ 9°1(x’.y’)d(i)|2 (4.43)
A. 8K 20 a

where a = case is 4505195135", and a = sine if OSIOIS45° or l35°SIthetalSl80° .
For the TE case, the radar cross section of a circular cylinder is
3

0(9) 2 IH’F
= 1: —
p H

= 2w $iI-‘iiflaz’tkm - 008(9 - milizltkpnltmdeiz
J fllo

148

r 404-3) -'(kp——-1)
= 2an — 1‘3- —£[ — 3—"‘-e ’ 4 + cos(9 - 6’)e ’ 4 2 ]l(6’)d9’l2
4] “(‘9 mo
[Caz amt I I I 2
= —| [ — + cos(9 - 9 )]I(9 )dB I (4.44)
4k l‘ 710

and the normalized radar cross section is

6(9) =

-k—aZ-|J[ -—'1- + cos(0 — 9’)]l(8’)dO’|2 (4.45a)
A. 81: c

o I
110
While the radar cross section of a rectangular cylinder is

A (— - a)ej"°'°°“°' “IO/maxi)? (4.45b)
81: a

7—0

o -£a—2ll[ omt

4.2.3 Some Properties of New Surface Integral Equations

The new surface integral equations for treating the scattering problem of a metal-
lic object with thin magnetic coating have been derived in section 4.2 as given in the
eq.(4.17) and eq.(4.19). The integral equations of (4.17) and (4.19) are of simple
forms. Compared to the integral equations for a pure metallic conductor, only one
extra integral, which is contributed by the magnetic thin layer, has to be evaluated.
The Green’s function in (4.17) and (4.19) remains the same as that in the free space
without the complex argument involved. The permeability used has been assumed to
be isotropic, but the integral equation can be easily extended to an anisotropic mag-
netic coating by introducing a dyadic permeability as J,,, = joxli’ - u0)H = 8,3.

Both the EFIE of (4.17) and the MFIE of (4.19) are the second kind of Fredhelm

integral equations which are normally diagonal dominant.

The strategies used in the derivation of (4.17) and (4.19) can also be applied to
solve a homogeneous lossy body with either thin magnetic or electric coating. This

will result in a set of coupled surface integral equations.

1 49.

4.3 Eigen-mode Expansion Solution to an Infinitely Long Circular Cylinder

It is well known that an exact solution for an infinitely long conducting circular
cylinder can be pursued by the eigen-mode expansion method. However, when a per-
fectly conducting circular cylinder is coated with a very thin layer of magnetic lossy
material, the eigen mode expansion solution needs to be modified. In this section, a
new approach with the eigen-mode expansion is developed for a circular cylinder with
a thin coating based on Taylor series expansion of Bessel functions. The electrically

large cylinders can be efficiently treated via this method.

The general solution to the wave equation in the cylindrical coordinate system
can be expressed by the sum of eigenmodes, which are cylindrical harmonic functions,
in different homogeneous regions. The coefficients in the expansion can be deter-
mined by matching the boundary conditions. As shown in Fig. 4.5, the space is

divided into a few regions.

First the TB polarization is considered. In free space,

H; = Egg—imam) cos (ne) (4.46a)
5', = fizcflmfilxkop) cos (n9) (4.46b)
H“; = zaflfiiz)"(k0p) cos (n6) (4.46c)
_ 17‘0 2»
5‘6 - —Eanfll. (kop) cos (n9) (446d)
0080 n
where

.4 m

Inside the film,
H; = Z[b,Hf,l)(kp) + walk/cm] cos (n9) (4.48a)

£5 = fiztbfiytm + ammo» cos (n9) (4.481»

Ei

150

Free Space

 

III.

k Perfect

 

\

 

 

Fig. 4.5 A circular cylinder coated with a magnetically lossy layer

151.

At surface of a perfect conductor p = a , the tangential component of electric field
f,(p = a) = 0 . Using this condition and the orthogonal property of the sinusoidal func-

tions yields:

 

 

HEM/ca)
b = — — 4.49
Substituting it back into Eq.(48) for the b, leads to:
Hm k
H; = m .. £14k: Hf,”(kp) + H§2)(kp)l cos (n9) (4.50a)
jko HE.’>'<ka)
35 = (Tag: c. I - Hf,1>’(ka)H"2)l(kp) + 119%ka cos (n9) (4.50b)

At the interface between the film and the free space, the tangential components of
fields are related by

Eé(a + £1) + E30: + F) = E§(a + t‘) (4.51a)

H§(a + t+) + H§(a + t+) = H§(a + t‘) (4.51b)

Note that the thickness t is very thin, Bessel functions can be approximated by

their first order Taylor expansions:

HS,‘)(a + z) = Hf,”(a) + 119%) t (4.52a)
115,2)(a + t) = Hf,2)(a) + Hfiz)’(a) t (4.52b)
119% + z) = HS,‘>'(a) + HS,””(a) t (4.52c)
Hf,2>’(a + t) = llama) + Hf)"(a) t (4.52d)
1,,(a + t) = 1,,(a) + Jn’(a) t (4.526)
J,,'(a + t) = ma) + J,,"(a) t (4.520

Using the approximations of Bessel’s functions and matching the boundary condi-

tion of O-component electric field lead to:

ko _ ko
—— "1,: +1,” k t +— ,HS,” +HS,‘>"k
Cal [ (koa) (oa)kol at (koa) (oa)kot]

MEN/ca)
“i - Haw/co
k2: HS.2*<ka)H§.‘*'<ka) - max/confine)
we "[ MSW/ca) ]

[HM/ca) + ”Mm k2] + HSEWka) + H‘.”'(ka) kn}

 

(4.53)

152

It can be shown by using Wronskin’s relations that

H9)’(ka)Hf,‘)”(ka) .. HS,”’(ka)Hf,2)”(ka) = in — <-”->21 (4.54)
rtka ka

Thus

_"9_ _.. ' ~ _"_0__ I» I)»
(0804.; [1,. (koa) +J. (koa) Icon + moms. (koa>+HS. thou) Icon

 

=__4&_ __ _n_2 1
“mon t [can Hf,‘)’(ka) (4.55)

As t approaches zero, k0: approaches zero, but quu — no): remains a constant.

Equation (4.55) also holds in the limit of t —) O :

417i [1 — (f5)? Hm
— Mme Cn ”9),(ka) = “Olga! Jr; “(00) + an n (koaH (456)

 

If A is defined as:

 

. [1 - <—” >12
4} ”rt ka
nave, Hf,”’(ka) ( )
Then the coefficient a,, is given by:
A ,, - ""J,’ k
n = C CH.) ( 061) (4.58)

 

Hf,2)’(k0a)

Following the same arguments and Using the boundary condition of the tangential
magnetic field result in another equation, which can be used to determine the

coefficients a, and c, in (4.53) and (4.54)

HS?" ka .
[Ca—"Mm + aJIffikkoan = c.[ - H(1)’:ka:H9)(k0) + Ira-lawn

 

_ 41c. _
_ nkaHf.‘)’(ka) -

 

86,. (4.59)

where B = .. From (4.59), the coefficient on is given

__4J__
nkaHS."’(ka)

c. = 'fi'lCJ—"JJ/(oa) + amt/coon (4.60)

153—

Finally the solution to the coefficient a, is obtained as

[ko’llr/nU‘oa) + Jn’(koa)]C,./”"
Hf?”(koa) + kotuer.”(koa)

Subsequently the magnetic field Hz on the surface of the film is represented in

 

(4.61)

n-*

terms of the coefficients { a, ] as

H,(a + F) = H;(a + 2+) + H:(a + F)
Romp/Alma) + J.’(koa)]C,J"'
HSY<koa> + ko‘lerizz)(koa)

= -" '3 1 e 4 62
$4 "koa H£2>'<koa> + WEEK/coo) COS“ ) ( ' )

where the Wronskin’s relation is used. The surface electric current is represented by:

 

= — ZCJ’"[J.(koa) - H9)<koa>1cos<ne>

 

19 = 'Hz

= __ ..,, ~12 l
33”“ Wt H:Z>'<koa) + kotllrl'lftz)(koa)

The radar cross section can be evaluated by using the asymptotic expressions of

 

cos(n8) (4.63)

Bessel functions when the arguments are large.

Hz 2
O = anlTi-‘rl (4.64)

I"!

= 3k;123a,,ej-2—cosme):2

fl = -2-I)Eauf'cosmen2 (4.65)
A 1: ":0

The same procedure can be followed to find the solution to the case of TM polar-

ization. In the case of TM polarization,

E=z‘E,, H=éH9+pHp

In free space:

E = noZCnI”"J.(/cop) cos (n9) (4.6621)

154

112', = — jXCf’UnTkop) cos (n9) (4.66b)

E; = nozaflfflkom cos (no) (4.66c)
3 = - Jxaflzmop) cos (n9) (4.66d)
Inside the film,
3; = n2[b,,H§,2)(kp) + c,H§,‘>(/cp)1 cos (n9) (4.67a)
Hz; = — 1212.4153ka + cwa‘rtkon cos (n9) (4.67b)

By approximating Bessel funcrions with their first order Taylor expansions and

matching the boundary conditions, a, , b, , and c,, are determined as

 

 

b Hf,‘)(ka) (4 68)
= — —c,I .
" Hfizkka)
—l——4' = ""1 ’ k W 4.69
Cu TtkaHfizkka) Cm] n( on) + 61an (koa) ( )
with
... komrln’(koa) - Jn(k0a)
a,, = " , (4.70)
C” HEM/coo) - komflffwkoa)
The surface current is determined by
lz=H9(p = a) (471)
. .. kOturJn’UcOa) — Jn(k0a) 2)’ r
= — 12;) "r , Hf. (koa) + J.(koa>1cos(ne)
{3 Hall/coo) — 1441.142) (koa)
and the radar cross section is by
2 = A -n-l 2
A 1: I?” cos(n6)l (4.72)

The numerical results based on this method will be shown in Section 4.5 to serve

as a comparison with the results by the MoM or the FD-TD method.

4.4 Finite Difference Time Domain Method

In the proceeding sections, the integral equations and the eigen mode expansion

have been used to solve the scattering problems of a perfectly conducting cylinder

155

coated with a thin film. When the body is electrically large, the integral equation tech-
nique becomes inefficient due to the matrix inversion required by the moment method.
In this section, the FD-TD method is introduced to solve the same problem and to pro-
vide an alternative approach to validate the integral equations derived. As an example
of the FD—TD method, only two dimensional metallic objects coated with layers of thin

lossy magnetic materials are studied in this section.

The basic used algorithm of the FD-TD method is developed from that discussed
in Chapter 3. As shown in Fig. 4.6, an infinite two dimensional space is truncated into
a finite region, and Mur’s second order radiation boundary condition is applied to the
truncated surfaces. The truncated region is then divided into the scattered field region
and the total field region to ease the application of the radiation boundary condition in
the presence of an incident wave. Two dimensional Yee’s model is used for interior

points.

The thin magnetic coating is taken into account by using the integral form of
Maxwell’s curl equations, which leads to a modification of Yee’s difference scheme.
It will be discussed in detail for the cases of TM polarization and TE polarization

respectively.

The Farady’s law can be represented by an integral form as:

ind? = _ (cognac - I (o - u0)§t-H-d3 _ l 3,43 (4.73)

where J", is magnetic conducting current. Hereafter the magnetic current is related to
the magnetic field by J," = omH as the electric current is related to electric field by
J, = 0E . In the time domain, 6,, represents the magnetic loss of magnetic material
which is equivalent to the imaginary part of complex permeability (upui) in the fre-
quency domain by a relation of o,,, = (nu,- when the excitation wave is harmonic in

time.

156

Y
(r

Radiation Boundary Condition

oooooooooooooooooooooooooooooooooooooo

 

 

 

Fig. 4.6 Truncation of an infinite two-dimensional space into a finite
region

">5:

157.

Shown in Fig. 4.7 is one cell of Yee’s model in the case of TM polarization near
the coating of a square cylinder. To update the x component of magnetic field H, , for
example, Farady’s law is applied. According to the model in Fig.4.7, the left hand
side of Farady’s law in (4.73) can be approximated as:

1’ Ed? = (Ego; j+1) — £30.13] Az (4.74a)
C

The third term on the right side in (4.73) is discretized as

n+ 1 n- I
j Jm-dK = 0.5(o,,,z)Az[H, 7(1, 14%) + H, 7(1, 14%)] (4.74b)
S

where t is the thickness of the thin layer. The first and the second terms on the right

hand side of (4.73) are approximated as:

8 . - _3.. _
juognfi+j<u 1618:1143- (4.74c)

n+21- 1 n-zl- 1
uOAsz/At[H, (i, j+—2-) — H, (i. 143)]

"+11. . 1 1e}. . 1
By rearranging (473-474), we can obtain a modified Yee’s difference scheme

suitable for the 2D scattering problem with TM polarization.

"+1 1 CA -’ 1 Com
Hf.’._=__H2'.,.__ E"',' —E"’°, .
x (1 1+2) CB . (1 1+2) —ZOCB[ ,(1 1+1) ,(1 D] (4 75)
where CA and CB are defined as
CA = 0.50,,tc0At/zo - (u, — 1)t — Ay (4.76a)
C3 = 0.50mtcoAt/zo + (u, — Dr «1» Ay (4.7 6b)

This modified formulation is for updating H, in the cells adjacent to the square
cylinder walls which cross the y-axis at y = -nldy and y = nldy . Similarly, a formula-
tion for updating Hy in the cells adjacent to the walls which cross the x-axis at

x = —n1dx and x = nldx can also be derived by using the Farady’s law.

158

 

.................................

 

 

ll

 

 

 

 

 

 

 

E2 E2

 

Fig. 4.7 A basic cell of Yee’s model adjacent to the coating
in the case of the TM excitation

159

Next, consider the case of TE polarization, and focus on the cell as shown in Fig.

4.8. To update the 2 component of magnetic field, Farady’s law is used again. On left
hand side of the Farady’s law, the contour integral yields:

l[rt-:71 = - 115:6. 14%) - 520.1%)14: + [E;(i+%.1) — Ext—g. moy (477a)

On the right hand side of the Farady’s law, the flux integral of the magnetic con-
ducting current and the magnetic polarization current are those included by the contour

integral path. From the third integral of (4.73):

11+ 1 n_ 1
j Jm-dis’ = 0.5mm, fa, j) + H2 2r(i, 1)]pr (4.77b)

where B is equal to 0.5 if Hz is on the corners of square cylinder, otherwise B = 1. The

second integral of (4.73) can be written as

..1 ..1
j (11 — uo>-§;H«d‘$ = (11 — 1101441111. : (i. 1) - Hz 7031311341: (4.77c)

where [3 is the same as defined in (4.77b). The first term of (73) is written as

 

M 1 ”_1
I 1103371143 = quxAy/MH. I (1. J) — H. 76'. 1)] (MM)
By adding (4.77) together and rearranging terms, a formulation for updating Hz
becomes:
n+2l- CA n-zl-
H .9 = - __.”; .9
z (1 1) CB (1 J)
+ COM [15"(1', j+-1—) — 5"(1', j—inAx - [E"(i+i. D - E"(i—i. mm (478)
zOCB ‘ 2 x 2 y 2 y 2
with
CA = 0.56,,tc0At/zo — (u, - 1): - Ay (4.79a)
. C8 = 0.50,,tcoAt/zo + (u, — 1): + Ay (4.7%)

The same strategies can be followed to treat the cylinders with an arbitrary cross

section. The irregular boundaries can also be modeled by using Farady’s law and

160

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Y 11
2 a
1523;: ..... “ ...... _“....-.. r:;,;,;,;,;_;,;;,_-;;. ;....: ............. A
2 b
; X
Ex
............ v
E E
l y . Hz 4) y
TEX

 

Fig. 4.8 A basic cell of Yee’s model adjacent to the coating of a
retangular cylinder in the case of the TB excitation

161.

Ampere’s law to conform the integral paths.

Numerical results have been obtained by using the schemes discussed in this sec-
tion for the TM or the TB case. They will be compared with those generated by the

integral equation method and the eigen mode expansion method in the next section.

4.5 Numerical Results

Now we have three different methods for solving the scattering problems involv-
ing a metallic object coated with thin layers of magnetic materials: the new surface
integral equation has been developed for an arbitrarily shaped metallic object coated
with a thin layer of magnetic material; the solution with eigen mode expansion for a
circular cylinder has also derived for a circular cylinder case; a modified scheme of
the FD-TD method has been developed for a cylinder with an arbitrarily shaped cross
section. Based on these three indepedent methods, numerical results are obtained,
respectively. To validate the new surface integral equations, the current distribution
and the radar cross section of a circular cylinder are calculated and compared with that
from the eigen mode expansion method. To study the singularity behavior of sharp
corners, square cylinders are chosen and the results from the integral equation tech-
nique and the FD-TD method are compared. The excellent agreement is obtained.
The effects of both complete coating and partial coating on radar cross section are

investigated for a circular cylinder and a square cylinder.

4.5.1 Perfectly Conducting Circular Cylinder Coated with Magnetic Thin Layer

a. Completely Coated Cylinder

Consider an infinitely long, perfectly conducting circular cylinder with koa = 51:
where a is the radius of the circular cylinder. The perfectly conducting cylinder is

completely coated with a magnetic thin layer which has a parameter of

162

u,:/a= 0.01 — j0.03. A plane wave is incident on the cylinder along the direction of
increasing x. For the TE excitation, the current distributions and the radar cross sec-
tions by the MoM and the eigen mode expansion method are plotted in Fig. 4.9-4.10.
As shown in Fig. 4.9-4.10, the results by the MoM based on the new surface integral
equations for both the current distribution and the radar cross section have excellent
agreement with that of the exact solutions of eigenmode expansion. The comparison
of the current and the radar cross section between the cases of coated cylinder and
non-coated cylinder is also plotted in Fig. 4.9-4.10, and significant reduction in the
radar cross section is obtained if a magnetic thin layer is coated on the cylinder. The

back scattered filed is reduced by more than 10 db.

Figures 4.11 and 4.12 show the currents and the radar cross section of a coated
and uncoated cylinder with koa = 21: under a TM excitation. Both results from the
MoM and the exact solution also have an excellent agreement. The back scattered
filed is reduced by 3 db. If the size of the cylinder is increased to koa = 51: , the back
scattered field is decreased by 7 db as shown in Fig. 4.13-4.14.

Comparing the results with the TM and the TE excitation as shown in Fig. 4.9-
4.14, we see that magnetic coating on the reduction of radar cross section is more

effective for the TB excitation than for the TM excitation.
b. Partially Coated Cylinder

A perfectly conducting circular cylinder of koa = 21: is partially coated with a thin
layer of u,t/a = 0.01 — j0.03 which covers 25% of the circumference within

180° - 4505951800 + 45° .

With the 'IM excitation, the radar cross section and the current on the fully
coated, partially coated or bare cylinder are plotted in Figs. 4.15—4.16. The current

distribution of partially coated cylinder exhibits a singular behavior at the edges of the

163

coating. On the coated surface the current is very close to that of fully coated cylinder
except for at edges. On the remaining portion without coating, the current is almost
the same as that of the uncoated cylinder. The radar cross section of a partially
coated cylinder is some value between that of fully coated and bare cylinders except

for around 180° degree.

Figures 4.17 and 4.18 show the currents and the radar cross sections of a cylinder
with koa = 2n under the TE excitation for fully coated, partially coated and uncoated
cylinders. The current distribution in the TE case does not have the strong singularity

as observed in the TM case.
4.5.2 Perfectly Conducting Square Cylinder Coated with Magnetic Thin Layer

Next we examine an infinitely long square cylinder with koa = 21: where a is the
length of one side. The plane wave is propagating along the x-axis at 0° degree and
the parameter of magnetic coating is u, = t/a = 0.01 - j0.03 . The curves in Figs. 4.19,
4.20 and 4.21 are the numerical results obtained by the MoM, based on the integral
equations, and the FD-TD method.

Figure 4.19 shows current distributions of a perfectly conducting square cylinder
without coating based on the MoM and the FD-TD method. Fig. 4.20 exhibits current
distributions on the fully coated square cylinder. The current distributions on the par-
tially coated square cylinder are also plotted in Fig. 4.21. The magnetic films with one
half of side length are symmetrically coated on the two comers at 1350 and 225°
degrees. The numerical results have shown good consistence between the MoM and

the FD—TD method, even including singularity behaviors.

As the last example, consider a plane wave propagating along the direction of 45°
degree. Figs. 4.22-4.25 depict the variations of current distributions on the square

cylinder for the TM and the TB cases. Figs. 4.23 and 4.25 are the radar cross sections

164

associated with the current distributions in Figs. 4.22 and 4.24.

4.6 Extension to Three Dimensional Case

The integral equation method and the FD-TD method used in this chapter can be
extended to the three dimensional problems. To use the integral equation method for a
3D problem, numerical algorithm developed in Chapter 2 can be applied. The moment
method with vector basis functions would be an appropriate choice. The eigen mode
expansion method can be used to solve the problem of a perfectly conducting sphere
with thin magnetic coating. The FD-TD scheme for 3D problems can be developd in

a similar procedure as described in this chapter.

165

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

° Exact solu. f
18" nuance MOM $011.1.
1 ...... No coating 1
1 ii
12 —
5-
Cl}. 0..
(dB)
—6-
-12 -1
-13 _
o ’ °
-24 - 1 1 1 1 1
0 60 120 180 240 300 360

6(indegree)

Figure 4.9 Radar cross sections of an infinitely long conducting circular cylinder
with or without a magnetic coating in the case of TE excitation (
koa = 51:, u,t/a = 0.01 - j0.03 )

166

 

2_4

lJeI/lHil 1.2—

0.8 —

      
   
      

—— Exact solu.
...... MoM solu.
...... N o coating

0.4 -

 

 

 

  

- - l l 1 °
0 60 120 180 240 300 360

9 ( in degree)

Figure 4.10 Amplitude distribution of the G-component current on the surface of an
infinitely long conducting circular cylinder with or without a magnetic coating in the
case of TE excitation ( koa = 51:, tut/a = 0.01 —- 1 0.03 )

167

 

 

 

 

 

   

 

 

 

 

14 _-
12 -

— Exact solu.

10 -— ........ MoM solu.

- - - — No coating
all). 3 -‘

(dB)
5 .4
4 _.
2 _
o ._
1 1 1 I I
0 60 120 180 240 300 360
9( in degree )

Figure 4.11 Radar cross sections of an infinitely long conducting circular cylinder

with or without a magnetic coating in the case of TM excitation (
koa=21t,tt,tla=0.01-j0.03)

168

 

 

 

 

/ \
1.7 — / \
/ \
o I \
/ \
1.4 —- ’ ‘
l \
/ \
11,1/1H11 ’ ‘
1.1 — / \
/ \
/ \
l \
0.8 — l \
/ \

0.5 -«

— Exact solu.
0.2 __ ........ MoM solu.

— - - - No coating

~ I L l l L
0 60 120 180 240 300 360
9 ( in degree )

Figure 4.12 Amplitude distribution of the z-component current on the surface of an
infinitely long conducting circular cylinder with or without a magnetic coating in the
case of TM excitation ( koa = 2n, u,t/a = 0.01 - j 0.03 )

169

 

—— Exact solu.
l8 — Y -
...... lNo coating
14 _
o/k
(dB)
10 —
5 _

 

 

 

 

 

 

 

 

 

 

 

0 60 120 180 240 300 360
9( in degree )

Figure 4.13 Radar cross sections of an infinitely long conducting circular cylinder
with or without a magnetic coating in the case of TM excitation (
koa = 51:, 11.,t/a = 0.01- 10.03 )

170

 

   
   
    

  

 

 

 

2 — IA).
1.6 -,
IJ, l/lHi I 1.2-
0.8 —
0.4 ‘ Exact solu.
...... No coating
0 - - 1 -
0 60 120 180 240 300 360

9 ( in degree )

Figure 4.14 Amplitude distribution of the z-component current on the surface of an
infinitely long conducting circular cylinder with or without a magnetic coating in the
case of TM excitation ( koa = 511:, tut/a = 0.01 — j0.03 )

171

 

 

 

 

 

 

 

 

 

 

 

14 .-
12 — _
__. Fully coated
10 __ ........ Partially coated
- - - .. No coating
o/71. 3"
(dB)
6 _
4 _.
2 __.
o _.
l 1 l I l'
o 60 120 180 240 300 360
6( in degree )

Figure 4.15 Radar cross sections of an infinitely long conducting circular cylinder
partially coated with a magnetically lossy thin layer in the case of TM excitation (
koa = 21:, tut/a = 0.01 - j0.03 )

172

 

 

 

 

/ \
1.7 — -"- ’ \ 2‘

3'2/ \:"-.

37 ‘1':

fl: "'-..-\§

1 E: '-E t
1.4—1 ' ' '

- _1 - \
1 \.
1' .I \.
lJ,l/|H l 1.1.. I \_
l \.
l \
/ \
0.8 — / \
l
0.5 -
—- Fully coated
0.2 _ ........ Partially coated
- .. - - No coating
1 1 ! l
0 120 180 240 300 360
9 (in degree)

Figure 4.16 Amplitude distribution of the z-component current on the surface of an
infinitely long conducting circular cylinder partially coated with a magneticallt lossy
thin layer in the case of 'IM excitation ( koa = 21:. pvt/a = 0.01 — j 0.03 )

173

 

 

 

 

— Fully coated
........ Partially coated
104 - - - - No coating
1 I
9 I
l I
4.. ’1 A f\ /\/V\/\ /\ A "
i l ..1 '5 V .. __ V..‘ {3'1 1‘
1.211. -°
4 . - '-
o/‘A. 1 ,1 311.: 3115 l 1
(dB) -2- l .1 EV: 1. I
1 j : .' | t
, . .. . ,
. I
-‘ 1.
3: :1
1', it
-3... .13 it
1.1 1t
:13 U
3:! [I
I1 ti
II I!
ll i!
I I
-14- : :
L l l l l
0 60 120 180 240 300 360

9(indegree)

Figure 4.17 Radar cross sections of an infinitely long conducting circular cylinder
partially coated with a magnetically lossy thin layer in the case of TE excitation (
Icoa = 21:, u,t/a = 0.01 - j0.03 )

174-

 

 

 

 

 

1.7-—
1.4—
~ 1.1-
119 l/IH‘ 1
0.8—1
0.5 —
\ I I
\ 1' 1 I
1 l _— Fully coated . i
02 .2 '\ ........ Partially coated j
. ---- No coating ’
1 1 1 ! l
0 60 120 180 240 300 360
9 ( in degree)

Figure 4.18 Amplitude distribution of the G-component current on the surface of an
infinitely long conducting circular cylinder partially coated with a magneticallt lossy
thin layer in the case of TE excrtation ( koa = 21:, u,t/a = 0.01 - j0.03 )

175

 

3.4 «—

2.6 -

2.2— ' ‘

 

IJ, I/IH‘I
1.8—-

 

1.4— ’

 

0.6 —
FDTD solution

02 _ —— MoM solution
0 60 120 180 240 300 360

 

 

9(indegree)

Figure 4.19 Amplitude distribution of the z-component current on the surface of an

infinitely long conducting rectangular cylinder in the case of TM excitation (
Icoa = 21:, tut/a = 0.0)

176

 

252-“

1.8—1 1W»

1J1-
q 1

‘

 

 

 

Ile/IH’I

  
 

...... FDTD solution

— MoM solution

 

l 1 1
180 240 300 360

 

 

 

9(in degree)

Figure 4.20 Amplitude distribution of the z-component current on the surface of an
infinitely long conducting rectangular cylinder coated with a magnetically lossy thin

layer in the case of TM excitation ( koa = 21:. 11,2/a = 0.01 - j 0.03 )

177

 

21$—

22%—

 

. 1.8—
IJzI/IH‘I

 

 

L4- ‘

 

 

 

   

 

FDTD solution

 

 

 

 

—— MoM solution
1 1 - l
180 240 300 360
9 ( in degree )

Figure 4.21 Amplitude distribution of the 2 -component current on the surface of an
infinitely long conducting rectangular cylinder partially coated with a magnetically
lossy thin layer in the case of TM excitation (Icon = 211:, u,t/a = 0.01 - j0.03 )

178

 

   

 

 

 

12
— Fully coated
8 _ -------- Partially coated
_ - - - No coating
4 _
/
. . '/_.-'
cl}. o-t/ '-._1 If
. h I .
(dB) ':\ I 5
'. I I .°
'1 \ I
'. \ I .'
5.. /\ 1;
3 V \/ '
-3 _
-12...
1 1 1 1 1
0 60 120 180 240 300 360
9 ( in degree )

Figure 4.22 Radar cross sections of an infinitely long conducting rectangular cylinder
partially coated with a magnetically lossy thin layer in the case of TM excitation (
koa = 21:, u,t/a = 0.01 — j0.03 )

179

 

3.8 -«

3.4 —

2.6— I!

2.2— :'
lJ,I/IH"|
1.8—

1‘-—
3'--

I l

1.4—-

.—-
.0"-
-—-

- --

0.6 — —— Fully coated

........ Partially coated
024 _ ---- Nocoating _
“ ~ - 1 L ! l
0 so 120 180 240 300 360

 

 

 

9(indegree)

Figure 4.23 Amplitude distribution of the z-component current on the surface of an
infinitely long conducting rectangular cylinder partially coated with a magnetically
lossy’thin layer in the case of TM excitation ( koa = 21:, tut/a = 0.01 - j0.03 )”

180

 

10-

 
 
   

__ Fully coated
........ Partially coated

---- Nocoating
5—

601. ‘5‘
(dB)

-10—

-15..

 

 

 

 

 

 

 

 

O 60 120 180 240 300 360
9 ( in degree )

Figure 4.24 Radar cross sections of an infinitely long conducting rectangular cylinder
partially coated with a magnetically lossy thin layer in the case of TE excitation (
koa = 21:, u,t/a = 0.01 — 10.03 )

181

 

 

 

 

 

/\ A
2.3 — I ‘ ’ \
1:"; \ I .-".‘
if '.I\ I 3“,.
2 — - ’2‘ I .5 '-.
g .. \ ,5. .3
E \\/\ Wig: ]
1.7 — , : '
t t
. 1' ‘1 1’ -- -'
Ila l/IH‘ | 1.4—1 . -._\/:' 1
'1
1
1
xix,
—— Fully coated
........ Partially coated
03- ---- Nocoating
! ! ! ! !
0 60 120 180 240 300 360
9 ( in degree )

Figure 4.25 Amplitude distribution of the 0~component current on the surface of an
infinitely long conducting rectangular cylinder partially coated with a magnetically
lossy thin layer in the case of TE excitation ( Icon = 211:, u,t/a = 0.01 — j0.03 )

CHAPTER v

EFFECTS OF AN IMPEDANCE SHEET ON THE
CHARACTERISTIC PROPERTIES OF CAVITY BACKED ANTENNA
BY FINITE DIFFERENCE TIME DOMAIN METHOD

5.1 Introduction

In some military applications, it is desirable to hide an airplane from the detection
of radar systems. To achieve this goal, the scattered electromagnetic fields from the
airplane need to be substantially reduced. Since an antenna on an airplane is an
efficient scatterer, it is necessary to cover the antenna with a lossy layer to reduce its
radar cross section. However, by doing so, the receiving characteristics of the antenna
may be hampered. In this chapter, we will study the effects of an impedance sheet,
covering a cavity backed antenna, on the scattering and receiving characteristics of the

antenna.

The antenna system to be analyzed is a cavity-backed antenna as depicted in Fig.
5.1. Its receiving and scattering characteristics will be investigated by the finite
difference time domain method. As shown in Fig 5.1, an open rectangular cavity is
situated on an infinite ground plane, and an impedance sheet covers the aperture of
the cavity. Intuitively, it is expected that the impedance sheet will attenuate the scat-
tered fields more than it will do to the incident field of the antenna, because the back-
scattered wave crosses through the impedance sheet twice but the incident wave to the
antenna crosses through the sheet only once. Since the interaction of an EM wave
with a cavity backed antenna covered by an impedance sheet involves complicated
phenomena, quantitative. information on the characteristics of this antenna can be

obtained after a complete solution is obtained.

182

183.

E

H
9

a V

 

—_

//////// //

 

 

 

Fig. 5.1 Cavity backed antenna with an impedance sheet on the aperture

184

The integral equation technique is difficult to apply to this problem due to its
complex geometry and the infinite conducting surface involved. Therefore, the finite
difference time domain method is applied due to its simplicity and efficiency in treat-
ing complex problems. Based on the fundamental theory described in Chapter III, the
Yee’s model is modified for treating an infinitely thin impedance sheet. The radiation
boundary condition is changed to adapt the infinite structure. The backseattered field

and the current distribution on the antenna are obtained by the FD-TD method.

5.2 Basic Finite Difference Time Domain Scheme

The basic algorithm used in this chapter is similar to that discussed in Chapter III.
The Yee’s model is used to discretize the Maxwell’s equations for the interior points
as shown in Fig 5.2. The Mur’s second order radiation boundary condition is used to
truncate the infinite space into a finite one. The total field region and the scattered
field region are also defined. However, several modifications have to be made when
the basic algorithms are applied to this specific problem. First, the Yee’s model is
reconstructed from the integral form of Maxwell’s curl equations to facilitate the deal-
ing of the infinite thin sheet. Second, the Mur’s second order radiation boundary con-
dition is applied to the infinite ground plane by using the image theory. And last, to
reserve the applicability of the radiation boundary condition in the scattered field
region, the total field is decomposed into a part generated with the aperture absent and

another part contributed by the equivalent sources on the aperture.

5.3 Integral Interpolation of Maxwell’s Curl Equation

The differential forms of Maxwell’s curl equations are easily discretized into rec-
tangular or cubic cells of finite differences. In practice, a complex geometry or con-

stituents of materials need a better modeling of the curvature of smooth surface or

185

 

 

 

 

 

 

 

 

 

 

 

 

Radiation Boundary

 

Scattered Field Region

 

Total Field Region

Total Field Region / /

 

 

 

 

 

 

 
   

 

/////

Fig. 5.2 Radiation Boundary And Boundary Seperating
Total Field and Scattering Field

186

irregular inhomogeneity than the staircase modeling. The integral interpolation of
Maxwell’s equations provides a simple but efficient method to overcome the shortcom-
ing of the classic finite difference method. By using the integral interpolation, curva-
ture can be handled by conforming the paths of integration and inhomogeneity which
need not to be some units of the cells, such as the case of thin layers, thin wires etc.
In this section, as an important application, integral interpolation is used to treat an

infinitely thin sheet covering a cavity.

A Yee’s model with an impedance sheet is shown in Fig 5.3. The idea of han-
dling the thin sheet is based on the integral forms of Maxwell’s equations which are

stated as follows:

Farady’s Law is

112-117: — {rm-43 - 111%1163‘ (5.3.1a)
Ampere’s Law is
£11.17 = 11,113 + £c%r:d‘s (5.3.1b)

where J,,, = 6,11. J, = oE , C is the contour path of integration and S is the surface
area surrounded by the contour C . When a magne..c *‘ield is to be updated, Farady’s
law is employed. While Ampere’s law is used to update an electric field. As an
example, to update H, , the contour integrations of electric fields E, and E), which sur-

round the H, as shown in Fig.5.3, are conducted.

The integral equations of (5.3.1a) and (5.3.1b) can be equivalently written into the

forms of

tiEJI = -- Pm-d-S - I(u — u0)-aat-H-d3’ — £flo§£fifg (5.3.2a)

£11.17: pas,» I‘E‘WoiE'dg“ loo-5°;E-d3 (5.3.21.1

187

., 11.111.11.111"

 

 

 

 

Fig. 5.3 Yee’s Model Near The Irnpedence Sheet

188

The equivalent sources can be defined as:

1;? = 1,, + (11 — “0’68?“ (5.3.3a)
1:4 = J, + (e — eo)—§?E (5.3.3b)
0r
1:: = omH + (11 — 110%11 (5.3.4a)
Ji" = GB + (e - %)§?E (5.3.46)

The equivalent magnetic current consists of two parts: magnetic conducting
current J,,, = omH comparable to the electric conducting current, and magnetic polariza-
tion current. 0,, represents the magnetic loss of material which is represented by the
imaginary part of complex permeability in the frequency domain. 0,, can be related to

the complex 11* = (u, — 111,-) by on, = (1111,- when the excitation wave is time harmonic.

In the thin film structure, the sheet current manifests if the thickness t approaches
zero but cm: or (0(11 - no): and or or 01(e — co): remain finite. According to the Yee’s
model as shown in Fig 5.3, the integration of (5.3.2) is performed. For simplicity,
we first assume that on, = 0 and u = us so that only the electric current exists on the
sheet. By using (5.3.2b) and integrating along the contour shown in Fig 5.4(a), the
integral equation (5.3.2b) can be interpolated as:

1H" 21-(i+-1-, j+-1—, k) - H"+ Ami, j—l, k)]Az
’ 2 2 ‘ 2 2
2' . 1 . 1 2' . 1 . 1
- [1'1y (Ii—5,], [(+3) - 1'1y (l+—2',j, k--2-)]Ay
= orAy0.5[E§”*’)(i~t-% ,i, k) + Efi")(i+% J, k)]

+ (s — eo>rAy/44£§"*”<1+—;—.i. k) — 121:1»; .1. k)]

+ coAyAz/Aztci"+”<i+§ .1. k) — arm; .1. k)] (5.3.5)

189

 

 

-

A .

7 .

I

I

C

O

O

O

l

I

O

O

O

I

O

O

O

C

I

y.» 6.: ..o . '.‘.' ._. e - u o n. n u e e o . 'u‘a' - . n .' . .'- o - --------- g ________ .3,- -.»_4 ._._. -_-,-_- - - J ,- - o - I I
:-.-:-.-:-.-:-:-.-.-:-..:-:-'-:-:-.=:-:-:-:-:«- .!:-:-.r. ,.::,.:.;.:.';,;.:.,z..:., . .-.-':¢.-.n 1 -.-:’:~ --------- -' '3‘- -.>. ....... ’5 ----- :~ .-.-.~:.-.-.-'=:'a ..... ~. :'. . > <‘:- -' . . . :-: --:~:-'

0

O

O

C

O

I

I

C

O

I

I

I

I

O

I

C

I

O

O

O

I

‘ :
r

 

 

(1+1/2 ,j-1/2 , -1/2) (1+1/2 ,j+1/2 , -1/2)

Fig. 5. 4.a Integral Path For Ex

 

 

 

 

(iaja'l) (i9j+19'1)

Fig. 5.4.b Integral Path For Hx

190

where the first and the second terms on the right hand side are contributed by the con-
ducting and the polarization sheet currents. Equation (5.3.5) can be further simplified

to

“1.1.1 n+‘.1.1
1H, 7(1+-2-.1+-2-, k)-H, I(1+3, )+—2-, k)]Az

:- . l . 1 2' . 1 . l
- [H, (z+-2-, j, [(+5) — H, (1+3, j, k—-2-)]Ay
= [otAy0.5 + (e — eo)tDELTAy/At + eoAyAz/At]E§,"+l)(i+-% ,1, k)
+ [CtAij — (c — coy/mo: — cooyAz/A:]E§,">(i+—; J, k) (5.3.6)

The modified Yee’s scheme is then obtained as

"+1)"_1. =__§_B_ n)..l
E; (1,}+2,k) CAE§(z)+-2-,k)

 

CoAtz ”$1.1 “$1.1
+ H —, —,k-H, +—, —,kAz
AyAZCA {[ 2 (1+2 1+2 ) (1 2 1-2 )1
+1 1
" 2' . 1 . 1 "+2- . l . l
- [Hy (1+3. 1. k+-2-) - H, (l+-2-. j. k--2-)]Ay} (5.3.7)
where coefficients CA and CB are defined as

CA = [coAtzoot/A205 + (e, — 1)t/Az + 1] (5.3.8a)
C8 = [coAtzoot/Az05 - (e, — 1)t/Az - 1] (5.3.8b)

Note that E is continuous across the sheet, only averaged value of E was used to

approximate the fields on the area closed by the integral contour.

In a similar way, a modified scheme for E, can also be derived:

n l . . 1 _ CB ’1 . . 1
5f, + )(IJ'I'E'. k) - " ETA-E; )(lJ‘I—Z'. 1‘)

 

_CoAtz MI.1.1_~+%.1._1_
AxAzCA {[H‘ (”2’92” H’ ° 2’fl2’knAZ

_ z. -_1. .1. _ s--1 -1

[Hz (1.1+2.k+2) Hz (LJ‘I’ZJC 2)]AX} (539)

where

191

= [coAtzoot/A205 +‘(e, - l)t/Az + 1] (5.3.10a)

= [COAIZOO‘t/Az05 — (e, - 1)t/Az — 1] (5.3.10b)

As discussed in Chapter IV, the scattered fields can be significantly modified
when a magnetic material is introduced in the sheet. If 0 = 0 and e = so, we have a
magnetically lossy sheet, which is characterized by the parameter of on, or u, - 1, over
the cavity on k= 0 plane as shown in Fig 5.4(b). Then the finite difference schemes
for magnetic fields are created in a similar way as we did previously for electric fields
with an electrically lossy sheet. Since equivalent electric current Jfi" is zero when
G = 0 and e: co , the magnetic fields are continuous across the sheet. By using
Farady’ law of (5.3.2a) and integral along the contour shown in Fig 5.4(b), a

difference scheme for the magnetic field H, can be constructed as follows:

15:11. j+1. 14%1— 15:11.1. 14%)]112 — 15311.11; 1+1) - 15:0 11% k114y=

(n+}) 111—}1
— omtAy0.5[H, (i, j+2, k+— 2) + H, (i,j+— ,k+—)]

2

—-(11 110)tAy/At[H, 7(1', j+- ,k+—)—H:" 7(1, j+— ,k+—)]

2’ 2’
1 l n-
- uOAsz/A1[H,n+7)(i, 14%, 14%) — H: 7 (1', 14%, [tr-1%)] (5.3.11)

Physically, the first and the second terms on the right hand side of (5.3.11) are contri-

buted by the equivalent magnetic currents. Rearranging (5.3.11) yields

H5711. 11%. 14%) = -3311 x 7 (1 1+; 1.1%)- 6121:2152",
{[520 j+1, k+—)-E:(1,j,k+—)1Az-[E’,’(i 1+5, k+1)—E’,’(i, 14%, k)]/5y} (5.3.12)
where
CA = OmtcoA105/(20Az) + (u, — l)t/Az + 1 (5.3.13a)
CB = omtcoAtOj/(ZOAZ) - (11, - 1)t/Az — 1 (5.3.136)

Similarly, the y component of magnetic field can be updated by

192

( 1 (1.. 1) c t
H“r’(1+i,,~, k+i) = — 9111, 7 (1+i.j.k+i)1— ——£——
7 2 2 CA 2 2 CAAzszo
{$3141. j, k+%)—E;‘(i, j. k+-;-)]Az—[E,‘(i+-;-, j, k+1}—E;’(i+-;-, j. k)]Ax} (5.3.14)

Both the electrically lossy sheet and the magnetically lossy sheet have been con-
sidered separately, and the Yee’s difference schemes have been modified for these
two special cases. However, most practical materials have both electrically and mag-
netically lossy properties. Thus, the electric current and the magnetic current are
presented simultaneously. In this general case, the formulations for the two special
cases discussed before can be combined and modified. It is informative to note that
the tangential components of the electric fields and the magnetic fields are no longer
continuous across the covered sheet when both the magnetic and the electric losses
exist. Therefore, the difference equations for updating both the electric fields and the

magnetic fields comprise both electric and magnetic sheet currents.

By the Ampere’s Law, we have
£11.17: 3[1.11:1 + 111 _ 1539:2421 + lac-(387E113

The discontinuity of electric fields across the sheet can be compensated by a magnetic
sheet current of 15,31 . Inside the sheet, electric fields can be approximated by the aver-

age value of the electric fields on both sides of the sheet. As shown in Fig 5.4(c), if

E (1+l ' 0*)

x 2 J.
or
. . 1 +
Ey(lJ+—s 0 )
2

is above the sheet, then the electric field inside the sheet is approximated by

. I . +
E +- ,0 +0.5!“

Ez

    

Ez

    
 

193

By

 

 

 

W

Ez
(ot,et,omt,ut)

 

 

 

 

1 IOCOOOIIOOOOIIOMOORR \\

V

Hy

(ot,et,omt,p.t)

   

 

Fig. 5.4

Hy

(c) Integral Paths

194

or
5,0,)»; 0*) - 0515,31,
where
r“? = c mal-1'. 0) + (u — 110—3- H<i+-‘— 1'. 0)
”' "' 2 a1 2
While the electric fields are below the sheet, they can be evaluated by a difference of
magnetic sheet currents as
E (143-) 0') = E (1+l ,1, 0*) + 0 1H (i+—1-J 0) + (11 — 110)1-a-H (1+1) 0)
x 2 ’ x 2 "' y 2 ’ a1 7 2 '
or
502141. on = 5 (ml. 0*) - o :11 (mi. 0) - (u - 11011311 (114—1. 0) .
7 2 7 2 "' " 2 a1 ‘ 2
It is seen that the right hand side of the Ampere’s law can be completed by finite
differences. Unfortunately, it is found that those approximations will fail if any time

derivative is involved. The difficulty can be viewed from the example of the third

integral on the right hand side of the Ampere’s law as shown in Fig.(5.4c).

The third integral ( i eo-EaTE-dts’ ) is approximated by

31 <~+%>.1., <n+§>.1.,
— —AS E +— ,0 + E, +— ,0 +

(n ‘1 (n ‘)
o,,,1Hy +7 (1+7: ,1“, 0) + (11 — gong-H, +7 (14% ,1, 0)) } (5.3.15)
and the derivative in time domain leads to:

ASEO/A"{E§"H)(i+%1f1 0*) - atria-$— .1'. 0*) + 0.56mrtH§"+”(i+-§ .1'. 0) — H;~>(1+% 1'. 0)]

M l
+ (11 — 110)1/A17[H§"+‘>(1+% ,1, 0) - 2H: 711%.), 0) + H§")(i+% J, 0))} (5.3.16)

However, the left hand side of the Ampere’s law l H-JI yields:

195

( 1) (n 1)
i 1111? = [H,"+7 (14%, 14%10) — Hz +7 (14%, j—%-,O)]Az —

<»+)>.1.1 (“9.1.1
[Hy (1+3! 1! '5') _ Hy (1+?! 19 _EHAy (5.3.17)
The iteration of Hy needs its information at time steps of n+1 and n on the right

hand side of the Ampere’s law, but time step of 114—3- on the left hand side. To over-

come this difficulty, the electric field E,(i+-;-,j, 0) or Ey(i,j+%, 0) is assumed to be the

field inside the sheet. Then, the integration over the area as shown in Fig 5.4(c) is
equivalent to taking the electric field inside the film as an averaged value. With this
assumption, the representation for the Ampere’s law remains the same, but that for the

Falady’s law is changed as shown in Fig 5.4(c). By the Farady’s law
_ 3 a
E-d‘l _ — Jm-d's‘ - (11 - 11,)_n.ds - 110—1111?
31 31
Consider the electric field inside the film, the total current closed by the integral con-

tour is -;-o,,,tdyH, or é—omtdyH, . After some manipulations, the updating formulas for

magnetic fields become

(M111). .1 1 CB <-}>. .1 1 6013!

H , —,k+—=——H, ,+—,k+— ——

7 (”72 2) CA ('1 2 2) CAAszzo

[52(1, j+1, k+-;—)-E’,‘(i, j, k+%)]Az—-[E;(i, j+%, 1+1)—E;(1', 14%, k)]Ay} (5.3.18)
and

(Mi-1.1 . 1 CB (rt-i). 1 . 1 coAt

H +—,,k+—=-—H —,,k+— __—

7 (l 2’ 2) CA ’ (”2’ 2)] 01111sz0

{1520“, 1'. k+%>—E:(i. 1‘. “firm—[Bani 1'. k+1)-EZ(i+%. 1'. 1mm} (5.3.19)

where

CA = omtcoA10.25/(20Az) + (u, - 1)t/Az + 1 (5.3.20a)

196
C8 = 0,,1c0A1025/(20Az) — (11, - 1)t/Az - 1 (5.3.20b)
Equations (5.3.18) and (5.3.19) hold also for k = -1 and k = 0 where the magnetic

. ( ‘)
field components are Just above or below the sheet. For example, Hyfiz (14%, j, -%)

l

(n ) . . .
or H,”- (14%, j, %) rs below or above the sheet Wthh rs located at plane k = 0 . The

formulas of (5.3.18) and (5.3.19) evaluated at k= -1 and k= 0 are almost the same

except for some changes in the coefficients of CA and CB .

5.4 Modification of Radiation Boundary Condition

For the geometry shown in Fig 5.1, a radiation boundary condition needs to be
applied to truncate the half space into a finite one. The outer surface of the truncated
space is depicted in Fig 5.2 where the infinite ground plane is included. Because both
the second order and the first order radiation boundary conditions use the points below
the ground plane, it is necessary to use the image theory. By the image theory, the
second order radiation boundary condition can be kept for the corner points on the

ground plane while the stability criterion is unchanged.

As discussed in Chapter III, with the Yee’s model only the components of electric
fields are determined from the field values at points outside the considered region if
the outmost surface is placed as in Fig.5.6. Thus, the radiation boundary condition is
only applied to the electric field components on the truncation surface. The Mur’s
second order radiation boundary condition for a three dimensional problem as stated in
Eq. (3.3.35a) of Chapter 111 can be illustrated in Fig. 5.5. If the outmost truncated sur-
face is assumed to be on the Jr: 0 plane, then x>0 region belongs to the interior
region. The field at point (0, j, k) is determined by all its adjacent points as shown in

Fig.5.5.

197

W(O , j, k+l )

 

 

 

 

W(l, j, k+1
W(O, j-l, 1:) Wm. j. k ) W(O. j+1, k)
W(OoS’ j’ k)
. w 1, ' 1, 11
W0. r-l. k) W(l, j, k) ( 1+ )
W(O. i. k-1 )
W(lv j, k'1 )
X

Fig.5 .5 Points Used in the Mur’s Second Order Approximation

198.

Figure 5.6 illustrates a radiation boundary condition for a structure with an
infinite ground plane. W represents any components of electric fields at a proper posi-
tion in the Yee’s model. It is assumed that the ground plane is located at k = 0 plane.
Since E, and Ey at k: 0 are tangential components which lead to zero, the Mur’s
second order approximation and the first order approxirnationcan be directly applied to

them.
To update the normal component of electric field above the ground plane, E, com-

ponent at k =-% , the knowledge of field below the ground plane is necessary. For-

tunately, the image theory can be used and the fields produced by the sources and their
images are intuitively illustrated in the Fig. 5.7. The real sources above the ground
and their images provides the actual solution in the space over the ground plane. The
tangential components of electric field produced by the sources and their images are
antisymmetric w.r.t. the ground plane, while the normal components of electric field

are symmetric w.r.t. the ground plane.

In our specified problem, the aperture fields can be regarded as equivalent sources

on the ground plane, so the same image theory can be applied to analyze their sym-

me ro erties. For exam 1e, to update E,(z', j, -1-), the normal com onent of electric
try P P P 2 P

field E,(1’, j, — g? in the radiation boundary condition is replaced by the field E,(i, j, %)

according to the field symmetry from the image theory. This technique is applicable
to the second and the first order approxiamtions on the points at the intersection of the

ground plane with the outermost plane. It can be shown that the resulted scheme is

stable. As an exmaple, the radiation boundary condition for E2 at points (1', j, k+-;-) on
the plane x = 0 is represented by:

coAt-Ax
coAt+Ax

 

' 1 1 n— . 1 .
E',‘+l(0,1,k+-2-)= —E, 1(1, j, 1+3) + [E’,'+1(1,j, k+-;-) +

199

 

 

 

Ex Ex
11
E2 1 Hz
Ex Ex
11 . 11 E2
(0 ,j,1+1/2) E2 (1 ,J,1+1/2)
Ex
Ex
11 . 11
(O,j,1/2) Ez (1,1,1/2) Ez
Ex

 

//

 

//://.

//7/

 

////

 

 

, 1i Ez , 11 E2
(011,-1/2) (1.11-1/2)
X=0 X: A x
Fig. 5.6 Mur’s Difference Scheme Near Ground Plane

‘
v

 

£[///////{//////
1 C7: ..... 0""

Fig. 5.7 Sources and Their Images

201

2A1

E... . 1
1 9 9 k — + _—
z (0 J + 2 )] coAt+Ax

[£10. 1'. 141%) + 52(0. 1'. keg-)1

+ (comm
2Ay7(coA1+Ax)

 

15:10. 1+1. 1%) - 251(11): 1%) + 51(1). 1—1. 11%)

+ 15*;(1, j+1, k+—;-) — 2153(1, j. 1%) + £20.14. 11%))

(c0A1)7Ax
2Az2(coAt+Ax)

 

[52(0, j, k+1+-;-) — 252(0, j, 141—é) + 53(0, j, k—%) +

52(1, j, k+1+%)- 252(1. j. 11%;) + 53(1, j, k-%)] (5.4.1)

Equation (5.4.1) is the Mur’s second order formulation at x = 0 for the wave pro-
pagating against the x-axis direction. It is exactly the same as (3.3.35a) if E, is inter-

changed with W .
When k = o, E,(O, j, +-;-) is the field just above the ground plane, but 53(1), j, —-;—)
or E’,‘(1, j, —-;-) in (5.4.1) are the fields below the ground plane which should be substi-

tuted by 52(0, j, +7}? or E’,‘(O, j, +%). Therefore (5.4.1) becomes

. 1 ".1 . 1 Com-Ax
E"+‘0.,+- =—E 1..— +
7 ( J 2) 7 ( J 2) coAt+Ax

 

152%. 1'. -;-) + E'r‘to. 1'. ~50]
2Ax
coAt+Ax
+ «0130711:
2Ay7(c0A1+Ax)

l
2

 

[53(1, j, 2') + E,'(0. j. )1

 

[152(0, j+1, 4%) — 253(0, j, 1%) + 13*;(0, j—l, —5—-)

. 1 .1 . 1
El, 19— -2: 09— z 9_9—
+ z( 1+ 2) 5(1) 2)+E"(111 2)]

 

(coAt)2Ax , 1 . 1
E09 ’1 — -5: ‘l ’—
2Az7(coA1+Ax)[’( 1 +2) (01 2))+
520.1. 1%) - 520.1; -;—) )1 (5.4.2)

202

5.5 Fields in the Scattered Field Region

In the scattered field region, only outgoing waves are considered so that one can
apply the radiation boundary condition. As shown in Fig.5.8, the total fields can be
decomposed into incident field E‘, reflected field E’ with the aperture absent, and scat-
tered field produced by the aperture field E’ . It is easy to see that the field produced
by aperture field is an outgoing wave exiting the outermost truncated surface, and the
incident field and the reflected field are known. The total fields then can be

represented by:

E‘“"’(r) = E‘(r) + E'(r) + E‘(r) (5.5.1)
where E‘(r) + E'(r) are known functions and E" is an outgoing wave exiting the trun-
cated region. A radiation boundary condition now can be applied to the E’ in the

scattered field region.

5.6 Backseattered Fields

When an excitation or driving force is a sinusoidal function, a steady state is
reached after a certain period of time. The fields at any points in the truncated region
can be obtained after sufficient iterations. To calculate the far zone scattered fields or
radiation fields, it is easer to use equivalent surface sources based on the equivalence

principle as shown in Fig 5.9(a).

On the aperture S, , WE at 0 and 1b<H a: 0 . On the ground screen So , ME = 0 but

fixH #3 0. The scattered fields are evaluated as

133(1) = .4L 1 [(fi’xE)xV’<I> + (r‘i’-E)V’<l> - jtou(r’i’xH)<I>] dS’ (5.6.1a)

7‘ Sa+5p

H‘(r) = 311; I [(ii’xH)xV’<D + (fi’-H)V’<b + jwu(fi’xE)<b] dS’ (5.6.1b)
SO+SP

203

 

/ k / Incident Fields

Scattered Fields

 

Incident Fields

Reflected Fields

 

 

Scattered Fields From Aperture And Cavity :

\(Total Field - Incident - Reflected Field )

 

 

\'\

 

Fig. 5.8 Decomposion of Total Fields in the Scattered Field Region

\\\\\\
\“\
\
\“‘
.\
‘\
‘\
\
~
‘\
\
‘\
\
‘\
\
‘s
‘\
u

s
\
I
u
\
u
h
s
I
‘i’ ‘
§
\
\
\

s
s

\\\:

 

Excluded Source Region

 

a
'o
C
l
I
’I
o
a
'a
'0
.0
'0
'1
’I
0"
0",
'''''
llllllll

llllll

+ A
Electric Images 3: p-e (n . E )

K=/r\1xH

b

——>I

=3)

+++

\\ \ .\ \¥\\

3+
;: p,=u(?1 - H)

Magnetic Images

 

A
Km=-(n x E)

(b) 3 so

Fig. 5.9 Aperture Fields and Their Images

205.

where the integration is over the aperture and the entire ground plane.

Since

equivalent surface currents are only known on the finite aperture region, it is desirable

to eliminate the contribution of rod! and fioE over the ground plane. Because the sur-

face fields maintain the null field within the excluded region, then nothing is changed

if that region is backed by a perfect conductor. The resulted image system is indicated

in Fig.5.9(b):

By the image theory, electric sources are annulled by their images while the mag-

netic source and their images are additive in effect. This leads to the conclusion that

r‘ixE is the only effective surface current such that:

133(1) = J—J' (fi’xE)xV’<l> 115'
21: P

1

H’(r) = “27:. [(rY-HW’QJ + jcott(n“’xE)<I>] dS’

..ij

where (b = 573-— and R = Ir — r’l . Under far zone approximation:

vb .
V'cb = 171—‘3 e77"
r
when Irl>lr’| and kR>1 . The radiated fields are:

E .1: i 8-170
21:

 

1x! (fi’xE) e17” 115’
P

H = (fXE)/T|

and radar cross section is given by
c = 211r714-1E or 12
E
As r approaches oo (5.6.5) becomes:
0 = 1‘3-1 (rb<—E-)xf e177” 11512
215 E0
0
or

.9.

k4 E 'kf'f‘ 2
= —| W— X? 5" (15,1
12 811:3 (f E0)

(5.6.2a)

(5.6.2b)

(5.6.3)

(5.6.4a)

(5.6.4b)

(5.6.5)

(5.6.6)

(5.6.7)

206

In the y-z plane, the unit vector r“: cosef+ sine)? , ?’=2(5c‘+ y’y‘ and

fixE = Ej — Er)? . Radar cross section in y-z plane becomes:
k4 . ,.
8—;IJ(E,£C038 + Eyy‘cose — E,z‘sin9)e”‘y smedS’l2 (5.6.8)
TC 0

Similarly in the x-z plane, f = c0502" + sinef , r" = £16 + y’)? and m<E = Ej — E): .

Radar cross section in y-z plane is given by

4 . .
filflEficose + ch056 — E,.i‘sin9)e”“"““6dS'l2 (5.69)
n 0

5.7 Numerical Results

5.7.1 Open Cavity situated in the Ground Plane

Use the algorithms discussed in proceeding sections, fields at any points in the
truncated region can be calculated. First consider an open cavity situated in the
ground plane without an antenna inside. The field distributions inside the cavity and
on the aperture, and the radar cross section are calculated when a plane wave is
incident on the cavity. The effects of an impedance sheet on the radar cross section

and field distributions are also studied.

As shown in Fig.(5.10), each side of the aperture is one wavelength long, and the
depth of the cavity is also one wavelength long. A plane wave is normally incident on
the ground screen. In the FD-TD computer code used for all the following numerical
data, each side of the cavity is divided into twenty-one cells and the radiation boun-
dary is set about fifteen cells away from the aperture. Two cases are considered which
are an open cavity with empty aperture and an open cavity with its aperture covered
by an electrically lossy film. Figs. (5.10a) and (5.11a) show the distribution of the
tangential components of electric field E’ on the aperture without film covered, and

Figs. (5.10b) and (5.11b) show that when the aperture is covered by an infinitely thin,

207

electrically lossy sheet of or: 0.01. By comparing these two cases, the difference
between field distributions can be seen. Fig. (5.12) shows the radar cross section in
the X-Z plane and the Y-Z plane for the two cases. The radar cross section is reduced

by 10 db in Fig.(5.13) when a resistive film is presented.

When aperture is empty, Figs.(5.10a), (5.11a) and Fig. (5.14) exhibit three com-
ponents of the electric field E, on the aperture and E, is plotted half cell above, E2, and
half cell below, 53W, the plane of aperture according to the location of Yee’s model.
It can be seen that the x component is much greater than the y component of electric

field and E, component is changed greatly across the aperture plane.

When the aperture is empty, the tOtal fields inside the cavity at the plane of
z = —d+5*dz which is 5 cells away from the bottom are plotted in Figs. (5.15-5.16). It
is also interesting to see the variation of the total field distributions along the z direc-
tion in both the x-z plane and the y-z plane as shown in Figs. (5.17-5.21). The fields
of E, and Hy are two dominant components which are also plotted for the case of an
electrically lossy impedance sheet with or = 0.01 covering over the aperture as shown
in Figs. (5.40-5.44). Comparison of Figs. (5.20-5.21) with Figs. (5.17-5.18) indicates
that the structures of field distributions are the same while the amplitudes of the fields

are reduced as an impedance sheet is covered over the aperture.

5.7.2 Cavity Backed Antenna

Receiving Characteristics

Figure (5.22) shows a cavity backed antenna where cavity size is the same as
described above. The antenna is placed parallel to the incident electric field. The
plane wave is normally incident on the cavity. The antenna is about half wavelength

long and one cell size thick. It is centrally fed by 50 ohm. The antenna is located

208.

17/20 wavelength away from the bottom of the cavity where the x component of elec-
tric field reaches maximum in the absence of the antennna as shown in Fig.(5.20).

The loaded impedance is accounted by taking the conductivity of central cell of the

 

antenna to be c = where Z", = 500hm .

AszZ,,,

The aperture fields are shown in Figs. (522-523) and the radar cross sections are
plotted in Fig. (5.24). We didn’t see a significant difference in the field distributions
for the case of an open cavity and the case of a cavity-backed antenna for this size of
cavity. The radar cross section is reduced more than 10db if an infinite thin film with

or = 0.01 is placed over the aperture.

As shown in Figs. (5.25-5.26), the field distributions inside the cavity are greatly
changed when the film is on the aperture. The discontinuities of the fields at both
ends and the center of the antenna are observed. The boundary conditions on both the

walls and the antenna are seen to be satisfied.

A comparison of current distributions on the antenna is given in Fig. (5.27) which
shows the effect of an electrically lossy sheet on the antenna. The received power of
the antenna without a film covered is -54.4 db, while it is reduced to -65.6 db when
the film is covered. Unfortunately, the received power is also reduced about 10 db,
which is almost the same as the radar cross section. This is not what we expected.
We expect to reduce the radar cross section of an airborne antenna but not to degrade

its receiving property by covering the antenna with an impedance sheet.

The comparison between the open cavity and the cavity backed antenna on the
ground plane shows that the reductions in their radar cross sections when the film is on

the aperture] are almost the same, which is about 10 db, for the geometry chosen.

The unsatisfactory reduction in the radar cross section may be due to the fact that

the scattering properties are strong functions of cavity parameters. In our numerical

209.

example, the cavity is at resonance. The scattered fields produced by the aperture and

the cavity are dominant, which is much stronger than that produced by the antenna.

210

h

 

 

 

 

 

 

(a)
36523
X 18261
1.11
0 0000
Ex w
11;
fl "hfllV/
7' 1.
(b)
(OI-0.01)
X
LA

 

Figure 5.10 (a) :1 -component of scattered electric field on the empty aperture of an
opened cavity in the ground screen; (b) 1: —component of scattered electric field on the
aperture covered with an impedance sheet.

211

R:

 

 

(a)

 
   

M9951

R)" 20972 i

 

/ If; 7‘ I.
(b) /;
(at-0.01)
\“\ 4’

 

Figure 5.11 (a) y-component of scattered electric field on the empty aperture of an
opened cavity in the ground screen; (b) y -component of scattered electric field on the
aperture covered with an impedance sheet.

212

 

15

 

 

 

 

 

 

 

( a )
-15 -J '
all2
((18)
-75 _ -
. In Y-Z Plane
-35 - — In X—Z Plane
45 ..
-55 1 1 1 1 L 1 1 1 1 L 1 1 1 1 1 1 1 i
-% —80 -70 -& -50 -40 -30 -20 ~10 0 10 20 30 40 50 60 70 80 90
0( in degree )
0
( b ) -1o —
-20 _.
-3o ..1 ’ \
I \‘
I \
01112 . .
(dB) ‘40 “ ‘1 1‘
' '.
i .
—50-Ju 1
i 1
40-5 ----- In Y-Z Plum '1
E In X-Z Plane 1
-70 .3 ‘
- 1 1 1 L 1 1 1 1 1 1 1 1 1

 

 

 

-90-80-70-60-50—40-30-20-10 O 10 20 30 40 510 610 710 80 90
0( indearee)

Figure 5.12 Radar cross section of an opened cavity ( koa =21t, dll=l ) when a

plane wave is normally incident on it ( 6 = 0.0 ): (a) an opened cavity with its aperture

empty; (b) an opened cavity with its aperture covered by a thin film ( or = 0.01 ).

213

 

  
     
     

 

 

15
5—
_5_ .' °.
O/xz .°. .0
(dB) 100......0. 0........00
-15...
. . Without Thin ilm
—25— _ Thin Film or = 0. 1
-35 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

 

-9O -80 -70 -60 -50 —-4O -30 —20 -10 0 10 20 3O 40 50 60 70 80 90
0( in degree )

Figure 5.13 Comparison of radar cross sections in the X-Z plane of an opened cavity
( koa = 21:, d/k = l ) when a plane wave is normally incident on it ( 6 = 0.0 ).

214

 

 

 

(87

1\’

El

 

 

 

 

//
(b)

   
     

Q \ \\\\\\\\
1‘0‘.’ “ \

 
 
        

,.‘\\\\\ \\ 2 35
3;;"‘ \\\» 57
\ \\\)". \\
/ ’0. ““§\‘§§‘;“‘" 19995
arm‘s/111') .~
71 [1 ":“:“‘}’7I“"’ \ 142‘s \ett‘
1 / "Oo‘e\§o,51u.. - ,e
3593 ] I ’ :‘éofgl”’.‘ 03
‘J 0 1”, bztilllllllA we" ,0

o,

vity
ed ca w
ope“ be1°
f an cell
0
half
t
8W5. a
mPty . (b)
thc C my
field on the ape
f 6 cc
t 9 31f
11°“ at h
'CO’ . (a) E‘
4 7' no
- d
Flgure gro‘m
. the
m

the aperture

215'

I

53235

 

132‘: “

‘ ‘LDD -
. .._-

nnnnn
ddddd

 

22:29 s H/ ' :\‘,';1E , ‘(ixX?<\
‘ , I .‘\‘\\'zl”'I. \\\ I I I, ( .\\. " ‘1:
- :1 “i /. .:‘In‘h§?blh~‘0\7umllbfi.b \\1\
. ,. I. ‘ I i V i 1 ‘“.. O.‘ ._ \\\: ”~96
..oo , . 1
' I 1 ““ we.“

’-~z‘E)cz

Ez :omoooeot at 1‘9 2:.

I
- 31.10:. 15039
I ' l ‘
33393 ' l \‘v V‘
. ' \, - 13929 «or
27739 ,1 ‘19
5' ,n 5323: - u7u35 ‘,39
k“9:49) \O
1 Shh ,3
1th

54 ccrooneot at the z=<-~zeS):z

Figm'e 5.15 Total electric field distribution at the plane 2 = - d + Sdz which is 5 cells
away from the bottom of the cavity ( koa = 21:, 111/11 = 1 ) when a plane wave is nor-
mally incident on it ( 0 = 0.0 ) and no film is covered on the aperture.

216

11!"

.44--

 

 

 

  

   

~z :ceccceot :“e z=*-“z-5>:z

‘9 :croocect at the z=<-*ze§>:z

Figure 5.16 Total magnetic field distribution at the plane 2 = - d + 51!: which is 5
cells away from the bottom of the cavity ( koa = 21:, all}. = 1 ) when a plane wave is
normally incident on it ( 0 = 0.0 ) and no film is covered on the aperture.

217.

 

 

 

//
\
.‘\\
( a ) ':O.:§.;\\
1:111 3377‘
. 'o.-f~t:"£”l: ’II””I;””’ '
3 E19221I ’ I. ‘c~..... I ’ ”’ 1 l
1 [I’l'”'3;;$.::0011 ”I .. l 1
x 7““6 1 / /1%"’”’:";;::.\”.££’ ' I: l
7 ./’/"' ”11%,0,” "'0' '1'?
“1’3 // //I/ If] l1/ ’/1: I, / /"” l ‘i l: f :. I ’ - U-‘r‘e
22360 / f I” ,. , ‘ 1| 1 I 36‘06
0993'- 1 5 E 6806“ (69:0
~209 - /

7 m ””95 ' - 99721 ~37

 

Figure 5.17 (a) x-component of total electric field in the X-Z plane inside the cavity

with its aperture empty; (b) x-component of total electric field in the Y-Z plane inside
the cavity with its aperture empty.

218

 

(a) i217 77'

 

0
/’~
i
ii .
“ea
Q
~‘
Q 9
.9
1‘:
1,:
/

Ht
.
II
I!
§
0’

- 34798

: 1.‘:
331.21 /1:tE
‘1“

99150 I 1 I ' ‘ ' I 35-706
“ l
i p - 680611 0Q
29809 2\
*7 1” *7986 _ 9979 \ABQ
e<Jg \0
19,, c
A
y C r7

(b)

  

2 8065 '1

Figure 5.18 (a) y-component of total magnetic field in the X-Z plane inside the cav-
ity with its aperture empty; (b) y -component of total magnetic field in the Y-Z plane
inside the cavity with its aperture empty.

219

 

 

r”
(a) 477EW

141:}. 6
2 t" :51 71.0 .9471; ’3‘“
WW I: 22 A? “2' "696
I?“ «’17?— 01,11“ d‘xw". 4,:
: H 7083 VII/// “.\:\‘ '::Z!\;
S i / I ”‘.' - ‘..7~3
_\_4 3300 a /,."/ ‘F”: 9.. .
t - 22360 ’ 7"”;lr _ 36906
0098‘ - 68064
3“204
Jr Neal _99731Qaue\e
n/‘JSUQ
/ A
S~965 V
C)
>-
\
:1

- 0H7u8

 

 

Fi 5 19 (a) 2 2component of total magnetic field in the X-Z plane inside the cav-
itygmw'cith . empty aperture; (b) 2 2component of total magnetic field in the Y-Z plane

inside the cavity with its aperture empty.

(a)

 

 

 

///

 

 

Figure 5.20 (a) 1: -component of total electric field in the X-Z plane inside the cavity
with its apertutre covered by an impedance film; (b) 1: -component of total electric field
in the Y-Z plane inside the cavity with its aperture covered by an impedance film.

221

 

 

 

 

(a)
fl M. 7/
///:;z;z;7,
(at-001)

th/YU

2789‘.

Diane with tfiin ‘11

Hg/YO

 

Figure 5.21 (a) y -component of total magnetic field in the X-Z plane inside the cav-
ity with its aperture covered by an impedance film; (b) y -component of total magnetic
field in the Y-Z plane inside the cavity with its aperture covered by an impedance film.

(a)

”Lb—77'

 

 

 

Lx

C )

 

5" ”‘39
It:

(at-0.01)
//fl*_._//
W

(Zm-SOohm.l.engthofDipole-M2)

(b)

 

 

 

 

Figure 5.22 (a) 1: -component of scattered electric field on the empty aperture of a
cavity-backed antenna; (b) 1: component of electric field on the impedance-film
covered aperture of a cavity-backed antenna.

223

a. __9u

 

 

(a) //-/7x_+_;77

 

I
196 ,
s 1 ,f’ifsz‘2“ 9 ..
‘ /\ .2'771’Q..‘\“‘ 943:“.
0 ‘9‘ ..‘\ 410 Q Q. “
L ° « ' “\‘“ b2‘:il”o":“:‘
no / l ‘5 V :fi‘QQO. .“$ 5 9 o
l !3 -.‘ ' I \\\\\\\ ’ .$ 5 799
\\\).
J95 ' E | I \\“\‘\v‘ o
‘“9 9 2693 1 ‘ (\\\ s 2799 2‘» Ct“
‘ 7 2‘? 215 a 0 (sad
Us 1909:”) g,

( (It-0.01)

// A //
/}.7

(Zm-SOohm.LengthofDipole-KIZ)

(b)

 

 

 

 

 

 

Figure 5.23 (a) y -component of scattered electric field on the empty aperture of a
cavity-backed antenna; (b) y -component of electric field on the impedance-film
covered aperture of a cavity-backed antenna.

224

15

 

(a) . j

5-4

all}
(d3)

-15-<

 

-15-

 

 

-35114lllLllll llll
-90-80-70-w-50-40-30-20-10 0 10 20 30 ‘0 $0 60 70 80 9O

9(indegree)

r—

Radar Cross Secrion in Ill-Z Plane

Cavity Backed Antenna koa = 21: . Depth of Cavity m. = 1
Normal Incident 6 = 0.0

Location of Antenna hi}. =17/21. Zm = 500hm:

(b)

 

-10—J

00.2 .
(dB) ‘3“

 

..so- . . mammal-um
_ ThinFilmOt =0.01

 

 

 

llllJJLLlllllllll
m40-70-d3-50-40-30-20-100102030405060708090

9(indegree)

 

Figure 5.24 Comparison of radar cross sections of a cavity-backed antenna (
koa = 21:, d0. = 1, location of the antenna kl}. = 17/21, Z", = 50 ohms ) when a plane
wave is normally incident on it ( 9 = 0.0 ): (a) Radar cross section in the X-Z plane;
(b) Radar cross section in the Y-Z plane.

225

(a) ”7. 77’

 

 

 

 

 

 

 

 

 

( (It-0.01)
(b) fl 1 / /
/ /

(Zm-SOohm.LengthofDipole-MZ)

 

produced by a cavity-backed antenna with its aperture empty; (b) x -component of total
electric field in the X-Z plane inside the cavity produced by a cavity-backed antenna
with its aperture covered by an impedance film.

226

(a) " ELFV

(Zm-SOohm.Lengtho{Dipole-1/2)

 

 

 

6 6905 *

HI y’YU

 

 

 

 

 

 

( at-0.0l)
(b) // k 7/
/

(Zm-SOohm.LengthofDipole- x12)

1

liq/Y0

 

Figure 5.26 (a) y-component of total magnetic field in the X-Z plane inside the cav-
ity produced by a cavity-backed antenna with its aperture empty; (b) y -component of
total magnetic field in the X-Z plane inside the cavity produced by a cavity-backed
antenna with its aperture covered by an impedance film.

II“

227

The Comparison of Total Currents on the Antenna

 

 

 

 

0.0005
00 nofilmcovered
0.0004— AA lossyfilmOt = 0.01
I"--‘.---‘\
1’ ‘9
I \
I \
I, \\
I \
0.0003— 9
\\
\
\
\
\
\
V
\
\
0.0002— I ‘\
I \
I \
I \
l \
¢ b
0.0001— &,-2—*"'*“'+--_$
z” “t\
\‘\
3” \“
I l I I l
—0.2 —0.1 0 01 0.2
x/K

Figure 5.27 Comparison of the total current density distribution on the cavity-backed

antenna 2

( centra fed impedance Z, = 50 ohms; the antenna is (I). = 17/21 from the bottom
of the cavity; receiving power in the case of no film covered P = - 54.4085 db;
receiving power in the case of an impedance covered aperture P = - 65.6098 db )

CHAPTER VI

SUMMARY

A new set of coupled surface integral equations, which are capable of handling
arbitrarily shaped heterogeneous bodies either with or without a perfectly conducting
body inside, has been derived based on the equivalent electric and magnetic currents
on the interfaces of the body. Its numerical solution for the equivalent surface currents
on the interfaces via the method of moments with vector basis functions has been con-
ducted. The numerical procedures developed have been tested by applying them to
concentric spheres and the results are in very good agreement with the exact solutions.
The absence of numerical anomalies in the procedures used is atributed primarily to
the use of basis functions which are free of ficticous line or point charges and in which
the expansion coefficients are not associated with current flowing parallel to boundary

edges.

The effects of a thin magnetic layer coating a metallic object on the radar cross
section of the object have been investigated. An efficient integral equation formulation
for this extreme case of thin coating has been proposed. To validate the new integral
equations, an alternative method based on the eigenmode expansion for an infinite cir-
cular cylinder coated with a thin, lossy magnetic layer has also been derived. The

agreement between the two approaches is excellent.

A new algorithm for the finite difference time-domain method has been developed
based on the integral forms of Maxwell’s equations for treating the problems of
scattering by a cylinder coated with thin magnetic films. In addition, the effect of the
impedance sheet, covering a cavity backed antenna, on the scattering and receiving
characteristics of the antenna have been studied by the finite difference time domain

method. Yee’s model and the radiation boundary condition have been modified and

228

229

then applied to the infinitely thin impedance sheet and the infinite ground plane struc-
ture. A systematic approach of analyzing the stability of the FD-TD method has also
been introduced and illustrated through examples of several different FD-TD schemes.
Conclusions about the stability of several popular FD-TD schemes have been made,

which are the same as numerical empiricism.

APPENDIX A

In this Appendix, we will prove the following relations:

“a __. J—V-(rbdl) (A1)
(08

an = el—V-(WE) (A2)
qu

By the Maxwell’s equation, we have:
VxH = jcer (A3)

VxE = - jmuH (A4)

Taking the inner product of the unit normal vector r? of the surface with both

sides of Maxwell’s curl equation yields

mm = jcoer’i-E (A5)
riVxE = — jtotui-H (A6)
Using the vector identity:
V-(AXB) = B-(VXA) — A-(VXB) (A7)
we have:
V-(nXE) = E-(eri) — fi-(VXE) (A8)
V'(r’b<H) = H-(er’i) — fi-(VXH) (A9)

Next, on any smooth surfaces, we want to show that Vxfi = 0 :

er’i = (V, + fiainm = Vyxri (A10)

230

231

where V,x denotes surface curl and V,xii= 0 from reference [35]. V90? = 0 can be
easily proved by the defination of surface curl [35].

Since stri = 0 , equations (A8) & (A9) become:
V'(rb<E) = — fi-(VxE) (All)
V°(r‘b<H) = — ri-(VXH) (A12)
Putting (A11) & ( A12) back into (A5) & (A6), we have:
— V-WH) = jwefi-E (A16)
V-(fixE) =11»th (A17)
Finally, we have the relations:

ri-E = -LV-(rb<H)
(D8

fi-H = J—V-(M)
1th

APPENDIX B

In oder to evaluate the various elements of matrix equations in Chapter 2, calcula-

tion of the following three integrals (2.3.16) is required.

 

j r,, . A" d5: j dS r,,,(r)- j f,,(r’) if (15’ (B1)
7’ 7’ 7'
2-j rmvvasr. jrm- jv;-r,,(r)v'<r> dS’ (132)
7’ 7’ 7"
Pi” = _ZEi—Hpfil. [”fn(r’)xV’<D(rP,r')dS’] dS (B3)
m 7’ 7"

.1”;
with R” = IrP - r’l and ch = _6_R_.

The first two integrals are easy to be performed in a local coordinate system

First, we consider integration of Eq. (82).

[rm- J‘ V,'-r,,(r')V'¢ d 5' = — jrm- v IV,’-f,,(r’)<b d 5' (B4)
7’ 7" 7’ 7"

= j d5 V - r,,,- jV,'-r,,(r')<r> d S’
7’ 7"

 

 

1m I,I -ijP I
=7J;i(;) ;i(?) e
l l—n l l-n' _j
naming g dam! I 9R, dé’a‘n’

ll—n

“41.11.! I dédn 11

232

233

Consider Eq. (B1), we first attempt to evaluate AW,
1 l-n
Am= 1! gone

ll-n

- iii}! { (r’ — r‘) e

= :tlfi,[r11€g + rzl€% + r31€g_réqu]

     

      

= i]; [0'1 — rfllqg + (r2 - r3)l’1’% + (r3 - ri)l’l)q]

(BS)

The two integrals have be represented as the functions of the following three sin-

gle integrals:

ll-n

=I J

n=0 §= 0

 

ll-n

=I I

11205,: 0

        

      

1 _
11,1 Zia; I211

and these three integrals have been referred in the Appendix C of [10].

Next, the integral I fm-Afi’m (15 is conducted as follows:
77

I fm-Afm d5 = I fm- ilfI [(r] — rfllfi’g + (r2 — r3)I€% + (r3 - r,)l€q]
17 7P

1 l-n

(B6)

(B7)

(138)

(B9)

= 2m.) I d§ emu... — r...) 2 [(r: - mitt + (r: - r9173. + (r: — ram]

234

After some manipulation, we can obtain:

= ignil;{(rml — rnB)-[(r1 - r3)!5 + (r; - r3)!8 + (r3 — r912] +
= (rmz — r,,,3)-[(rl — r3)!6 + (r2 — r3)l9 + (r3 - r913] +

= (rm3 — r,,,j)-[(r1 - r3)!“ + (r2 — r3)l7 + (r3 - r911] +

where l1 - - - 19 are defined as:

ll—n

11:! I[Ildtdn (BIO)
ll-n

12:; I tildédn (1311)
ll-n

13:! (I nlldtdn (812)
ll—n

14:! lllgdédn (1313)
11-11

15 = a lad: dn (814)
ll—n

16:! I fllrgdidn (815)
11-1]

17: [I g Imdacin (816)
ll—n

1%; l tilndédn , (1317)

19:! ( nlmdédn (318)

235

Now we consider the third integral

 

p1,, =

2:": H 935:“ [an(r’)><V’<D.(r”.r')dS’]
m 7’ 7"

in (B3) , where

 

6””?
3

V’CI>(r,r’) = (r — r’)(1 + ij) R

What we need to do is the double integrations over 7’ and 7" . First, we evaluate the

following Pf"? , the integral over T" with observation point at r in triangle 7’ .

 

 

 

P j 1" px(ri ')(1 + 'th) [M (15 (319)
mn = r' m - l’ j
1" 2A9 (R53
Then, we calculate P3,," by
1
P1” = M $13.53: B20
2A: I7]? ( )

Starting from Pfifi’, and Changing varibles by using area coordinate:
r’=§r1+nr2+(1—§—n)r3 result in p,=r’-r‘=§r,+nr2+(1 —§-n)r3—r‘
Then

1 H.
PM = fit}! g [En + nr2 + (1—§—n)r3 — r’]x (B21)

 

_‘kRP
[r32 -(§r1+ Tlrz +<1-§- 21111212191, at m

= in.“ o [(érr + nrz + (1 - é - n)r3)x(r: - r‘) + (r:><r‘)](1 +jk.-R*) 6;, dfi dn

 

= i1;[12§r1x(r,*,, — .9) + lzflquXO’i ; r.) + 12£4r3x(r,*,. - r‘) + rfilxr‘lzpq]

236

           

               

             

Where
ll-n
"=11 d"
ll-n
’25"! 1’
rr-n
”“11 d"

125 = [ZN—[£74254

Note that these three integrals 12"',,12§q 12,1 P" are independent integrals

Second, substituting Eq. (821) into Eq. (820) for PM, we have:

 

 

               

1" ( w as
' —— 'in
(:qu RH)
pm=i(rrin—r£n)

1 l-n

Pm, = 1:11,! i d g d n (r3, — rmJ)-i1,[12§r1x(r§, — r‘)

+ Izfiqrzxafi, — n) + 12€4r3x(rfn - n) + mergzpq]

l l-n
= 111,, l (E d: d n i1£,[12§r£,,-(r1><r1) + Izfi’rmj-(rzxn) + ’2€q'm,-

('3er + rmj'(rixrrtn)l2pq - lzgrfn'fl'txri) — lzn "I: (rzxrr) — 12?":

”I

(’3er " lzgrfllj.(rlxr§l) " 1251‘" mj'(r2><ri) - szqrmj'flsn’IiJ]

(822)

(823)

(824)

(825)

(826)

(827)

(828)

237

since:
r..,-(r.—><r*..) = #0.,er (829)
1 1—11
P... = 21.21:.) g d a d n [lzfiqri-(m — r.,)><r1) + lzfiqri‘flr. — r.)
Xl'fl + szqri-«ri - rm)xr3) + rin-(rllzg + r212,” + r312§)>(r,- + rmj-(rixrina12]
l l-n
= 21.21:. g i d ad n [Izzfi-«r. — r.)x<12g"r1 + Izfiqrz + Izi’qr.) +
rmj-(rllzé + r212n + r312c)xr,- + rfi'(rmjxr,-)IZP4] (830)
where
r3. = (r..,§ + r..,n + r..,§) (B31)

Finally, we have an expression for PM:

I l-n
P... = $11.11;! 1 d &d n [mg + rim + r..,<1 - é - n))-((r.- — r.,>>< (1332)

(liqu + 1251qu + szqr3)
+ rmj'(r112§ + "2’2n + r302?" - Izfi" - 125:»er + (rmfi + rm,” 4’ rm3(1 - é - Tl))'(rmjxri)lzpq]

ll-n

= il’mil‘n! I! d g d n [[(r,,,l -— rm)§ + (rm2 - rm3)n + rm3]-(r,- - rm)x

(Izgqcrl - r3) + 125%. — r.) + 12%.) + raj-(agar: - r3) + 125402 - r.) +

Izpqr3)xrr° + [(1.011 _ rrn3)§ + (rm? — rm)” + rm]-(rfi,,xré)lzpq]

238

l l-n

= finial I! d 5'; d ll {03", - rm)'((ri " rm/)x(E,12€q(r 1 " r3)

+ 521251402 - r3) + 5.12” ql’3) + (I'm: - rm3)-((ri - rm)><(n12€"(n - r3) + n12fiq(r2 - r3) + nlzpqr3)

+ rm,'(rt — r m)x(1 2501 ‘ l'3) + 1251402 - I'3) + 12p qrs) + rmj'(12£q(r1 - l'3) +

[25,402 — r3) + 12Wr3)xr,- + [(rm1 — rm3)§ + (rm2 — rm3)n + rm3]-(rmxri)lzpq}

If the observation point r and field point r’ are within the same triangle, then R
will be zero for some value of §,n,§’,n’ and the integrals lzquzfiqufiq will be singular.
Fortunately, PM equals to zero when r and r’ are in the same triangle. The above

integral can be numerically calculated by seven-point quadrature rule.

APPENDIX C

In Appendix C, we demonstrate that E on the left hand side of Eq. (4.13) is the
electric field on the out surface of a lossy layer, if the thickness of the layer
approaches zero. The proper boundary condition posed for the Eq. (4.13a) is
r’z‘xE = -J,,,t where t is the thickness of the layer. To justify this statement, it is nec-

cessary to check the singular integral around a field point.

Using a vector identity and Maxwell’s equations [37] leads to

l [— jqud) - Bold/me’CD + -%V’<I> (18’ (C1)

 

= i [— jwu(rb<H)<l> + (an)wi + (92B)V'<1>]d5’

Equation (Cl) holds in any singularity-excluded region . Refer to Fig. C1, the
observed point is excluded by a semisphere Sc and S is composed by 5., + 51+ Sa .

The integration on 5., is zero by the Sommonfeld radiation boundary condition.

Performing the integration around the observed point on S, yields:

J [— jmu(fixH)<D + (WE)XV’<D + (fi'E)V’<D] dS’ (C2)

The first term of the surface integral in (C2) is:
IformSa [- jmu(rb<H)<l>]dS’ (C3)

. . e-jBR 2
= —}(0|.L(ID(H)llm[ R 41ER ] = 0

 

The second and the third terms of the surface integral can be combined:

239

240

J [(r’ixE(r))xV’<b + (r’i-E(r))V’<I>] (15' (C4)

 

_ I: ,. . 1 . e‘VfiR
— J (nXE(r))xn + (n-E(r))](; +18) R dS’

1 . e 1BR I
= J E(r)(-R-+JB) R as

 

This procedure for calculating the surface integrals around singular point is the
same as that described in [37]. However, special attention should be paid to Eq.(C4)
when a thin coating on a metallic object is concerned. In [37], it is assumed that E is
continuous and finite on So, so that an average electric field can be taken out of the
integral. But when the observed point is close the surface of a perfectly conducting
body coated with thin magnetically lossy layer, the tangential component of electric
field is not continuous across the magnetically lossy layer. Since the thickness of the
layer is very thin, the contribution from the integral inside thin layer can be ignored
and the significant contribution roots in the integral of E on the out surface of the
coated layer over S, . Therefore, (C4) is calculated as

J [(fixE(r))xV’<I> + (fi-E(r))V’<I>]dS’ (C5)

e’JBR

dS’
R

 

= E(rlon the our-surface)! (_é- + jB)

. 1 . e-JBR 2
= E(rlon the out-surface) hm (E +JB)—R- 27d?

= ZREO'IM the our-swface)

Thus we see that when thegenral integral formulas is used to represent the fields
for an object with a magnetically lossy coating, the electric field E resulted from the

singular integral should be understood as the field on the out surface of a magnetically

241

lossy layer.

242

BIBLIOGRAPHY

[4]

[5]

[6]

[7]

[8]

[9]

BIBLIOGRAPHY

T. K. Wu and L. L. Tsai. " Scattering from arbitrarily-shaped lossy dielectric
bodies of revolution, " Radio Sci. vol. 12, pp. 709-718, Sept. 1977.

P. W. Barber and C. Yeh, " Scattering of EM waves by arbitrarily shaped dielec-
tric bodies," Appl. Opt., vol. 14, pp. 2864-2872, Dec. 1975

M. A. Morgan and K. K. Mei, " Finite element computation of scattering by
inhomogeneous penetrable bodies of revolution, " IEEE Trans. Antennas Pro-
pagat., vol. AP-27, pp. 202-214, Mar. 1979.

D. E. Livesay and K-M Chen, " Electromagnetic fields induced inside arbitrarily
shaped biological bodies, " IEEE Trans. Microwave Theory Tech., vol. MIT-22,
pp. 1273-1280

J. J. H. Wang, " numerical analysis of three-dimensional arbitrarily-shaped con-
ducting scatterers by trilateral surface cell modelling," Radio Sci., vol. 13 , no. 6,
pp. 947-952, Nov.-Dec. 1978.

J. R. Mautz and R. F. Harrington, " Electromagnetic scattering from a homogene-
ous material body of revolution, " Arch. Elec. Ubertragung., vol. 33, pp. 71-80,
Feb. 1979.

R. F. Harrington, Field Computation by Moment Methods. New York: Macmil-
lan, 1968.

S. M. Rao, D. R. Wilton, and A. W. Glisson, " Electromagnetic scattering by sur-
faces of arbitrary shape, " IEEE Trans. Antennas Propagat., vol. AP-30, pp. 409-
418, May. 1982.

A. W. Glisson and D. R. Wilton, " Simple and efficient numerical methods for
problems of electromagnetic radiation and scattering from surfaces, " IEEE Trans.
Antennas Propagat., vol. AP-28, no. 5, pp. 593-603, Sept. 1980.

[10] S. M. Rao, " Electromagnetic scattering and radiation of arbitrary-shaped surfaces

by triangular patch modeling, " PHD. dissertation, Univ. Mississippi, University,
MS, Aug. 1980.

[11] K. R. Umashankar and A. Tafiove, " Analytical models for electromagnetic

scattering, " Part-I, Final Rep., Contract F 19628-82-C-Ol40 to RADC/ESD, Han-
scom Air Force Base, MA, June 1984.

[12] K. R. Umashankar and A. Taflove, " electromagnetic scattering by arbitrary

shaped three-dimensional homogeneous lossy dielectric objects, " IEEE Trans.

Antennas Propagat., vol. AP-34, no. 6, pp. 758-766, June. 1986.

243

244

[13] D. H. Schaubert, D. R. Wilton and A. W. Glisson, " A tetrahedral modeling for
electromagnetic scattering by arbitrary shaped homogeneous dielectric bodies,"
IEEE Trans. Antennas Propagat., vol. AP-32, no. 1, pp. 77-85, Jan. 1984.

[14] O. C. Zienkiewicz, The Finite Element Method in Engineering Science. New
York: McGraw-Hill,197l

[15] P. C. Hammer, 0. P. Marlowe, and A. H. Stroud, " numerical integration over
sinmplexes and cones, " Math. Tables Aids Comp., 10,

[16] K. D. Paulsen, D. R. Lynch, and J. W. Strohbehn, " Three-dimensional finite,
bondary, and hybrid element solutions of the Maxwell equations for lossy dielec-
tric media, " IEEE Trans. Microwave Theory Tech., vol. MIT-36, pp. 682-693,
April. 1988

[17] A. Taflove, K. R. Umashankar, " Finite-difference time-domain ( FT-TD) model-
ing of electromagnetic wave scattering and interaction problems, " IEEE Trans.
Antennas Propagat., Newsletter, vol. 30, no. 2, pp. 5-20, April. 1988. IP [18]
Ke-Li Wu, Gilles Y. Delisle, " Coupled finite element and boundary element
methods in elecreomagnetic scattering problems, " 1988 IEEE AP—S International
Sympsium.

[19] Dennis M. Sullivan, David T. Borup and OM P. Gandhi, " Use of the finite-
diffemce time-domain method in calculating EM absorption in human tissues, "
IEEE Trans. Biomedical Engineering, Vol. BME-34, No. 2, ppl48-157, Feb.
1987.

[20] Dennis M. Sullivan, OM P. Gandhi and Allen Taflove, " Use of the finite-
diffemce time-domain method for calculating EM absorption in man models, "
IEEE Trans. Biomedical Engineering, Vol. BME-35, No. 3, pp179-185, Mar.
1987.

[21] A. Taflove and KR. Umashankar, "Finite Difference Time-Domain (FD-TD)
Modeling of Electromagnetic Wave Scattering and Interaction problems," IEEE
Antennas and Propagat. Newsletter, Apr. 1988.

[22] K. S. Yee, " Numerical Solution of Initial Boundary Value Problems Involving
Maxwell’s Equations in Isotropic Media", IEEE Trans. on Antennas and Pro-
pagat., vol. AP-14, pp. 302-307, May 1966.

[23] A. Taflove and M. E. Brodwin, " Numerical Solution of Steady-State Electromag-
netic Scattering Problrms Using the Time-Dependent Maxwell’s Equations", IEEE
Trans. on Microwave Theory and Techniques., vol. MIT-23, pp. 623-630,
August. 1975. '

[24] B. Enguist and A. Majda, " Absorbing Boundary Conditions for the Numerical
Simulation of Waves" Mathematics of Computation, Vol. 31. no. 139, pp. 629-

245 ‘

651, July 1977
[25] or). Smith "Finite Diffrence Method,"

[26] A. Taflove and KR. Umashankar, "The Finite Difference Time-Domain Method
for Numerical Modeling of Electromagnetic Wave Interactions with Arbitrary
Structures," Chapter 8 in "....", to be published.

[27] J. Fang and K. K. Mei, " A super absorbing boundary algorithm for solving elec-
tromagnetic problems by time domain finite difference method," 1988 IEEE AP-S
International Symposium, Vol. II, pp. 472-475

[28] G. Mur, " Absorbing condition for the finite difference approximation of the time
electromagnetic field equations, " IEEE Trans, Electromagn. Compat., Vol.
EMC-23. pp377-382, Nov. 1981

[29] A. Taflove, K.R. Umashankar, B. Beker, and F. Harfoush, " Detailed FD-TD
anlysis of electromagnetic fields penetrating narrow slots and lapped joints in

thick conducting screens," IEEE Trans. on Antennas and Propagat., vol. AP-36,
Feb. 1988

[30] A. Taflove and KR. Umashankar, " Advanced numerical modeling of microwave
penetration and coupling for complex structures," Final Report, Contract No.
6599805, Lawrence Livermore National Laboratory, Sep. 1987

[31] C. D. Talor, D. H. Lam and T. H. Shumpert, " Electromagnetic pulse scattering in
time varying inhomogeneous media," IEEE Trans. on Antennas and Propagat.,
vol. AP-17. pp. 586-589, Sep. 1969

[32] D. E. Merewether, " Transient currents induced on a metallic body of revolution
by an electromagnetic pulse media,"IEEE Trans. Electromagn. Compat., vol.
EMC-13. PP. 41-44, May 1971.

[33] K. S. Kunz and K. M. Lee, " A three dimensional finite difference solution of the
external response of an aircraft to a complex transient EM environment: Part 2 -
The method and its implementation," IEEE Trans. Electromagn. Compat., vol.
EMC-20, pp. 328-333, May 1978.

[34] K.R. Umashankar and A. Taflove, " A novel method to analyze electromagnetic
scattering of complex objects," IEEE Trans. Electromagn. Compat., vol. EMC-24,
pp. 397-405, Nov. 1982.