How does the presence of dielectric interfaces affect the computation of scattered fields in electromagnetic applications?

Dielectric interfaces significantly influence the computation of scattered fields by necessitating methods such as the finite element method (FEM) to tackle the complexity introduced . These interfaces create challenges in accuracy due to discontinuities at boundaries that require precise treatment, often employing boundary integral techniques or absorbing boundary conditions (ABCs). Furthermore, the FEM methods have to account for the singularity of the integral equations associated with boundaries, which can be reduced using potential formulations or hybrid finite element-boundary integral approaches . The presence of dielectric materials also leads to shifts in the wave impedance and reflections, impacting the overall field behavior and requiring meticulous domain discretization to model these effects accurately .

What role do higher order basis functions play in finite element analysis of complex geometries, and how are they constructed?

Higher order basis functions provide improved accuracy for complex geometries by allowing for more precise interpolation within finite elements. They are constructed by increasing polynomial order or by adding more nodes within elements, such as higher degree polynomials in Lagrange functions or increased nodes per edge in quadrilateral elements. They are particularly important in providing flexibility and reducing errors in finite element approximations .

What are the advantages and limitations of using edge-based basis functions in finite element analysis?

Edge-based basis functions offer several advantages in finite element analysis, particularly for three-dimensional vector problems in electromagnetics. They provide a robust solution to issues like spurious modes and difficulties in maintaining tangential continuity across material interfaces, which are typical pitfalls of node-based elements . Edge-based elements associate degrees of freedom with edges rather than nodes, which helps eliminate non-physical solutions, also known as spurious modes . They enforce tangential continuity and handle complex geometries and material inhomogeneities efficiently . Additionally, edge bases are curl conforming, which leads to more accurate representation of electromagnetic fields . However, edge-based functions have some limitations. Their integrals do not lend themselves to easy evaluation, and for certain element shapes, such as quadrilaterals, they may not be divergence-free . They also require more complex formulations compared to node-based elements and might involve higher computational costs due to the association of unknowns with edges and faces, potentially increasing the number of unknowns to handle . Despite these limitations, their ability to model complicated shapes with fewer unknowns than alternatives like nodal bases remains a significant advantage .

What techniques are crucial for parallelizing global matrix assembly in large computational problems?

Parallelizing global matrix assembly involves distributing elements across processors, allowing each to compute and update the global sparse matrix. The use of hardware locking mechanisms in shared-memory architectures, efficient communication patterns during matrix-vector operations, and minimizing inter-processor communication are critical for optimization. The balance of computation and communication significantly affects speedup and overall performance .

How does the storage complexity differ between partial differential equation techniques and integral equation approaches in computational electromagnetics?

Partial differential equation techniques are favored in computational electromagnetics due to their O(N) storage and solution times, offering more efficient handling of large problem sizes compared to integral equation approaches that require O(N^2) storage. As the problem size escalates, integral and hybrid methods become unmanageable due to their dense matrices and higher storage requirements .

In what way do variational techniques apply to solving wave equations in electromagnetic problems?

Variational techniques are employed to find extrema of functionals which describe the entire domain of functions of interest, such as the total electric field in electromagnetics. By manipulating wave equations into functionals, derived under the calculus of variations, one can apply principles like the divergence theorem to formulate the problem in terms of minimizing, for instance, residuals in the electromagnetic field .

How do one-dimensional and two-dimensional basis functions differ in terms of their application in finite elements?

One-dimensional basis functions are primarily used for problems involving curves or contours and are expressed in terms of Lagrange polynomials, whereas two-dimensional basis functions are used for structures where the third dimension can be disregarded. Two-dimensional elements allow for simpler computational solutions compared to three-dimensional ones and include structures like rectangular and quadrilateral elements, where the shape functions have unit magnitudes at nodes .

Why is it important to manage the number of unknowns effectively in three-dimensional electromagnetic problems, and what strategies help in this task?

Effective management of the number of unknowns in three-dimensional electromagnetic problems is crucial because the computational complexity, storage, and solution time increase significantly with the number of unknowns. This is a limiting factor in solving large-scale problems as the matrices involved can become dense, leading to O(N²) storage requirements in traditional integral equation methods, whereas partial differential equation techniques like the finite element method (FEM) favor sparse matrices that require O(N) storage . Strategies to manage unknowns include employing sparse matrix storage schemes to efficiently handle the sparse nature of matrices in FEM and using iterative solvers like GMRES, which are better suited for large sparse systems . Preconditioning techniques are also vital, reducing the condition number of a system and enhancing convergence, thereby managing computational resources more effectively . Parallel computing and efficient data structures further optimize the handling of unknowns in large electromagnetics problems ."}]}]}]

What is the jagged diagonal storage scheme and how does it improve matrix processing in vector machines?

The jagged diagonal storage scheme is used for storing sparse matrices efficiently on vector machines by arranging matrix elements into jagged diagonals instead of full rows or columns. Initially, matrix rows are sorted by increasing sparsity, allowing the first element of each row to form the first jagged diagonal. Subsequent diagonals are constructed similarly, producing a collection of subvectors with decreasing lengths. This organization enables entire lengths of jagged diagonals to be traversed during matrix-vector multiplications, substantially enhancing vectorization . This approach also allows similar storage requirements to the CSR format but improves performance by reducing zero-padding, which is crucial for vector processors that otherwise suffer from performance degradation due to short vector lengths associated with traditional schemes .

What are the boundary conditions derived from the integral form of Maxwell's equations for material interfaces?

The boundary conditions derived from the integral form of Maxwell's equations at material interfaces include: 1. **Continuity of the Tangential Components of the Electric Field**: The tangential component of the electric field (E) must be continuous across the boundary between two materials. This means that there is no abrupt change in the direction of the electric field as it crosses the interface . 2. **Continuity of the Tangential Components of the Magnetic Field**: Similar to the electric field, the tangential component of the magnetic field (H) must also be continuous. This condition is crucial because abrupt changes would imply the presence of surface currents, which are not present in ideal material interfaces . 3. **Normal Components of the Electric Flux Density**: The normal component of the electric flux density (D) can be discontinuous at the interface if there is a surface charge present . However, if there is no surface charge, the normal component of D must be continuous. 4. **Normal Components of the Magnetic Flux Density**: The normal component of the magnetic flux density (B) is continuous across the interface if there is no surface magnetic monopole . 5. **Impedance Boundary Conditions**: At interfaces involving lossy or complex surfaces, impedance boundary conditions (also called Leontovich boundary conditions) may be applied. These conditions relate the tangential components of the electric and magnetic fields through the surface impedance, providing a way to account for lossy surfaces or thin dielectric layers . These boundary conditions ensure accurate modeling of electromagnetic fields at interfaces and are essential for solving Maxwell's equations across different materials.

Open navigation menu

Upload

0% found this document useful (0 votes)

12K views360 pages

Finite Element Method For Electromagnetics

This document provides an overview and summary of the book "Finite Element Method for Electromagnetics" by John L. Volakis, Arindam Chatterjee, and Leo C. Kempel. The book covers the theory, development, implementation, and applications of the finite element method for solving electromagnetic problems. It begins with the basic theory and variations of the finite element method, and provides modern applications to both open and closed-domain problems in 2D and 3D. The book is intended for graduate students, engineers, and scientists working with computational electromagnetics.

Uploaded by

Felipe

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as PDF, TXT or read online on Scribd

0% found this document useful (0 votes)

12K views360 pages

Finite Element Method For Electromagnetics

Uploaded by

Felipe

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as PDF, TXT or read online on Scribd

Finite Element

Method for

Electromagnetics

ANTENNAS,
MICROWAVE CIRCUITS,

AND SCATTERING APPLICATIONS

John L Volakis
Arindam Chattcrjet
Leo \320\241
Kempel

The IEEE/OUP Series

IEEE
on Electromagnetic Wave Theory
Donald Q. Series Editor
PRESS
Dudley.
ANTENNAS, MICROWAVE CIRCUITS, AND
SCATTERINGAPPLICATIONS
A volume \\n the IEEE/OUP Series on Electromagnetic Wave Theory

Donald G. Dudley, Series Editor

Employed in a large number of commercialelectromagneticsimulation packages, the finite clement method is

one of the popular and well-established numerical
most techniques in engineering. This book coversthe theory,
development, implementation, and application of the finite element method and its hybrid versions to
electromagnetics.
Finite Element Method for Electromagnetics begins with a step-by-step presentation of the finite element method
and its variations, and then provides up-to-date coverageof three-dimensional formulations and modern
applications to open- and closed-domain [Link] include:
and Ritz methods
Calerkin's MATLAB sample codes
One- and two-dimensional theory and applications Efficient implementation of the finite element
Thrce-dlmeDSicmal development of the method method, sparse matrix storage schemes,popular
using edge elements and applications iterative solvers, eigenvalue solutions
Mesh truncation schemes Experiences on code porting to parallel computers

Integral algorithms forfastimplementation'of the boundary integral matrix-vector products. Written by

experts who have extensive experience in both teaching and implementing this method to many applications,
Finite Element Method for Electromagnetics can be used as a textbook tor lirst-year graduate students, as well as a handy
reference Jor engineers and scientists interested in computational electromagnetics.

John L. Volakis is professor

Department of Electrical Engineering and Computer
at the Science at the University
of Michigan. He has published
than 140 refcreed journal
more articles and more than 140 conference papers
on numerical and analytical techniques in electromagnetics. Dr. Volakis is also coauthor of Approximate Boundary
Conditions in Electromagnetics (IEE Press, 1995) and several book chapters.

Arindam Chaterjce has developed three-dimensional computer simulation of electromagneticfields for scattering
and microwave circuits, and is currently a member of the finite element development group for the HFSS finite
element commercialpackageat Hewlett-Packard.
Leo C. Kempel developed three-dimensional antenna simulation packages using the finite element-boundary
integral method and has extensive experience with all popular numerical techniques in electromagnetics. He
is currently at Mission Research Corporation, Florida, conducting research and development on all aspects of
electromagnetics.

Formerly tbv E Press

I\320\230 Series on Electromagnetic V ave< ttiis joi it series between IH ?,Press and
Oxford Univeisny Press, offers outstanding coverage ol the field, with new titles as well 3s \320\263\320\265\321\200\320\277\320\277\
and revision? (>f'n \302\246 clones
>*\342\200\242.,ncd that mat -tain 1o ig-tcrrn archival srgnifitaiKC in [Link]
wav^S and applications. Designedspecilic.i!!y
for graduate student*-, piactiung engineer, and rcvcarchcis,
ibis\" series provides' all'\302\[Link] volumes tbit explore electromagnetic waves* and applications beyond ihc
iindcryi level
\302\273d\302\273ate

\320\263

Published by John Wiley and Sons. Inc.

0780334 56
IEEE Press 90000
445 Hoes Lane
P.O.Box1331
Piscataway, NJ 08855-1331
1-800-678-IEEE (Toll-Free [Link] Canada)
or 1-732-981-0060

9||78078\320\262||33425\320\262
Order
\320\250\320\225\320\225 No. PC5698
FINITE ELEMENT METHOD
FOR ELECTROMAGNETICS
lEF.E/Ol'P SERIES ON ELECTROMAGNETIC WAVE THEORY

Tlie IEEE/OUPSerieson Electromagnetic Wave Theory consists of new title* as well as

reprimings and revisions of recognized classics that maintain long-term archival significance
in electromagnetic waves und applications.

Series Editor Associate Editors

Donald G, Dudlev Electronuigneth- Theory. Sivtiuring. and
University of Ehud
\320\220\320\273\321\202\320\276\320\277\321\217 Hcyinan
Tel-Aviv University

Advisory Board Equation Methods

\320\251\320\246\321\202\320\265\321\202\320\253

Robert K. Collin Andreas C. Cangellans

Ciise Western Reserve University University at Arizona

Akiru Ishimaru Iitfvgrid Equation Methods

University of Washington Donald R. Wilton
University of Houston

D. S. Jones
University of Dundee Anwrnm, Propagation* uittl
David R. Jackson
University of Houston

BOOKS IN THE IEEE/OUP SERIESON ELEC'I'ROMAGNETIC WAVE THEORY

Christopoulos. \320\241The Tranmihxhn-Lwi1 Modeling Mvthods: Tl.M

Cltimmow. P. C. Tin' Plane ll'avv Sfietlrwu Rf/nvxt'iiiuihn <*( ffleiirvnwftriinfc Fields
C'nllin. R. E,. Fh'lil Thtfttrv of Gukkd Waves. Second Edition
Dudley. D, G.. MalfnttKltifttl Ftnmtiiilhnxjhr Elevtwnuignrtit Tfwory
Elliot, R. S.. Kh'ftrontttfpwtirs1 History. Theory, and Applhttlltms
helsen. L. \320\222.. and Murcuvitz. N.. Radiation tmd SaJilfring i\\t Wun:t

Harrington. R. F.. Field Computation hy Munwtii Meiliud.\\

Joiics\302\27315. S.. Methods in b'/t'itrontagHeilr Wove Propagation. Second lidition

Lindell. 1. V.. Method* for KfartnunagfU'lk Field Anulysis

Peterson c\\ \321\217\\\342\200\236
I'tmtpuwthmat Melhmk
far Electrmnufowiirs
Tai. \320\241, Gi'in'rii/Pi'd
\320\242.. iiml Dyadic Analysis'
I't't'tnr Applied Maduwoik v in Field Theorr
Tui. \320\241.
\320\242..
Dyadic (iiei'n Ftmctitmx in Merirminigiit'tic Theory. Second Edition
Van Hlailel. J.. Singular Eh'cinmmglii'tic Fields and Source*
Wail. J.. Klecinmwpicih' Waves in Sivaiifwd Media
FINITE ELEMENT METHOD
FOR ELECTROMAGNETICS
Antennas, Microwave Circuits, and Scattering Applications

1EEE/OUP Series on
Electromagnetic Wave Theory

John L. Volakis
University of Michigan
Arindam Chatterjee
Hewlett-Packard
LeoC. Kempel

Mission Research Corp.

IEEE Antennas & Propagation Society. Sponsor

IEEE
PRESS

The Institute ol Plcctricul Oxford University Pns

and Electronics Engineers. Inc.. [Link].

New York Melbourne
This book and other books may be purchased at a discount
from the publisher when ordered in bulk quantities. Contact:

IEEE PressMarketing
Atln: Special Sales
445 Hoes Laae, P.O. Box 1331
Piscataway, NJ 0\320\2505-1331
Fax: G32) 981-9334

For more information about PRESS

\320\223\320\225\320\225\320\225
products,
visit the IEEE Home Page:hltp:/'/[Link],org/

f 1998 by the Institute of Electrical attd Kettrunies Engineers, Inc.

345 East 47th Street. New York, NY 10017-2394

All rights reserved. No part of this bonk may he reproduced in any form,
{\321\216\320\263
may it be stwed in a retrieval system or transmitted in any form,
without written permission front the puhlbber.

Printed in the United States of America

\320\256 987654321

ISBN0-7803-3425-6
IEEE Order Number: PC5698

OUP iSBN 0 19 $504799

Library, of Congress Catatoging-in-PubHcation Data

Volakis. John Leonidas,1956-

Finite element method for electromagnetics : with applications to
antennas, microwave circuits, and scattering / John L. Volakis.
Arindiim Chaliffrjse. Leo C. Kenipel.
pi cm.

Includes bibliographical references, .and mdes.

ISBN0-7803-3425*6 (alk. paper)
1. Electromagnetic fields\342\200\224Mathematical models. 2. Finite element
method. 3. Antennas (Electranies) 4. Microwave circuits.
5. Electrons\342\200\224Scattering. 1. Chalterjee. A. (Arindam) II, Kempel,
LeoC.
TK7867.1.V65 \\\320\250 97-4876\320\232

530,. 14'I\342\200\224dc21 CIP

IEEE Press
445 Hoes Lane. P.O. Box 133!
Pjseataway. NJ 08855-1331

Editorial Board
Roger f\\ Hoyt. Editor In Chief

J, B. Anderson A. H. Haddad M. Padgett

P. M. Anderson R, Herrick W. D. Ree\\4

M. Eden S. Kartalopoulos G, Zobrist

M E. El-Hawary D. Kirk
S. Furui P. Laplantc

Kenneth Moore. Director oj IEEE Press

John Griffin, Senior Acquisitions Editor
Linda Matarazzo. Assisumt Editor

Antennas
\320\250\320\225\320\225 & Propagation Socicly, Sponsor
AP-S Liaisonto IEEE [Link] Maiiloux

Cover design: William T. Donnelly. WT Design

Technical Reviewers

James T. Aberle. Arizona Suite University

Jin-Fa Lee. Wmvster Pntyierhnw Imtintfr

Andreas Cangellaris, University of Arizona

D. R. Wilton. University of Houston

Daniel T. [Link] Fane Research Lahorutorv

Oxford University Press

Walton Street, Oxford OX2 6DP

Oxford New York

Athens Auckland Bangkok Bombay

Calcutta Cape Town Diir cs Salaam Delhi
Florence Hong Kong Istanbul Karachi
Kuala Lumpur Madras\302\273 Madrid Melbourne
Mexico City Nairobi Paris Singapore
Taipei Tokyo Toronto

and associated companies in

Berlin Ibadan

Oxford is a trademark of Oxford University Press

Contents

PREFACE xiii

ACKNOWLEDGMENTS xv

CHAPTER 1 FUNDAMENTAL CONCEPTS 1

1.1 Time-Harmonic Maxwell'sEquations I

1.2 Wave Equation 5

13 Electrostaiicsand Magnetoslaties 6
1.3.1 Electrostatics 6
1.3.2 Magnctostatics9
1.4 Surface Equivalence 9

1.5 Natural Boundary Conditions 14

1.6 Approximate Boundary Conditions 17
1.6.1 Impedunce Boundary Conditions 17
1.6.2 Sheet Transition Conditions 14
1.7 Pointing's Theorem 20
1.8 Uniqueness Theorem 22
1.9 Superposition Theorem 23
1 10 Duality Thedrem 23
1.11 Numerical Techniques24
1.11.1 The RiU Method 24
I.I 1.2 Functionals for Anisotropic Media 27
I.I I 3 Method of Weighted Residuals 28
1.11.4Vector and Mairix Norms in Linear Space 2\320\243

vii
viii Contents

1.5
I\320\233 Some Matrix Definitions 31
1,11.6 Comparisonof Solution Methods and Their
Convergence 32
I.I 1.7 Held Formulation Issues 34

CHAPTER 2 SHAPE FUNCTIONS FOR SCALAR AND

VECTOR FINITE ELEMENTS 37

2.1 Introduction 37

2.2 Features of Finite Element ShapeFunctions 3JJ

2.2.1 Spatial Locality 38

2.2.2 Approximation Order 38
2.2.3 Continuity 38
2.3 Node-BasedElements 39

2.3.1 One-Dimensional Basis Functions 39

2.3.2 Two-Dirnensroml Basis Fimcrtons 40
2.3.3 Three-Dimensional Basis Functions 45
2.4 Edge-Based Elements 48
2.4.1 Two-Dimensional Basis Functions 49
2.4.2 Three-Dimensional Basis Functions 53

CHAPTER 3 OVERVIEW OF THE FINITEELEMENT

METHOD: ONE-DIMENSIONAL
EXAMPLES 65

3.1 Introduction 65

3.2 Overview of the Finite Element Method 66

3.3 Examples of \320\236\320\273\320\265-DimensionalProblems in
Electromagnetics 69
3.4 The Weighted Residual Method 71
3.5 Discretisation of the \"Weak\" Differentia! Equation 73
3.6 Assembly of the Element Equations 76
3.7 Enforcement of Boundary Conditions 79
3.7.1 Neumann Boundary Conditions (Homogeneous)80
3.7.2 Dirichlct Boundary Conditions (Homogeneous) 80
3.7.3 Nonzero Boundary Constraints (Inhomogenuous) 81
3.7.4 Impedance Boundary Conditions82
3.K Examples 83

Appendix 1: Sample One-Dimensional \320\274\320\273FEM Analysis

\320\270\320\264\320\270

Program 89
Appendix 2: Useful Integration Formulae for One-Dimensional
FEM Analysis 91
Contents ix

CHAPTER 4 TWO-DIMENSIONAL APPLICATIONS 93

4.1 Introduction 93

4.2 Two-Dimensional Wave Equations 94

4.2.1 Transmission Lines 94
42.2 Two-Dimensional Scattering 95
4.2.3 Waveguide Propagation (Homogeneous Cross
Section) 97
4.2.4 Waveguide Propagation (lahomogeneous Cross
Section) 98
4.3 Discretization of the Two-Dimensional Wave Equation 100
4J.I Weak Form of the Wave Equation 101
4.3.2 Discretization of the Weak Wave Equation 102
4.3.3 Assembly of Element Equations 105
4.3.4 Assembly Example: Waveguide Eigenvalues 108
4.4 Two-Dimensional Scattering 11$
4.4.1 Treatment of MetaLlic Boundaries 119
4.4.2 Absorbing Boundary Conditions 121
4.4.3 Scattered Field Computation 124
4.4.4 Scattering Example Using ABCs 127
4.4.5 Artificial Absorbers for Mesh Truncation 130
4.4.6 Boundary Integral Mesh Truncation 134
4.5 Edge Elements 137
4.5.1 Example I: Propagation Constants of a Homogeneously
Filled Waveguide 144
4.5.2 Example 2: Scattering by a Square-Shaped Material
Coated Cylinder 145
Appendix 1: Element Matrix for Node-Based Bilinear
Rectangles 149
Appendix 2: Sample matlab Code for Implementing the Matrix
Assembly 150

CHAPTER 5 THREE-DIMENSTONAL PROBLEMS:

CLOSED DOMAIN 157

5.1 Introduction 157

5.2 Formulation 158

5.2.1 Field Formulation 159

5.2.2 Potential Formulation 162
5.3 Origin of Spurious Solutions 163
5.4 Matrix Generation and Assembly 164
5.5 Source Modeling 168
5.6 Applications 171

5.6.1 Cavity Resonators 171

5.6.2 Circuit Applications 173

Appendix: Edge-Based Right Triangular Prisms 176

Fundamental

Concepts

The material book is generally considered to be at the level of a

presented in this
graduate student engineer. Thus, the reader is assumed
or a practising to have been
exposed to electromagnetics either through a suitable graduate course or practical
[Link] this chapter, fundamental electromagnetic concepts and theorems are
presentedalongwith the notation used throughout this text so thai a common base is
available to all readers.
Many good texts on general electromagnetics principles and techniques are
currently in print. Some tire considered classical treatises suchas [\\]~[2] while others
are more recent vintage such as [3] [4].Although this chapter presents the minimal
introductory material necessary for the study of the finite element method as applied
to electromagnetics, the interested reader is encouraged to consult these references
for a more completetreatment of electromagnetics.
Upon assumption of a harmonic field, the phasor or time-harmonic form of
Maxwell'sequations are presented along wilh complex material definitions which

permit of loss [Link] natural

the incorporation boundary conditions are
derived followed by fictitious, though useful, approximate resistive and impedance
conditions. Electrostaticand niagnelostatic formulations are discussed for use in
later
examples of element method. Several useful
the finite electromagnetic \321\201\320

are presented
\321\201\320\265\321\200\321\214 including the Poynting, uniqueness, superposition, and duality

theorem*. Time-harmonic Maxwell's equations wiJl fee covered firs!,

1.1 TIME-HARMONIC MAXWELL'S EQUATIONS

Maxwell's equations were originally written as a set of coupled, time-dependent

integral equations. However, of primary interest in this text is the study of har-

harmonically varying fields (i.c, frequency domain) with an angular frequency of

1
Contents

PROBLEMS:
CHAPTER 6 THREE-DIMENSIONAL
RADIATION AND SCATTERING 183

6.1 Introduction 183

6.2 Survey of Vector ABC's 1 \302\2534

6.2.1 Three-Dimensional Vector ABCs 184

6.2.2 Artificial Absorbers 194
6.3 Formulation 201
6.3.1 Scatteredund Total Field Formulations 201
6.4 Applications 204
6.4.1 ScatteringExamples 205

6.4.2 Antenna and Circuit Examples 216

Appendix: Derivation of Some Vector Identities 221

FE-BI
CHAPTER 7 THREE-DIMENSIONAL
METHOD 227

7.1 Introduction 22T

7.2 General Formulation 228

7.2.1 Derivation of the FE-BI Equations 229

7.2.2 Solution of the FE-BI Equations 233
7.2.3. Commentson the General FE-BI Formulation 236
7.3 Excitation and Feed Modeling 238

7.3.1 Plane Wave 238

7.3.2 Probe Feed 239
7.3.3 Voltage Gap Feed 240
7.3.4 Coaxial Cable Feed 241
7.3.5 Aperture-Coupled Microstrip Line 242
7.3.6 Mode Matched Feed 243
7.4 Cavity Recessed in a Ground Plane 245
7.4.1 Formulation 246

7.4.2 Solution Using Brick Elements 247

7.4.3 FIT-BasedMatrix-Vector Multiply Scheme 249
7.4.4 Examples 252
7.4.5 Aperture in a Thick Metallic Plane 255
7.5 Cavity-Backed Antennas on a Circular Cylinder 257
7.5.1 Examples 260

7.6 Recent Advances in the FE-BI Method 262

7.6.1 Finite Element-Periodic Methodof Moments 262

7.6.2 Finite Element Surface of Revolution Method 264

7.6.3 Fait Integral Solution Methods 266
Appendix 1: Explicit Formulas for Brick Elements 267
Appendix 2: Brick Finite Element-Boundary integml Computer
Program 272
Contents xi

CHAPTER 8 FAST INTEGRAL METHODS 277

by S. Buidiganavale and J. L. Volakis

\302\253.I The AdapUvc Integral Method 277

8.2 Fast MullipolcMethod 279

8.2.1 Boundary Integral Equation 279

8.2.2 Exact FMM 280
8.2.3 Windowed FMM 2K3
8.2.4 Fast Far Field Algorithm 284
8.3 Logic Flow 287
8.4 Results 294

CHAPTER 9 NUMERICAL ISSUES 299

9.1 Introduction 299

9.2 Sparse Storage Schemes 300

9.3 DirectEquation Solver 303
9J. 1 Factorization Schemes 303
9.3.2 Error Control 304
9.3.3 Matrix Ordering Strategies 305
9.4 Iterative Equation Solvers307
9.5 Preconditioning 313

9.5.1 Diagonal Preconditioner 313

9.5.2 IncompleteLU (ILU) Preconditioner 315
9.5.3 Approximate Inverse Precondilioner318
9.5.4 Flexible GMRES with Preconditioning 320
9.6 Eigenanalysis 320
9.6.1 Directand Inverse Iteration 322
9.6.2 Simultaneous Iteration 324
9.6.3 Lanc/osAlgorithm 325

9.7 PurallelLzation 327

9.7.1 Analysis of Communication 330

INDEX 337

ABOUT THE AUTHORS 343

Preface

The finite element method (FEM) and its hybrid versions (finite element-boundary
integral, finite clement-absorbing boundary match-
condition, finite element-mode
matching,etc.) is oue successfulfrequency
of the most domain computational methods for
electromagnetic simulations. It combines geometrical adaptability and material gen-
generality for modeling arbitrary geometries and materials of any composition. The
latter is particularly important in electromagnetics since nearly most applications
dealing with antennas, microwave circuits, scatterers, motor and generator model-
modeling,etc. require the simulation of nonmetallie/cumpositematerials. Also, the hybri-
hybridization of the finite element method with integral equation techniques leads \321\202\320\276
fully

rigorous approaches which combine the best aspectsof volume and surface formula-
formulation
techniques.
Because unique features, the finite element
of its method is becoming the work-
workhorse for
electromagnetic modeling and simulations. Many research and develop-
development codes are now available from universities and industry, and these have
demonstrated the utility and capability of the method. Also, a number of commercial

finite element analysis packages are currently available. Typically, these packages do
not yet incorporate the more rigoroushybrid versions of the FEM. However, they
are rapidly evolving to more sophisticated and capablepackages which incorporate
new technologies in geometrical modeling, simulation engines, and solvers.
With the increasing importance of electromagnetics simulation packages using
the FEM, this book should serve as a valuable text for students, practicing engineers,
and researchers in electromagnetics. The original goal of writing the book was to
serve as a text for beginning graduate students Interested in the application of the
finite element method and its hybrid versions to electromagnetics. However, the
authors also recognized a need to report (in a coherent manner) the many recent

advances in applying the method!*) to traditional and new problems in electromag-

The
electromagnetics. result is a book that can serve both beginning students and more
advanced practitioners, The first half of the book has already bee\302\273used in the

xiii
xiv Preface

classroom as part of a courseon numerical electromagnetics at the University of

Michigan. The secondhalf of the book covers primarily work on three-dimensional
CD) developments and applicationswhich have primarily appeared in the literature
over the past 5 \321\203\320\265\320\260\320\263\321\207.
The book assumes that the reader is a first-year graduate student who has
likely taken one advanced course in electromagnetics beyond the standard under-
undergraduate courses. For practicing engineers, it is assumed that the reader is familiar

with concepts of electromagnetic radiation and has an understanding of Maxwell's

equations and their implications. No previous experience in numerical methods is
necessary, but such experience will, of course, help the reader in understanding
the procedure of casting analytical equations into discrete systems for numerical
solution.
For [Link] is expected that the lirst four chapters will be thoroughly
covered with the exception of Chapter 2, which describes of basis/expansion
a variety
functions. At the introductory stage, only the initial sections of Chapter 2 need be
[Link] readermay then return to Chapter 2 as [Link] 3 and 4 (one-

dimensional (ID) and two-dimensional BD) formulations and applications] are

written in a step-by-step process with the assumption that this is the first exposure
of the reader to numerical methods and the finite element method ia general.
Chapters through 5 7 introduce the finite element method and its hybrid versions
for 3D simulations with applications to microwave circuits, scattering, and confor-
rrial antennas. These chapters are written at a more advanced level and cover the
latest applications and successes of the method in electromagnetics\302\273. Chapter 5
(closed-domain 3D applications) is a straightfoward extension of the two-dimen-
two-dimensional
development in Chapter 4 and can he part of a quarter or semester coiirse
which includes Chapters 1-5. FEM implementations with absorbing boundary con-
conditions and the Unite element-boundary integral method for 3D applications are
described in Chapiers 6 and 7. respectively. These arc practical simulations
realistic
and should he of particular interest to practicing engineers and researchers in the
field. Their 2D counterparts are describedin Chapter 4 at a significant level of detail

along with explicit formulas for developingcomputer codes.

Chapter 8 describes some lecent developments on the implementation of the
boundary integral methods for mesh truncation in conjunction with fast integral
methods. Fast integral methods have shown dramatic reductions in CPU and mem-
memory. They are currently the subject of research and will impact the utility and devel-
development of the finite element-boundary integral method
Finally. Chapter 9 presents an overview of storage techniques for sparse
systems, iterative solvers, preconditioning, parallelization, and a variety of details
pertinent to the development of finite element codes. These items were not mixed
with the earlier chapters which discuss the mathematics and applicationsfor the
FEM. Thus, the reader can refer to Chapter 4 at different stages, and as needed,
when developing finite element codes.

/, L. I'olukis
A. Chanerjee
L. \320\241
Krmpe!
June \342\204\22

Ann Arbor, Ml
Acknowledgments

Interest in Ihe finite clement method (FEM) at the Radiation Laboratory of the
University of Michigan began in 1987 by the first author and his graduate students.
The motivation was to model large domains without restrictions in geometry and
material composition. At thut time, two graduate students. Timothy J. Peters and
Kasra Barkeshli,had completed successful implementations of boundary integral
solutions using A-space methods. This O(.VlogiV) approach paved the way for a
fully OiN) linile element-boundary integral algorithm which combined the rigor of
the boundary integral for mesh truncation and the generality of the FEM for
volume/domain modeling. The first of these hybrid implementations was developed
by Dr. Jian-Ming Jin. a graduate ussistant of John Volakis, resulting in \321\217 highly
successful finite element-boundary integral computer program. Versionsof this code
are still in use by government, industrial, and academic researchersin the United
Stales. Another graduate student of ProfessorVolakis, Dt. Jeffery D. Collins, furth-
furthered this work to a body-of-revolution with an integral mesh enclosure. His later
students among them the co-authors. Dr. Chattcrjce and Dr. [Link] Dr,
Daniel C. Ross, Dr. Jian Gong, and Dr. Tayfun 6?demir\342\200\224made significant con-
contributions toward the understanding of 3D problems in antennas, scattering, and
microwave circuits.
The authors are indebted to the entire research group of ProfessorVolakis for

graciously helping in the preparation of the manuscript. We would particularly like

to mention Hristos Anastassia, Lars Andersen. Youssry Botros. Arik Brown, and

Tayfun Oxdemir for proofreading and in providing data and figures. The authors are
grateful to Dr. Sunfl Bindiganavale for co-authoring the section on fast integral
methods in Chapter 8. He also helped in proofreading various sections of Lhe manu-
manuscript. The acknowledgments would be incomplete without mentioning Mr. Richard
Carnes. whose expertise in LATEX made the typesetting of the book a much easier

task. Ruby Sowards typed some sections of the book, and Patti Wolfe helped in
preparing several figures. The authors are thankful to each of them.

xv
xvi Acknowledgments

A lot of encouragement was received from several people throughout the pro-
project.
The authors would like to particularly acknowledge Professor Donald G.

Dudley, Series Editor of LEEE Press, who was instrumental in publishing this

book with the Press and was supportive throughout the preparation of the manu-
manuscript. The comments and constructive criticisms of the early chapters and the final
manuscript by Professor Andreas Cangellaris. Professor Jin-Fa Lee. Dr. Daniel T.

McGrath, and the anonymous reviewersare very much appreciated. The authors are
also thankful to the entire IEEE Press staff {John E. Griffin, Linda C. Matarazzo,

Christy Coleman. and Savoula Amanalidis) for their help.

On the personal front, our acknowledgments would not be complete without
mentioning the support of oxir families. John Volakis would like to express gratitude
to his wife. Maria, and children. Leonithas and [Link] their patience, sacri-
lice. and understanding during the preparation of the manuscript. Arindam
Chatterjee would like to express deep appreciation for the constant encouragement
he received from his parents. He is also grateful to his wife for her insightful criti-

criticisms. Leo Kempel would like to express his appreciation for the support and
patience provided by his wife, Cathy.

John L. Volakrs

Arindam Chatter fife

Leu \320\241
Kempi'l

August 1997
Ann Arbor, Ml
Fundamental Concepts \302\246
Chapter I

eo = 2nf rad/sec since the finite element method for electromagnetics utilizes time-
harmonic fields. The interested reader is referred to one of the excellent general
electromagneticstexts cited in this chapter's introduction for a discussionof the
time-dependent form of Maxwell's equations. For our purposes,we begin with the
time-harmonic form of Maxwell'sequations.
The time-harmonic electric field is related to the lime-dependent electric field
that \342\200\224
(assuming j V\342\200\224T)
by

= + cos{<t>t +
xExq cas{cot + \321\203\320\225\321\203\321\212
\321\204\321\205) cos(w? + \321\204\321\203)
+ zE20 \321\204.) A.1)

where the complex vector

z)
\320\251\321\205.
\321\203.
=
xExOeJ*' + yEytf* + \320\263\320\225\320\273\320\265^ A.2)

is referred to as the field phasor, and similar representations can be employed for the
other field quantities. Introducing these into the time-dependent Maxwell'sequa-
equations [5], we obtain a simplified set

VxH = J+./WE A.3)

V x E =-M-\321\203\302\253\320\264\320\235 A.4)

a,, A.5)
p A.6)
where the corresponding vector field and current phasors are
E = electricfield intensity in volts/meter (V/'m)
H = magnetic field intensity in amperes/meter (A/m)
J = electriccurrent density in amperes/meter2 (A/m2)
M= magnetic current density in volts/meter2 (V/m2)
and the two scalar charge phasors are

p = electricchargedensity in coulombs/meter3 (C/m3)

pm = magnetic charge density in webers/meter3 (Wb/m3)

Both the magnetic current density (M) and the magnetic charge density (pm) are
fictitious quantities introduced for convenience.
Implied in these time-harmonic equations are constitutive relations for an

isotropic medium

D = eE = eoerE A.7)

\320\222
=/uH = \320\264\342\200\236\320\264\320\263\320\235 A.8)

J=crE A.9)
M = crmH A.10)

where two additional phasors

D = electricflux density in coulombs/meter2 (C/ra2)
=
\320\222 magnetic field density in webers/meter2 (Wb/m2)
Section i.I \302\246 MaxweJfs
\320\242\320\263\320\277\320\272\320\263-\320\235\320\260\320\263\321\210\321\201\321\210\320\265
Equations \320\2

are related and H. [Link]

to E equations are an important
link between the
original time-dependent form of Maxwell'sequations and the time-
harmonic form used in the finite element method. Similarly, the phasor forms of the

continuity equations [5] arc given by

(Ml)
= 0 A.12)
V.M+./\302\253?>*

In A.3)-( 1.12), the material constantsare given by

1:
f,, = free space permittivity
= #.854 x 10 farads/meter (F/m)
7
=
U\302\273 free space permeability = x
4\321\217 10 henrys/meter (H/m)

(.,
\342\200\224
medium's relative permittivity constant
=
\320\246, medium's relative permeability constant
G = electric current conductivity in mhos/m (jj/m)
am = magnetic current conductivity in ohms'm (fi/m)

The first two of these (e() and /x()) are fundamental constants while the others describe
the specific material. For example, er is a measure of the material's electric storage
capacity while a is a measure of the material's ability to conduct electric currents or
alternatively as an Ohmic lossmechanism. The relative permeability (ir and magnetic
conductivity am are the magnetic field analogues to e, and or. respectively. For the
purposes of the finite element method, all four of these material quantities may vary
spatially (inhomogeneous) and spectrally (dispersive).
Thecurrent densities J and M appearing in A.3) and A.4) do not include the

presence of impressed sources, In general. J and M can be written as a sum of

impressed (or excitation) and induced (or conduction) currents as

J = Jf + 4e = J, + aE A.13)
=
\320\234 = \320\234;
\320\234,+\320\234\320\263 + (\321\202\342\200\236\320\263\320\235 A.14)

where the subscript \"(\" denotes impressed currents while the subscript \"V refers to
conduction currents. When these are substituted into A.3) and A.4). the familiar
form of Maxwell's equations arc obtained

A.15)
A.16)
where

= - - =
iT j
\342\200\224
\342\202\254r = \302\253' /\302\253\" f'A -,/tanS) A.17)

and

A.1\302\273)

represent equivalent relative complex permittivity and permeability constants. For

notational convenience, the dot over the relative constitutive parameters will be
omitted with the understanding that these still represent all possible material losses.
Fundamental Concepts \302\246
Chapter I

Any one of the representations given in A.17) and A.18) are likely to be found in the
literature with the quantities

tan<5 = \342\200\224
A.19)
\320\270

tan im=?j A.20)

referred to as the material's electric and magnetic loss tangents, respectively.

To summarize, Maxwell's equations in phasor form for isotropicmedia are

V x H = J, +joj\342\202\254E A.21)

V x E = -M, -JtotiH A.22)

A.23)
A.24)

where the phasor form of A.11) and A.12)was employed to rewrite A.5) and A.6) as
given above.

The corresponding integral representations of A.21HI-24) are

H d\\ =
\342\200\242
(J, +>*E) \342\200\242
dS A.25)
1 \321\201 ih

E d\\ =
\342\200\242 - (M,
\302\246
+\320\234*\320\235) dS A.26)
f f

where is the
\320\241 contour bounding the open surface S illustrated in Fig. 1.1 and
rfS \342\200\224
ndS. The circle through the single integral indicates integration over a closed
contour, whereas the same symbol through the surface integral denotes integration
over the closedsurface 5,.which encloses the corresponding volume V, The surface S
associated with the integrals A.25) and A.26) is completely unrelated to Sr which
encloses the volume V.

Expressions A.21) and A.22) imply six scalar equations for the solution of the
six components associatedwith E and H. Thus, for time-harmonicfields. A.21)and

nds

Figure 1.1 Illustration of the differential

d-tctl clement ds and the contour \320\241
Section 1.2 Wave
\302\246 Equation 5

A.22) or A,25) and A.26)are sufficient for a solution of the electric and magnetic

fiejds. The divergence condjtwas A.23; and (J.24), \320\276\321\202

their integral counterparts,
A.27) and A.28). are [Link] fact, these two equations follow directly from
the first two upon taking their divergence and observing that {V x A) = 0 [6] for
V \342\200\242

any vector A. Equations A.21\320\2301 -28)can be easily modified for anisotropic material
fiH be replaced by f E and |f H. respectively,
\342\200\242 \302\246
as well. This requires that and
\302\253E

x
where ? and Ji represent 3 3 tensors[7].
In this text, open scattering and radiation problems will be considered.
Consequently, any valid and unique solution of the electric and'or magnetic fields
must also satisfy the Sommerfeld radiation condition, which describes the field
behavior at infinity

where is the wavenumber = = and A(I is the corre-

\320\272\320\276 free-space (An tOy/i^f)
2\320\273/\320\220,(,

corresponding free space wavelength. This simply states that the Held is outgoing and of
the form e~^\"r/f as r \342\200\224*\302\246
oo.

1.2 WAVE EQUATION

Ampere-Maxwell's Law A.2I) and Law A.22) arc independent first-order

Faraday's
vector equations, and noted
as> they lead to a unique solution subject to the
earlier,
specified boundary conditions. They may be combined together to yield a single
second-ordervector equation in terms of E or H known as the wave equation.
The finite element method is used to numerically approximate the solution of the

wave equation.
Specifically by takingthe curl of A.2I) or (L22) and making use of the other*
the following vector wave equationsare obtained:

kbn,H =
V x _ + V x HI
-Jkn\320\223\342\200\236\320\234
IYJLHJ

where the upper set of equations are for solution of the electric field while the lower

set is for solution of the magnetic field. In A.30). er denotes the relative permittivity
of the media and indicates
\321\206\320\263
the reUtive permeability of the media. For free space,
these two quantities are both unity. Also. Zu = I / >'u = is the free space
\320\243\320\264\320\276/\320\261\320\236 wave
impedance. In materials other than free space, the wave impedanceand wavenumber
are given by Z = 1/ Y -= yjuje and =
\320\272 respectively.
\321\210^/\320\265\320\264.

Utilizing the vector identity

=
\320\251\321\205
\320\247\321\205(\321\204\320\220) x A
\320\220+\321\204\320\247 A.31)

A.30) can be written in a more convenient form

Wa %,E + x v x = -A2\302\273J
- V x A.32)
\320\253^-) E] Ij-i|
Fundamental Concepts \302\246
Chapter 1

for the electric field and

H- +
*\302\247\320\264,\320\235 x V x Hi = -jka YQM + V x i A.33)
Wij j j
for the magnetic field. The significanceof this form of the wave equation is that for
homogeneous materials, the terms within the bracket are zero. Most implementa-
of the finite element method assume a homogeneous
implementations material within each finite
element and hence this bracketed term can be set to zero for those cases.
Another important version of A.30) for homogeneous media is obtained by
utilizing the vector identity

V x E = - V2E
\342\200\242
V x VV E A.34)
in A.32) to get

-V(V E) + V2E
\342\200\242
+ klfrfirE =jk{]ZonrJ + V x M A.35)

In a source-free region, A.35)simplifies to

2 2
= 0 A.36)
These equationsrepresent three vector field components each of which satisfies the
Helmholtz or scalar wave equation

vV + /cV = 0 A-37)

where V/ denotes Ev, E,. Similar

\320\225\342\200\236, partial differential equations can be formed for

the magnetic field from A.33).

1.3 ELECTROSTATICS AND MAGNETOSTATICS

Although for the majority of this text we are concerned with dynamic electromag-
some
electromagnetics, examples of static electromagnetics are included to illustrate basic finite
element principles. Therefore, we presentthe basic equations of electrostatics and

magnetostatics, namely equation and the potential relations.

Poisson's
In this section, we present a basic review of electrostaticand magnetostalic
expressions sufficient for this text. It is not a comprehensive review of either electro-
or magnetostatics, and the interested reader is encouraged to study [8] or [9] for
further information. We begin with electrostatics.

1.3.1 Electrostatics

The fundamental equations of electrostaticsare forms of Faraday's equation

A.4) and Gauss' Law A.6), namely
V x E= 0
V \302\246(<:?)
= A, (].38)
Section 1.3 \302\246
Electrostatics and Magnetostatics 7

The electric field can be expressed in terms of a gauge condition involving a

scalar quantity. \321\204(\320\263)

E=-V0(r) A.39)

With this field egression A.39), A.38) reduces to Poisson'i\302\273 equation (here given in
terms of general scalar fields iind sources since Poisson's equation is also used for
magnelostatics)

V \302\246
I\302\253V#r)| =/(r) A.40)

For electrostatics, the scalar quantity isa voltage = J'(r))

{\321\204{\321\202) and the sources are
volume =
(/'(\320\263)
\342\200\224
hence A.40) reduces to
charges fair)),

V-[fVK(r)] = -/>,,(r) A.41)

Solution of A.41). subject to the appropriate boundary conditions (Dirichlet.
Neumann, and/or impedance conditions), is equivalent to solving the original
Maxwell's equations.
Solution of A.41) is often accomplished either in closed form or numerically.

Closed form solutions are available for only a limited number of boundury condi-
conditions [V]. Hence, it is usually more practical to employ a numerical method such as
the finite element or boundary element methods.
The potential attributed lo a volume charge is given by the integral relation

K'(r) = f
f
\320\223
^p-
GiD{t. r') dV' <1.42)
where the three-dimensional static Green's function is given by

and volume charges are denoted by p,,. This representation is the solution of A.41) in
unbounded space and a pictorial relationship of the primed and unprimed par-
parameters is shown in Fig. 1.2.
J-'or two-dimensiona! situations (e.g., the sources of excitation run from
z = -oo to z = oc and are invariant with respect to r), the following potential
integral is appropriate

~ G-^ir.t1)^' A.44)

where /i, denotes surface charges and G2n is the two-dimensional static Green's
function

2n \\\\r-
Theseintegral
relations are used to determine the potential. V. at some point
in space due to an impressed charge [Link] are used to derive integral
equations for the lotal potential due to an impressed source subjectto boundary
conditions on surfaces within the domain. Such an integral equation (for surface
problems) is given by
Fundamental Concepts \302\246
Chapter 1

(a)

Figure 1.2 Illustration of the geometrical

parameters associated with field representa-
(a) near
representations; zone setup; (b) far zone setup.

rt ds,

where = n\302\246
VK(r),
\320\255\320\230(\320\263)/\320\255\302\253 r')/8n'
=
\320\255\320\241(\320\263, n'- V'G(r,r') = -\302\253'\302\246
VG(r, and
\320\263') n'
denotes the outward normal to the integration surface S. Note that n = w(r') implies
that the unit normal is a function of the integration variables, where = \320\277{\321\202)
\321\217 is a
functicm of the observation variables.
For perfect electric conductors, A.46) can be rewritten in terms of unknown
surface charges. Specifically,by making use of A.39), and relation pje = \302\253-E=
\342\200\242
\342\200\224n
VK(r), we obtain the usual expression

V(t) = V'(r) + \320\223\321\201(\320\263,

\320\263') dS' A.47)
j J ?{pl

where G represents either A.43) or A.45). as appropriate.

Section 1.4 \302\246
Surface Equivalence 9

1.3.2 Magnetostatics

The solution of Maxwell's equations for a stationary magnetic field is similar lo

the procedure given above for electrostatics. In this ease. Ampere's and Gauss'
Magnetic Laws are
VXH=J
V =
\342\200\242
\320\222 0
(..48)

where the static field density is assumed to be related to the magnetic field intensity

by the expression

(L49)

Note that in A.48), we have not assumed a fictitious magnetic charge density.
Rather, the fundamental sources of static magnetic fields are currents. J.
Wc can define a magnetic vector potential in terms of these currents

A(r) = /i f
\320\223
\320\223
J(r')G(r. V-\"
\320\263 U-50)

where the Green's function is the same three-dimensionalfunction used for electro-
electrostaticsA.43) or the two-dimensional function A.45). For the hitler cuse, the integral
must be reduced to a two-dimensional one over the domain of J. With the introduc-
of this vector
introduction potential. soJufiun of A.48) wilh ().49) yields the expression

J(r'}x V<7([Link]\" (LSI)

The derivation of <LSI) can be found in most introductory electromagnetic texts

[2, 3. 5].and similar expressions are possible for two-dimensional currents where the

integration over a surface rather

is taken than a volume.
Integral equations may be formed for magnetostatics in a simitar manner to
electrostatics and the interested reader is referred to [9].

1.4 SURFACE EQUIVALENCE

Surface equivalent currents are very useful in thv formulation and execution of a
numerical solution of Maxwell's equations. Their introduction can be readily justi-
justified in the context equivalenceprinciple,e.g..two sources thut produce
of the surface
the Mime field within a region are .said lo be equivalent within that region.
The surface equivalenceprincipleslatesthat the field exterior (or interior) to a
given (possibly fictitious) surface may be exactly representedby equivalent currents
placed on that surface and allowed lo radiate into the region external {or internal) to
that surface. For the exterior cast;, these equivalent currents are given in terms of the
total exterior (E. H) fields while the interior liekb\302\273are assumed to be ?cro (this is

Love's equivalence principle). The appropriate currents for representing the fields
exterior to the surface are given by
10 Fundamental Concepts \302\246
Chapter 1

n x H = J
A.52)

For the interior fields, the negative of A.52) are [Link] radiated fields due to these
equivalent currents are given by the integral expressions

E(r) = -J\302\261 V x G(R) \302\246

n' x E(r

\302\246
ti x H(r')dS' A.53)

H(r) = - \320\224 V x d(R) n' x

\302\246
H(r') dS'

fi' x E(r') dS' A.54)

where R = |r r'|,r and r'

\342\200\224
denote the observation and integration point, respec-
respectively, and h' is the outward directed unit normal at the point r'. The closedsurface
on which the equivalence theorem is applied is denoted by Sc. These geometrical
quantities are illustrated in Figs. L2 and 1.3. In A.53) and A.54),a dyadic Green's
function is required which at least satisfies the radiation condition A.29). When
J and M are radiating in free space, the dyadic Green's function is given in closed
form by

'II = A.55)

where / = xx + yy + zz is the unit dyad and the corresponding scalar Green's func-
function is given by

AnR
A.56)

Also,
V x G0(r. r') = -V x [lG(i(r,r')] = -VG0(r.r') x / A.57)

implying V x Go(r, r') \342\200\242

M(r') = -VC?0(r, r') x M(r').
When A.55)-A.57) are introduced into A.53) and A.54),and after the use of
common vector and dyadic identities, we obtain the representations

E,H E,H

Js=rixH
Sources
(a) (b)

Figure 1.3 Illustration I he application

ol\302\260 o( the surface equivalence principle.
Section 1.4 Surface
\302\246 Equivalence II

E(r) = 11 \320\223M(r') x VGu(r. r') -,/fcDZoJ(r')C0(r. r')

A.58)

H(r) = \320\231
[-J(r')
x VGn(r. r') -i ^ M(r')&'0(r. r')

I
A.59)
'kitZa

More explicit expressionsfor E and H can be obtained by introducing the identities

Grfr. r') V* = \320\2410{

=^~ -(/ft,, +1)

in which R = (r- r'Vlr \342\200\224

r'|. as depicted in Fig. 1.2. Specifically. A.58) becomes

[M(r') x

'
-ji-\342\200\224
lj(r')
\302\253\320\276\320\273
(A-(l\302\253)-J

and a corresponding expression for H can be obtained by [Link] (E-\302\273H.

H -E.
\342\200\224 J -> M. M -> -J. ,u -> e, f -> /./, Zo \342\200\224 Ko -\302\273
K\302\253, Z,().
We cun rewrite A.58) in a. less singular form by noting the identities.
VG = -V'a,

A.61)

= VJ(r') - + J(r') \342\200\242

V(V
-
J(r'K/0(r.T')l VG',,(r, r') VVG,((r. r')
= J(rVWC,i(r.r') 1.62)

to deduce that

J(r') \342\200\242
WC/,,(r. r')dS' = -vli J(r') -
V'<7,,(r. r')(fS'

V .J(r')VGnU.r')clS'
12 Fundamental Concepts \302\246
Chapter I

Introducing these into A.58) yields the expression

E =\320\224 x VC0(r, r') -jk0Z0J(r')G0(r, r')

JM(r')

which is most commonly used for integral equation numerical solutions and is also
valid for open surfaces since the normal components of J to the perimeter edges
of the surface vanish. The correspondingH field expression is again obtained by

duality.
For far zone computations (r -* oo>, the Green's function A.56) can be sim-
simplified as

**'*> A.63)
^
Using this in A.58) and A.59), carrying out the vector derivative operations, and
retaining only the terms that decay1 as C?(I/r), we get

E(r) ft\302\273
>0
e-jL
\320\257
[+
f x M(r') + Zof x(rx J(r'))]e**'* dS' A.64)

H(r) *\320\2240
\320\246
\342\200\224^
[-? x J(r') + Yof x (r x M(r'))] ****'*dS' A.65)

Theseare referred to as the fai zone field expressions and are typically used for the
evaluation of antenna radiated fieldsor for the calculation of the radar cross section
(RCS)ofa target. An acceptable criterion for using A.64) and A.65) is

? \342\200\236.\3

where D is the largest antenna or target dimension. In this case, the phase error in the
intervening approximations is maintained at less than f.
A typical setup for comput-
the
computing radiation from volume sources at points in the near and far zones is depicted
in Fig. 1.2.
Figure 1.4 also shows the spherical angles commonly used in electro-
electromagnetics and to be used throughout this text.

Note that use of A.55) requiresboth equivalent currents A.52). Alternatively,

when J and M radiate in the presence of certain bodies such as an infinite metallic

plane, a cylinder or a half-plane, a dyadic Green's function can be chosen to satisfy

the boundary conditions appropriate for that surface, hence avoiding the need for
currents to be placed on that surface. For example, the Dirichlet condition

\320\271\321\205\320\261|=0 onS A.67)

can be imposed for relating the electric or magnetic fieldsin the exterior region to the
magnetic or electriccurrents,respectively. The radiated field expressions, A .Si) and
A.54), now simplify to

'The notation O( 1 /r) is road as \"order of 1 jr.\"

Section 1.4 Surface
\302\246 Equivalence 13

Figure 1.4 Illustration of an infinitesimal

source J = aJa dt and the associatedcoordi-
and
coordinates angles.

E(r) = -<\320\253> V x G{{R) \302\246

n' x E(r')dS' A.68)

H(r) = - \320\224 V x n x H(t')dS' A.69)

implying that only a single current is requiredfor the representation of each field
quantity.
One use of this Green's function is in the calculation of the scattering by a
perfect magnetic conductor (a fictitious material) through the use of equivalent
electric currents.
If the Neumann boundary condition is imposed

nxVxG2 =0 onS A.70)

the resulting field integral expressions are again given in terms of only one current

E(r) = -\320\254/ V2(R)\302\246 h' x H(r')dS' A.71)

H(r) = -jk0Y0<B G2(R) -n'x E(t')dS' A.72)

This Green'sfunction (G2) is useful for calculating the scattering by a perfect electric
conductor using electric currents. Also, this Green's function will be used to relate
the electric field quantities of an interior finite element formulation to the magnetic
field of the
bounding surface.
_ _
The above dyadic Green'sfunctions, Gi and G2, are commonly termed the first

and the second kind dyadic Green's functions, respectively. A good discussion of
dyadic Green's functions used in electromagnetics is given in [10].
The above expressions are for three-dimensional fields. In the case of two-
dimensional fields (e.g., one dimension, such as r, is invariant), similar expressions
are used. These are scalar and typically written in terms of TM. and \320\242\320\225.
polariz-
14 Fundamental Concepts \302\246
Chapter I

ations. For TM- (#- = Ex = Ey

= 0), the electric field on or outsidea contour\320\241
due

to fields on that surface is given by

E. = i t')]dl' A. 73)
~'*<>Z\302\260I H,(r')[G2n(r,
?-<r')|\"^7 G2\302\260(t-r<)]dV

since dEJ'dn= +jk0Z0H,. The corresponding expression for \320\242\320\225.

(?. = Hx =
JJr = 0) polarization is

H. = l G2D(r,r')\\dl' +jkoYoi A.74)

H:(r')\\?;
In A.73) and A.74), the subscript \"/\" denotes the tangential component of the field
along the unil vector i = I, where nxt = z. Also, the two-dimensional Green's
function is given by

G2D(r. r') = -J- H\342\204\242(ko\\r

-
r'|) A.75)

In A.75), #o2>(-)denotesthe zeroth-order Hankel function of the second kind. We

observe that A.73) and A.74) can be reduced from A.46) the potential
by replacing
with
\320\243 E. or H, and making use of the equivalent current relations A.52) and
Maxwell's equations to relate H. For the special
E and cases of the TMr and \320\242\32

polarizations, the first two of Maxwell's equationsimply

H = 2 x V?r, A.76)
-\321\204- \320\225=^!\321\205\320\243\320\257;

the expressions presented in this section are given for currents radiating
All of
in free Fields within
space. a homogeneous media can be determined by replacing Ao
and
with \320\272 Zo with Z in all of these expressions. Throughout this text, will denote
\320\272

the wavenumber in the homogeneous media (k = je^.k0) and Z is the intrinsic

impedance of that material (Z =

IS NATURAL BOUNDARY CONDITIONS

Maxwell's cannot be solved without

equations the specification of the required
boundary conditions interfaces. The pertinent
at material boundary conditions can
be derived directly from the integral form of Maxwell's equations. Specifically, A.25)
is applied to the contour illustrated in Fig. l.5(a) with 5 being the area enclosedby
\320\241
Assuming At is small, 0 and
-\302\273\302\246
\320\224\320\233 At \302\273 A.25)
\320\224\320\233, gives

(H,
- H,) \302\246
/ = [J,
\342\200\242
(n, x /)] \320\224\320\220 A.77)

and in deriving this we set

lim + e2E2] \342\200\242nl=0

^\320\224\320\233[\320\265,\320\225!

which is valid provided eE is finite at the interface. When A.26) is applied to the

same contour in Fig. l.5(a) we find that

- E2) \342\200\242
i = -[M, \342\200\242
x /)] \320\224\320\220
(E, 0?! A.78)
Section 1.5 Natural
\302\246 Boundary Conditions IS

(a)

Infinitesimal volume
enclosed by Sc

Medium B
(E2.H2)

(b)

Figure 1.5 Geometries Torderiving the boundary conditions (a) for tangential
components, and (b) Cor normal components.

where we have again set

lim ^- = \320\236
&h\342\200\224n

implying that fiH is finite at the dielectric interlace.

The conditions A.77) and A.78) can be rewritten in vector form and more
compactly by introducing the definitions

A.79)
= \320\234,-\320\224\320\271
\320\234\342\200\236 A.80)

giving the conditions

\320\273,
\321\205(\320\235,-H2)
= Jit A.81)
16 Fundamental Concepts \302\246
Chapter I

w, x(E, = -\320\234\320\271
-\320\225\320\273) A.82)

The quantities Jft and M,> are referred to as the impressed electric and magnetic
surface current densities in A/m and V/m, respectively, at the interface. Nole that if

E2 and H; are zero, these conditionsare identical to A.52) except that in this case Ju
and \320\234/,refer to actual impressed currents rather than equivalent currents.
To generate the boundary conditions correspondingto A.27) and A.28), we
select St, to be the surface of a small pill box, shown in Fig. 1.5(b), enclosing the
volume V. The pill box is positioned at the dielectric interface so that half of its
volume is in medium I and the other half in medium 2. It is again assumed that
0 so that only its
\342\200\224\302\273\342\200\242
\320\224\320\220 flat surfaces need be considered in performing the integra-
integrations. Through direct integration of A.27) we obtain the interface conditions

= pm A.83)
\321\217,-(<=|E, -e2E2) = A A.84)

where p, denotesthe unbounded electric surface charge density in C/m\" at the inter-
interface and pnvt is the corresponding fictitious surface magnetic charge density in

Wb/m2.
The boundary/interface conditions A.81HL84). although derived for time-
harmonic fields, are applicable for instantaneous fields as well. In the time-harmonic
case, only A.81) and A.82) are required in conjunction with A.23) and A.24) fora

unique solution of the fields.

If we ignore the fictitious magnetic currents and charges appearing in A.81)-

A.84), the boundary conditions are

x - H2) =
\320\231, (H, Jfc A.85)

x (E,
\320\270,
- E2) = 0 A.86)

\321\217,-(/ijH,-\320\2642\320\2352)
= 0 A.87)
=
\320\273, -\342\202\2542\320\2252)
.(\320\265,\320\225, \321\200, A.88)

The first two of these state that the tangential electric fields are continuous acrossthe
interface whereas the tangential discontinuousat the
magnetic fieldsare same loca-
location by an amount equal to the impressed electriccurrent. Unlessa source (i.e., free
charge) is actually
placed at the interface. Jfc is also zero and in that case, the
tangential magnetic fields will be continuous across the media as well.
When medium 2 is a perfect electric conductorthen E2 = H2 = 0. In addition,

Mix and p,m vanish and A.81)\342\200\224A.84) reduce to

=Jft
\320\233|\320\245\320\235, A.89)

/51 xE| =0 A.90)

= 0 A.91)
\320\231,
\342\200\242(/i,H,)

\320\273,
=
.(\320\265|\320\225,) \321\200\320\273 A.92)

The first two of these now imply that Ihe tangential electric field vanishes on the
surface of the perfect electric conductor whereas the langential magnetic field is

equal to the impressed electric surface current on the conductor.

Section 1.6 \302\246
Approximate Boundary Conditions 17

1.6 APPROXIMATE BOUNDARY CONDITIONS

In the previous section, the boundary conditions which must be imposed at the
interface of different dielectrics were [Link],it is difficult to utilize
these conditions sinceexcessive computational cost is required or the resulting for-

formulation is numerically unstable such as the case of a thin dielectric sheet. In many
cases, much simpler approximate boundary conditions that account for the presence
of an inhomogeneous medium, coaled metallic surface, or a thin dielectric layer can
be employed to simulate the actual surface. Below we discuss two types of such
approximate conditions:impedance boundary and sheet transition conditions. The
interested reader is directed to [II] for a general treatment of approximate condi-
conditions.

1.6.1 Impedance Boundary Conditions

The most widely applied approximate conditions are referred to as the

Standard Impedance Boundary Conditions (SIBC) or Leontovich Boundary condi-

conditions. It is derived by considering a plane wave impinging upon a material half-space.
Consider a material-air interface which corresponds to the = 0
\321\203 plane. The SIBC
takes the form

E. = -i)Z0Hx, Ex = rjZaH; A.93)

where the free space impedance is given by Z()
= vTW^o and the normalized ma-
material characteristic impedance is rj. An important concept to understand is that these
conditions are appliedslightly above the interface (assuming a plane wave originat-
in
originating the upper half-space) at =
\321\203 0+. Combining these two conditions, the vector
form of the SIBC is given by

\321\205\320\235 A.94)

where the outward directed unit in Fig. 1.6. This form

normal, is
\302\253,of the shown
SIBC is not restricted to a particular [as is the case interface
with A.93)] and is
commonly applied to convex surfaces such as a sphere, cylinder, etc.
All of the quantities used in A.94) are familiar and well defined except for the
normalized impedance, rj. One means of deriving this quantity is to demand that the

reflected field attributed to A.94) is identical to that due to the natural boundary
conditions. Then,
nrr

A.95)

This is exact for an infinite planar interface while it is approximate for a curved

boundary provided that

\\lm{yfcift\\kt,r,-\302\273l A.96)

where Im(-) denotes the imaginary part of the complex argument and the principle
radii of curvature, is associated
/\342\200\242\342\200\236 with the surface at a [Link] condition assures
that the material is sufficiently lossy so that the fields which penetrate into the
material does not re-emerge at some other point.
18 Fundamental Concepts \302\246
Chapter 1

Impedance
surface

(a)

Impedance
surface

(\320\253

Impedance
Dielectric surface
(e,M) \320\273
coating

Perfect conductor

(c)

Figure 1.6 Simulation of dielectric boundaries and coatings with SIBCs.

For a coated conductor,the choice of typically

\321\206 is found by considering a

shorted transmission line model with length corresponding to the coating thick-
thickness, t:

A.97)

However, this condition is derived at normal incidence and deteriorates at oblique

angles with increasing inaccuracy for thicker coatings.

Section 1.6 \302\246
Approximate Boundary Conditions 19

The SlBCs can be applied for modeling surfaces whose material properties vary
slowly in the transverse plane. For a planar interface, the coating can have a varying
composition in the normal dimension, and Rytov [12] found the following
impedance
I \320\257 1

is where N = y/]Z^~r
useful is the index of refraction and the normal derivative is
applied at the surface.
More accurate approximate conditionscan be developed by incorporating
higher order derivatives in their constructions. These are referred to as Generalized
Impedance Boundary Conditions (GIBCs), and these are discussedin [11].

1.6.2 Sheet Transition Conditions

A thin dielectric layer may be replaced with an equivalent sheet model to

simplify the analysis. Consider a thin dielectric layer with thickness r, as shown in
Fig. 1.7. This layer has conductivity a and will support a volume current density

J = <tE A.99)

where E is the electric field within the layer. For r A.,

\302\253: this volume current may be
replacedwith an equivalent sheet current (having units of A/m)

J., = rJ A.100)
and from A.99)

E=^ = Z0ReJ, A.101)

zcr

This condition is referred to as a resistive sheet transition condition which only

supports a single surface current, J,. The parameter /?,, is the normalized resistivity

of the sheet and is measured in Ohms per square.

The electric field was assumedto be tangential to the sheet in the above deriva-
derivation. A more general expression for the resistivesheet condition is given by
ft x (n xE) = -Z0/i,J, A.102)
Furthermore, it is desirable to utilize fields just outside or inside the sheet surface.
Sincethe tangential electric field is continuous across the sheet

n x x + =
[\320\231 (E+ \320\225-)] -22\320\223\342\200\236/\320\263\320\233

Figure 1.7 Simulation of thin dielectric layer ///////////^r.

with a sheet condition.
20 Fundamental Concepts \302\246
Chapter I

The superscripts denote

\302\261 the above and below the sheet, and the intro-
fields just
introduction of the second equation in is necessary
A.103) to maintain equivalencywith
A.102). The natural boundary conditions may be used to rewrite A.103) as

x [n
\302\253 x (E+ + E\]")= -2ZaR,.nx (H+
-
H)
, A.104)
n x (E+ - E~) = 0
As long as the loss in the layer is sufficient to assure that no multiple field penetra-
penetrationswill occur, these resistive transition conditions may be used for curved layers.
The dual to the resistive sheet condition is the conductive sheet condition which

supports only a surface magnetic current. It is given by

+ H\]") = 2 Y0Rmn x (E+

- E~)
\320\273\321\205[\320\271\321\205(\320\235+
(LI\302\2605>

The normalized conductivity of this sheet is denoted by Rm with units Mhos per
square. This condition is requiredfor the simulation of materials which have non-
trivial permeability. Also, a specialcombination of coincident resistive and conduc-
sheets
conductive with respective resistivity and conductivity

yields the as an impedance sheet with

same result impedance ij and A.106) implies
that 4ReRm This set of sheetsare
= 1. useful in simplifying the analysis of a planar
impedance sheet since coplanar resistive and conductive sheets are uncoupled.
The resistivity of a dielectriclayer can be determined by considering the equiva-
equivalentpolarization current

J=jk0Y0(er-\342\204\226 A.1071
and A.100). It follows that the tangential components of the field are given by

E, = Z0ReJ,, A.108)
with

R,.= . ^ n A.109)
koz(er- I)
A dual conductive sheet is given by

More accurate representations can be formed using Generalized Sheet Transition

Conditions (GSTCs) [11] which incorporate higher order derivatives in their con-
construction.

1.7 POYNTING'S THEOREM

The quantity (\"*\" indicates complex conjugation)

Section 1.7 \302\246
Poynling's Theorem 21

is known as the complex Poynting vector and has units of Watts/m2. It represents the
complex power density of the wave, and it is therefore important to understand the
source and nature of this power. To do so, we refer to A.21) and A.22), where by

dotting each equation with E and H\",we have

E V x H* = j; E -/WE* E =
\302\246
j; \342\200\242
E -./W|E|2 A-\320\2302)

H* \342\200\242
V x E = -M/ \342\200\242
H\" -jcoixH
\302\246
H* = -M' \342\200\242
H* ->\320\274|\320\235|2 A.113)

From the vector identity [6]

V.(ExH) = H\"VxE-E-VxH A.114)

we then obtain

V \342\200\242
(E x H*) =yW|E|2 ->m|H|2 - j; E - M,
\302\246
H* A.115)

which is an identity valid everywhere in space. Integrating both sides of this over a
volume V containing all sources, and invoking the divergence theorem yields

A.116)

which is commonly referred to as Poynting's theorem. Since Sc is closed, based on

energy conservation one deduces that the right hand side of A.116) must represent
the sum the power stored or radiated,
of escaping, i.e.. out of the volume V. Each
term of the volume integral of A.116) is associatedwith a specific type of power but
before proceeding with their identification, it is instructive that e* be first replaced by
eoer+Jw- Equation A.116) can then be rewritten as

(E x \320\251 ds =
\342\200\242
Pei + Pmi -Pd A.117)
s,
X- 1\321\202<\320\246
(E x H*) \302\246
ds
-
lo^W, - Wm\\
-i Im f f [j;
\342\200\242
E + Mr H*]dv (I.I 18)
f

where
1 f f f
Pej
\342\200\224
\342\200\224r\\ Re(J*
\342\200\242
E) dv \342\200\224
averaging outgoing power due to
J J Jv current J
the impressed A.119)

^mi
\342\200\224
Re(M,
\342\200\242
H*) dv \342\200\224
average outgoing power due to the
~~Z)\\\\\\
\342\200\242'\302\246'\342\200\242'v
impressed current M, A.120)

= - <r|E|2rfi)= in \320\232 A-121)

Pfl average power dissipated
2 v

Wc =
- I I I ener|E|2 dv = average electric energy stored in \320\243 A.122)

Wm =2 /^(i^rlHI2^ = average magnetic energy stored in V A.123)

22 Fundamental Concepts \302\246
Chapter 1

The time-averaged power delivered to the electromagnetic field outside V is clearly

the sum of />,., and Pml, whereas Plt is that dissipated in V due to conductor losses.
Thus, we may consider

(ExH')-rfe A.124)

to be the average or radiated power outside V if a is zero in V. Expression A.118)

gives the reactive power, i.e., that which is stored within V and is not allowed to
escape outside the boundary of S,..

1.8 UNIQUENESS THEOREM

Whenever one pursues a solution to a set of equationsit is important to know a

priori whether this solution is unique and if not, what are the required conditions for
a unique solution. This is important because depending on the application, different
analytical or numerical methods will likely be used for the solution of Maxwell's
[Link] that Maxwell's equations (subject to the appropriate boundary
conditions) yield a unique solution, one is then comforted to know that any con-

convenient method of analysis will yield the correct solution to the problem.
The most common form of the uniqueness theorem is: In a region \320\243
completely
occupied with dissipative media, a harmonic field (E. H) is uniquely determined hy the
impressed currents in that region plus the tangential components of the electric or
magnetic fields on the closed surface Sc bounding V. This theorem may be proved

by assuming for the moment that two solutions exist, denoted by (E],H|) and
(Ei, Hi). Both fields must, of course,satisfy Maxwell's equations A.21) and A.22)
with the same impressed currents (J,, M,). We have

V x H, = J, +>eE,. V x H2 = J, +jaxE2
V x Ei = -M/ \342\200\224j<afiHlt V x E2 = -M, -jto(iH2

and when these are subtracted we obtain

A.126)
A.127)

where E' = E| H' = H|

\342\200\224
E2 and
H2. To prove the theorem it is then necessary
\342\200\224
to
show that
(E', H') are zero or equivalently. if no sources are enclosed by a volume V,
the fields in that volume are zero for a given set of tangential electric and magnetic
fields on Sc.
As a corollary to the
uniqueness theorem, it can be shown that if a harmonic

field has a zero tangential or magnetic field on a surface enclosing

electric a source-free

region V occupied by dissipative media, the field vanishes everywhere within V.

The usual proof of the uniqueness theorem can be found in many electromag-
texts
electromagnetics (see, for example, [5]).
Section 1.10 \302\246
Duality Theorem 23

1.9 SUPERPOSITION
THEOREM

The superposition theorem states that for a linear medium, the total field intensity

due to two or more sources is equal intensities attributed to

to the sum of the field
each individual source radiating independent of the others. In particular, let us
consider two electric sources J, and J2. On the basis of the superposition theorem,
to find the total field caused by the simultaneous presence of both sources, we can
considerthe field due to each individual source in isolation. The fields (E|, H|) due to

Ji satisfy the equations

VxH, = J, +./WE, A.128)

A.129)

and the fields corresponding to J2 satisfy

A.130)
A.131)
By adding these two sets of equations, it is clear that the total field due to both
sources combinedis given by

E=E,+E2, H = H, +H, A.132)

where Hi)
(\320\221;,
and (Ej, are
\342\204\226>) obtained by solving separately A.128\320\235\320\25329) and

A.1 1.131),
\320\227\320\236\320\235 respectively.

1.10 DUALITY THEOREM

The duality theorem relates to the interchangeability of the electric and magnetic
fields, currents, charges, or material properties. We observe from A.3) and A.4) that
the first can be obtained from the second via the interchanges

M->- -J

HE:HE

Similarly, A.4) can be obtained from A.3) via the interchanges

J-> M

The duality theorem can reduce formulation and computational effort when

one is able to invoke it for a particular application.

24 Fundamental Concepts \302\246
Chapter 1

1.11 NUMERICALTECHNIQUES

For numerical solutions, all governing equationscan be written in operator form as

?u-f =0 A.135)
subject to appropriate boundary or transition conditions

B(u) = 0 A.136)
within the domain (fi) and on its boundary (Sr = dQ). In these, the operator ? is
based on oneof the following: an integral representation of the fieldssuch as (I.52)-
A.53), on the vector wave equation A.30), or the Helmholtz equation A.36) for
scalar fields. It is understood that \320\270
must be replaced by a vector field u when dealing
with the vector wave equationsA.30) or A.35). The forcing function/ is a known

excitation function while or u

\320\270 is the unknown quantity. Throughout this text, or u
\320\270

will denote a field or current density.

Unfortunately, very few analytical solutions for A.135) are available in elec-
electromagnetics. One such solution, the fields due to a magnetic dipole in the presence
of an infinite metallic plane or cylinder, will be used in Chapter 7 to form the
appropriate dyadic Green's function for those [Link], most useful
electromagnetic scattering and radiation problemscannotbe solved using analytical
methods. Rather, an approximate numerical solution is sought which in some way
closely resembles the exact solution. Two methods of formulating such an approxi-
approximatesolution are: the Ritz method and the method of weighted residuals.

1.11.1The Rib Method

The Ritz or Rayleigh-Ritz method2 [13, pp. 74-78],[14,pp. 13-63] seeks a

stationary point of a variational functional. For operators which are self-adjoint and
positive-definite (see later subsection for definitions), the stationary point of the

following functional

F(u) = \\{?u,u)-(f.u) A.137)

is an approximate solution of A.135). In A.137), the inner product over the domain \320\277
(volume, surface, or contour) of the two functions is defined as

{a.b)= ahdQ. A.138)

\\
Jn

or for vector functions

<a.b)= f
abdu A.39)
.In

The choiceof this inner product extends the validity of the varialional expressions to
vectorial fields. When the operator Cit and/ in A.137) are chosen as

2The method was originally introduced by Ritylcigh in 1877 and was extended by Ritz in 1909.
Section 1.11 Numerical
\302\246 Techniques 25

i = A.140)
Vxl-l-=-=|-*fcru
Mr J
/\\\321\217\\

A.141)

it can be shown that setting the first variation of F(u) to zero is equivalent to

satisfying the vector wave equation A.30) over the computational domain ?2.
Similarly, when
A.142)
setting the first variation of F(u) to zero is equivalent to satisfying the inhomo-
geneous Helmholtz wave equation

V2U + k20u=/ A.143)

To show that setting the first variation of F(u) to zero is equivalent to satisfying
the Helmholtz equation A.143), we begin by rewriting the functional (we assume a
two-dimensional domain)

+ klu]udn- f f fudu A.144)

n[V2tl
q J Jn
as

F(u) = \320\250 (-V\302\253

\342\200\242
Vm + k2ou2)dQ
+\\
\320\257
\302\253|j
dl-U fudQ A.145)

in which Vm \342\200\242
=
\321\217 denotes
\320\241
\320\254\320\270/\320\264\320\277, the contour enclosing the region ?2 (see Fig. 1.8)
and h is the unit normal vector to Note
\320\241 that in deriving A.145) we used the
identity

=
V\302\253) -Vm \342\200\242
V^ + V \342\200\242
(ifVu) A.146)

and the divergence theorem

=| Vu-nds A.147)
a ic
Next we proceed to evaluate the first variation of F(u) given by

SF = -
F(u + \320\220\320\270) F(u) A.148)

Figure 1.8 Illustration of (he region \320\277

and
the enclosing contour \320\241
26 Fundamental Concepts \302\246
Chapter 1

where 0 is
-*\302\246
\320\224 a scalar quantity. The evaluation of SF involves the quantities

= \320\2702
+ \320\220\320\270I
(\320\274 + (\320\224\320\270J
+ 2(\320\224\320\270)\320\275 A.149)

[Vh +
\342\200\242
V( \320\224\320\270)]
[Vm + V( Aw)]
= Vm \342\200\242
Vm + 2V( \342\200\242 \342\200\242
Vw + V(\320\224\320\270)
\320\224\320\274) V(\320\233\320\270)A150)

\320\264\320\270
3(\320\224\320\274) \320\264\320\270 \320\255(\320\224\321\213) .
\320\264 \320\220

These can be simplified by neglecting the last tenn of each expansionwhich is of

order \320\2242.
Doing so yields the approximation

F(u + % \342\200\242
Vm + klu2] dU - fudU
\\ f
\320\224\320\270) [-V\302\253 f [
2JJq JJq
9m
\342\200\236

\302\246
Vm
+ \320\224 f [-mVm + katr]dQ
I

When this expression is comparedwith that for F(u) in A.145), we have

F(u + \321\212
F(u)
\320\220\320\270) + A
[ [
m[V
\342\200\242
Vk + k\\u -f] </?2

-\320\233 \320\270*\320\250 +
+ A.153)
Jc dn 2 Jcl
\302\2611\\\320\270? \320\270?\\\320\260
\320\255\302\273
BnJ

where we also used the divergencetheorem A.147)and the identity A.146) to obtain

the second and third terms. Clearly, the last two terms in A.153) cancel each other
leading to

SF = F{u + -
\320\224\320\260)F(u)

= \320\224
f f h[V2m + klu -JidQ A.154)

Thus, setting SF = 0 implies that V2m + klu =f provided \320\224\320\274

is nonzero. That is,
from A.154) we conclude that the extremization of F obtained by setting
8F
SF = 0 or -r- = 0
\320\260\320\270

is equivalent to enforcing the Helmholtz wave equation over the domain In

\320\257.

practice, the condition SF \342\200\224

0 is enforced by setting
-
dF
~
F{u + \320\220\320\270)F(u) = 0 A.155)
\320\264\320\270 \320\220\320\270 \320\224\320\274
\320\273_\320\277 \320\273\342\200\224\320\276

i.e.. by settingto zero the derivative of the functional with respect to u.

Having established the equivalence between SF = 0 and the Helmholtz equa-
equation we can proceed with the discretizationof F(u) and SF to obtain a discrete system
of equations. The discretization begins with the trial function, ii, expanded in terms

of N basis functions
Section 1.11 Numerical
\302\246 Techniques 27

= <\320\230\302\273\320\223<\320\235'} A.156)
\320\233\">
./=1

where iv, are the basis functions and are

\302\253,
the unknown expansion coefficients. In
A.156), column data vectors are denoted with while
(\342\200\242} row data vectors involve a
transposition {}r. Substituting A.156) into A.137), the functional becomes

[w)fda A.157)
JJ
where we used the innerproduct definition

=
}. (\302\253}) [u)T{v]

for discrete data vectors. This functional is extremized by allowing all partial deri-
derivatives with respect to the coefficients, to
{\320\274}, vanish

A.158)

A single equation is obtained by differentiating with respect to each For

\321\211.

/=1,2 , N we obtain N equationswhich can be written as a matrix system

[A][u) = [b\\ A.159)

The elements of the matrix [A] and excitation vector \\b\\ are given by

A.160)
\\\\
b, = w,fdQ

A word of caution: electromagnetics differs from other branchesof engineering in

that no physical significance can be attached to the stationary point of the functional
A.137). In mechanical systems, for example, minimizing this functional represents
minimization of the total potential energy of the system. However, since electromag-
involves
electromagnetics complex quantities, such a statement may not be asserted.

1.11.2 Functionate for Anlsotropic Media

In three dimensionswith anisotropic media [15, 16] the appropriate operator is

of the form

U = E
?(u)=(Vx(^ \302\246Vx\302\273)-|1\302\260' A.161)

with the corresponding source function given by

- x \321\204.~[\342\200\242 u = E
I ->^oJ
V M),
f = A ]62)
+ Vx(?;'.J), u = H
I ->e0M
28 Fundamental Concepts \302\246
Chapter 1

The associated functional to be extremizedis then of the form (for u = H)

F(H) = lf [VxH-f;'-VxH-^H-^H]f/r- f
H(dV A.163)

where

\320\241\321\203\321\205
\342\202\254\320\243\320\243

are the permittivity and permeability tensors of the media and \320\257
represents a
volume. In general, for arbitrary anisotropy, this functional will lead to an asym-
asymmetric (non-Hermitian) system. One way to obtain a symmetric system is to use the
functional

= <?u, ua)
- (u. fa) -
(u0, f) A.165)
where ua and f(l satisfy the partial differential equation

A,ua =./\342\200\236 A.166)

in which Ca is the adjoint operator to ?. That is,

)= <v,?eii> A.167)

1.11.3Method of Weighted Residuals

The method of weighted residuals [17], [18] begins with the residual

A.168)

and seeks a solution for = \320\270

\320\271 by satisfying the condition H = 0 within the domain
In
\320\277. general, such a solution cannot be achieved at all points in \320\257.
Instead, it is
more practical to find a solution which satisfies the residual condition in some

average or weighted sense over jV subdomains of ?2, viz.

[t,C[w)r{u)
-
tif] d?l=0, i= 1.2,3 N A.169)
f
Jo

In general any testing function, r, (also referred to as trial or weighting functions),

may be used; however, since these functions modify the enforcement of the boundary
conditions throughout the domain, the choice of testing functions affects the quality
of the solution A.169). One popular testing procedure is called collocation or point
matching. In the testing or weighting
this, function is a Dirac delta function,
= S(.x implies enforcement of the boundary conditionsonly
\342\200\224which at discrete
t/ .X/).

points (e.g., xh i = 1,2,3 Another

\320\2330- popular choice is termed Gaterkin'spro-
When
Gaterkin'sprocedure. employing the Galerkin's testing procedure, the testing function is

identical to the expansion function used in A.156), e.g., /,

= 117 and the weighted
residual equation is given by

I
f
w,C\\Vjdn\\ {u)=\\ w,f d& A.170)
Section 1.11 Numerical
\302\246 Techniques 29

which is identical to the Ritz procedure given above. Thus. Galerkin'smethod leads
to the same linear system A.159)as the Ritz method.
As a generalization, when F(u) is chosen as

= i
F(M) (?\302\273.\302\253*)-(/,\302\253*) A.171)

where the \"*\"

indicates complex conjugation, the extremization of F{u) leadsto a
linear system that is identical to that obtained from the weighted residual method
with /,= w*.

1.11.4 Vector and Matrix Norms In Unear Space

A norm is a real valued function that of the size or \"length\"

provides a measure
of a multicomponent mathematical quantity
such as a vector or a matrix. It is usedin
numerical analysis to provide a measure of how well a given vector approximates the
exact [Link] matrices, norms provide a single value to quantify the \"size\" of

the matrix [A].They are often used in evaluating the numerical system's condition,
which in turn affects the stability of the solution. That is, of interest is how a small
change in the excitation or right hand side of the matrix system (I.I59) affects the

solution data vector.

[Link] Vector Norms

Euclidean Norm. The most popular form for a given discrete data vector

=
{\320\270} {\302\253|,Mi-U3,...,uN)T is the Euclidean norm. It is defined by3

llulh = ,J = !\"}>
((\302\253}.
= MT[u) A.172)

for real valued [u] and by

112
= = }r[u*
{u}r[u*\\ A.173)

for complex vector (\320\270).Here \\ut\\ implies the absolute value of the quantity.
Throughout the book, the notation will
f|\302\253|| imply the Euclidean norm of a vector
or data column unless otherwise noted.
Infinity Norm, The infinity norm of a data vector is defined by
=
||\320\270||,*, max of |m,| \\<i<N A.174)
This norm is also referred to as the uniform vector norm or maximum magnitude norm.
H6lder Norm. The Holder or p-norm is a generalization of the Euclidean

norm and is defined as

| {>p

(\320\25375)
X>
In
where denotes
|\321\213/[/* the /?th power of the quantity \\u,\\.

JHere u is a simpler notation for |\302\253|.

Fundamental Concepts \302\246
Chapter 1

[Link] Matrix Norms

Frobenius Norm. The Frobenius matrix norm is a generalization of the

Euclidean vector norm. For a square matrix [A], it is given by

\\\\A\\\\F
= (LI76)

the (l,j) of the [A] matrix.4 Note

where Al} denotes entry that A.176) can be general-
generalizedto nonsquare matrices and to the/?-matrix norms.

Matrix Infinity Norm. The infinity or uniform matrix norm is denned by

= \321\202\320\272
\320\276\320\223 (\320\25377)
\\Y,\\Aii\\\\

This specific norm is also referred to as the row-sum norm. Similarly the column-sum
matrix norm is given by

\\\\A\\\\X
= max of f
53\\\320\220-\320\233
A.178)

The infinity norm is referred to as the natural norm of [A]. It can be shown that

\\\\A\\\\X
= max(||M](u}||) = max of 0-179)
\\Y,\\*A

For Hermitian matrices A.178) is identical to A.177).

Matrix Condition Number. The condition of a matrix is related to the

natural norm of [A] as

Cond(^) = \\\\A\\\\\\\\A-11|
>
?dh1 A.180)

where \\\\A~l1| refers to the natural norm of the

of [A] (the Frobeiiius norm is inverse

not a natural norm). Here [Link] the maximum and minimum

and Amjn denoteeigen-
eigenvalues of the matrix, respectively. Since |[Link]| is a lower bound for the natural norm
of [A], the ratio |Amax|/|[Link]| gives a conservative estimate for the condition of [A].
As an example of the importance of Cond(.A),let us assume that due to
truncation or arithmetic errors. [A] is instead approximated by [\320\220\320\264].
A computer

with / decimal digits of accuracy gives

\\\\A-AA\\\\ .
= 10\"
1\320\230\320\230

''Here we use the notation \320\246.\320\224\320\246

= IIMlII, i.e.. the calligraphic capital letter implies a matrix. Such a

notation will be used later in Chapter 9.

Section [Link] Numerical
\302\246 Techniques 31

which is a measureof the normalized error in approximating [A]. If the condition of

the matrix is \321\201
= Cond(.4), the corresponding (normalized) error in the computed
solution [u] will then be [19, 20]

||u \342\200\224
uAll
,0-,

where

That is, s is always /, implying a larger error for the final solution
smaller than
vector. More specifically, in the seventh decimal place for the norm
an error of [A]
(i.e.. / = 7) translates to an error in the third decimal place for (the norm of) the

solution vector when the matrix condition number is 104. Alternatively, if

CondM) = 1, the solution vector error is of the sameorder the same decimal (i.e.,to
place) as the matrix error itself.

1.11.5SomeMatrix Definitions

1.1 we give the mathematical and descriptive definitions

In Table of matrices

often encountered in numerical analysis.

Additional definitions can be found in several numerical mathematics books
(e.g., see [19]and [20]). As can be understood, the operators which generated the
matrices in Table 1.1 also carry the same definition. That is, an operator is referred
to as Hermitian or self-adjoint if the resulting matrix is also Hermitian.

TABLE t.l Definitions of Matrices Often Encountered in

Numerical Analysis

Mathematical Statement Descriptive Statement

[A]r = [A] Symmetric

[A] = [A'f Herraiiian (self-adjoint)
\\A) = -]A*f Skew Hermitian
= \\A*\\T
[\320\233]\021 Unitary
[Ar][A] = [J] Orthogonal
[/] = identity matrix
[A')T[A] = [A][A'}' Normal
Ay
> 0 Positive-definite
A/j
> 0 Nonncgativc or positive scmidclinitc
= 0 for i ?\321\214
\320\220 / Diagonal
\320\273
= 0 for / >j Upper triangular
Asj
A,j
= 0 for i >j Strictly upper triangular

A,, > Y^,\\Ai,\\ for all i Diagonally dominant

> 0.
|u| 0
|\320\270)\320\263[\320\233]|\321\213|
\321\204 Positive-definite
>
0.
(u| 0
|\320\270}\320\263\320\230]{\320\270|
\320\244 Nonnegaiive
7'rMI\302\253)) <0 Indefinite
32 Fundamental Concepts \302\246
Chapter

TABLE 1.2 Some Common Relationships Between Operators

and Eigenvalues

Operator Type Eigenvalue (X.) Properties

Hcrmitian Real eigenvalues and >0

Unitary Eigenvalues on unit circle

Skew Hcrmiiian Eigenvalues on imaginary axis
Positive semidefinite Eigenvalues >0
Positive-definite Eigenvalues >0
Indefinite Sonic eigenvalues are >0 and some arc <0

Given the natureof the malrix or operator, we can immediately make a state-
statement about the eigenvalues of that operator. Some of the most common relation-
relationshipsbetween operators and eigenvalues are given in Table 1.2.

1.11.6 Comparison of SolutionMethodsand Their Convergence

The Rayleigh-Ritz and Galerkin's methods are standard solution approaches

for solving differential equations arising in practical engineering problems. Both
methods project a continuous space onto a finite separable Hilbert space.5 The
mathematical problem is then rephrased to seek a discrete solution set whoseentries
are the coefficients of the expansion. The premiseof each method is summarized in
Table 1.3. The third entry in the table is referred to as the Least Squares Method and
is generally more robust than the other two, as will be noted later.
It was shown in Section l.H.3 that Galerkin's method is equivalent to the
Rayleigh-Ritz approach when the first variation is set to zero. At the heart of this
statement is the assumption that the operator \320\241
(or the resultant discrete matrix \320\233) is

positive-definite. Unfortunately, in most practical problems in electromagnetics,

particularly as k0 becomes larger, the operators

= Vh2 - klu and
\320\241\320\270 = V x
\320\241\320\270

(/li^'v x u)
\342\200\224
kleru do not guarantee positive-definiteness. That is, if the operator
is not positive-definite, the Rayleigh-Ritz method fails to ensure minimization of the
functional since a global stationary point may not exist. However, the application of
Galerkin's method to yield a discrete system does not require that the operator is

positive-definite or even symmetric.6 The resultant solution is simply a stationary

point which is not guaranteed to be a minimum.

TABLE 1.3 Mathematical Statement of Solution Methods

Method Mathematical Statement (subject to boundary conditions)

minimize -
Rayleigh Ritz \320\270)(/'.\302\253)
\\ {\320\241\320\270,
Galerkio solve - (/.
(\302\253\320\241\320\270. = 0
ttj) \302\253/>
\342\200\224
Least Squares minimize Q{u) = {\320\241\320\270
\320\241\320\270
\342\200\224f. f)

'Hilbert space refers to a linear space where a given intcrproduei has been defined and which is
complete with respect to this interproduct.
6Because of the complex t, and y.,, in electromagnetics, the operators may be symmetric but not
Hcrmitian (i.e., self-adjoint).From Section 1.11.5, Hermitian operators have positiveand real eigenvalues
and are therefore positive.
Section 1.11 Numerical
\302\246 Techniques 33

In some cases, the problem statement is that of minimizing a functional and

consequently the Rayleigh-Ritz procedure is the natural method of choice. Examples
include problems associated with system energy minimizations and resonance
computations. However, in solving the Helmholtz or vector wave equations (subject
to given boundary conditions), the minimization of the functional must simul-
simultaneously imply a solution of those equations. For those cases, Galerkin's method
is the appropriate approach for constructing the linear system. However, introduc-
of
introduction an appropriate variational functional can simplify the problem statement
when dealing with boundary conditions other than Neumann-type (also referred
to as natural conditions). As seen from the derivation given in Section 1.11.1, the
presence of the boundary integral term provides a direct means for imposing the
Dirichlet, impedance, or other type of boundary conditions. In the case of Galerkin's
method, the boundary are introduced after application of the divergence
terms
theorem. FinaJJy, we note some cases a variattonaJ (minimization)
that in statement
of the boundary value problem may not be possible. In those cases, Galerkin's
method is the only approach for constructing a linear system of equations.
When dealing with operators that are not positive-definite, apart from the
breakdown of the functional minimization process, most iterative linear system
solvers also break down. Specifically, convergence of the approximate solution \320\27
to the exact solution w cannot be proven [21], [22] without invoking positive-definite-
ness. When in doubt, a positive-definite operator (and a correspondingpositive-
definite matrix) can be generated by instead solving the differential equation

CaCu = C\"f A.182)

where satisfies
\320\241

v) =
{\320\241\320\270, {u, Cv) A.183)
and is referred to as the adjoint operator. Clearly, A.182)is obtained by multiplying
the right and left hand =f by C. The new
sides of \320\241\320\270 operator V, where

Vu = g A.184)
(g = Cf) is now positive-definite and self-adjoint. The corresponding matrix system
resulting from A.184) is of the form

) A.185)
or
[B]M = {g}
where [B] It is thus seen that the desired property of positive-definiteness
= [A*]T[A].
comes at the
price of squaring the matrix condition number. As is well known, large
matrix conditions lead to less accurate solutions and slowerconvergence when an
iterative solver is used.
It should be remarked that minimization of the functional for the Least
Squares Method is equivalent to solving the differential equation A.182).
Consequently,the Least Squares Method leads to positive-definitesystemsat the
expense of squaring the matrix condition. Also,the Least Squares Method minimizes
the square of the norm (as -*\302\246
\320\277 u), viz.
Fundamental Concepts \302\246
Chapter I

Jim ||?a-/||2->0 A.186)

whereas Galerkin's method minimizes the norm

lim -*0 A.187)

N-*oo ||?\320\271-/'||

Again, it is important to note that nothing can be said about convergenceunless the

operator is positive-definite.

1.11.7Field Formulation Issues

The finite element method can assumevarious forms depending on the desired
field quantity. Many applications prefer either a total or secondary electric field
formulation. Other applications desire a result in terms of either the total or second-
magnetic
secondary field. Some applications can utilize a potential formulation. Thus, even
though Maxwell's equations relate these various quantities, an accurate field com-
computation often demands a particular formulation. The advantages and disadvantages
of each of these formulations are discussed below.
The total electric field formulation very popular choice. This is because
is a
enforcement of the
boundary conditions associated with perfect electric conductors
(pec) is particularly easy. Since the tangential electric fields on a pec surface must

vanish, the edges of the mesh associatedwith those surfaces are a priori set to zero.
Three methods are commonly used in practice to enforce this condition. The first is
accomplished by forcing a null field condition to zero out all entries of the matrix
associated with that edge (except for the self-term which is set to unity), and by also
setting the excitation entry to zero. Thus, as the unknown fields are solved, the edges
lying on pec surfaces are forced to zero. The second method involvesa preprocessing
step where the edges associated with a pec surface are removed from the list of
unknowns. Thus, the number of edges greateris than the number of unknowns
and matrix entries for these pec edgesare never computed. This approach has the
advantage of reducing the order of the matrix and therefore reducing memory and
compute cycle demands. The third method, useful when an iterative matrix solver is

employed, involves forcing the unknowns associated with the pec edges to zero

during each iteration.

Thus, use of a total electric field formulation and edge-based elements(see
Chapter 2 for a discussion on edge-elements) reduces the order of the matrix and

computational burden. However, this is not the only valid formulation and in certain

circumstances, a scattered field or a magnetic field formulation may be preferred

Scattered (or secondary) field formulations are usedto simplify the use of absorbing

boundary conditions (see Chapters 4 and 6). However, they also have an added
advantage when a boundary integral is used for mesh closure. Experiencehas
shown phase errors in the computed
that interior field tend to increase within the

mesh locations
at distant from boundaries on which boundary conditions arc

imposed. This is due to unavoidable numerical inaccuracies t hat increase as the effect

of the boundary conditions propagate throughout the mesh. That is, previous errors
in the adjoining field are incorporated and magnified as the field is evaluated at a

more distant field point. Since boundary conditionsalways are enforced with total

fields, the total field formulation enforces such conditionsonly on the boundaries of
References 35

the mesh and pec surfaces (for E-field formulations). For very large computational
domains,significant distance can lie between a field point within the interior of the
mesh and the mesh [Link],the potential for error propagation throughout
the mesh. A scattered field formulation enforces the boundary conditions on the
mesh boundary, pec surfaces, and dissimilar material interfaces. Therefore, the dis-
distance between a boundary condition and any interior field point is reduced and

accordingly the phase error throughout the mesh may be also [Link] scattered
field formulation has the disadvantage of higher matrix order (i.e., more unknowns
and equations) and explicit enforcement of the boundary conditions associated with
pec surfaces and material discontinuities.

Magnetic field formulations are also possibleand can be obtained by applying

duality to the corresponding electric field equations. A total magnetic field formula-
has
formulation the same advantages for a fictitious perfect magnetic conductor (pmc) as the
total electric field has for pec surfaces. A scattered magnetic field formulation also
can reduce phase error propagationin the same manner as the scattered electricfield
formulation.

Magnetic field formulations are preferred for applications where the desired
result is the magnetic field within the computational domain. This is due to the fact
that although Maxwell's equations relate the electric and magnetic fields, in practice
one quantity cannot be accurately obtained from the other by numerical differentia-
This
differentiation. is due error occurring when
to the inherent continuous derivatives are
replaced with discrete Rather, a computationally expensive integral
differences.
expression is necessary for accurate field differentiation provided a suitable
Green's function is available. Hence, an accurate solution demandsa formulation
consistent with the desired result.
Finally, some finite element practitioners utilize a potential formulation which
employs the scalar or vector potentials as the unknown quantities (see Chapter 5).
The use of this approach is related to the hybrid finite element-boundary integral
method where the singularity of the integral equation associated with the boundary
can be reduced with a potential formulation. However, if this reduction is present,
we note that a numerical differentiation operation may be requiredto obtain the
desired field quantity and this operation may lead to inaccuracy.

REFERENCES

[1] J. A. Stratton. Electromagnetic Theory. McGraw-Hill,New York. 1941.

[2] R. E. Collin. Field Theory of Guided Waves. IEEE Press, New York, 1991.
[3] C. A. Balanis. Advanced Engineering Electromagnetics. McGraw-Hill,New
York, 1989.

[4] D. S. Jones. Methodsin Electromagnetic Wave Propagation. IEEE Press, New

York, 1994.
[5] R. F. Harrington. Time-Harmonic Electromagnetic Fields. McGraw-Hill,New
York. 1961.

[6] \320\241. Tai.

\320\242. Generalized Vector and Dyadic Analysis. IEEE Press, New York,
1992.
36 Fundamental Concepts \302\246
Chapter I

[7] J. A. Kong. Theoryof Electromagnetic Waves. Wiley InterScience, New York,

1975.
[8] J. D. Jackson.
Classical Electrodynamics. Wiley InterScience, New York, 1975.
[9] J. van Bladel. Electromagnetic Fields. HemispherePublishing Corp.,New York,
1985.

[10] \320\241. Tai.

\320\242. Dyadic Green's Functions in Electromagnetic Theory. IEEE Press,
New York. 1994.
[11]\320\242.\320\222.A Senior and J. L. Volakis. Approximate Boundary Conditions in
Electromagnetics. IEE Press, London, 1995.
[12]S. M. Rytov. Computation of the skin effect by the perturbation method.
/. Exp. [Link]
Theor. Phys., 10:180-189, by V. Kerdemelidis and
\320\232. Mitzner,
\320\234. Northrop Navair, Hawthorne. CA 90250.
[13] S. G. Mikhlin. Variational Methods in Mathematical Phvsics. Macmillan, New
York, 1964.

[14] J. N. Reddy. An Introduction to the Finite Element Method. McGraw-Hill,New

York, 1984.

[15] A. Konrad. Vector variational formulation of electromagnetic fields in aniso-

tropic media. IEEE Trans. Microwave Theory Tech., MTT-24:553-559,

September 1976.
[16] \320\241H. Chen and \320\241D. Lieu. The varialional principle for non-self-adjoint
electromagnetic problems. IEEE Trans. Microwave Theorv Tech.. MTT-
28:878-886, August 1980.
[17] R. F. [Link] Computation by Moment Methods. Macmillan. New

York, 1968.

[18] J. J. H. Wang. Generalized Moment Methods in Electromagnetics. John Wiley &

Sons, New York, 1991.

[19] R. L. Burden and J. D. Faires. Numerical Analysis. PWS Pub. Co., Boston, fifth

edition, 1993.

[20] G. H. Golub and F.

\320\241 Van Loan. Matrix Computations. Johns HopkinsUniv.
Press, Baltimore. MD, 1983.
[21] G. Strang and G. J. Fix. An Analysis of the Finite Element Method. Prentice
Hall, Inc., EnglewoodCliffs. NJ. 1973.

[22] D. G. Dudley. Mathematical Foundations for Electromagnetic Theory. IEEE

Press, New York, 1994.
Shape Functions
for Scalar and
Vector Finite
Elements
2.1 INTRODUCTION

The finite element method is used for modeling a wide class of problems by breaking

up the computational domain into elementsof [Link] interpolation

polynomials (commonly referred to as shape or basis functions) are used to approxi-

approximate the unknown function within each element. Once the shape functions are
chosen, it is possible to program the computer to solve complicatedgeometries by

solely specifying the basis functions. The element choice, however, needs human
intervention and intelligence to ensurea reliable solution of the problem at hand.
As will be shown later in this chapter, the development of a specialclassof elements
which mimic the character of electric/magnetic fields has proved to be the key in
obtaining robust solutions to three-dimensional problems electromagnetics.
in

In this chapter, we will discuss the derivation of node-based and edge-based

shape functions for one-, two- and three-dimensional finite elements. Node-based

shape functions have been used extensively in civil and mechanical engineering
applications as well as in scalar electromagnetic field problems. However, a full
three-dimensional vector formulation brings out numerous deficienciesin these tra-
traditional element shape functions [1], [2]. Edge-basedvector basis functions with
unknowns associated with element edges have thus been derived overcome
to the

problems related to nodal basis and theseare now extensively used for solving three-
[Link] will also describe the hierarchical nature
of the edge-based functions and their possible applicability in \321\203\321\202-based refinement

techniques.

37
38 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

2.2 FEATURES OF FINITE ELEMENT SHAPE FUNCTIONS

The polynomials used to interpolate finite element solutions on specific element

shapes have some distinct features over the wide variety of basis functions used in

other partial differential equation (PDE) or integral equation (IE) techniques.

2.2.1 Spatial Locality

Finite element shape have compact support within

functions each element. i.e.,
their scope of influence is limited only to the immediate neighboring [Link]
feature plays a pivotal role in the viability of finite elements over integral equation
(IE) methods. The limited scope of influence for the basis functions is a distinguish-

property
distinguishing of PDE techniques and leads to very sparse matrices in finite elements,
whereas IE techniques give rise to full, dense matrices resulting in poor scalability as
problem size increases.

2.2.2 Approximation Order

The of the approximation depends on the completeness

order of the poly-
polynomialsmaking up the finite element basis functions. Moreover, the form of the
polynomial function must remain unchanged under a linear transformation from
one Cartesian coordinate system to another. This requirement is satisfied if the
polynomials are completeto a specific order such as

m(.y, v) = ct + c2x + c3y + c-tx2 + csxy + cby2 B.1)

or when the extra terms are symmetric with respect to one another, as in the follow-

incomplete
following third-order polynomial
=
\302\253(\320\273-,
\321\203) C[ + c2x + W + +
\320\263**2 Cixy + c6y2 -f c7a:2.v -f- t^xy1 B.2)
Such approximation functions have the characteristic that, for fixed x or y, they are

always complete polynomials in the other variable. In general, we seek expansion

polynomials that will yield the highest order of approximation for a minimum num-
number of unknowns associated with that element shape. The two examples shown above
apply to two dimensions, but their extension to three-dimensional elements is
straightforward. Typically, the higher the order of the approximating polynomial,
the lower the error in the final solution if element size remains [Link] usual,
there is a trade-off here between the desired accuracy and the degrees of freedom
required to solve the problem.

2.2.3 Continuity

The order of the differential equation to be solved determines the order of

shape function to be employed. Functions with continuous derivatives up to the

nth order are said to be C\" continuous. For elliptic PDEs of order 2k (k = 1,2),
the continuity requirement is C*\021 for Galerkm methods. In most electromagnetic
problems,functions which exhibit C\302\260
continuity (i.e., function continuity) are used
since the discontinuous first derivatives are integrable. However, it is difficult to
Section 2.3 Node-Based
\302\246 Elements 39

impose continuity of order 1 and higher since the determination of suitable shape
functions is very complicated. For example, ninth-order polynomials are required to
obtain C1 continuity for tetrahedral elements. In electromagnetics, Wong and

Cendes [3] used C1 node-based triangles to avoid the problem of spurious modes
in the determination of cavity resonances. AH shape functions derived in the follow-
sections
following impose function continuity or C\302\260
continuity (not derivative continuity)
between elements.

2.3 NODE-BASED ELEMENTS

In node-based finite elements, the form of the sought function in the element is
controlled by the function values at its nodes. The approximating function can

then be expressed as a linear combination of basis functions weighted by the

nodal coefficients. If the function values \320\270/

at the nodes are taken as nodal variables,
then the approximating function for a two-dimensional element e with p nodes has
the form

p
v,.v) B.3)

Since the expression B.3) must be valid for any nodal variable uj\\ the basis function
N\"(x.y) must be unity at node / and zero for all remaining nodes within the element.

Shape functions can be derived either by inspection {Serendipity family) or

through simple products of appropriate polynomials {Lagrangefamily). It is easier

and more systematic to construct higher order bases in the Lagrange family while

progression to higher orders is difficult in the Serendipity family. However, Lagrange

shape functions have undesirable interior nodes and more unknowns than

Serendipity shape functions of the same order.

2.3.1 One-Dimensional Basis Functions

One-dimensional are employed for solving

finite elements problems where the
discretization domain involves a curve or a contour around a two-dimensional

structure; for example, the bounding curve of the cross section of an infinite cylinder.
These basis functions can also be used in conjunction with higher-dimensional finite

elements when the modeled structure can be decomposedinto a single dimension

without loss of [Link] is most convenient to derive one-dimensional basis
functions in terms of Lagrange polynomials [4]. Let us considera straight line

with endpoints .V| and xi. The basis functions for element e are then defined as
\342\200\224
x
x\\

vf _ ,-e

The basis functions have unit magnitude at one node and vanish at all others with
linear variation between the nodes. Higher order basis functions can be constructed
40 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

by inserting nodes between the endpoints of the finite element. If \320\273?,

/ = 1,2 \320\273,

are the nodes of the one-dimensional element, and we are interested in finding the
basis function for the /ah node x%. then the corresponding Lagrange polynomial
describing this basis function is given by [4]

{x-x4)(x-x'2)...(x-4-\\){x-xUi)---{x-x>tt)

The basis function defined above is of (n \342\200\224

1 )th order and correspondingly passes
zero n \342\200\224
1 times.
through

2.3.2 ltoo-Dimensional Basis Functions

Two-dimensionalfinite elements have found widespread use in modeling struc-

structures whose third dimension is significantly larger or smaller than the cross section,
thus ensuring little variation in the unknown parameters in this third direction. Two-
dimensional finite elements used to obtain reliableestimatesof three-
have also been
dimensional problems since the computational cost for obtaining two-dimensional
solutions is vastly less expensive than for three dimensions.

[Link] Rectangular and Quadrilateral Elements. The simple shape of the rec-
rectangular element permits its shapefunctions to be written down merely by inspec-
inspection. On examining the element shape given in Fig. 2.1. the shape functions can be
cast in the form

Figure 2.1 Reeiangular element.

Section 2.3 Node-Based
\302\246 Elements 41

where and
\320\273? y* denote the coordinates of the midpoints of the element edges, l(x and

hy represent the edge lengths and Ae denotes the area of the element. Each basis
function NJ has unit magnitude at the /th node, vanishes at the remaining three
nodes,and varies linearly. Hence it can be used m B.3) to represent u\". Higher order
rectangular elements include the eight node (three equispacednodes per edge)and
the twelve node (four equispaced nodes per edge)quadrilateral element discussed in
[4]. However, these elementscan model only regular geometries and decline in accu-
accuracy with excessive shape distortion. Thus, often they are not very useful in practice.

Irregular geometries can be modeled by using quadrilateral elements which can

also be viewed as distorted rectangles. To construct basis functions for a quadrilat-
quadrilateral
element, we need to use a transformation that maps a quadrilateral element in
the xy plane to a square element in the ?r\\ plane (Fig. 2.2). Such a transformation can
be found by satisfying the following relation at the four nodes of the quadrilateral
element
a- = a + />? + tJ7 -f dt-\321\202) \321\203
= a' + />'$ + -f d'
\321\201'\321\202) B.5)

The unknown coefficients\342\200\224a,b,c,d and a',b',c'.d'\342\200\224are solved by mapping the

four corners of the quadrilateral in the xy plane to the corner points of the unit

square in the tq plane. The eight equations thus obtained are

\342\200\224
)'\\ a'-b'-c'-\\-d'

y2 = a' + b'-c'-d'
+ b' + c' + d'
B.6)
=
.*\302\246, a + b + eH-d. yy
= a'

xA=a-h + c- d, v4
= -h'
\302\253' + c' - d'

On solving for the unknown coefficients in the above equation, the basis functions
can be cast in the following form

/= 1. B.7)
where ?0
= ff, and = and
\321\211 \321\211,

=
\321\203 B.8)

The variables (&, rj() denote the coordinates of the /th node in the (?, rj) coordinate
system. The linear quadrilateral is also known as an isoparametric element since the
shape functions defining the geometry and the nodal values are the same.

4 3

= -1 5=1

Figure 2.2 Transformation of a quadrilateral ele- X 1 1 = -1 2

\320\273
elementin the .v,r plane to a unit square in the f?j
plane.
Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

Higher order quadrilateral elements include the eight-node element (four cor-
corner nodes and four midside nodes) and the twelve-node element (four equispaced
nodes per edge).The basis functions for such elements can be found in [5].

Due to the irregular shape of the quadrilateral element, it is not easy to inte-
integrate the basis functions in the xy plane. To facilitate and generalize the integration
process, the conceptof the Jacobian is introduced. The Jacobian matrix, or more
specifically its determinant, transforms the infinitesimal area or volume element from
one coordinate system to [Link] we consider the above example, we are essen-
essentially transforming between the global (x, y) coordinates and the local (f, 17) coordi-
coordinates. By the chain-rule of partial differentiation, we can express the f derivative of
Nf as

3N[_-dN[0x \320\250[\320\264\321\203
\321\215\320\273-
\320\267\302\247 a? By a?

Similarly, taking the derivative

\321\211 and combining the expressions, we can write the

result in matrix form as

1 \320\223 BNf/Bx If I
= \320\255\321\205/*
\320\252\342\204\226\320\233
j l
\320\264^/\320\264\321\205

8N4/dr) J J
\\ dNf/\320\255\321\203
[\320\264\321\205/\320\264\320\263]
\320\264\321\203/\320\222\321\206] J| BNf/By J

Since (x, y) are known explicitly in terms of the local coordinates (?, rj), the Jacobian
matrix can be found explicitly in terms of local coordinates. Care must be taken in
the choice of the local coordinate system such that the Jacobian matrix is non-
singular. To find the derivatives with respect to x and v, we merely needto invert
the Jacobian matrix to yield

{2\320\251
\\dN!/<>y\\-[S] \\BNf/Bn)
The technique can easily be generalized for n-dimensional transformations if necess-
[Link] infinitesimal area element <IA can now be written as
ctxdy
= det[J]d^di] B.11)

Thus the integration of the irregular quadrilateral element is simplified considerably

by performing it over the local coordinate system instead of the global Cartesian
coordinate system.

[Link] Triangular Elements. Triangular elements are popular because they

can model arbitrary geometries. We will determine the shape functions of triangular
elements by using Lagrange interpolation polynomials. In their final expression, the
shape functions will be expresed in terms of the so-calledarea coordinates. Let us
consider a point P within a triangular element (Fig. 2.3) located at (.v. r). where
(.V;,.\302\273f) denote the coordinates of the /th triangle node. The area of the smaller
triangle formed by points P. 2. and 3 is given by

I \321\205 \321\203
I
=\302\246 l
\321\207 4 /2 B.1-

The area coordinate\320\246is then given by

Section2.3 Node-Based
\302\246 Elements 43

2J
Rgm\302\273 Triangular element.

Area/\302\27323
41 I '

is the area of the whole triangle and can be found

where \320\224 from B.12) by replacing x
and \321\203
in the first row with .v, and yt. Similarly, the two remaining area coordinates
L2 and Lj are given by

A2_AreaP31
~
*~
\320\224 Area 123
._ \320\2243 Areafl2
3 ~ ~_
~A Area 123

The values for x and v inside the triangular element reduce to

(=1
> J>
i=l
?>,
/=l
B.14)

where the latter condition 52'=1L,- = 1 a result

is of the area identity +
\320\224|

+ \320\2243
\320\224\320\263
= \320\224.
Alternatively, \320\246,IX,
and L\\ can be obtained in terms of x, y, and
the vertex coordinatesby solving the system of equations B.14).
The coordinate\320\246is zero on the edge opposite to vertex 1 and unity at vertex 1.
Its variation along the height of the triangle is displayed in Fig. 2.4. The remaining
two area coordinates associated with the other two vertices behave similarly, vanish-
on
vanishing edge opposite to the corresponding vertex
the and having unit magnitude at
the vertex it belongs to. This feature combined with spatial locality and C\302\260
continuity

qualifies the area coordinates as suitable basis functions N]' for a triangle when the

interpolation order is linear. That is,

Nf^L'l B.15)
Higher order basis functions for triangles can be derived using the procedure
given in [4] and [6]. In general, the shape function /V/ for node /, labeled as (/, /. K),

is given by

= n B.16)
44 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

t,=0

=
\320\2460.25

\\/ \\ \\ \\

\\ \\ \\ \\
Figure 2.4 Area coordinates of a triangle.

Figure 2.5 Six-noded triangular element

2 supporting quadratic bases.

where i
\320\246,
= I e are area cordinates defined previously and is the
\320\240\"(\320\246) poly-
polynomial

j -
.t=o

B.17)

Similar definitions apply for P\",{L\\) and

Using the formulae given above, the basis functions for a quadratic triangle

(see Fig. 2.5) can be conveniently defined as

N* = -
\320\246B\320\246 I), i = I. 2, 3 CORNER NODES
B.18)
M1DS1DE NODES

Besides simplifying the treatment of basis functions for a triangle, simplex

coordinates greatly facilitate the integration of an arbitrary function over a tri-
triangular region. A very useful integration formula in terms of the area coordi-
Li,
s^\342\200\224Li, and L$\342\200\224over a triangular domain is given by
Section 2.3 Node-Based
\302\246 Elements 45

where a, b, and are

\321\201 integers and is
\320\224 the area of the triangular region.

2.3.3 Three-Dimensional Basis Functions

Shape functions for three-dimensionalelementscan be described in a precisely

analogous way to their two-dimensional counterpart. However, the simple rules for
inter-element continuity given previously must be modified. The nodal field values
should now interpolate to give continuous fields across the face of each element.

2.3J.I Rectangular Bricks. The simplest polynomial approximation to a rec-

rectangular brick element is the trilinear function
tf{x, y4 z) = <f + b'x + c\"y + d*z + e\"xy +fyz + gezx + ifxyz B.20)
whose eight parameters are uniquely determined by matching z) to the field
uc'(x, y,
values u\" at the eight corners of the brick. This results in eight equations that are
solved to determine the coefficients <f, h\",..., h\". The final expression is in the form

8
ue(x, y. z) = J^ u4Nf(x, y,z) B.21)

However, this approach of formulating B.21) is cumbersome and can be easily

avoided on writing down the required basis functions by mere inspection. Since
the basis function N\302\260
must be unity at node i and zero at the remaining nodes,
the eight interpolation functions can be written down as

where jc'!. y*r, and z\"c denote the coordinates of the center of the element, hex, hey, and he:
represent the edge lengths of the element and V is the element volume.
Brick elements of the Serendipity family are derived in [4]. To obtain higher
order basis functions, progressively larger number of nodesare uniformly placed on
46 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

the element edges. Bricks with nth-order interpolation functions (with \320\270
+ 1 per

edge) require 8 + 12(n- I) degrees of freedom. Higher order bricks are rarely
used in electromagnetics applications since the regular shape requirement and

decline of accuracy with excessive shape distortion place severe limitations on the

generality of the geometry to be modeled.

Shapefunctions for hexahedral elements or distorted bricks can be derived by

mapping the element in the xyz coordinate system onto a standard cubein a new ^in-
^incoordinate system. We proceed along lines similar to the derivation of bases for the
quadrilateral element in Section 2.3.2. We express the Cartesian coordinates (x,y,z)
m terms of (?, rj, ?) as follows:

x = o,
=
\321\203 a2 + *rf + e2n+ && + e,Jv +\320\250 + &ft + Arf 4f B-22)
2=

The unknown coefficients\342\200\224a,. A;, c,, dh / = 1,2,3--can be obtainedby a one-to-one

mapping of the corners of the hexahedral element to the corner points of the unit
cube. The desired transformation thus yields the basis functions Nf of the hexahedral
finite element

= + + WXI + M) B-23)
Nf id \320\2501

with (?,-, ?//, f,-) denoting the coordinates of the ith node in the ?>jf coordinate system.
As before, the relationship between the (\320\273;, z)
\321\203, and (?, ?;. f) coordinates is given by

8 8 8

=
\342\200\242v
?

As in the case of the two-dimensional quadrilateral element, the Jacobian

comes to our aid for calculating the volume integral over the arbitrary hexahedral

domain. The Jacobian matrix [J] is of orderthree and is given by

f ajvf/at 1
fax/af \321\215\321\203/ag Wax
I dN4 = \320\255\321\205/cit]
\320\264/\320\264
= \320\243]
{ BNf/By
\\
) [ HNf/d:

The volume element transformation from the global xyz to the local ?qf coordinate
system is expressed as

dx dy dz = delM dt- di) d% B.25)

23.3.2 Tetrahedral Elements. The three-dimensional analogue of a two-

dimensional triangle is a tetrahedron (four-faced element).Onceagain, we can intro-

introduce special coordinates, called volume coordinates or simplex coordinates, to sim-

simplify the derivation of shape functions. If P is a point within the eth tetrahedron
shown in Fig. 2.6, the four volume coordinates are given by
Section 2.3 Node-Based
\302\246 Elements 47

Figure 2.6 Teirahedral clement. ' z

Volume P234
~
1
Volume 1234
Volume P341
2
\"Volume 1234

Volume P412
3~ Volume 1234

_ Volume PI 23
4 ~
Volume 1234 B.26)
and any position within the element is specifiedby
4 4 4

with (xhyhZi) being the coordinates of node /. As for the two-dimensional case
(triangular elements), the basis functions ,Vf are equal to the volume coordinates,
i.e.,

Nf
= Li, / = 1 4 B.27)

Quadratic shape functions for a tetrahedron necessitates the use of ten node
points: the four corner nodes and the remaining six at the midpoints of the edges.
The shape functions for the quadratic tetrahedron are given by

Nf = -
\320\246B\320\246 1), 1=1 4 CORNER NODES
B.28)
N) = 4L}'Z.?, 1 = 5 , 10, MIDS1DE NODES

j and \320\272
are endpoints of each edge
Similarly to triangular elements, volume coordinates greatly simplify integra-
over
integration tetrahedral elements. A useful formula for integrating over the volume of a
tetrahedron is

= 6V B.29)
volume
dxdydz
(a + h ^f
+ C+ d + 3)!
where a, b, and
\321\201 d are integers and V is the volume of the tetrahedron.
48 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

[Link] Triangular Prism Elements. Other three-dimensional elements that

have simple shapes include the triangular prism and isoparametric elements. To
ensure that a small number of elements can model a relatively complex region,
distorted prisms can be used in conjunction with rectangular bricks. The shape
functions for the first-order prism element (shown in Fig. 2.7)is given by

/=1.4

'= 2,5 B.30)

/ = 3.6

Figure 2.7 Linear triangular prism clement.

where {;u
\342\200\224\342\200\224
2(: zc)/heighi varies linearly from -1 to +1 over the height of the
prism and is zero at the midpoint ze of the vertical edge (joining nodes 1-4, 2-5.
or 3-6 in Fig. 2.7). Here \320\246
refer to the area coordinates of the triangle that forms the
cross section of the prism (i.e., the triangle formed by nodes 123 or 456).
The quadratic node-based triangular prism has 15 nodes\342\200\224one each at the
corners and at the midpoints of each of the nine edges. Shape functions for the
quadratic and cubic triangular prisms can be found in [4]. For a more efficient
discretization using the fewest unknowns, prisms and bricks can be combined.
This is easily done because prisms and brickscan be readily connected by sharing
the same nodes and edgesat their boundaries.

2.4 EDGE-BASED ELEMENTS

In electromagnetics, we encounter serious problems when node-based elements are

employed to represent vector electric or magnetic fields. First, spurious modes are
observed when modeling cavity problems using node-based elements [7]. Nodal basis

functions impose continuity in all three spatial components whereas edge bases
Section 2.4 \302\246
Edge-Based Elements 49

guarantee continuity the tangential component. This feature

only along mimics the
behavior of field components along discontinuousmaterial boundaries, and its
importance in resolving the problem of spurious solutions will be discussed in detail
in Chapter 5. Second, nodal bases require specialcarefor enforcing boundary con-
conditions at material interfaces, conducting surfaces, and geometry corners [8]. The
first limitation can alsojeopardizethe near-field results of a scattering problem, the
far-field typically escapes contamination since spurious modesdo not radiate.
Edge-based finite elements, whose degrees of freedom are associated with the
edges and the faces of the finite element mesh, have been shown to be free of the
above shortcomings. The detailsof how edge bases avoid the pitfalls of nodal basis
functions will be discussed in Chapter 5. Edge basis functions were described by
Whitney [9] over 35 years ago and have been revived by Bossavit and Verite [1],
Nedelec [10], and Hano [11]in the recent past. It was Nedelec's landmark paper [10]
that laid down the guidelinesfor constructing finite element basis functions that span
his curl conforming space with degrees of freedom associated with the edges, faces,
and elements of a finite element mesh. Mur and de Hoop [12]:van Welij [13]; Barton
and Cendes [14];Jin and Volakis [15],[16];and Lee et al. [17] among severalothers
have extended their applicability to various two- and three-dimensional shapes and
even constructed higher order elements for a more accurate approximation of the
field values. More recently, Monk [18] provided error estimates and convergence
proofs for edge bases. The derivations leading up to the convergence proof is beyond
the scope of this book but the final result is stated below. If we denoteby \320\257\320\233'(\320\257
the

standard Sobolev space of functions of order ,y in ft and by || \342\200\242 ||, the norm on this
space, we can then define the space

H(curl;ft) = e (L2(ft)K|V
(\320\270 x \320\270
\320\265
(L2(ft)K)

and its corresponding norm

\\W\\h- =(ll\302\253llo + IIVl/2

If E is the exact solution of Maxwell'sequation in ft and EA is the finite-dimensional

approximation to it (in essence, the edge basis functions), ihen assuming that
E e (//*+l(ft))\\ the following result holds

\320\246\320\225-\320\225'-\321\203,,,
<\320\241\320\233\320\220-||\320\225||\320\264.+1 B.31)

provided is not
\321\201\320\276 an interior eigenvalue and h is sufficiently small. In B.31), A is the

discretizationparameter, \320\272 is the order of the basis function and is a constant

\320\241

independent of h. Thus higher order bases lead to lower errors in the solution when
the sampling size is sufficiently small. The convergence is also optimal in \320\233.

2.4.1 IVvo-Dlmenslonal Basis Functions

[Link] Rectangular Elements. We consider the rectangular element first since

its vector basis function is usually the easiest to formulate. For the elementshown in
Fig. 2.1. we can find its edge-based finite element basis function merely by inspection.
If the edges are numbered accordingto Table 2.1the vector basis functions can be
written as
50 Shape Functions for Scalarand Vector Finite Elements \302\246
Chapter 2

TABLE 2\320\233Edge Numbering for Rectangular Element

Edge No. h k
I 1 2
2 4 3
3 1 4
4 2 3

where x, y, and z are the unit vectors in the Cartesian coordinate [Link] above
basis functions have unity value along one edge and zero over all others, i.e..

where is the Kronecker

\320\246
delta and ij is the unit vector along the yth edge.
A graphical illustration of the W\\ vector basis function is given in Fig. 2.8. In
this figure, the largest value of IH^I corresponds to the largest vector length.
Using the above Wf the electric field within the finite element can be repre-
represented as

? B.32)

where now ?f denotesthe average tangential field along the /th edge. The has\302\273

functions W'( guarantee tangential continuity across inter-element boundaries sinoc

they have a tangential component only along the /th edge and none along the other
edges. They are also divergencelesswithin the element and possess a constant non-
nonzero curl. It should be noted that by taking the cross-product of z with W,, we obtain

Figure 2.8 Illustration for rectangular

of \320\251
clement.
Section 2.4 \302\246
Edge-Based Elements 51

basis functions which possess normal continuity across element boundaries, have zero
curl and non-zero divergence. The latter are ideal for representing surface current
densitiesand are known as rooftop basis functions in electromagnetics. They have
found extensive use in the solution of integral equations [19] and hybrid finite ele-
element-boundary integral implementations.
Edge-based vector elements can be derived
bases for quadrilateral by carrying
out the transformation detailed in of nodal basis for quadrilaterals
the derivation in

the previous section and then taking the gradient of the resulting expression for each
[Link] two shortcomings. First, the integrals associated with edge-
based quadrilateral elements do not lend themselves to easy evaluation. Second, they
may not be divergence free. However, their ability to model complicated shapes with
a lesser number of unknowns than tetrahedra and the inherent property of enforcing

tangential continuity across elements makes them attractive for use in two-dimen-
vector
two-dimensional formulations.

[Link] Elements. We consider

Triangular again the triangular element
depicted in the edges of an arbitrary
Fig. 2.3. Since triangular element are not
parallel to the x- or ^-axis,it is not easy to guess the form of the vector basis function
by inspection. Therefore, the vector basisfor a triangular element will be expressed
in terms of its area coordinates,L\\, Lj, and LLj. These are the Whitney elements. If
the local edge numbers are denned according to Table 2.2, then edge bases for a
triangular element are defined as

W{ = N'u
=
\320\246(\320\246VL/
-
LJjVUj), ij = 1, 2.3 B.33)
where W\\ denotes the basis function for the A'th edge of the eth element and ly
= ?*
is the length of the edge formed by nodes i and j of the triangle. The vector field

inside the triangular element can, therefore, be expanded as

B.34)

where ?? denotesthe tangential field along the /cth edge. It can be easily shown that

the edge-based functions defined in B.33) have the following properties within the

element

VxJVj
= x
VLJ
\320\230\321\206\320\247\320\246

If is
\321\221| the unit vector pointing from node 1 to node 2 in Fig. 2.3, then
\342\200\242=
VL|
\321\221|
\342\200\224
\\jt\\ and \302\246
VL2
\321\221|
= \\jl\\. Since L\\ is a linear function that varies

TABLE 2.2 Edge Numbering for Triangular Element

Edge No. Node i\\ Node t;

1 I 2
2 2 3
3 3 1
52 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

from unity at node 1 and zero at node 2 and L$ is unity at node 2 and zero at node 1,
we have

e, Afo
=
\342\200\242
L? + Z.1 = 1 B.35)

along the entire length of edge 1. This implies that \320\233^2has a constant tangential

component along edge 1. Moreover, since U\\ vanishes along edge 2, VL\\ is normal to
edge 2, U, vanishes along edge3, and VL$ is normal to edge 3, N\\2 has no tangential
component along these edges. Similar observations apply to N13 and [Link],
tangential continuity is preserved across inter-element boundariesbut normal con-
continuity is not. Fig. 2.9 shows the actual variation of the basis function for the edgeof

a right triangle that is opposite to the node associated with the right angle. A
different method of constructing edge bases for triangular elements is given in [20],
[21].
Higher order vector basis functions involve a node at the midpoint of
adding
each edge and including the contribution of facet to the approximating
elements

function. Unknowns in the triangular element are assigned as shown in Fig. 2.10 [17].

Figure 2.9 Variation of edge basis function

for the edge opposite to the right angle.

_2 Figure 2.10 Triangular edge element with

3 two unknowns each per edge and per face.
Section 2.4 \302\246
Edge-Based Elements 53

The tangential projection of the vector field along edge {i,j} is determined by two
unknowns E\\ and E)and two facet unknowns\342\200\224-F\\ and F2\342\200\224areprovided to allow a
quadratic approximation of the normal component along two of the three edges.
Only two facet unknowns are required to make the range space of the curl operator
complete to first order. Therefore, there are eight degrees of freedom for each trian-
triangular element. Since the edge variables provide common unknowns across element
boundaries, tangential continuity of the field over the boundary is [Link],
an obvious disadvantage of these elements is that the two-facet variables cannot be
symmetrically assigned. This disadvantage can be avoided by employing third-order
edge bases [22]. The higher order approximation to the vector field within the ele-

element is given by

B.36)

where we have arbitrarily chosen the facet variables to lie on edges 1 and 2. These
variables are local unknowns associated with each separate triangular element and
are included to provide a linear approximation for V, x E,, where the subscript t
denotes the tangential components of the operator. This property turns out to be
very important in the selection of the order of the basis function to be used in the
modeling process. The basis described by B.36) can be classified as belongingto the
Hl(curl) space. The Hk(curl) space consistsof those vectors whose inner products
are square integrable and whosecurl consists of complete polynomials of order k.
The basis given in B.33) thus belongs to the H\302\260(curl) space, since its curl is merely
constant within the finite element. The basis generated by excluding the facial con-
contribution would result in six unknowns\342\200\224 two per edge\342\200\224but the order of the approxi-
approximation would still be H\302\260(enr[) and does not add to the accuracy of modeling the H-
field while doubling the unknown count. It should also be noted that the form of the
facet bases in B.36) are different in the original paper [17].This is due to recent
analysis [22] that shows that the Nedelec constraints [10]are met by B.36), resulting
in smaller dispersion error and better conditionedmatrices.

2.4.2 Three-Dlmenslonal Basis Functions

Edge-based elements have facilitated to a great degree the finite element ana-
analysis of three-dimensional structures in electromagnetics. Linear nodal bases with
their problem of spurious modes and difficulty in maintaining only tangential con-
continuity across material interfaces are not as convenient for electromagnetic field
simulations in three dimensions. On the other hand, the introduction of edge-
basedshape functions provides a robust way of treating general three-dimensional
problems having material inhomogeneities and structural irregularities like sharp
edges and corners.
In the following section, we will consider first the simple rectangular bricks and
will proceed to present edge-based shape functions for more complicated finite ele-
elements such as tetrahedrals and curvilinear hexahedrals. The chapter is concluded
with a brief discussion on hierarchical edge elements.
54 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter]

[Link] Rectangular Hexahedrals. As in the two-dimensional

Bricks and case,
we derive the edge-based shape for a rectangular brick (seeFig.2.11)
function bj

simple inspection. Since a constant tangential field component must be assigned to

each edge of the element, we can expressthe shape function along each edge of the
element as [15]

Figure 2.11 Rcciangular hrick element, Tfe

numbers denote the local node numbensj
scheme.
Section 2.4 \302\246
Edge-Based Elements 55

where tfx, denote the edge lengths in

hev. If. the x, y, and z directions, respectively, and
the center coordinates of the brick are given by zj.'). If the local edge numbers
v?,
(\320\273?,

are defined as in Table 2.3, the vector field within the element can be expressed as

B.37)
k=\\

where El represents the value of the electric field along the A*th edge of the <?th

element. The vector bases W? defined for the rectangular brick element have zero
divergence and a nonzero curl. Furthermore, the expansionB.37)guarantees tan-

tangential continuity of the electric field across the surfaces of the elements.
A rectangular brick element has limitations in the sense that it is unable to
model irregular geometries. For this reason, the analog of the two-dimensional

quadrilateral (the hexahedral element) is more attractive for modeling practical

three-dimensional problems. As in the case of the quadrilateral element in two
dimensions, a hexahedral element in Cartesian coordinates can be seen as the
image of a unit cube under a trilinear mapping to the t/tf coordinate system (see
Fig. 2.12).
Let us consider those faces for which ? = constant. Therefore, V? must then

possess a normal component

only on that face. Since ? varies linearly along Ihe edges
that are parallel to the ?-axis, the vector function has
\320\251 nonzero tangential compo-
components only along those edges that are parallel to the ?-axis. Using the node-based

TABLE 2.3 Edge Definition for Rectangular Brick

Edge No. Node Node

1 1 2
2 4 3
3 5 6
4 8 7
5 1 4
ft 5 8
7 2 3
8 6 7
9 1 5
10 2 6
11 4 8
12 3 7

4
Figure 2.12 Mapping of a hexahedral de-
dement\321\216
unit
\320\260 cube.
56 ShapeFunctions for Scalar and Vector Finite Elements \302\246
Chapter 2

expression for the shape function in a hexahedral element given in B.23). we may
write the corresponding edge bases as

W%
= A + n)(
\321\211 1 + UO V? edges || to $-axis B.38)
^

\\
= -\302\261
A + \320\2501 + M) V^7 edges || to rj-am B.39)

= T <' + **#<'+ W V< edges || to f-axis B.401

\320\276

where (&, &)

\321\211, denote the coordinates at the kth edge and hek is the length of the

klh edge belonging to the eth element.

The vector bases derived above possess all the desired continuity properties of

edge elements and generally result in about half the number of unknowns generated

by telrahedral gridding. The difficulty in generating a finite element mesh of \320\260\3

arbitrary structure using hexahedra can be a seriouslimitation. In practice, often,

a combination of hexahedra and telrahedra in a finite element mesh is used and

continuity is imposed across the different element interfaces to solve the problem
However, this may result in ill-conditioned matrices.

[Link] Elements.
Tetrahedral Tetrahedra are. by far, the most popular el-

element shapes employed for three-dimensional [Link] is because

to be
the tetrahedral element is the simplest tessellation shapecapableof modeling arbi-

arbitrary three-dimensional geometries and is also well suited for automatic mesh gen-
generation. The derivation of shape functions for these elements follow the same pattern
as that for triangular vector basis functions. If we consider the tetrahedron shown in
Fig. 2.6 and define the edge numbers accordingto Table 2.4, we have

= = - 4
W\\ N4; \320\246\320\247\320\246)./.y=l
\320\225\302\253(\320\246\320\247\320\246 B.4li

where again ty =
tk denotes the length of the edge between nodes / and./, which in
turn define the A-th edge. The vector field within the element can then be expandedas

B.42:
ft=l

TABLE 2.4 Edge Definition for Tetrahedron

Edge No. Node if Node /j

1 1 2
2 1 3
3 1 4
4 2 3
5 4 2
(y 3 4
Section 2.4 \302\246
Edge-Based Elements 57

where the coefficients represent

\320\225\320\246 the average value of the field along the Ath edge of
the fth element.
An explanation of the physical character of the edge-based
elegant interpola-
function
interpolation by Bossavit [23]. Let us consideredge number
is given 1 connecting

nodes 1 and 2 in Fig. 2.6. Since V/\320\233is orthogonal to facet A34) and VLf is orth-

orthogonal to facet B34}, the field turns around the axis 3-4 and is normal to planes
containing nodes 3 and 4. The field thus has only tangential continuity across el-
element faces. Edge elements can also be describedas Whitney elements of degree one
and can be broadly classified as belonging to the ^(curl) space.
Whitney elements of the second degree are calledfacet elements becausethey
are constant over the face of the tetrahedron. The vector function for the facet
elementcan be written as

=
2(UtVLe] x
\342\204\226\321\202 VL? + L]VLck
x \321\205
\320\247\320\246+\320\246\320\247\320\246
VL'j), ij.k=l 4
B.43)
As explained in [23], we now have a central field (as if emanating from node 4 in Fig.
2.6) on each of the two tetrahedra that share the face A,2,3). The field can be
imagined as coming from the 'source' 4, growing, crossing the facet, and vanishing

into the 'weir 4'. the fourth vertex of the other tetrahedron. Thus, this field has
normal continuity and the flux across the facet forms the degree of freedom for the
element.
Alternative expressions for linear basisinside a tetrahedron have been derived
in [14]. They are given by

w _ 1 f7-< + x f. r in tne tetrahedron _

^
Wl~i ~
87w
otherwise {'m>
10,
with

f7-/ = B.45)
^r/lxr,1

^ =
'4f B.46)

in which /=1,2 6, Vt. is the volume of the tetrahedral element, e,- = (r,, \342\200\224
r(|)//\302\273,
is the unit vector of the /th edge, and //, = |r,^ \342\200\224
r,, |
is the length of the /th edge with
and
\320\263,, r,, denoting the position vector of the /| and i2 nodes. It can be shown that

B.41) is identically equal to B.44) when simplified. Therefore,

g7_,
=

where i| and i2 are given in Table 2.4. The basisfunctions given in B.44) have zero
divergenceand constant curl (VxHf = 2g,). The form of the basis functions given
in B.44) is similar to the zeroth-order edge elements postulated by Nedelec [10].
order of the polynomial approximation for the first-order
The edge element
given in B.41) or B.44) can be taken as 0.5. This is because the value of the basis
function is constant, i.e., 0A), along the edge it supports and is linear everywhere
else within the element. Mur and de Hoop [12]presented edge elements which are
consistently linear, yielding a linear approximation of the field both inside each
tetrahedron and along its edges and faces. However, the curl of the basis is still
58 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

0A). Since this requires two unknowns per edge, there are twelve degrees of freedom
per element. The basis functions in [12] are derived by first defining the outwardly
directed vectorial areas of the faces as

where r/, / = I..... 4 denote the position vectors of the vertices of the tetrahedron
and i,j, k, I are [Link] the edge-based vectorial expansion function is defined \320\253

,4. '' \320\244.1 B.48)

where V is the volume of the tetrahedron and is a

\321\204(\320\263) linear scalar function of

position given by

in which 17, position vector of the centroid of the tetrahedron.

is the We observe thai
equals
\321\204{(\320\263) unity when r = and zero for
\302\273\342\200\242,\302\246 the remaining verticesof the tetrahedral
element. In that sense, they are very similar to the simplex or volume coordinates
mentioned earlier. They also satisfy the following equalities:

1=1 1=1

The edge basis function

/V'} is
a linear vector function of position inside the tetra-
tetrahedral element, and its tangential component vanishes on all edges of the elemeni

except the one joining vertices i and j. jVJ/

varies linearly along the edge formed bt

j such that = 0 while

nodes i and \342\200\242
/V? r,

These basis functions of divergence and curl.

have nonzero values
An
inspection expressions for the vectorial areasA; reveals that the form
of the
is identical to that obtained by taking the gradient of one of the simplex or volumt
coordinates mentioned earlier, in other words,the three components of the vector \320\233,
have the same functional dependence as that obtained by VL^'

V 1
~
X-> 1 \\h
\"
x det -y det V 1
+ 2 det x3 1 \320\233
V I *\"
\320\2434 1 \024 1
\342\200\242V4 >'4

where the volume coordinate

is
\320\246 for a tetrahedron defined in B.26), \"det\" indicates
the of the determinant of the matrix
value and (Xj,yt,2() denote the coordinates of
the ith vertex. This is only to be expectedsincethe gradient of the shape function is

normal to its corresponding edge in two dimensions and normal to its corresponding
face in three dimensions. The basisfunctions with consistently linear interpolation in
the tetrahedron can thus be rewritten in a more convenient notation as
= \"(I= /.7=1. ,4. A49)

where the normalization factor ty

was introduced. Note how the first order edge
basis is similar in form to the zeroth-order edgebasis in B.41).
Section 2.4 \302\246
Edge-Based Elements 59

Figure 2.13 Tetrahedral clement.

Still higher order basis functions are sometimes necessary for rapidly varying
fields. The second-order edge basis @{rim5)) for a tetrahedral element was first
presented by Lee, Sun, and Cendes [24].We need 20 degrees of freedom to achieve
a quadratic approximation of the vector field inside a tetrahedron (see Fig. 2.13).
Accordingly, the field within a tetrahedron can be written as

fa! fal

- +
-
'IUKUjVLi UkVL'j) FiLJj(Lek\320\243\320\246 \320\254'\320\243\320\246) B.50)

where i,j\\ form

\320\272 cyclic indices. The facet variables F[ and Fi are common unknowns
for two letrahedra that share the same face. Even higher order edge-based elements

complete up to polynomial order two can be [Link] tetrahedral element

now has 30 unknowns\342\200\224three along each edge and three on each face.

[Link] Triangular Prism Elements. The primary attraction of triangular

prisms lies in the fact that they yield fewer unknowns than tetrahedrals while retain-
retainingthe ability to mesh arbitrary geometries unlike hexahedrals. Moreover, it is
sometimes possible to extrude the volume mesh out of an existing surface mesh
using triangular prisms. This feature, however, may not always lead to good quality
elements, especially when the geometry is non-planar with sharp corners. Finally, it
is not easy to construct edgebasis functions for such elements. 6zdemir and Volakis
[25] proposed edge-based shape functions for right-angled and distorted triangular
prisms. The vector basis functions derived in [25] are a combination of edge basis
over the triangular cross section and a linear variation over the height of the prism.
A sketch of the basis function over the triangular and quadrilateral faces is presented
in Fig. 2.14. One of the shortcomings of these bases is the lack of tangential con-
continuity across element faces when the prisms are distorted, i.e.. the vertical arms are
not at right angles to the plane of the triangular faces.
60 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

(x,y,z)

(a) (b) (c)

Figure 2.14 Sketch of edge basis function over triangular and quadrilateral faces of
prism element. [Courtesy o/T. Ozdemir.]

Figure 2.15 Di&toclcd triangular pjisa

The geometry of an arbitrary triangular prism is shown in Fig. 2.15. The edge

basis functions for the top triangle edges are given by

= /,./=1,2.3; =
\320\233 1.2,3 B.5\320\246
\320\251^\320\251 \320\254^\320\246\320\247\320\246-\320\246\320\243\320\246)\320\260,

and those for the bottom edges are

W%
=
\320\251
= - 1
\320\254\321\206{\320\246\320\247\320\246
\320\246V?D(
-
s), l,j = 4. 5, 6; =
\320\272 4. 5, 6 B.5:,
and the vertical edges are

/=1.2,3; = 7,8,9 B.531

M). \320\272
Section 2.4 \302\246
Edge-Based Elements 61

In the above equations, \320\246are the node-based shape functions (area coordinatesof
the triangle) defined earlier and ,v is a normalized parameter which is zero at the
bottom face and unity at the top face of the prism. It should be noted that the basis
functions for the top and the bottom edges are exactly similar to that of a triangular
basis scaled by a dimensionless parameter. For the vertical edges, the vector v is a
linear weighting of the unit vectors v,. vj, v3 associated with the vertical arms and is
defined as

This particular choice of v minimizes tangential discontinuity across inter-element

faces.
It is worthwhile to note that the edge basis functions for the triangular prism
stated above are not Nedelec-type [10] elements since the first-order bases are not of
the form

\320\251
= = 1
\320\272 number of edges B.55)
\320\260\320\272+\320\240\320\272\321\205\320\263,

where a*. fik are constants and r is the position vector insidethe finite element. The
edge bases for the top and bottom faces fall into the Nedelec form but the bases for
the vertical forms do not, hencethe loss of tangential continuity across elementfaces.

[Link] Curvilinear Elements. Wang and Ida proposed a systematic method

for the construction of curvilinear elements in [26]. The vector shape function is

expressed in the following form:

= = 1 \320\234
**(*\342\200\242). \320\272
\320\251(\320\263)\321\204*($.ij. \320\236 B.56)

where rj, f)are completely defined in the local coordinate

\321\204\320\272(%, system, vk contains
the edge and facet information and M denotesthe number of degrees of freedom in
the element. These basis functions differ from the bases described earlier in the
chapter in that they are constructed in the local coordinate [Link] direc-
direction vectors are defined in global coordinates, the bases are uniquely defined. The
joint effect of \321\204\320\272
and vk ensures that W*k is unity at node and
\320\272 zero elsewhere. These
basis functions usually lead to a symmetric system of equations. Antilla and

Alexopoulos also proposed a curved brick superparametric element with 8, 27, and
64 nodes in [27] for solving scattering problems. Superparametric elements attach

unknowns to a lesser number of points than required to define the geometry. The
advantage of curvilinear elements lies in the fact that they can model curved surfaces
with more accuracy and lesser number of unknowns than rectilinear elements.
Analytical surfaces and even complicated non-planar surface features can thus be
modeled exactly at low computational cost. However, many mesh generation
packages cannot construct curvilinear elements for arbitrary geometries.
62 ShapeFunctions for Scalar and Vector Finite Elements \302\246
Chapter 2

[Link] Hierarchical Vector Elements. Finite elementsare saidto be hierarch-

when
hierarchical the basis functions for an element are a subset of the basis functions for

any element of higher order [4]. Hierarchical elements find use in a class of adapiivc
finite elements (called /^-refinement) where the order of approximation is improved

by refining the order of the polynomial basis functions instead of refining the mesh

density. However, there is usually a trade-off when higher order basis functions
extract a heavy price m terms of computer resources.A major problem with going
to higher order bases is the increased density of the finite element matrix and the

likely worsening of the matrix condition [Link], some structures may

have features for which a lower order approximation is sufficient to model the field
variations. This is especially the case where field variations are either uniform or
constant and a higher order of interpolation may actually degrade the solution. It
turns out that the vector finite elements defined by Nedelec [10] and subsequently
derived by Barton and Cendes [14] and Leeet al. [17]have a hierarchical structure.
Vector elements complete up to polynomial order two are available, and basis func-
functions of a given order are fully compatible with basis functions of lower or higher
orders. Thus elements of different orders could be used in the same mesh.
Specifically, lower order elements could be used in regions where field variation is

slowly varying and higher order elements in regions where the field varies rapidly.
The implementation of hierarchical vector elementscan be difficult, especially
at the transition boundaries where elementsof one order mergeinto the elements of

higher or lower order. If several vector elements share an edge, the field tangent to
the edge must be made identical in each of the tetrahedra. This is done by carefully
matching the coefficients of the vector basisfunction corresponding to that edge. For

tangential continuity across a face, the same equality must be enforced between the
coefficients of all the edge and facet functions associated with the face. Table 2.5
given in [28] shows the basis functions for hierarchical vector finite elements, li
should be mentioned that for the zeroth-order edgeelement,the described polyno-
polynomialapproximation to be of order 0.5 is somewhat of a misnomer. It should be taken
to mean that the field variation along the edge is constant, i.e., 0(ru), and the
variation normal to the edge is O(r'). Averaging the orders, albeit a mathematically
dubious procedure, yields the described polynomial [Link] the plus side, the table
offers a concise view of the hierarchical nature of these edge elements. Higher order
basis functions are constructed by systematically adding the extra terms up to the
desired order. It should be noted that the bases for the tetrahedron with six and 20

unknowns shown in Table 2.5 is identical to the \320\2570(\321\201\321\2107)

and Hl(curl) edge basis
given in B.41) and B.49), respectively.

TABLE 2.5 Hierarchical Basis Functions for a Tctrahcdral Element

Unknowns per
Element Type Polynomial Order Element Basis Function

Edge 0.5 6 -
L,VLj LjVL,
Edge 1 12 V(L,L,)
Face 1.5 20 Ul.,VLk-UVL))
Face -
Lj(LkVL, L,V^)
Edge 2 30 V\\L,L,a,-Lt)]
Face 4L,L,Lk]
References 63

REFERENCES

[1] A. Bossavit and J. \320\241Verite.

A mixed FEM-BIEM method to solve 3D eddy
current problems. Trans. Magnetics, 18:431-5,March 1982.
IEEE

[2] J. P. Webb. Edge elements and what they can do for you. IEEE Trans.
Magnetics,29:1460-1465. 1993.
[3] S. H. Wong and Z. J. [Link] finite element-modal of three-
solution
dimensional eddy current problems. IEEE Trans. Magnetics, 24F), November
1988.
[4] Zienkiewicz.
\320\241
\320\236. The Finite Element Method. McGraw-Hill, New York, Third
edition, 1979.
[5] K. Tuncer, D. Norrie, and F. Brezzi. Finite Element Handbook. McGraw-Hill,
New York, 1987.

[6] J. M. Jin. The Finite Element Method in Electromagnetics. John Wiley & Sons,
New York. 1993.
[7] Z. J. Cendes and P. Silvester. Numerical solution ofdielectric loaded wave-
waveguides: 1\342\200\224Finite element analysis. IEEE Tram. Microwave Theory Tech.,
1970.
118:1124-1131,
[8] X. Yuan, D. R. Lynch, and K. [Link] of normal field continuity
in inhomogeneous scattering calculations. IEEE Trans. Microwave Theory
Tech., 39:638-642, April 1991.
[9] H. Whitney. Geometric Integration Theory. Princeton Univ. Press, NJ, 1957.
[10]J. Nedelec.
\320\241 Mixed finite elements in r\\ Numer. Math., 35:315-41, 1980.
[11] M. [Link] element of dielectric-loaded
analysis waveguides. IEEE Trans.

Microwave Theory 32:1275-1279, October 1984.

Tech.,
[12] G. Mur and A. T. de Hoop. A finite element method for computing three-
dimensional electromagnetic fields in inhomogeneous media. IEEE Trans.
Magnetics, 21:2188-2191, November 1985.

[13] J. S. van Welij. Calculation of eddycurrents in terms of H on hexahedra. IEEE

Trans. Magnetics, 21:2239-2241. November 1985.
[14] M. L. Barton and Z. J. Cendes. New vector finite elements for three-dimen-
three-dimensional
magnetic field computation. J. Appl. Phys., 61(8):3919-3921, April

1987.

[15] J. M. Jin and J. L. Volakis. Electromagnetic scattering by and transmission

through a three-dimensional slot in a thick conducting plane. IEEE Trans.
Antennas Propagat., 39D):543-550, April 1991.
[16] J. M. Jin and J. L. Volakis. Scattering and radiation from microslrip patch
antennas and arrays residing in a cavity. IEEE Trans. Antennas Propagat.,
39:1598-1604, November 1991.
[17] J. F. Lee,D.K. Sun, and Z. J. Cendes. Full-wave analysis of dielectric wave-
waveguides using tangential vector finite elements. IEEE Trans. Microwave Theory
Tech., MTT-39(8): 1262-1271, August 1991.

[18] P. Monk. A finite element method for approximating the time-harmonic

Maxwell equations. Numer. Math., 63:243-261, 1992.
64 Shape Functions for Scalar and Vector Finite Elements \302\246
Chapter 2

[19] D. H. Schaubert. D. R. Wilton, and A. W. Glisson. A tetrahedral modeling

method for electromagnetic scattering by arbitrarily shaped inhomogeneous
dielectric bodies. IEEE Trans. Antennas Propagat., pp. 77-85, January 1984.
[20] D. R. Tanner and A. F. Peterson. Vector expansion functions for the numerical
solution of Maxwell's equations. Microwave Opt. Tech. Lett., 2B):331-334.
1989.
[21] R. D. Graglia, D. R. F. [Link]}
Wilton, and A.
vector bases for computational electromagnetics. IEEE Trans. Antenm

Propagat., pp. 329-342, March 1997.

[22] J. S. Savage and A. F. [Link]-ordervector finite elements for tetra-
tetrahedral cells. IEEE Trans. Microwave Theory Tech., 44F):874-879,June 1996.

[23] A. Bossavit. Whitney forms: A class of finite elements for three-dimensional

computations in electromagnetism. IEEE Proceedings. 135, pt. A(8), November
1988.

[24] J. F. Lee, D. K. Sun, and Z. J. Cendes. Tangential vector finite elements for

electromagnetic field computation. IEEE Trans. Magnetics, 27E):4032-4035.

September1991.
[25] T. Ozdemir and J. L. Volakis. Triangular prisms for edge-based vector finite

element antenna analysis. IEEE Trans. Antennas Propagat., pp. 788-797, Ma\\

1997.

[26] J. S. Wang and N. Ida. Curvilinear and higher order 'edge' finite elements in

electromagnetic field computation. IEEE Trans. Magnetics,29B):1491-1494,

March 1993.

[27] G. E. Antilla and N. G. Alexopoulos. Scattering from complex three-dimen-

three-dimensional
geometries by a curvilinear hybrid finite element-integral equation
approach. /. Opt. [Link]. A, 11D): 1445-1457, April 1994.
[28] J. P. Webb and B. Forghani. Hierarchal scalar and vector tetrahedra. IEEE
Trans. Magnetics, 29B): 1495-1498,March 1993.
Overview of the Finite
Element Method:
One-Dimensional
Examples

3.1 INTRODUCTION

The finite element method (FEM) belongsto the class of partial differential equation
(PDE) [Link] origin is frequently traced to Courant [1] who in the 1940s first
discussed piecewise approximations in the appendix of his paper. In the 1950s,
Argyris [2] began putting together the many mathematical ideas(domain partitioning,
assembly, boundary conditions, etc.) that comprise the FEM for aircraft structural
analysis. The introduction of FEM to the engineering community occurred in the
1960s, and some feeJ that the conferences on finite elements held in 1965, 2968, and
1970at the Wright Patterson Air Force Base in Dayton, Ohio, U.S. played an
important role in advancing the method. Finite element activity in electrical engin-
engineering also began in the late 1960s with the papers by Silvester [3] (see also the
reprints volume [4] and Arlett, Bahrani and Zienkiewicz[5])addressingapplications
to waveguide and cavity analysis. Later developments on absorbingboundary con-

conditions, perfectly matched absorbers and hybridizations with boundary integral

methods have led to the successful application of the FEM to opendomain problems
in scattering, microwaves circuits, and [Link]'smain advantage is its
capability to treat any type of geometry and material inhomogeneity without a need
to alter the formulation or the computer code. That is, it provides geometrical
fidelity and unrestricted material treatment. Moreover, the application of the
FEM leads to sparsematrix systems which can be stored with low memory require-
requirements when iterative solvers are employed for the solution of these systems. We
typically state that the FEM systems have 0{N) storage requirements, implying

65
66 Overview of lhe Finite Element Method: One-Dimensional Examples \302\246
Chapter 3

that the memory needed for a solution of an FEM system is proportional to the
number of unknowns N. For most casesthese memory requirements may range from
ION to 40iV depending on the type of problemconsidered and the employed basis or
expansion functions approximating the field within the computation domain. This is
in contrastto boundary integral solutionswhich lead to fully populated systems

having O(N2) storage and OiN*) CPU requirements. However, it should be pointed
out that the number of unknowns for boundary integral equations are generally
much less than those of FEM for the same problem. Nevertheless, when dealing
with nonmetallic structures, the FEM and its hybrid versions is the most attractive
choice.

3.2 OVERVIEW OF THE FINITE ELEMENT METHOD

The geometrical adaptability and low memory requirements of the FEM have made
it one of the most popular numerical methods in all branches of [Link]

application to boundary value problems [6] involves the subdivision of the computa-
computationaldomain (region where the fields are to be determined) into smaller elements [7].
[8]. For two-dimensional these elementsare typically
problems, triangles or quadri-
quadrilaterals as discussed in Chapter 2 and illustrated in Fig. 3.1. Additional example
meshes are given in Figs. 3.2 and 3.3 with the latter referring to a three-dimensional
mesh around a sphere.
The subdivisionof the domain into small elements is referred to as meshing or

discretization of the geometry and is an important part of the FEM solution pro-

procedure. By keeping the elements small enough (typically less than 1/10 of a wave-

wavelength per side), the field interior to the elementcan be safely approximated by some
linear or, if necessary, higher order expansion. The collection of these elements and

their associated expansion or shape function is therefore capable of modeling arbi-

arbitrary and rather complex fields in terms of unknown coefficients which may repre-
represent the field values at the nodes (node-based basis)or the average field values over
the edges (edge-based basis).
In the FEM, the equations for the unknown
context of the coefficients of the
expansions are
by enforcing the wave
constructed equation in a weighted (average)
sense over each element.A subsequent step involves the application of the boundary
conditions leading to a matrix system of the form

[A][x) = [b] C.1)

where \\b) matrix and is determined on the basis of the boundary con-
is a column
conditions forced excitation (current source,incident
or the field, etc.). The matrix [A] is
square of size N x N. very sparse and typically symmetric unless nonreciprocal
material existsin the computational domain. Its nonzero entries provide the relation-
among
relationship field or voltage of adjacent elements within the computational domain,
and its specific form is a characteristic of the problem geometry and discretization of
the domain. Once the system C.1) has been constructed, its solution proceeds with
the application of an iterative or direct [Link] solvers are primarily used for
large systems (i.e., large numbers of unknowns, N) sincethese solvers avoid explicit
storage of the entire matrix. That is, only the nonzero entries need be stored by
Section 3.2 Overview
\302\246 of the Finite Element Method 67

(*!./!>

Quadrilaterals

(four-sided elements)

V=Q

(8)

(b)

Figure 3.1 Example illustrations of finite clement meshes: (a) shielded strip-
conductor transmission line problem; (b) shielded circular conductor
transmission line problem.

Figure 3.2 Finite element mesh around an airfoil for scattering computations.
[Courtesy of Daniel C. Ross.]
68 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter 3

Figure 3.3 Structured tetrahedral mesh

around a metallic sphere.

employing established storage schemes such as the compressed row and ITPACK
formats discussed in Chapter 9. Direct solvers such as LU decomposition are still

better suited for smaller size systemssincethey require storage of the entire matrix
including its nonzero entries.
The steps involved in the generation and solution of an FEM system can be

summarized as follows:

\302\246
Define the problem's computational domain
\302\246
Choose mesh truncation schemes (in the case of open domain problems)
\302\246
Choose discrete elements and shape functions
\302\246
Generate mesh (prepocessing)
\302\246
Enforce the wave equation over each element (or Laplace's/Poisson's equa-
equation for statics) to generate the clement matrices
\302\246
Apply boundary conditions and assemble element matrices to form the over-
overallsparse system C.1)
\302\246
Ensure matrix symmetry (for domains with reciprocal materials)
\302\246
Choose solver and solve matrix system
\302\246
Postprocess field data to extract parameters of interest (suchas eigenvalues,
capacitance, impedance, insertion loss, scattering matrix, radar crosssection
and so on)

In this chapter we will present these steps for one-dimensionaJ problems before

discussing them for two-dimensional applications in Chapter 4. Although many of

the two-dimensional problems are approximations of three-dimensional ones, they

are nevertheless attractive becauseof their simplicity. Consequently, they can be used
to illustrate the solution procedure without burdens due to geometrical and formu-
Section 3.3 \302\246
Examples of One-Dimensional Problems in Electromagnetics 69

lation complexity. Of course, one-dimensionalproblemsprovide the simplest way to

illustrate the computational approach and belowwe begin with a brief illustration of
the FEM for solving the classic Sturm-Liouville differential equation.

3.3 EXAMPLES OF ONE-DIMENSIONAL PROBLEMS IN ELECTROMAGNETICS

Consider the solution of the differential equation (one-dimensionalSturm-Liouville

problem) [6]

d ( dU\\
-
-r />(.v)-j- + <i(x) U(x) =j\\x) Q<x<xa, C.2)

where p(x), q{x),and f(x) are known functions and U(x) is the unknown field or
voltage quantity. Depending on the interpretation of U(x), this equation can repre-
represent any of the following problems illustrated in Fig. 3.4.

Parallel Plate Capacitor

U{x)= V(x):potential between plates

Boundary Conditions: V@) = 0, V{xa)
= Va,
p{x) = -I, g(x)= O.f(x)= -p/t

Differential equ.: -^V(x) = -- (orV2V = -p/e) C.3)

Wave Between Parallel Plates

U(x) = Ey(x): electric field between plates

Boundary Conditions: Ey@) = = 0.
?,.(.\302\253\342\200\236) p(x) = -\\/fir,
q{x)= kl(n f{x)
= source function

Differential equ.: ~
?,\320\224+ C.4)
(\342\200\224 =/(*)
-^ ^\320\265\320\263?\342\200\236

Reflection From a CoatedMetallic Conductor

pe\321\204endicular polarization
or

H.{x):
IE.(x): parallel polarization

Boundary Conditions:

EAx
= 0) = 0, \342\204\226
+ jlcQE:\\ I - 2jk0e*\302\260*\"
Overview of the Finite Element Method: One-Dimensional
Examples \302\246
Chapter 3

v=va

(a) infinite capacitor problem

4-
(b) Parallel plate

QQ
*
x-t \321\205~\321\205\320\271

Figure 3.4 Problems represented by the one-

dimensional differential equation. Here
Reflection
(\321\201) from a coated metallic conductor = a in (a) and (b).
\320\273\320\263\342\200\236

\320\255\321\205

= 0 or = \320\236
[Link].: J*(l*f)+kfcrE: -f f1 \320\250 + /\321\201^\320\263\320\257. C.5)

In thislatter case, the boundary condition at .v = xu is a consequenceof the fact that

the reflected field must be of the form

C.6)
Section 3.4 The
\302\246 Weighted Residual Method 71

where R is the reflection coefficient of the coated ground plane and is not known
until the FEM solution is completed. Thus, the total field

E. = ?t\"c + or
\320\257\320\223, Hz = + \320\257\320\223
\320\257!\320\277\321\201 C.7)

satisfies the stated boundary condition. As indicated in Fig. 3.4, ?lnc, and H\342\204\242
represent the z components of the plane wave incident upon the metal backed
dielectric slab.

3.4 THE WEIGHTED RESIDUAL METHOD

We proceed now with the solution of C.2) on the basis of the finite element method.
As a first step we introduce the residual

R{x)= - \342\200\224
(p{x) ^-\\ + a(x) U(x) -f(x) C.8)
which must be zero in accordance
of course with the state problem. However, it is
impractical enforce R(x) = 0 at every
to point in the domain from x = 0 to

x = xa = a. Since U(x) is not expected to vary substantially over a small distance,

say \320\224\320\264\320\263,
we subdivide the domain into small segments (as shown in Fig. 3.5) and
instead enforce the condition

I Wm(x)R(x)dx = 0 C.9)
J Domain of Wm

over of the segments. Remarking that

each Wm(x) is some weighting function to be
defined later, C.9) enforces the differential equation on an average senseover the mth
segment. By changing the integration or testing interval (and weighting function Wm)
from m = I up to m \342\200\224
N, we can construct N equations for the solution of the
discretized potential or field values. Before proceeding to do so, wemake the follow-
observations
following about tVm(x):

\"' 2= = N\021
Unknown Uz = U% = U,e=2 Unknown UN-,
= U|=N\" \320\246\"

XN-1

(e=2)

Element or Segment (N -1)

segment #1

Figure 3.5 Tessellation the line

\320\276\320\223 segment 0 < \320\273-
< x0.
72 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter 3

\302\246
If Wm(x) = S(x - xm) or WJx)
- (xm+]
= 8[x + xm)/2], the resulting

weighted residual procedure is referred to as point matching and leads to

a form of the finite difference method.
\302\246
If Wm(x) is set equal to the basis functions used for the representation of

U(x), the procedure is referred to as Galerkin'[Link] is the most

popular testing/weighting method for casting the differential equation to a

linear system.
\302\246
The choice of Wm{x) is not completely arbitrary. For the mathematical steps
in the FEM procedure to hold rigorously, Wm{x) and its derivative must be

at least square integrable over the domain. Specifically, for the problem al

hand, it must satisfy the condition

< oo

In addition, Wm(x) must satisfy conditions at the boundary nodes (end-

points at .y = 0 and x = xu for the one-dimensional problem) which are

compatible with the imposed boundary [Link] smoothness con-

conditions on Wm(x) and dWm(x)ldx may also need to be imposed.
Before generating a linear system of equations from C.9) subject to the bound-
conditions,
boundary it is first necessary to cast it in a more suitable form by following the
steps:
Step 1 Take advantage of the weighting function Wm(x) to reduce the order
of the derivatives contained in R(x). To do so we employ integration by parts, giving

\320\223
duy
\"Jo J*-.v=O

The first right hand side (RHS) term can be evaluated by enforcing the known

boundary conditions at the endpoints. Its effect on the overall system will be con-
considered later.

Step 2 Derive the weak form of the differential equation. The weak form of
the differential equation is most appropriatefor numerical solution and is obtained
by substituting C.10) into C.9). We have

+ \320\257{\320\245) ~ ~
Wm(x) = \302\260
Wm{x) U(X)
Jo\"\\PiX) ^T ? W/'\302\273(V)-/(A')](lx [p{x) Sf]
C.11)

which holds provided C.2) is valid. However, becauseof the integral, the weak form
C.11) enforces the differential equation on an average (and therefore weaker) sense.
Equation C.11) is often referred to as a varialional statement of the problem. What \320\27
remarkable about C.11) is that it incorporates in a single mathematical statement the

requirements imposed by the differential equation and the boundary conditions at the

endpoints. That is, upon substitution of the boundary conditionson U(x) and
Section 3.5 \302\246
Discretization of the \"Weak\" Differential Equation 73

dU(x)/dx, C.11) is not only an alternative statement of the differential equation

C.2), but also includes information about the boundary conditions which are essen-
essentialfor the uniqueness of the solution. This is at the heart of the FEM, and later in
this chapter we will observe how the boundary conditions impact the discrete form of
C.11).

3.5 DISCRETIZATION OF THE \"WEAK\" DIFFERENTIALEQUATION

The discretization of C.11) to a linear set of equations is done by introducing an

expansion for U(x) and then making appropriate choices for the weighting functions

Wm(x). This is an essential step in all numerical solution procedures, and the pre-
previous chapter served to introduce the various classes of basis functions used in a
discrete representation of the unknown function. We choose the linear representa-
representation

C.12)
\320\225

where are
\320\251 the unknown coefficients of the expansion and (seeFig.3.6)
x-x\\
X*\\ < X < \320\233'2

0 otherwise
otherwise 10 otherwise

are the shape functions discussed in Chapter 2. That is.

C.14)

Linear approximation between nodes

Global node # eto segment Local Global (e- 1)th elb element
n=e node# node# element

Figure 3.6 Illusiralion of lhc clh segmentor element and the linear shape functions:
(a) the nodal expansion functions: (h) overlay of the nodal expansion
functions.
74 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter 3

where U\\ and are

\320\251 the unknown values at each node. On the basis of this expan-
expansion, the field or potential over the domain 0 < x < xa is determined by a linear

combination of the field or potential values at each node. Although not necessary for
this one-dimensional example, two types of node-numbering schemesare typically
used to facilitate the programming and implementation of the finite element solution.

Local Node Numbers. These are assigned node numbers unique only within a

single element. For the line segments in Fig. segment is formed by two nodes
3.5 each

having the local node numbers 1 and 2. as illustrated in Fig. 3.6. That is, the notation

x\\ refers to the location of the local node 1 of the <>th element. Similarly, refers
\320\251 to
the field or potential at node I of the eih element. Using this type of notation, we can

develop formulations and equations for a single element which can then be incorpo-
incorporated into the overall solution by attaching the superscripte to all local or elemem
variables. In this manner, the uniquenessof the equations is maintained when com-
combined with those from the other [Link] elements (two-nodesegments for one-

dimensional problems) in the computational domain are assumed to have a unique

number from e = 1 to e = Nc.

Node
Global Numbers. Each node of the discretized domain is also given a
unique number from 1 to N, as shown for example [Link] 3.5, The assignment of
these global numbers is necessary since eventually all unknowns from each element
must be collected (a process referred to as element assembly) into a matrix system

such as that in C.1). At that stage, it is necessary to maintain a single subscript or

array dimension. This necessitatesthe correspondence between the notation \320\273)'
and

where
\321\205\342\200\236, first refers to the local numbering
the notation and hereon the single
subscript be
will reserved for the global numbering notation. As can be realized,
since every node in Fig. 3.5 belongs to two elements, multiple local notations can
refer to the same global node or field value. As an example (see Fig. 3.6),the node
location is identical
\320\273\342\200\236 to the locations implied by the notations x\\ and .v>\"'.
Likewise, for the held values we can state that ?/,,
= V* = ?/-Tl a\"d so on.
When the expansion C.14) is substituted into C.11) we get
f
\320\233', 2 [Link]

,.=1 [ ,=j J-n

- -fix)
\302\273\302\246\342\200\236(*>
/\321\204\320\263> Wln(x)
= 0 C.15)
^ ^ v=0

The latter terms in the brackets are due to contributions from the endpoints of the

domain and their evaluation is subject to the specific boundary conditions. This

equation now explicitly shows how the boundary conditionsenter into the construc-
of
construction linear system. Hereon, we will refer to their contributions
the as [endpoints]
since we have not yet specified the type of boundary condition to be imposed.
We are now ready to make different choices for the weighting function to

generate a system of linear equations for the solution of \\Un). As stated earlier,
this step is also referred to as testing and Galerkin's method is usually employed
Section 3.5 Discretization
\302\246 of the \"Weak\" Differential Equation 75

in the finite element method. Specifically, we choose Wm(x) =

Nf(x) and for each of
these testing or weighting functions a single linear equation is generated. From C.15)
we have

+ [endpoints] = 0 C.16a)

where

[endpoints]= endpoints | +endpoints2

.v=0
*=*.J

~p{x) C.16b)
~a\\
.\321\202=0 x=xd

Since Nj(x) is nonzero only over the eth element, the summation over the elements
can be eliminated at this stage. In other words, integration is carried out only over
the nonzero portion of the integrand. We can rewrite C.16a) in matrix form as

+ fendpoints] = C.17a)
[A'ji\\{Uf)

which can be considered as the weighted discrete form of the differential equation.
This matrix system provides a relationship only between the two nodes forming the
eth element a localizedrelationship among the
and is therefore node fields/potentials.
The endpoint contributions appear only when e = 1 or e = Ne and vanish when the

Neumann boundary conditions^

~
=^ = 0 are imposed. OtherwiseC.17a)
becomes

This is commonly referred to as the elementalmatrix system. Also, the matrix

C.18)

is referred to as the element matrix and its entries are given by

C.19)

The excitation vector entries are specifiedby

=
ri{x)f(x)dx C.20)

Since Nf(x) are linear functions, the evaluation of can be carried

\320\220\320\265\321\203
out in closed
form provided p(x) and q{x) are taken as constant over the integration of the e\\h
element. Specifically, setting p(x) \302\253*
//' and #(.v) = q\" for x\\ < x < x'i, we find that
76 Overview of the Finite Element Method; One-Dimensional Examples \302\246
Chapter 3

Also, on approximating f(x) by f(x)

\302\253=
f over the eth element, we obtain

The above testing procedure will result in 2Ne equations obtained by letting
e = 1,2 Ne in C.17). Since only Ne + 1 unique unknowns exist, it is necessary \321\216
condense or combine the 2Ne equations down + 1. The additional set of
to Ne

equations is a result of testing at the same from the left using the testing
nth node
function N2~\\x) and from the right using the testing function N'(x). Their reduction
to Ne + 1 equations is referred to as assembly of the element equations and is a
standard step in all finite element solutions.

3.6 ASSEMBLY OF THE ELEMENT EQUATIONS

The essence of the assembly procedure is to take the average of the test equations
from the left and right of the node.
\302\253th That is, we consider the weighted average.

dx + N\\{x) R(x)
\320\251-1(x)R(x) dx\\=0

\"'*'
\320\223J Tm(x) R{x) dx = 0, m = 2,3 N - 1 C.24)

where Tm(x) is illustrated in Fig. 3.6(b) and is given by

= ~x
Tm(x) x <x < C-25)

0 otherwise

From C.24), the test equation at the with node has the explicit form

= *m. m = 2,3 N-\\ C.26)

?/Ut/\302\253
Section 3.6 \302\246
Assembly of the Element Equations 77

where

1
/ ,dTmdTn
= \320\223\"+'\320\223 -\342\200\224\302\246
p(x)
\320\233\342\200\236\321\210 +
\342\200\224r- q{x) Tm{x) \320\242\342\200\236(\321\205)
dx
Jxa-X dx dx J
A^\\\" + &\320\221.
^
\320\277
\342\200\242 w

\320\273
12 . n = m \342\200\224
1 l->^',)

/4^\"' . n = m + 1

'
f-v\302\253+i

*m
= = \320\271\320\223\321\210\"'
\320\233
\320\243\320\266\320\270)\320\233\321\205) + *\320\223\" C.28)
I
Jtm-I

We have temporarily excludedthe cases of m = 1 and m = N since testing at the

boundary nodes 1 and N requires the inclusion of the endpoint contributions. These
must be dealt with after the specification of the boundary conditions. Thus, C.24)
gives a set of N \342\200\224
2 equations with the other two equations to be supplied later.

When the Dirichlet boundary conditions U{x= 0) = U(x = = \320\236

are
\321\205\342\200\236) imposed at
the endpoints of the domain, the implication is that U\\ (or Ux) = UN (or ?/^')
= 0
and thus the number of unknowns is reduced from N down to N - 2. Consequently,
when U(x) satisfies the Dirichlet boundary conditionson both ends of the domain,
the system C.24) is sufficient for the solution of the unknown vector [Un] provided
the node fields or potentials are set to zero a priori.
In practice(i.e.,in a computer implementation of the FEM), the construction
of the linear system is done by manipulating the element matrix system C.17). This
procedureappears rather tedious on paper because of its repetitiveness but is most
suitable for computer implementation. We will illustrate it below for the one-
dimensional formulation and in connection with the three-element segment given
in Fig. 3.7.
For convenience, let us assume that the Neumann boundary condition is satis-

satisfied at nodes 1 and 2, i.e., \320\251 =\320\251

=0. Thus, the element matrix equations

C.17b) even when testing at elements 1 and

are applicable 3. More specifically, for
e=\\, C.17b) becomes

C.29)
'22jl^2j l\022j

Also, for e = 2 and t' = 3, we have

'*!! (\320\260\320\264

\320\270,

1 4, 2 \\ 3 \\ 4

Figure3.7 Three-element tessellation of a ' e=1 e=2 e=3 '

line scgmenl.
x~ 0 x- xa
78 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter i

and

C.31)
b\\

These are six equations for the four unknown coefficients U\\. U2, U$< and C/4. To
reduce C.29)-C.31) down to four add those which
equations, we simply
correspond
to testing from the left and right This addition of equationsamounts
of the node. to

simply performing the first sum in C.15) and is not an arbitrary decision in the

process. For nodes 1 and 4, there is only testing from one side and therefore the
first equation of C.29) and the second equation of C.31) are left unchanged. For
node 2, we add the second equation of C.29) and the first equation of C.30).
Likewise for node 3, we add the second equation of C.30) with the first equation
of C.31). The resulting (i.e., assembled) system of four equations is

\320\220\\\320\263 0 0 t/,

22 + \320\2202\320\270\320\220\\\320\263 0 \320\2702
C.32)
\320\276 Ah + \320\220]
\320\220222 \320\260]2

\320\276 0 42- \320\2704 b\\

where we employed the global node notation for the {U} vector. In placing this into
the compact notation

[A][V\\
= C.33)

we note that [A] is a tridiagonal matrix regardlessof the number of elements used for
the tessellation of the line segment 0 < x < xa. That is, a maximum of three nonzero
entries appear in each row of [A] and except for the top and bottom row, the
diagonal entry and its adjacent elements are the only nonzero entries. We can
state that the bandwidth of the matrix is three regardless of the number of elements.
Simply put. as the number of nodes/elementsincreases, matrix takes the generic
t he

sparse form

X X 0 0 . ..
X X .V 0 0 . \342\200\242
.

0 X .V x 0 0 ...
0 0 JC X X 0 0
C.34)
0 0 0 X X X 0 0

0 0 X X X 0

0 0 .v .v

\" \"
where the \"x\" symbols imply a nonzero entry and the . . . denote a continuation
of the zero entries, in applied mechanics [A] is referred to as the stiffness matrix
because of the similarity of C.33) to the equation Kx =f for the deflection* of a
linear spring with stiffness under
\320\232 an applied force/. In electromagnetics. [A] can
be interpreted in several ways depending on the physical quantity represented by

U{x). For example, if U(x) = voltage or electric field, then [A] can represent an
admittance matrix with {/>} being the electric current excitation. Alternatively, if

U = electric current or magnetic field, then [A] may be compared to an impedance

matrix.
Section 3.7 Enforcement
\302\246 of Boundary Conditions 79

Becauseof the sparsity of [A], only its nonzero entries are stored when solving
the matrix system C.33). Also, global numbering as used in C.26)-C.28) must be
employed in defining the assembled matrix system. Clearly,the matrix system in
C.32) is identical to that in C.26H3.28) except for the first and last row of C.32)
which correspond to the omitted m = 1 and m = N equations in C.26). In practice,
the assembly of [A] is done by employing C.27) and C.28)directly or by implement-
a double
implementing loop and keeping track of the correspondence between the local and
global numbering. A possible double loop which generates the nonzero entries of
[A] is

Initialize
\320\241 the [A] matrix to zero
DO Wm = l,N
DO 10 \320\270
= I./V

10 A(m,n) = 0.0
\320\241
Loop through all elements and construct [A]

DO20e= [Link]-1
\320\241
Compute element matrix [A*]
DO 30/= 1,2
DO 30./= 1.2
30 Compute AE(iJ) from equations C.21) and C.22)
Assemble
\320\241 [A'] into global matrix

A(e+ l,e + l) = A(e+\\.e+l) + AEB,2)

A(e.e+l) = AE(l,2)
20 A(e +l,e) = AEB, I)

3.7 ENFORCEMENT OF BOUNDARY CONDITIONS

So far we have postponeda comprehensive discussion on the imposition of boundary

[Link] is the case with any differential equation, a unique solution can be
obtained only after the specification of boundary conditions which constrain the
values of the field/potential at the boundaries (endpoints for the one-dimensional

case) of the domain. These boundary conditions(also referred to as boundary con-

constraints) come in various forms. In their most typical form they provide a specifica-
of
specification the field at the end nodes 1 and N or a specification on the value of the
normal derivative of the field. However, the boundary condition may simply provide
a relationship between the normal derivative and the field at the boundary node(s).
The derivative of C/(.v) at the boundary nodes shouldnot be approximated by using
the expansion C.12). This representation is an interpolation function between the

boundary nodes, but the fields/potentials and their derivatives at the boundary nodes
must be independently specified for a unique solution of the differential equation. As
a reference, we remark that if the spatial variable .v in C.2) was replaced by the time
variable i, the boundary conditions would become the initial conditions of the tem-
temporal response U(i).
80 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter 3

The various conditions to be encounteredin the solution

boundary of differ-
differential equations are with a well-establishednomenclature. Belowis a list
associated
of the boundary conditions typically imposed at the boundary nodes. The type of
boundary condition to be used will depend on the physics of the problem, as dis-
discussed in Section 3.2.

3.7.1 Neumann Boundary Conditions (Homogeneous)

For one-dimensional problems, this boundary condition states that

at the left or right endpoint of the domain. In two- and three-dimensional problems,
it can be stated as

n Vf/ =
\342\200\242
^oft
=0 on S or \320\241 C.36)

where n denotes the outgoing unit normal vector of the domain boundary, as
illustrated in Fig. 3.8. In acoustics this is referred to as the hard boundary con-

condition, and in electromagnetics the magnetic field obeys this condition on metallic
boundaries.
The Neumann boundary condition is the easiest to be numerically enforced in
FEM solutions. In this case, the [endpoint] contributions C.16b) vanish and the
elemental equations C.17b) lead to the global system C.32) without any special
considerations.

3.7.2 Dlrlchlet Boundary Conditions(Homogeneous)

The Dirichlet boundary condition specifies a vanishing field or potential at the

endpoints of the computational domain, i.e.,

U(x)= 0 at endpoints C.37)

In acoustics this is referred to as the soft boundary [Link] two- and three-
dimensional electromagnetic problems, the Dirichlet boundary condition is satisfied

by the tangential electric fields on all metallic surfaces.

When the Dirichlet boundary condition must be satisfied in connection with
the elemental equations C.17a), the [endpoints] are not zero and must be considered
in the assembly of the final system. Since the [endpoints] contribute only to those

Figure 3.8 Illustration of the enclosures for

\320\241
(enclosing contour) S (enclosing surface) tw0. ;ind three-dimensional domains.
Section 3.7 Enforcement
\302\246 of Boundary Conditions 81

elementequationsresulting from testing at nodes 1 and N, the global system C.32)

for the three-element tessellation (see Fig. 3.7)takes the form

a\\2 0 0 U\\ endpoint, h\\

Av
A A-,-, + A\"\\\\ Ah 0 U2 0 b2 H- b\\
4-
0 \"Ah An + A11 Uj 0 b2 -f- b'\\
0 0 A2l fAn. 4 endpoinh
C.38)

However, U\\ = C/4 = 0, as dictated by C.37), and thus N 2 =

\342\200\224 2 unknowns (f/2 and
l/j) remain to be determined of C.38).This
from the solution allows us to discard the
first and last equations of C.38) provided we set U\\ = f/4 = 0 wherever they occur.
In doing so, we have
1'2
\320\233 C.39)

which can be inverted to yield

U2\\_ 1

Thus, even though the [endpoint] contributions are not computable, we can still
solve for the node fields or potentials. We remark that if C.38) was a system of N
equations, the reduced system C.39) will consist of N \342\200\224
2 equations after the en-
enforcement of the Dirichlet boundary conditions.

3.7.3Nonzero Boundary Constraints (inhomogeneous)

We may also encounter situations when the field or potential at the end node is
assigned a specifiedvalue. An example of this situation is the parallel plate capacitor
problem where the upper plate has a potential equal to Vu. As a more general
example, let us consider the situation where in reference to the three-element example
in Fig. 3.7, we set

= 2o. ^7 C.40)
which are typically referred to as inhomogeneous Dirichlet and Neumann boundary
conditions, respectively.
Substituting these values into the system C.38) gives

/4J, A\\-> 0 0 1 Qo 1 endpoint,

4i + \320\260\320\252
\320\273\320\252 A]2 0,1/2. 0
\302\246\342\200\242*
0 Ah A222 + A*n A]2 U} 0

0 O\" A\\i /422 J t74 J Q\302\273

Q,, = = (/, is already specified,the

where \342\200\224 since
Q,,p(xa) N2\\xa) \342\200\224Q,,p(xa). Clearly,
first of the equations can be [Link] doing so and rearranging, we obtain the
system

A 12 0
A\\\\ \320\220\\\320\263 Al C.41)
0 \302\246L b]
82 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter]

which can be solved for U2, V-i and U4. We again point out that, when dealing with
N nodes, C.41) would be a system of N - 1 equations and with the exception of the
first and last equations, the test will have three nonzero elements as illustrated by tk
matrix C.34).
The procedure of reducing the system C.38) to C.39) or C.41) is often referred
to as condensation of boundary [Link] reduction is typically performed
during the assembly process by eliminating for example the rows which test at a

boundary node assigned a specified value. Finally, we remark that the condensation

process modifies the excitation column implying that the specification of a potential

or a field value on a boundary is equivalent to a source excitation. Thus, the excita-

excitation of the domain can either be specified through a nonzero /(a) or through the

enforcement of nontrivial boundary constraints.

3.7.4 ImpedanceBoundary Conditions

The impedance boundary condition provides a relationship between the field

and its normal derivative. Referring to Fig. 3.8, it is typically stated as [9]

+ at/ = 0
^\320\264\320\277 on S or \320\241 {\320\252\320\2

where or is a constant. This boundary condition has been found very useful in

modeling the
presence of thin dielectric coatings without a need to tessellate tk

region interior to the dielectric. In finite element simulations, C.42) also plays the
role of the radiation condition or a first order absorbing condition (to be discussed
later). These boundary conditions will be discussed extensively in later chapters and
are used for truncating the computational domain of open domain problems as in
the case of scattering by an airfoil (see Fig. 3.2). They basically provide a statement
on the field behavior at the boundary [Link] need to mesh the region beyond the

boundary enclosure is therefore eliminated provided the absorbing boundary con-

condition gives the proper field behavior beyond the boundary enclosure.
A generalization of C.42) is

+ aU = f) on 5 or \320\241 C.431

where or and j8 are constants and we can refer to this as the inhomogeneous boundun
condition. The treatment of C.43) is no different than that for C.42). For the one-

dimensional problem considered here, C.43) reducesto

\342\200\224
+ aaU = x =
\321\200\342\200\236.
= a
\321\205\342\200\236 C.44|
ox

When these conditions are used in the finite element solution of the three-elemeni

segment example (see Fig. 3.7), the resulting system is again of the same form as

C.38). From C.16b). we have

[endpoint], =/?@)(A,-or0f/|)
[endpoint],
= -p(x,,)(fiu - a<tU4)
Section 3.8 \302\246
Examples 83

and when these results are incorporatedinto C.38), after rearranging, we obtain the
system

0 0
0
0 V,

-\320\233/\320\232\320\236)
\320\276
\302\246
+ C.45)
\320\276
b\\ + b\\

bl +A,/K-va)

This systemcan now be solved for [U] and we remark that the middle two equations
are unchanged by the imposition of the impedance boundary [Link] is clear
that when iV equations are involved, the imposition of the impedance boundary
condition will only alter the first and last equationsof the overall system.

3.8 EXAMPLES

To illustrate the application of the various steps in constructing and assemblingan

FEM system, let us consider two simple examples.
EXAMPLE 3.1

Solve numerically the differential equation

2W + n2 ?,.(*) = sin
2\320\2732 nx, 0 < x C.46)

using ten segments and subject to the boundary conditions

?v@) = ?,,(!) = 0 C-47)

Compare the numerical results with the exact solution [10, Chapter 11]
Fv.(.v)
= sin;r.v C.48)
Solution
This problemcorrespondsto a source in a parallel plate waveguide as illustrated in Fig. 3.4(b).
By comparison to the general form of the differential equation in C.2), we note that

p(x) = 1. q(x)= ?r, j\\x)

= 2jC sin jr.v

Also, from Fig. 3.9 and the formulae C.21) and C.22), we find that

Aeu = 10.328987
=^2=-^-

and

. = ,*,=\342\200\224L +\321\217*\320\224*=_9.\320\2305507
\320\224\320\264- 6
84 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapier!

(x=0.

\\ e

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

\320\232\342\200\224H

Ax=0.1,N=#ofnodes = 11
Figure 3.9 Tessellation of a line segmra;
Ne = # of elements = 10
into ten equal length dements.

Thus, the elemental equations are given by

10.328987
\320\223 -9.835507 II ?J, 1 _ I b\\ \\
C.49t
[
-9.835507 10.328987 J | \320\251\320\263
\\~ \\ b\\ \\

where from C.20)

\302\246>
\320\223*/**-.
lit1 = 2\320\2732 x.
f ~
x\\

with .r' = (e - 1) and x$

\320\224\321\205
- e Ax for e = 1,2,..., 10.
After assembling the elemental equations C.49) subject to the Dirichlet boundary con-

conditions, we get the 9 x 9 system

20.6580 -9.8355 0 0 0 0 0

-9.8355 20.6580 -9.8355 0 0 0 0 0

0 -9.8355 20.6580 -9.8355 0 0 0 0

0 0 0 -9.8355 20.6580 -9.8355

0 0 0 0 -9.8355 20.6580

Ey2
Eyi h
EyA b,
=

\320\272

with

{/>}r = {0.605, 1.1507, 1.5838,1.8619,1.9577,1.8619,

1.5838,1.1507.0.605)
This system in form to that in C.26).
is identical The node fields Eyl and ?,.|, were noi
included they are zero as dictated
since by the boundary conditions. Also, we remark that
the symmetry of the system (i.e., Aq = Aj,) is not [Link] problems wiik

reciprocal permittivity and permeability tensors are inherently reciprocal, and this properly\302\273
exhibited in the Hermitian form of the matrix.
Section 3.8 \302\246
Examples 85

Solving the above system via matrix inversion, for example, gives the following results:

Nodc# fi\342\204\2424 ??\"\342\200\242 Error = |?fEM - ??\"ct|

1 0 0 0
2 0.3103 0.3090 0.0013
3 0.5902 0.5878 0.0024
4 0.8123 0.8090 0.0033
5 0.9550 0.95II 0.0039
6 1.0041 1.0000 0.0041
7 0.9550 0.9511 0.0039
g 0.8123 0.80\320\255\320\224 0.0033
9 0.5902 0.5878 0.0024
10 0.3103 0.3090 0.0013
I! 0 0 0

As seen, the difference between the exact and numerical solution is in the third decimal place,
indicating that the employed number of elements are sufficient for anaccurate representation
of the field distribution. To reducethe solution error, more the elements can be used for
tessellation of the line segment. However, as \320\224\321\217 0, the system condition \\\\.s/\\\\ ||.\302\253/~'||,
\342\200\224>

where is a natural norm of the

||.\302\253/|| matrix [A]. Eventually the error would not reduce with

increasing N. unless the machine precision is also increasedwith the inclusion of additional
decimal places in the calculation of the matrix entries and in carrying out the system solution.
Numerical precision is particularly important for solving problems with many thousands of
unknowns (the case with many problems).
practical
Having the node fields, the field is found from C.12) or C.14). Specifically, since

x\\ = 0.1(e [Link], we have

\342\200\224
1) and x% =

Ey(x)
= 10 -
\321\201
?[?\342\200\242;,@.1 x) + - [Link]
\320\225*.\320\263(\321\205 + 0.1)]P^Jx -2e + 0.05)
9

0.1

where

x < A.v/2
otherwise

is a pulse function.
EXAMPLE 3.2

The field reflected by a metal-backed dielectric slab due to a plane wave excitation E'J\" = e'*\"*
is given by

Er. = Re'\021\"

where R is the unknown reflection coefficient to be determined.

The dielectric slab is of thickness t, as shown in Fig. 3.4(c). Consider the case of
f = 0.25A0, t, = 4 \342\200\2247/3,=
\320\224, 1 and employ the finite element method to compute the reflection
coefficient as a function
\320\257 of the loss parameter fi. Compare the result with the analytical
reflection coefficient given by
86 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter}

where

and

are the wave impedances in free space and in the dielectric medium, respectively.

Solution
As discussed al the beginning of Section3.2,the pertinent differential equation is

d2E.,ir. _

subject to the boundary conditions

We will set xa = t + At, where Al is the chosen element length, as illustrated in Fig. 3.10. The

choice of \320\224\320\263
(the element size) should be somewhere between Ad/10 and Ad/20, where
Ad = 2jr/Re(ftd) is the wavelength in the material. For our case, Ad = and
B\320\267\320\263/\320\2200)\321\20

Ad = ko/y/Re(trttr), and we will choose \320\224/

= 0.025A0 to ensure that this criterion is satisfied
Since t = 0.25A0, this choice of At leads to Nv = 11 layers or elements, including the air-filled
element occupying the region t < x < t + At. We also observe that for \320\224\320\263= iIk
\320\236.\320\23625\320

sampling rate is 20 elements per wavelength in the dielectric since ReFr) = 4. From |3.2ii
and C22) the entries of the elemental equation are given by (with / = 1,qc = -k\\tK)

= Ah = -L - = 40.- 0.32899e

At -
J_ = -40. 0.164493<n
At

-?,= 0

\302\246

Figure 3.10 Tessellation of the dielectric sb

using N,.
= N - I layers (elements).
Section 3.8 \302\246
Examples 87

where e,c denotes ihe relative of the eih element and we have suppressed
permittivity the
presence of the factor Ao the free-space wavelength
since cancels out in the final result.
After assembly, we will obtain 11 equations since ?-t = E:(x = 0) = 0. Except for the last
and first, all other equations of the assembled system will be of the form

[0.0 Aw.i)m. Amm, Am{m+l), 0 0]|?.} = 0

where {?,lT = |?.:, ?., ?.12(is the unknown node Gelds vector. Anm
= A\",\\+ A'S\021'1.
=
-4-j'f .
\320\220\321\210-1\\\321\202 and = The
\320\220,\321\204\342\200\236+\\)
\320\220\"^- last equation of the system is obtained by recogniz-
that
recognizing the boundary condition at * = t + Al is of the form C.44) with aa=jkQ and
ft,
= 2/*0<J A'\"-Thus.

[endpoint], = +7*(,?,,i- 2jk0eJk\"u+Al)

and the tenth equation of the assembled system will be

+Jko)E:l2 =
That is, b\\2 = 2jk()e'*|>1'+\320\224|> js the only nonzero entry of the excitation column and is a result of
the assumed plane wave incidence.
Upon solution of the assembled FEM system, the reflection coefficient is extracted from
the node field values as

E2U
\342\200\236hem
- E^(x =t+ At)
?Inc(.v = I + At)
A plot of R as a function of the loss parameter /3
= Im(er) is given in Fig. 3.11 (courtesy of
S. Legaull).
It thai the computed values
is seen of RFEM are in good agreement with the exact result
even when the discrete layers in the dielectric are reduced from ten down to only N,. - 1 = 2
(corresponding lo a sampling rate of eight elements per \\d). In this case the decay of the field as it
propagates in the dielectric is very rapid and the employed discretization of the dielectric slab is
not sufficiently fine to pick up the large changes in the field values from one discreteslab
(element) to the next. This is a good example of the fundamental assumptions made when
diseretizing the computational domain. For large jS, the field decrease is slower and thus less
samples are neededto approximate the field within the region while maintaining the same level
of accuracy.

\342\200\224
Exact
\342\200\224
Elements in Dial. = 10
-
\342\200\224
Elements in Dlei. = 5
Elements in Diel. = 2

DC 0.7

[Link] Plot of ihc reflection coefficient

normal
(\302\253I incidence) from a grounded dielcc-
incblab of thickness / = = 4 -iff, 0.4
0.25X,, <<=,
it, = I): comparison of the numerically com-
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Rfm with the exact result
computed \320\233\"\".
88 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter 3

EXAMPLE 3.3

Repeat Example 3.2 assuming that and with

\320\265\320\263=\320\264\320\263=4-\321\203/* all other parameters fcli

unchanged.

Solution

By selecting \321\206,
= er, the impedance of the wave inside and outside the dielectric is

^ = 120jt

and thus the wave does not exhibit any reflection at x = t (i.e., at the dielectric interfacci

The interface is then referred to in the literature as reflectionless. Also, as it enters ilw

dielectric, the wave is absorbed due to the nonzero imaginary components of the conslituuvi;

parameters. If p is chosen so that the wave has decayed to negligible levels by the time it

exits into the air medium, the dielectric layer can be considered as \"perfectly absorbing'
since no reflected field is returned. Such layers have been proposedrecently [11], and it ha>
been shown that certain anisotropic layers can lead to perfectly matched interfaces for
incidences away from normal. These types of layers are important in finite element simula-
simulationsbecause they can be used for simulating a nonrenecting surface. The latter is essential
in solving open domain problems, as is the case with scattering and antenna radiation

problems considered later.

For the one-dimensional problem considered here, we will restrict ourselves to ,t
\342\200\224
simple isotropic dielectric =
layer having \302\253r/i, 4 \342\200\224jfi.
The numerical results for a layn
of thickness / = 0.25\320\233.\320\276
are given in Fig. 3.12 (results are courtesy of S. Legault) as a
function and for three
of p different sampling rates {Nt = 6.11,21). It is again seen ihii

as fi increases, the numerically computed reflection coefficient begins to deviate from ihc

analytical/exact result given by R = -\342\200\242\">\"\302\246

[12]. This is expected and is due to the rapid

decay of the field within the absorber as fi is increased. Thus, higher sampling is needed \321\213
belter model the field values from one discrete layer to the other.

-10

-20
Elements In Diet. =5

solution |
[\342\200\242\342\200\224Exact 4^ Elements In piel.= 10
-40

\302\246\\Elements in Diel. = 20
\342\200\22450

Figure 3.12 Plot of the reflection coeffiacm

-60 n_ 1
. (at normal incidence) from a grounded die!\302\273
trie slab of thickness / = 0.25A,, having
Appendix 1 \302\246
Sample One-Dimensional FEM
\320\274\320\260\320\277\320\273\320\262
Analysis Program 89

APPENDIX 1: SAMPLE ONE-DIMENSIONAL matlab FEM ANALYSIS PROGRAM

%This MATLAB code can be used to reproduce the data in Fig. 3-11
\\ CoUTtesy Of LARS ANDERSEN

%t=thickness of slab
%kO=free space propagation constant
%xa=location of the left
computational domain endpoint
%alphaa=alpha coeff. to be used in the boundary condition at xa; see C.43)
%betaa=betaa coeff. to be used in the boundary condition at xa; see C.43)
%p=p coefficient appearing in the differential equation C.2)
%epr=relative permittivity of the slab
%beta=as defined in Fig. 3-11
%q=coefficient appearing in the differential equation C.2)
%N=number of nodes (N-l=number of layers)
%R_FEM=computed reflection coefficient (to be plotted)

% Initialization

clear;
N=7;
t=0.25;
Dx=t/(N-2)(
xa=t+Dx;

P=-l;
kO=2*pi;
f=0;
alphaa=j*KO;
betaa=2*j*kO*exp(j*kO*xa) ;
%values to generate plot in Fig. 3-11
Z=26;

4epsr=4-j*beta;
for 2=1: 7,

beta=(z-l)/(Z-l)*5;
epsr=4-j* beta;

% Initialize global matrix and vector

for m=l:N,

b(m,l)=0;
for n=l:N,
A(m,n)=Os
end

end

% Compute element matrices and assemble global matrix

for n\302\253l:N-l,

xel=(n-l)*xa/(N-l);
xe2=n*xa/(N-l) ;
if n=\302\273N-lf

eps=l;
90 Overview of the Finite Element Method: One-Dimensional Examples \302\246
Chapter]

else
eps=epsr;
end
q=eps*k0'2;
Ael(l,l)=p/abs(xe2-xelLq*abs(xe2-xel)/3;
AelB,2>=Ael(l,l>;
Ael(l,2)=-p/abs(xe2-xel)+q*abs(xe2-xel)/6j
AelB,l)=Ael(l,2);
belB>=bel(l);
A(n,n}=A(n,n)+Ael(l,l)!

A(n,n+l)=A(n,n+l)+Ael(l,2>;
A(n+l,n)=A(n+l,n)+AelB,l)\320\263
A(n+l,n+l)=A(n+l,n+l)+AelB,2)
b(n,l)=b(n,l)+bel(l) ,-

end

% Enforce boundary conditions

A(l,l)=l?
for n=2:N,
A(n,l)=0;
A(l,n)=0;
end

A(N,N)=A(N,N)+alphaa*p;
b(N)=b(N)+betaa*p;

% Solve the resulting NxN system

x=inv(A)*b;

% Compute reflection coefficient

R_FEM=abs{ (it(N)-exp( jkOxa>)/exp ( jkOxa) ) ;

zeta0=377;
zetal=sqrt(l/epsr)*zeta0j

R_exact=abs((zetal*tanh(j*kl*t)-zetaO)/(zetal*tanh(j*kl*t)+zetaO));
Rf {z,l)=beta,-
Rf(Z,2)=R_FEM;
Re(z,l)=beta;
Re(z,2)=R_exact

end

%Plot Reflection Coefficient

elf;
plotfRf(:,1),Ret:,2>'r')(
hold;
plot(Re(:,l),Re(:,2));
References 91

xlabeK 'beta' );
ylabeK ' | ft i \342\200\242)
s

% End of program

APPENDIX 2: USEFULINTEGRATION FORMULAE FOR ONE-DIMENSIONAL

FEM ANALYSIS

- -
xn)q{x)dx * + q(xn+l)]
p'*\342\200\236+,
-v)(.v
^ [q{xn)

/2 = (x - \321\205\342\200\236_,J
<?(\302\246*>rfx
(j-J p
- .vJ
(xn+1 q(x) dx % -^
[iq{xn) + <jf(.vtt+1)]

=\302\273 + /K*,,-i')]
p{x)dx \\[\321\200(\321\205\342\200\236)

= -L \320\223\" -
.v,,_,)/(.v) *
(.v
^- [2/(-vn) +/(*\342\200\236_

=
\320\263
\320\223
.v,,

In all cases, hn
= -
xn\\
\\\321\205\342\200\236+\\and hn_\\
= \342\200\224 are
xn-i
\\\321\205\342\200\236
I assumed to be small. These
are derived by introducing a linear approximation for the functions g(.v), j\\x), and
p(x). For example,in deriving the approximation for /4, p{x) was approximated as

\"n-l \302\253n-1

REFERENCES

[1] R. [Link] methods for a solution of problems of equilibrium and

vibrations. Math. Soc, 49:1-23,1943.

Bull. Amer.
[2] J. H. Argyris. Energy theorems and structural analysis. Aircraft Engineering,
26:347-356, 1954.
[3] P. Silvester. Finite element solution of homogeneous [Link]
Frequenza, 38:313-317, 1969.
[4] P. Silvesterand G. Pelosi,editors.
Finite Elements for Wave Electromagnetics:
Methods and Techniques. IEEE Press, New York, 1994.
[5] P. [Link], A. K. Bahrani, and \320\236. Zienkiewicz.
\320\241 Application of finite elements
to the solution of Helmholtz's 1968.
equation, Proc. IEE. 115:1762-1766,
[6] 1. Stakgold. Boundary-Value Problems of Mathematical [Link],
New York, 1968. Volumes I and II.
92 Overview of the Finite Element Method: One-Dimensional
Examples \302\246
Chapter ]

[7] S. \320\241
Charpa and R. P. Canale. Numerical Methodsfor Engineers. McGraw-

Hill, New York, Second edition, 1988.

[8] R. L Burden and J. Faires. Numerical Analysis. PWS Pub. Co., Boston. Fifth
edition, 1993.

[9] T. B. A. Senior and J. L. Volakis. Approximate Boundary Condition,*; in

Electromagnetics. IEE Press, London, 1995.

[10] E. Kreyszig. Advanced Engineering Mathematics. John Wiley & [Link]
York, Fifth edition, 1983.
[11] D. M. Kingsland, J. Gong, Volakis, and
J. L. J.-F. [Link] of an

anisotropic artificial absorber for truncating finite element meshes. IEEE Tram

Antennas Propagat., 44:975-982. July 1996.

[12] S. R. Legault, T. B. A. Senior, and J. L. Volakis. Design of planar absorbmc
layers for domain truncation in FEM applications. Electromagnetics, 16D):451
464, July-August. 1996.
Two-Dimensional

Applications

41 INTRODUCTION

Having discussed one-dimensional examples, we now proceed with the application of

the FEM to two-dimensional problems. Although these represent simplifications of
the real-world three-dimensional situations, they are much simpler to formulate and
solve. Thus, they are appropriate for illustrating the FEM procedure and all
aspects
of the three-dimensional analysis can be discussed in the context of two dimensions.
Specifically, matrix assembly, absorbing boundary conditions,and various hybrid
formulations of the FEM can be discussed in two dimensions without loss of general-
This
generality. chapter is therefore important in understanding the machinery needed to
carry out FEM analysis and in illustrating its capabilities.
Two-dimensional models can used also be
to generate useful results for a
number of
practical problems in electromagnetics. For example, \320\242\320\225
and TM
mode analysis for straight waveguides of arbitrary cross section can be carried out

using a two-dimensional formulation [l]-[2] or the capacitance of transmission lines,

such as those in Fig. 3.1, can be found by solving Laplace's equation in two dimen-
dimensions. Similarly, the can be computed by carrying
inductance out a two-dimensional
analysis of Maxwell's equations. Also, an understanding of three-dimensional prob-
problems can be obtained by performing an analysis of closelyrelated two-dimensional

problems. Since two-dimensional analysis is much simpler, it is often a quick way to

obtain results for many practical open domain [Link] latter are among the
most difficult to solve and include those of radiation and radar scattering. Scattering
involves the computation of fields returned from a given structure due to plane wave
excitation or. more generally, the impinging radar wave. Two-dimensional analysis
has been extensively used in scattering and has provided engineers with an under-
understanding of scattering phenomena in the absence of three-dimensional simulation
which may be prohibitive due to their greater CPU and storage requirements.

93
94 Two-Dimensional Applications \302\246
Chapic

In this chapter we cover many aspects of the finite element method (and hybr

versions) at sufficient detail to provide the reader with a comfortable level of uui
standing its implementation. Such an understanding is essentialfor a three-dim,
sional analysis where many of the steps must be discussed symbolically due loi
size of the matrices even for very small problem sizes. We begin by first perform;
reduction of Maxwell'sequations to two-dimensional wave equations and \321\200\320
with the solution of the latter by following the same FEM steps discussed in i
previous chapter. That is, we generate the weak form of the wave equation and cur
out its discretization with the introduction of linear shape [Link] as\302\253

bly and boundary conditions are then discussed for determining the propagaii
constants in waveguides, and we give examples of this type of [Link] proceed
with the solution of open domain problems, we first discuss absorbing bound;
conditions and material absorbersfor truncating the finite element mesh. The \320\277\32

ary modifications of the FEM system are presentedand example calculations,!

given for scattering by cylindrical [Link] use of the boundary integral!
truncating the FEM mesh is presented in Section 4.4, and an example application
scattering by two-dimensional recessed cavities is given. Finally, we discuss I
implementation of two-dimensional solutions using edge-based basis funeiiu
since these are essential in three dimensions and are used extensively in subsequc

chapters.

4.2 TWO-DIMENSIONAL
WAVE EQUATIONS

Throughout this we shall assume that

chapter the fields or potentials are till
independent ofcoordinate
the z or have a known z dependence as is the case
propagation in waveguides. That is, it is assumed that the geometry's cross section
any xy cut remain invariant for all z values. The appropriateassumption depends
the type of physical problem being consideredas discussedbelow.

4.2.1 Transmission Lines

The characteristic impedance Zc and phase velocity vp

of the transmission!
can be computedby solving Laplace's equation

where V? denotes the Laplacian in two dimensions. The characteristic impedt

and capacitance of the transmission line is found by following these steps (seef
4.1):

1. Carry out the FEM solution and find V(x, y) in the absence of dielecir.

2. Determine the charge per unit length of one of the conductors by carry
out the integration

eh Well |j
Contour
Section 4.2 Two-Dimensional
\302\246 Wave Equations 95

Finite
element
Integration
contour for
evaluating
charge on
enclosed
conductor

Figure 4.1 Shielded microsirip line.

where n denotes the outward directed unit normal and the contour is shown
in Fig. 4.1.
3. Evaluate the capacitance per unit length of the free space filled transmission

line from

D.3)
\320\264\320\272

where potential difference between the two conductors.

AV is the
4. Repeat steps 1 to 3 for the original dielectrically filled transmission line to
obtain the per unit length capacitance of the line \320\241
=
-$p.
5. Determine the line's characteristic impedanceusing

}
z ... .. D.4)

where v0 denotes the wave velocity of the air-filled transmission line. We

also note that vp = vo/y/E^ where ee is referred to as the effective dielectric
constant of the substrate and, from D.4)

=
?\342\200\236 D.5)

4.2.2 TWo-Dimenslonal Scattering

In scattering analysis, the excitation is a plane wave of the form

D.6)
for incidence
\320\242\320\225 (also referred to as H. polarization) or

E' = fg D.7)
96 Two-Dimensional Applications \302\246
Chapter 4

for TM incidence (also referred to as E: polarization), as illustrated in Fig. 4.2. Boih

of these plane waves are impinging upon the scatterer at an angle They
\321\2040. may bt

alternatively written as

where k1 = \342\200\224
(j?cos0o +.vsin0o) 'S tne direction of the incident wave and

r = xx + yy is the position vector at which the field is observed. More generally,a

magnetic current source zMiz{x, y) for \320\242\320\225
incidence or an electric current sourct

zJb(x, for TM incidence can be assumed

y) as the excitation. For these types oi

excitation, the fields scattered by the cylinder will also be z directed and thus la
the \320\242\320\225
case the vector wave equation becomes

V, X \320\223-
V, X - zklflrHz
= -ZJ(D8OMI: D.81
B#s)j

where we included the possibility of a magnetic source as the excitation

(k0 = 2jt/A0 = u\\/jup(j. =

\320\272 2jr/A = Using
\321\210^/\320\251).
the identity1

V, x \\- V, x (?#jl = V, x \\- (-z x V,#r)l = -fVr

-
\302\246
V,#:
L?r J Ler J ?r

the vector wave equation D.8) is now reduced to the scalar wave equation

D.41

Scattered <ft
field \\

Scattering
cylinder

Contour, \320\241

Artiftcfai mesh truncation

boundary

Figure 4.2 Two-dimensional scattering configuration.

'Derived from the vector identity [3] V x (A x B) = AV -

\342\200\242
\320\222 BV \342\200\242 V)A - (A \342\200\242
A + <B \342\200\242 V)B.
Section 4.2 Two-Dimensionai
\302\246 Wave Equations 97

for inhomogeneous [Link],H. denotes the z component of the total magnetic

field. In scattering, it is typically decomposed as

//.- = \320\257!\320\2501+\320\257: D.10)

where Hi is the component of the magnetic field generated by zMt. in isolation or is

simply given by D.6). component The Hfal is referred to as the scattered magnetic
field (z component only) and is the unknown quantity of interest. As will be dis-

discussed later, there are advantages in solving for tff01 directly and in this case the
pertinent wave equation is
\302\246 \\ + = -V, \342\200\242
VfH{\\
_ ^fr =f(x,y)
+j(OSoMt:
V,
(~ V,#fat klnrff?al
(y

D.\320\237)

obtained by substituting D.10) into D.9). From duality, the corresponding wave
equation for TM incidenceis

Vr
\302\246 + ftfer?f\302\253
= -V,
\342\200\242 - kierEl+>\320\274\320\276^- D\320\2332)
f-L V,?f\") (~ V,?i)

where denote
??\320\270\" .*ie z component of the scattered electricfield and El. is the \320\263
component of the incident electric field generated by zjt in isolation or is simply
given by D.7).

4.2.3 Waveguide Propagation (HomogeneousCrossSection)

In waveguide propagation, the field is assumed to be of the form

|U,(a-,v)
\302\246jU/X-v.y)

where /? is the propagation constant along z and V(x,y) is the field value over the
waveguides cross [Link] waveguide whose nonmetaUic region ?2 is filled with

a homogeneous material certain standard simplifications are afforded. In these wave-

waveguides, the configuration as well as the field behavior is independent of the third
variable. This would be the z direction in a Cartesian coordinate system. Thus, in the
absence of a source, the identity

Vxfi Vxu) =- Vx VxU=--V2U D.14)

holds. In D.14), U denotes the electric or magnetic field vector.
The identity D.14) can be used to rewrite the vector wave equations as (see
Section 1.2)

V2H + fc2H = 0 D.15)

V2E + *2E = 0 D.16)
with =
\320\272 Given that V2H = xV2Hx + + ?V2#:, these equations per-
ku>/JI^. yV2Hy
permit a decoupling of the differential equation among the field components. Thus for
modes
\320\242\320\225 \342\200\224
[\320\225- we
\320\236.\320\235.\321\204
\320\236), only need to solve the scalar wave equation
98 Two-Dimensional Applications \302\246
Chapter 4

= \320\236 D.17)

where V; denotes the surface gradient and we set 82H:/dz2 = -P2H:. This is also
called theHelmholtz equation and is basically the scalar equivalent of the curl-curl

vector formulation. We note that once #. is found from D.17), Maxwell's equations
can be used to obtain the other field components using the expressions
Ey
= Ex
= -Utofi/ytX8HI/8y),
\320\270<\320\276\321\206//)(\320\264\320\230:/\320\264\321\205), Hx = (-jfi/fWHJSx). =
\320\251

(-JP/Y KdH;/ty) (see also Section 1.4).The quantities

=2-
D,18)
we
denote the wavenumber and characteristic mode impedance, respectively.
Alternatively, we could opt to solve the vector wave equation

V, V, = 0 D.19)
x(\342\200\224 xE,)-(^,-^)E,

in which

D.20)

is the total transverse electric field in the guide. Again, D.19) is valid only for cases
where the field variation is independentof the third dimension.
Alternatively, for the TM modes (\320\257.= \320\236, 0), the appropriate
E. \321\204 scalar and
vector wave equations are simply the duals of D.17H4.19).However, the boundary
conditions are of the Dirichlet type when solving for Ez and of the Neumann type
when solving for H-. As was shown in Chapter 3, the relation = 0 serves as
\320\265\320\250_,/\320\255\320\273

the boundary condition for the case.

\320\242\320\225

4.2.4 Waveguide Propagation (Inhomogeneous Cross Section)

Let us again consider the in Fig. 4.3 consisting of a center conductor

waveguide
enclosed shell. The region between the outer boundary
by an outer metallic of the
inner conductor and inner boundary of the outer conductor is the pertinent compu-
computational domain. When the waveguide's material cross sectionis inhomogeneous, the
fields can exhibit variation in the third dimension. Therefore, the normal component
of the field is required for an accurate representation of waveguide propagation,

Also, the resulting differential equation is not easily separable,which leads to a

coupled differential equation and normal components of the field.
in the tangential
To simplify the vector wave equation, we begin by rewriting the del operator
using D.13) as

V = 1 e = x
V, + ? V, -m
| + 1-jfii
\321\203 D11,

implying that

V x E = (V, -jfiz) x \342\204\226,+??\342\200\242;)

= V, x E, + V,?, x z -jfi x E,
= V, x E, + (V,?r+\320\243\320\224\320\225\302\273)
x f

in which E= E, + zE: as [Link],

Section 4.2 \302\246
Two-Dimensional Wave Equations 99

Finite element
y* mesh

Metallic
<a) (b)
boundary

figure 4.3 Waveguide configuration: (a) cross section of waveguide; (b) three-
dimensional view.

V x i- V x E = (V, -jfiz) x \320\250

[V, x E, + (V,?: +J0E,) x f)

- V, x
\342\200\224 x
V, E,
- z x \342\200\224
[(V,?. \321\205
f]
-\320\254\321\203/\320\227\320\225,)
Me Mr

+ Vf x (V'?:+/)8E') x D.22)
\\j f]

Next we introduce D.22) into the vector waveequation V x (l/fir)V x E \342\200\224

?j)?rE
= 0
and set each of the vector components to zero. This permits the decomposition of the
original wave equation into a pair of differential equations [4, 5, 6]

V, x \342\200\224x
V, E,
- ^ (V,?. -
*ge,E,
+\320\243/\320\227\320\225,)
= 0 D.23)
Mr Mr

V, x \\\342\200\224
(V,E. +J0E,) x - Ager?:J= 0 D.24)
fj

for the transverse and z components of the wave equation. Clearly, D.23) and D.24)
a
represent pair coupled of differential equations which either needs to be decoupled
or solved as is. Decoupling them using the divergence property yields a nonsym-
metric generalized eigenvalue problem, the solution of which is numerically ineffi-
inefficient. However,
using simple variable transformations [2],

D.25)
100 Two-Dimensional Applications \302\246
Chapter 4

the coupled pair of differential equations D.23) and D.24) can be expressedas

V, x V, x +
<\320\233 p2
\342\200\224
(V,e, + e,) = klere,
(\342\200\224

r. 1= D.26)
/?4 x (V,e: + e,) x fJk2oere.J
\\\342\200\224 fJ
The coupled pair of differential equations D.26) can now be solved for f}2\342\200\224the

square of the propagation constant of the inhomogeneous waveguide- subject to

the following boundary conditions:

on PEC surfaces and

<V* +''>\342\200\242*
=

; D.28)
V, x e, = 0

on PMC [Link] differential equation to be discretized is obtained by adding

the two coupled differential equations and weighting them with the necessary weight-
functions.
weighting This fonnulation leads to a symmetric generalized eigenvalue problem
which is a direct result of the transformation D.25).

4.3 DISCRETIZATIONOF THETWO-DIMENSIONAL WAVE EQUATION

From the above presentation, a general form of the two-dimensional wave equa-
equation is

\342\200\242
V \320\253-v,y)VU(x, y)] + klq{x, y)U{x, =\320\224\321\205,
\321\203) \321\203) D.29)

and it is understood that here V = V, = x ^ + \321\203 This can be specialized to the

-jfc.
problems of scattering and waveguide propagation by choosing p{x,y) and ^(\320\273\\.|')

appropriately. When we consider \320\242\320\225

(or Hz polarization) we must choose

U{x,\321\203)
=
\321\203). p(x,
\320\235\320\263{\321\205, >') =
\342\200\224
. <l(x< = Pr
\320\243)

while for TM or Ez polarization

The steps to be followed for the solution of D.29) via the FEM parallel those
given in Sections 3.4 to 3.7 for the solution of the corresponding one-dimensional
(ordinary) differential equation. They involve
\302\246
casting of the original wave equation to its weak form to obtain a single
functional incorporating the conditions imposed by the wave equation and
the boundary conditions.
Section 4.3 Discretization
\302\246 of the Two-Dimensional Wave Equation 101

\302\246
tessellation of the computational domain allowing for a discretization of the
weak form to a linear system of equations element by element.
\302\246
assembly of the element equation and imposition of the boundary conditions
to obtain the final linear system of equations.

In this section we carry out the first two steps, and in the subsequent section we
considerthe assembly of elemental equations to solve for the fields and eigenvalues
associated with a metallic waveguide (closed domain problem).

4.3.1Weak Form of the Wave Equation

The pertinent residual of D.29) is

R(x,y) = V \302\246
/>(r)Vt/(r) + k20 q{r)U{r) -/(r) D.30)

where as usual r = xx + yy the position

denotes vector. To derive the weak form, we
multiply D.30) by a weighting or testing function W(x,y) and enforce \320\233(\320\263)
= \320\236
over
the domain of each element rather than point by point. This gives

mm of
W(r)R(r)dxdy= Q D.31)

which is the generalization of C.9) for the one-dimensional case. As discussed in

Section 3.4, the weighting function must again be compatible with the boundary
conditions and be square integrable over the domain. Its derivative must also be
square integrable.
To reduce the order of the derivatives in the residual and r*tcoducethe bound-

boundaryterms, we must make use of the identities

WV-pVU- 4{pWVU)-pVW 4U D.32)

f f V \302\246
(p WVU) ds = I pW(VU \342\200\242
n) dl D.33)

in which Q denotes the pertinent computational domain, is the contour

\320\241 enclosing
Q as illustrated in Fig. 4.4, and n refers to the outward directed unit normal vector to
the contour \320\241We note that the latter is the divergence theorem, which can be
considered as the generalization of integration by parts to two dimensions.

Figure 4.4 Parameters for the application of

ihe
divergence theorem.
102 Two-Dimensional Applications \302\246
Chapter 4

Making use of D.32) and D.33) into the weighted residual equation D.31)
yields

-Mr)V^(r) \342\200\242
V?/(r) + kfab) W(r)U{t) - ds
\320\251\321\202)\320\224\320\263)]

+1 p(r) \302\246
\320\251\321\204V?/(r)] dl = 0 D.34)
ic
This is the weak form of the two-dimensional scalar wave equation and should be
compared to the one-dimensional weak form in C.11). Again, we note the presence
of the integral over \320\241
boundary which boundary allows for the imposition of the
conditions. D.34) provides a single statement
Thus, incorporating the conditions
implied by the wave equation and the pertinent boundary [Link] noted
for the one-dimensional case, this is at the heart of the finite element method.

4.3.2 Discretization of the Weak Wave Equation

D.34) we proceedby introducing

To discretize a discrete representation of the
field U(x, making appropriate choicesfor the weighting
y) and function W(x, y). As
a first consider a tessellation of the computational
step, we domain ?2 in small
triangular elements (see Fig. 4.1 to Fig. 4.3). We next choose to approximate the
fields in each triangle as a linear function and in so doing, we can choose the
expansion or basis functions Nf(x.y) given in Section 2.3.2 to represent U[x,y)
over Q. Specifically, we expand ll(x,y) as

U(x,y) = YtYtUfNXx<y) D.35)

i=l
\302\253\342\200\242=1

where Uf are the unknown coefficients of the expansion and represent the field or
potential values at the nodes of each triangle. This representation is therefore
referred to as a node-based expansion. As usual, Ne denotes the number of elements
used for tessellating the domain. The procedure of tessellation is referred to as

meshing and typically each side of the triangle is chosen to be less than 1/10 of a
wavelength.
From Chapter 2, the explicit form of the shape functions Nf (x, y) is

0 \\!fl.
otherwise D.36)
v

where

Ar = ldet 1 \321\203\\ =^[(^-4)(>'\320\267-^1)-D-4)(>'5->'|)]

\321\205\320\247 D.37)
.1 A y\\_

is the triangle area. The coefficientsin D.36) are given by

with the indices (ij, k) following the cyclical rule. That is, (ij, k) = A,2,3),B,3.1),
orC,1,2) for the first, second, and third local nodes, The shape func-
respectively.
functions Nf(x,y) are pictorially illustrated in Fig. 4.5. They are equal to unity at the Ah
Section 4.3 Discretization
\302\246 of the Two-Dimensional Wave Equation 103

Linear field
approximation
over the eth
triangle

2 1
x, y)

Figure 4.5 Node coordinates for the eth triangle and illustration of the node-based
expansion functions Nf{x,y).

node of the eth element and taper

the other two nodes. linearly to zero at
Consequently, expansion
the eth element
J^=l UfN*(x,y) scales the shape functions
depending on the value of the field or potential at the nodes of the eth element.
Substituting the expansion D.35) for U(x,y) into the weak wave equation
D.34)yields
\320\232\320\227

E E
w

f f f D.39)
Jc J Jn-
Notice that we have not used the expansion D.35)to approximate VU
\320\273 =
\342\200\242 on the
\320\251
boundary since
\320\241 the behavior of \320\251
on \320\241
must be provided through the boundary
conditions.
A linear set of equations can now be obtained by employing Galerkin's method
where we choose the weighting function lV(x.y) equal to the expansion basis
Nf(x,y). i = 1.2,3- Doing so yields

Y,Vf\\\\ \302\246
VAff(r)
[-\321\200(\320\263)\320\2511\320\263) + *grfr)A(f(ryvftr)] dx dy

+
f p(r)NJ{t)n
\302\246
Vt/(r)rf/ = I
f /=
\320\233\320\223/(\320\263)/(\320\263)dxdy, 1,2,3 D.40)
where we have temporarily dropped the sum over the elements since Nj(x,y) is
nonzero only over the eth element. We will later perform the summation over all
104 Two-Dimensional Applications \302\246
Chapter 4

elements during the assembly process of the matrix system. Also, the presence of Ihe
boundary integral is required only if element has an edge bordering
the eth the
contour The
\320\241 contour segment C, refers to the edge of the eth triangle which is
part of \320\241

From D.40), we can obtain a 3 x3 system of equations by running through all

choices =
of \320\243 1,2,3. If we assume that on the
\320\241 field satisfies the Neumann bound-

boundarycondition h \302\246
Vt/ = = 0,
\320\251
we have

Ah Ah Ab~ U\\
Ah Ah Ah U{ =. bl D.41)
Aji Ah Ah. v\\.
or

which is the element matrix system. The explicit form of the matrix entries is

4 = -/ \302\246 dx dy
V/Vf (r)
jJ V/Vj(r)

Nt(r)Nj(r)dxdy, /=1.2.3, 7=1.2,3 D.42)

/,;=[[ Nf(r)/(r) dx dy, j = 1.2, 3 D.43)

J Jar
wherewe assumedp(r) \321\217= and q1' are constant
q(r) \302\253s
\321\200\320\265 over each element to permit a
closed form evaluation of the integrals. SinceN\"(r) are linear functions, the evalua-
evaluationof the [A*] matrix entries can be done in closed form. We have

where imply
((\320\270) column vectors)

=
[\320\232'] L\" =
f J NfNfdxdy] qe^JNf\\{Nj)Tdxdy
Evaluating the entries of the submatrices [K?] and [Ke] yields [note that the constants

b* below are those in D.38) and are not related to the excitation column in D.41)]

D.44)

2 1 \320\223

1 2 1 D.451
1 1 2
Section 4.3 Discretization
\302\246 of the Two-Dimensional Wave Equation 105

The latter is independent of the triangle coordinates (except for the area multiplier)
and is referred to as a universal element matrix [7]. By making use of D.44) and
D.45) in D.42), the [Ae] matrix entries can be more compactly written as

= :? + D.46)
^ (l s,j)
=
a% /Gt + *\302\247kj
+
\321\204\321\211
^\321\204
+ kfc

where

8=l
j
'' 0 otherwise
\\

Thus, linear triangular

for the chosen tessellation, the matrix entries are given in

closed be easily evaluated. Dependingon the form of the function

form and can
f(x, y), the excitation column entries
bj may need
to be computed numerically.

4.3.3 Assembly of Element Equations

As was the case for the one-dimensional solutions, the next step in the finite
element procedure is the assembly of the element equations D.41). This refers to the
procedure of carrying out the sum

? {/,') D.47)

as implied by the original discrete

form of the \"weak\" wave equation D.39). Since
several of the elements the same node, the sum or assembly D.47) con-
may share
consolidates the surrounding elements to yield a single [Link] is simply the sum
of the element equations from all of the surrounding elements. Specifically, let us
consider node 1 which is shared by five elements, as shown in Fig. 4.6.
For this case, the sum D.47) producesthe equation

+ Atf'US+l) = D.48)
W2+l }_;H!l+>
1=0 i=0

(e'+2)th

e'thelement

Figure 4.6 Illustration of a node shared by

five elements (triangles).
106 Two-DimensionalApplications \302\246
Chapter 4

for node I with the matrix entries as given in D.46). The correspondingequations for
the other nodes are very similar, except for differences in the
superscripts/subscripts
and the order of the sum. We should remark that the assembled equation D.48) for
node the same regardless of the number
1 is of elements/nodes contained in the
computational domain. That is, even if the entire domain contains thousands of
nodes, D.48) will still involve only six nodal fields, implying that the corresponding
row of the assembled matrix system [/4](t/} = lf>] W'H contain only six nonzero
entries even though the rank of [A] is in the thousands. Thus, the assembledfinite
element matrix is always very sparse and this is a major advantage and characteristic
of the FEM. bandwidthThe and structure of [A] is determined by the connectivity of
the a result of the tessellation [Link] can be understood,
nodes, the bandwidth
of the matrix [A] is strongly dependent on the node numbering scheme. We can
reduce the bandwidth by numbering the adjacent nodes using consecutive numbers.

This is difficult to achieve, but sparse matrix storage schemes as those presentedin
Chapter 8 can be used to maintain their efficiency for matrices of the same sparsitv
but different bandwidths. Parallel computing architectures take advantage of spar-
sity but in the case of vector processors, narrow matrix bandwidths must also be
maintained for substantial efficiencyimprovements [8].
A key issue in performing D.47) is the transforma-
the assembly as dictated by
from
transformation local to global nodes. This was discussedfor the one-dimensional analysis.
However, because of the easily predictableconnectivity of the elements (i.e., each
element was sequentially numbered and each node was shared by the two adjacent
segments) the local to global transformation was not an issue for the one-dimen-
case.
one-dimensional The issue of node numbering becomes apparent when we look at the
assembled equation D.48).
Since the unknowns Uf must be eventually put into a single column (with one

subscript), it is necessary to have a readily available mapping between the local and

global nodes which are associated with the <?th element* Thus, in addition to the node

geometry data provided to the finite element program, we must also provide infor-

information about the local and global node numbering schemes. Four tables may be
required before carrying out the matrix assembly routine:
\302\246
Node Location Table
A listing of all mesh nodes (interior and boundary nodes) using global
numbers and their corresponding (x,y) coordinates.
This table specifies
the geometry of the input configuration.
\302\246
Triangle Connectivity Table
The global nodes comprisingeach triangle are given by this table. For
example, by referring to Fig. 4.7 we observe that element #3 (e - 3) is
formed by nodes 3, 5, and 2, as given in line 3 of the table. Basically, the

table defines three arrays: n(l, e). nB,e), and \321\217C, The
\320\265). first of these pro-
provides the correspondence between local node 1 of the eth element and the

global nodes. The other two arrays provide the same correspondence infor-
information for the other two nodes of the triangle. Let'ssay for example that we
are working with the local nodes of element e = 4 and want to get the

corresponding global nodes of that element. Then the value of \302\253A.4) will
give the global number of local node 1, \302\253B.4) will provide the global
number of local node 2. and so on.
Section 4.3 \302\246
Discretization of the Two-Dimensional Wave Equation 107

Global node
numbers

Node location table Triangle connectivity table

Global coordinates
{x, \321\203) Element Local node arrays
rtode# x \321\203 e Mte) nB,e) nC,e)
1 *\\ \320\243\\ 1 2 4 1
2 2 5 4 2
3 *\320\267 \320\243\320\267 3 3 5 2
4 5 6 4
\342\200\242\302\246 \342\200\242

\342\200\242

Boundary element
connectivity table
Surface Local node arrays
edge
number
5 rtsd.S) /JsB. S)
1 6 4
2 4 1
3 1 2

Outer surface
of mesh

Figure 4.7 Geometry and connectivity data tables required for matrix assembly.

Boundary Element Table

To impose the boundary conditions it is necessary to identify the surface

edges (line elements) on the outer boundary of the mesh and their associated

nodes. This table can be generated using the data in the previous two tables.
The data manipulations required to generate the surface node and element
information is typically part of the data preprocessorand is an important
step before assembling the final system. An example of such a boundary
element table is given in Fig. 4.7. For the two-dimensional case, it suffices to

identify all segments which bound the mesh and the nodes which form those
108 Two-Dimensional Applications \302\246
Chapter 4

elements. Again, the listed arrays provide the correspondence between the
local surface element nodes and the global nodes.
\302\246
Material Group Table
This is a look-up table the material
for specification of each element. 1\320\273

practice, the same material covers blocks or sectionsof the domain and il is
therefore not necessary to specify the material parametersfor each indivi-
individual element. Instead, one may choose to attach a material code column to
the element connectivity table. That is, the material of each element is
specified through the \"code\" which is in turn associated with specific values
of er and fir.
The first two of the above tables are always required but the latter two may or

may not be needed depending on the application at hand. Also, it may be convenient
to introduce other tableswhen different elements are used or more complexgeom-
geometries are modeled.

4.3.4 Assembly Example: Waveguide Eigenvalues

To show how the assembly is performed, let us do a specific example. For

illustrative purposes, consider the 18-element rectangular waveguide cross section
shown in Fig. 4.8. The node location and element connectivities are given in the

accompanying tables. If we assume that U = H,, then ^ = n \302\246

VI/ = 0 at the bound-

boundarynodes of the rectangular metallic waveguide. Thus, the boundary integral in

D.40) vanishes and the FEM system for the node fields is obtained by performing
the assembly
18

D.49)

where the dimension of (?/) is 16 and {bL'} was set to zero since no excitation is
assumed. A procedure for carrying out the assembly is as follows:
Note the correspondence between the local and global nodes. For example,
C/f=l() t/6, where, as usual, the
= single subscript refers to the global node 6. Thus,
the element matrix for e \342\200\224
10 is

\320\27310
.19
Aw \320\220\\\320\263

\320\220\320\236
Uu = 0 D.50)
,10
\320\260\\1 \320\233
32

If we are interested in obtaining the assembled equation for global node6, we

must also consider all of the elements sharing this node. In addition to element 10.
node 6 is shared by elements 2, 4, 9, and 7. The correspondingelement equations for
each of these are

'\320\220], An An' Hi

Ah Ah Ah
= 0 D.51)
.Ah Ah Ah.
Section 4.3 Discretization
\302\246 of the Two-Dimensional Wave Equation 109

Local
node
numbers (\\

Global node numbers Element number

Triangle connectivity table Node location table

Element Local node arrays Global

number node
e \320\2261.e) nB,e) nC, e) number X \320\243

1 1 2 5 1 0 0
2 2 6 5 2 0.5 0
3 2 3 7 3 1.0 0
4 2 7 6 4 1.5 0
5 3 4 8 5 0 0.25
6 3 8 7 6 0.5 0.25
7 5 6 10 7 1.0 0.25
8 5 10 9 8 1.5 0.25

9 6 7 11 9 0 0.5
10 6 11 10 10 0.5 0.5
11 7 8 11 11 1.0 0.5
12 8 12 11 12 1.5 0.5
13 9 10 14 13 0 0.75
14 9 14 13 14 0.5 0.75
15 10 11 15 15 1.0 0.75
16 10 15 14 16 1.5 0.75
17 11 12 15
18 12 16 15

Figure 4.8 Geometry, node location, connectivity, and boundary node data tables
required for matrix assembly.

Ax \320\220\320\263 V2'
\320\220\320\267'

Ax A>
\320\220\320\263 u7 = 0 D.52)
Ax A*.
\320\220\320\263
110 Two-Dimensional Applications \302\246
Chapter 4

\320\2707
= 0 D.53)
A\\\\

U \320\270
.Ax

\320\220\\\\ -4,2 \320\220\321\206

\320\220]] \320\220722^23
= 0 D.54)

Adding the element equationsin D.50) to D.54) which refer to testing at a common
node, yields the system

A\\2 0 0
Ah An + An a]3 0
a\\+a\\{
A . .9
An + Au + a\\\\

.9
0 + '
\320\220\\\321\212 \320\220

0 All 4? + /ill 0

0 0 ^21 + ' \320\220\\\320\263 + -422 .

D.551

U \320\270

We can continue assembling element equations onto this system until all elements
have been accounted for. However, the third equation in D.55), referring to tesling
at node 6, will not change through the remaining assembly [Link], this

equation is no different than the generic equation D.48). Thus, we have established a

pattern for assembling the final system of equations. Specifically, note that in the
final assembled system, the entries will be given by the sum

all D.56)
The index (see Fig. 4.8) e^ of the sum must be kept identical in n(i,e{) and
nU, ?i) when carrying out the sum over et, i and^. For example, A(,b
= a\\\302\260\\
+ \320\233^+

\342\200\224 +
a\\2
\320\220$\320\263 -^3|. and so on. A computer (matlab) routine for carrying out the

assembly is given in the appendix where the coefficients are adjustedfor the chosen
polarization. For a wavenumberof kn = 2\320\273
(i.e., A. = 1), the numerical values of the
assembled [A] matrix are given on the next page and can be used in conjunction with
a given excitation vector [b] to obtain the waveguide fields across its cross section.
0000000000*000*

1
.- \342\200\224
?, v> <
(\320\233

\320\2634
\320\276 <
\320\276
\320\276\320\276

- \342\200\224\342\200\242
\321\216 <\320\233
\\\320\236

\320\235
\320\276 \342\200\224'
\320\276 \320\276

OOOOOOOOOOOOOOsDOO
1\320\236 00 (\320\236
\320\236 1^. 1\320\233

ON \302\253N

\320\236
\320\236
\320\236
\320\236 f-l \320\237
\320\241 \302\253N
\320\236\302\273\320\236
\320\223\320\247
CI \320\236
\320\236
\320\223\320\247
\320\236

*t 't ^ \302\273\320\233
\320\236;O ^
W \320\276
\320\276\320\274\320\265
\320\276

\320\223\320\247 N \320\236
\320\236 f*J \320\236

\342\200\224
\320\276 \342\200\224
\320\276

lN-OOfN\302\273COOOO
\321\201\320\247\321\207\320\236

\342\200\224 \320\276
\320\276 \320\276

\342\200\224
\342\200\224
J5

\320\263\321\207 r-i \320\276

\320\276 \320\263\321\207
\320\276

\342\200\224
^\320\236\320\223\320\247\320\236\320\236
\320\234\320\236\320\236^\320\236\320\223\320\247\320\236\320\236\320\236\320\236\320\236\320\236
\342\200\224 \302\253\320\273
\320\276 \342\200\224
\342\200\224 \320\276
\342\200\224
\320\236 \320\236
\320\236 \342\200\224*
\320\276

\320\236\320\236^\320\236\320\241\320\236\320\236\320\241\320\236\320\236-\320\241\320\236\320\236\320\241\320\236\320\236\
m
\302\273ri \320\276

\320\265 \320\263\320\273
\302\273\320\276
%\320\276\320\276\320\276
\320\277\320\277
\342\200\224
^ \302\253\320\233
1\320\233 \342\200\224
vi
'\320\273
\320\276
\342\200\224
\320\276\320\276 ri \320\276

\321\201 \320\2766w
\320\276 \320\276

\320\236\320\241
** \320\223-1
\320\236\320\241 \342\200\224
I

111
112 Two-Dimensional Applications \302\246
Chapter 4

[Link] (Neumann ModesBoundary

\320\242\320\225 Conditions). One way to verify the
correct and implementation
evaluation of the FEM matrix is to use it in
computing
the known capacitance of a given transmission line such as the shielded microsLrip
line in Fig. 4.1. In this case, the wavenumber k0 is set to zero and thus only the [\320\251
matrix is used to evaluate the capacitance Co, as in Fig. 4.1. A situation where both
submatrices [Kv] [K] [obtained from the assembly of [Ay] and [K'] given in
and
D.46)] are used is
of computing
that the waveguide cutoff wavenumbers. Upon
comparing the wave equations D.17) and D.29), we observe that the cutoff wave-
numbers are obtained by solving the generalized eigenvalue problem

-[KV][H:] = /\342\204\226\342\204\226 D.57)

where y1 are the eigenvalues and \321\203

refer to the cutoff wavenumbers. Consideringa
rectangular (see Fig. 4.8) with
waveguide dimensions a/b = 2, the computed eigen-
eigenvalues are given in Table 4.1 and the corresponding mode field distributions (eigen-
(eigenvectors) are shown in Fig. 4.9 for TE|0 and TE| |. These calculations were carriedout
by Reddy et al. [9] using 400 triangular elements over the waveguide cross [Link]
seen, they are within 0.4% of the exact [10] eigenvalues given by

for the TEmn (m \321\204\320\236\320\276\321\202\320\277\321\204\320\236)

and the TMm/, (m 0 and
\320\244 n 0)
\320\244 modes. We remark
that the accuracy of the calculated eigenvalues deteriorates for the higher order
modes since the latter require a finer tessellation due to their more complexmode
structure.
We can solve the eigenvalue problem D.57)with equal ease for a coaxial cable
of inner radius r\\ and outer radius r2, as shown in Fig. 4.10. The typical triangular

TABLE 4.1 Cutoff Wavenumbcrs for a Rectangular Waveguide

ya {a/b = 2)
\320\242\320\225 TM Analytical 110] FEM Calculation

10 3.142 3.144
20 6.285 6.308
01 6.285 6.308
11 II 7.027 7.027
12 12 12.958 13.201
21 21 8.889 8.993

. t t \342\200\242
t t t t t t t t t ,
. t t I t t t t
t t t t t t t t t 1 T t t
t t t i t t t t t t t 1 1
..'it t t I , ,
. t t i t t II t t t t t
t t t i t t t t t t I t t
t t t i t t t t t t t t t
t i i i I
j
1 t 1 t t t
i t \342\200\242
i t t t t ! t t
I

Figure 4.9 Calculated TE|0 and \320\242\320\225\321\206

mode electric fields in a rectangular wave-

waveguidewith a/b
= 2. [Courtesy of Reddy el al. [9J.\\
5\302\253lion
4.3 \302\246
Discretization of the Two-Dimensional Wave Equation 113

mesh for this type of cross section is illustrated in Fig. 4.3. Calculations were carried
out by Reddy et al. [9] using 340 elements to model the cross section between the
inner and outer conductor for \320\263\320\2631\320\263\321\205
= 4. The first few eigenvalues (cutoff wave
numbers) and eigenvectors are given in Table 4.2 and Fig. 4.11, respectively. For

the shown mode

\320\242\320\225 values, the agreement between the FEM calculated and ana-

analytical [11] wavenumbers is within 0.8%.

Figure 4.10 Coaxial waveguide geometry.

TABLE4.2 Cutoff Wavenumbers for the Three Lowest Order \320\242\320\225

Modes
and Two Lowest Order TM Modes of the Coaxial Cable(the \320\242\320\225\320\234
mode
is excluded)

= 4)

Mode Analytical [II] FEM Calculation

\320\242\320\225,, 0.411 0.412

TE2, 0.752 0.754

\320\242\320\225,, 1.048 1.055

TMo, 1.024 1.030
TM,i 1.112 1.120

Figure 4.11 Calculated and TE2,

\320\242\320\225,, mode electric fields in a coaxial line with

r,/)-) = 4. \\Vourivsy of Reiklyet til- 1^1)

114 Two-Dimensional Applications \302\246
Chapter -I

[Link] TM Modes (Dirkhlet Boundary Conditions). For the TM modes,

U = E. and ?. = 0 on the boundary. However, = h
BU/\320\264\320\277 V?. =
\342\200\242
[Link] is non-
nonzero on the boundary and it is therefore necessary to treat =
\320\244 dEJdn as an addi-
additional unknown on the boundary. The contribution of the boundary integral then

plays the same role as the [endpoints]for one-dimensional problems (see Chapter 31
As discussed in Section 3.7, since the field values at the boundary nodes are zero, ii is
not necessary to test at these nodes. Instead,we set the boundary node fields to zero
whenever they appear in the system. By avoiding testing (or weighting) at the bound-

boundarynodes (or elements), the integrals over Cs do not enter in the construction of the
final system and can thus be neglected altogether. We also remark that the choice of
omitting testing at the boundary nodes is equivalent to setting the weighting func-
functions to zero when testing at these nodes.
Although the above arguments are sufficient to proceed with the FEM maim
assembly while neglecting the presence of the boundary integral, they are neverthe-
difficult
nevertheless to visualize without going through some of the details. Therefore, below
we will (for a moment) proceed with the assumption that the boundary integral
contribution is needed. We begin with the discretization of the boundary integral

by introducing the expansion

where Ns denotes the number of boundary edges A2, for example, in Fig. 4.8) and f,
carries the usual notation. That is, \320\244? refers to the value of \320\255\320\225:/\320\264\320\277
at the /th local

node of the sth segment of the boundary (see Fig. 4.7). The expansion bases arc
linear interpolation functions between the node values of They
\320\244(\320\263). are of the same
form as those discussed in Chapter 2 (see Section 2.3.1) and Chapter 3. Fora
constant value of .v or y, the two-dimensional shape functions reduce to the same
linear one-dimensional expansion functions given in B.3). Thus, we can write

x\\ - x , .
,boundaries
-
\342\200\224
on v = constant
4 *
= v? - v
IJ(r)
,,s
on.v =
,,.)\302\246
constant boundaries
\320\263> 1
\320\243

\320\236 outside the .vth segment

L\\(t) = 1 - Z.f(r). e C,
\320\263 D.59|

These expansion bases are applicableto piecewise rectangular boundaries, as in Fiji

4.8. In the case of circular boundaries,as in Fig. 4.10, an appropriate set of linear
basis functions is

otherwise

where is
\321\204 the angular variable ranging from = 0
\321\204 to =
\321\204 2\320\273\\
Section 4.3 Discretization
\302\246 of the Two-Dimensional Wave Equation US

To obtain the element equations for the TM excitation, we substitute D.58)

along with D.35) into D.39) and proceed as done for the case.
\320\242\320\225 The resulting
element system is

[Ae]\\Ue] )
= {be] D.61)

where the entries of [Ae] and (//} are given by D.42) and D.43). Those of [B\"] are

computed from

Nf(t)L](i)dl D.62)

L](r)L)(r)dl D.63)

(dt= dx or
dy for rectangular boundaries) and are associatedwith the last left hand

side of D.39).The latter expression D.63) for By results from the identification that
the ith surface or boundary segment must be an edge of the eth element for a
nonzero value of By. It is also understood that the boundary matrix [Bx] will be
nonzero only when the eth element is associated with a pair of nodes on the outer
surface/boundary of the computational domain.
of the elemental equation D.61)again
The assembly amounts to summing the
equations from each element weighting at the same global node. For the specific

rectangular waveguide example shown in Fig. 4.8. the resulting assembled system
will be of the form

D.64)

The expandedversion of this matrix system is

\342\200\242''is.:
116 Two-Dimensional Applications \302\246
Chapier-i

\320\276 0
0
Bn
\320\262\302\273 0
0 0
0 \320\236 \320\276 0
0 \320\276 \320\276 0
0 \320\276 \320\276 0
0 \320\276 \320\276 0
\320\262\320\274 0 8.12 0
0 \320\276 \320\276 \320\276
\320\276 0
0 \320\276 0
0 \320\276 0
0 \320\276 0 BaM
0 \320\276 ^1.1,14 \320\236

0 \320\276 ^14.14 0

0 \320\276 ?|5.!4 #15.15

0 \320\276 0 \320\236\320\236
\320\236 0

\321\204,

*2 /\302\2732

\320\244,

\320\2444 ^4

\320\2445 \321\2145

\320\244? bi

*8 \321\214%

\320\244, h
b\\o

bu
\321\204\342\200\236

\320\244|2 bn

\320\244|4 bu

\320\261''

in which \320\244\342\200\236
denotes the outward directed normal derivative of E. at the nth node
The values of \320\244 at the interior nodes are irrelevant becausethey are associated with

all zero rows, included for the proper addition of the [-4] and [B] matrices. It is
understood that [A] is a very sparse matrix, as discussed earlier in the chapter.
We observe that the above system involves 12 nontrivial unknowns from the
column
{\320\244} plus 16 unknowns from the {(J) column for a total of 28 unknowns

Clearly, this number of unknowns is much greater than the available 16 equations
For a solution for {[/} and {\320\244} we must add 12 more equations or conditions on the

values of {{/} and {\320\244).This is done through the introduction of the Dirichlet bound-
boundaryconditions satisfied by \\U) = {?.) on the boundary nodes. The procedure i;

referred to as condensation of boundary conditions and was discussed in Chapter 3

To see how it can be done in a consistent manner, let us divide the [U] column into
two parts, one for the interior nodes and another for the boundary or surface node
That is, we write (U) as
Section 4.3 Discretization
\302\246 of the Two-Dimensional Wave Equation 117

\\{US)\\

where in the case of Fig. 4.8, (?/'}= \\Ub, G7, U\\0, Un]T are the interior node fields
and [Us] contains the boundary node [Link], we formally define the column
(*} as {\320\244}
= *4.
*\320\267>
(\320\244,,\320\2442, *5. *8. *e. *i2' *i3- *u. *is. *ie)r which excludes all
interior node values of \320\244.
Using this notation, we can rewrite the system D.64) as

\"] Oil
Oil \320\236
[A'S]]\\{U')\\.\\O \\O \\_\\{b')
\\_\\{b)\\ 4

in which [-4\"] refers to the submatrix of [A] containing the interactions among the

interior nodes.
[AIS] and [ASI] are associated with interactions among exterior and
interior nodes, and [Ass] refers to the interactions among the boundary or outer
surface nodes. Similarly to{U1} and \\US), the excitation subvectors {b1} [bs] are
a nd

associated with the interior and boundary nodes,respectively.

When we invoke the Dirichlet boundary condition

the system D.65) can be decomposed into two independent subsystems

= {b')
[A\"]{U') D.66)
and

s'l
{bs} D.67)
The interior node fields are now
decoupled from (\320\244}
and we can therefore proceed to
solve D.66) without consider the solution of (\320\244).Thus, the boundary
a need to
integral can be neglectedwhen assembling the FEM system provided U or its normal
derivative are zero on the boundary C. After [U1] is found, we can return to D.67)

and determine {\320\244),if necessary, by solving the system

= \\b) D.68)

where

{b} = [bs}-[As'][U')
In the case of an eigenvalue problem, the excitation column {/>} is set to zero
and D.66) is written as

~[K\"]\\E;) = ftK\"\\[Es\\ D.69)

where [K\"] and [K11] are identical to those in D.57) except that only the entries
associated with the interior nodes are kept and that pejq\" are defined for the TM
polarization case. Somevalues of for
\321\203 the TM modes, obtained from D.69), are
given in Tables 4.1 and 4.2 for the rectangular and coaxial waveguides. Also, Fig.
4.12displays the fields of the lowest order TM mode for each of these empty wave-
waveguides. The analysis can be carried out using the matlab program in the appendix.
For the mesh in Fig. 4.8, the corresponding [An]. [Ky] and [Ku] are
118 Two-Dimensional Applications \302\246
Chapter -I

fttttttftf...

1 1 t 1 I t
'
\342\200\242
1 \302\246
i ! | |

TM,, mode TM., mode

Figure 4.12 Calculated fields for the lowestTM modesof the rectangular (a/b = 2)
and coaxial (ri/rt = 4) waveguides. [Courtesy of Rcddy cl at. (9j.\\

-2.9438 0.9112 2.4112 0.4112

=
0.9112-2.9438 0 2.4112
[A\"]
2.4112 0 -2.5326 0.9112

0.4112 2.4112 0.9112 -2.5326

0.521 0.0104 0.01040.0104
0.0104 0.0521 0 0.0104
[K\"\\ =
0.0104 0 0.0625 0.0104
.0.0104 0.0104 0.01040.0625
J
5.0000-0.5000-2.0000 0

-0.5000 5.0000 0 -2.0000

-2.0000 0 5.0000 -0.5000
0 -2.0000 -0.5000 5.0000

The above TM mode analysis confirms that the boundary integral in {AM
could be neglected from the start when dealing with metallic domain enclosures. The
final system for the \320\242\320\225
and TM analysis is then obtained by enforcing the boundary
conditions t/(r) only. For the \320\242\320\225
on case, testing is imposed on the boundary and
interior without any specification
nodes for the boundar/ values of U(t). However,
for TM analysis, the boundary values of U(r)areset to zero a priori. Thus testing at
the boundary nodes is avoided and the final TM system is smaller than that corre-

corresponding to the case.

\320\242\320\225

4.4 TWO-DIMENSIONAL SCATTERING

As another example application of solving the two-dimensional wave equation, w<

consider the scattering by a cylinder shown in Fig. 4.2. We shall assume H. incident
as given by D.6) and as based on the discussion in Section 4.2.2, it is best to work
with the scattered rather than the total field. From D.11). D.29), and D.34),the
pertinent weak form of the wave equation is
Section 4.4 Two-Dimensional
\302\246 Scattering 119

T'M da
|

+ - f
n \302\246
[Link](r)] dl=\\l (is D.70)
\320\251\320\263)\320\224\320\263)
Fr J(.\302\253w+ J

where

/(r) = -V \342\200\242
\320\224 1 +ja>E0Mi:

is the excitation function and cinncroulcr represent the closure boundaries of the
computation domain, as depicted in Fig. 4.13. For Hi = and
\320\265\320\264\"(\0208*\"+1!\"\320\277\320\233)

Mt = 0. it follows that

\321\217; D.71)

in the eth element, where we approximated \302\253r and by their

\321\206, average values eer and
M^, respectively, over the eth element. Also, we note the presence of the boundary
integral over =
\320\241 Coulcr + Cmncr and from D.33) it is necessary that encloses
\320\241 the
entire computational domain (see Fig. 4.13for the directions of the normals on each
contour).

To discretize D.70), we expand the scattered field as

D-72)
\320\201 E
i i i

where Hi* denotes the unknown scattered field values at the nodes and Nf(x,y) is
given by D.36). Subsequently, on choosing Galerkin's testing (i.e., W = NJ), we
obtain the element equations
r -1
f f \320\2231 1
\320\257*' \302\246 + tiVrNiNJ dx
E E VNi
\342\200\224L
e<-
VN> dy
/=l
\342\204\226| J-)S2,L J

i \342\200\242 - \302\246
e
' Vtff\302\260l]dl
\320\233'\320\224\320\275 j '6
N}[n V//f\"]dl

D.73)
,.=i

It has been assumed that the interior domain (dielectricand free space)
bounded by c\"\"\"\" and has
\320\241\321\210\" been subdivided into triangles, whereas the con-
contours c\302\260\"lerand been subdivided into N,{ and Nx, line segments, respec-
c\"nner have
respectively. Note also the excitation function/*' will be nonzero
that only when er and \321\206\320
are not equal to unity, i.e.. only if the element is within the material region \320\257,/,
illustrated in Fig. 4.13.

4.4.1Treatmentof Metallic Boundaries

In the previous section, we eliminated the boundary integral over \320\241\321\202\321\210

since

h \302\246
V//;tal was zero on the metal boundary. However, this is not the case here because
120 Two-Dimensional Applications \302\246
Chapter 4

-:x/\321\217
'
. , \321\201-

outer

Inner
Q Figure 4.13 Illustration of the eonii\302\273
and
C\302\260ula as well as boundary
\320\241\342\204\242\" vcui.-r.
* for the scattering geometry.

Hf*{ represents the scattered and not the total magneticfield. SinceH. = H'.+ \320\257!6\".
it follows that

.
\320\246
(vhI + V//.scaI) = n \342\200\242
VH, = 0, re Cinner D.74i

and thus

n \302\246
VH-scal = -n \342\200\242
VHt, r e CinMr D.75i

which is an inhomogeneous Neumann boundary condition, where n \302\246

VH'. is a

known quantity everywhere on C\"[Link] alternative form of D.75) can be obtained

by noting that

(see also A.76)) or

where EUin refers to the total tangential electric field on cmn\". Similarly, for the

scattered field

\".? D.77,
Z,, \320\262\320\263

and since Ewn = /

\302\246
(E' 4- Escal) = 0 on conducting surfaces, it follows that

- VtfP\" n =
\302\246
E'ttn r e Cinntir D.7h
+\320\244

with EJan
= Zfl;-(.vsin0o-Pcos0ok/*(l(lfClls*e+1'sin'*ul. Thus, the boundary integral
over Cnncr can be moved to the right hand side of D.73) to be included as pan of
the excitation column. This type of detail in treating a boundary condition (or
constraint) demonstrates how a knowledge of the field over a boundary becoi\302\273

an equivalent source excitation. It was discussed in Section 3.7 for one-dimensional

applications and is referred to as condensation of boundary conditions.
Section 4.4 Two-Dimensional
\302\246 Scattering 121

4.4.2 Absorbing BoundaryConditions

Before casting D.73) onto a matrix system, we must also consider the boundary
conditions on So far,
\321\201\302\260\321\210'. no information has been given with regard to the
boundary condition that must be satisfied by #fal on [Link] scattering prob-
problems, the field continues to propagate to infinity and at large distances from the two-
dimensional scaUerer it has the form Ae~'kr j-Jr, where r denotes the radial distance
from the origin and A is a constant. Thus, since

1 e~

flfal satisfies the condition

00 D.79)
This is the well-known first-order Bayliss-Turkel [12]absorbingboimdary condition

(ABC) and, by its nature, it can only be enforced on circular boundariessuch as that

shown in Fig. 4.14. in this case, f coincides with the normal h, and we then rewrite
D.79) as

D.80)
\320\255\321\217
^\320\271\342\200\224\320\276

where ra is of the circular

the radius boundary Couler. We observe that D.80) is
identical to the impedance boundary
in form condition given by C.42) and its
treatment in the FEM solution is therefore straightforward.
The purposeof the ABC is to create a boundary which does not perturb the
field that is incident upon it\342\200\224in
effect, to simulate a surface which is actually not
there. A measure on how well an ABC simulates a nonreflecting boundary can be
obtained by examining the reflection field due to a plane wave impinging upon the
ABC surface (see Fig. 4.14). F orthe Bayliss-Turkel ABC D.80), it will be seen that

its approximation of a nonreflecting surface may be poor unless it is placed far A to 2

wavelengths) from the scatterer. The latter is highly undesirable because the

Nonreflecting or
ABC surface

Incident
plane wave Reflected
Figure 4.14. Illustration of a circular ABC wave

twfacc and how to measure its effectiveness. (to be suppressed)

122 Two-DimensionalApplications \302\246
Chaptcr4

computational domain is enlarged substantially leading to unnecessarilylarge maim

systems. In the 1980s (see Seniorand Volakis [13] for a review of ABCs)much work
was carried out. aimed at deriving ABCs which provide a better simulation of non-

reflecting surfaces even when placed at a fraction of a wavelength from the scatterer
[12]. [14], [15], [16].These improved ABCsare associated with higher order tangen-
tangentialderivatives. For
example, D.79) is referred to as a first-order ABC because ii
involves a single derivative (with respect to the tangent) of the [Link] general form
of the second-order ABC is

where t and n are shown in Fig. 4.14. For the second-order Bayliss-TurkelABC [12],
the coefficients a and j8 are given by

a = \\
-Ao(

where again r0 is interpreted as the radius of the boundary contour c\302\[Link] on

their derivation, c\302\260ulcrmust actually be circular for the application of the second-
order Bayliss-Turkel ABC. However, it can be specialized to the case of planar
boundaries by letting oo.
rQ \342\200\224> On constant x boundaries, as in Fig. 4.15, D.81)
becomes

D.83i

with

and a similar expressionis obtained on constant boundaries.

\321\203 This simplified

second-order ABC D.83) is the Engquist and Majda [14]ABC and is extensively
used in connection with Finite Difference-Time Domain solutions.

Figure 4.15 Piecewise planar boundary to:

mesh truncation.
Section 4.4 Two-Dimensional
\302\246 Scattering 123

Regardless of the type of ABCbeingused to simulate a nonreflecting surface at

C*uW, mathematically the ABC replaces the normal field derivativeswith tangential
ones. In contrast to normal
derivatives, tangential derivatives can be carried out
using the approximate since there is no requirement
field expansion for information
outside the computational domain. Specifically, we have

f 1 f I / \320\257'\320\235\320\226\320\2331\\
-
Nf[n \342\200\242
Vtffdl] dt = - Nf I <5\320\257^\"
' 4- \320\254
-^~ dt
J (\320\274\320\275\320\270\321\202
Er J (-outer Er \\ at' I

' \" D.84)

](\321\216\320\270\321\216?\320\263 Jf-omer Br at at

where we used integration by parts to transfer one of the derivativesfrom the field to
the testing function as was done by obtaining the weak form of the wave equation.
To proceed with the discretization of D.84) we must introduce an expansion for the
field on the boundary C\"utcr. Choosing the linear expansion basisD.58)or D.59), we
have

and upon substitution into D.84) we obtain

- \302\246
f Nf[n Vtff'l]dl

4 4- *]
As noted before, L-'(r) = \320\233^(\320\263)
when e Cmcr
\320\263 is an edge on C\302\260ulcr
and \320\273-, belonging
to the eth element, as depictedin Fig. 4.7.
With the evaluation of the boundary integral over C\"nner as given by D.78) and
the discretization of the other over C\302\260ul\"as given by D.85), we are ready to cast

D.73) into a linear system of equations. The element equations are

N, \302\253\302\253, \320\232
- D.86)
?\320\270\321\207(\320\2421}'+\320\264^\320\275/\320\271\321\202 J2ibe)
C=l ,v,=l f=l

The entries of the [Ae] matrix are again given by [see D.46)]

4 = - + + tidr A + fy) D-87)

^\320\2647
<$\321\204
\320\244\320\251
^
and

roe\021

in which the angular cylindrical variable \321\204

has been used in the integral to emphasize
that the employed ABC is only applicable to circular enclosures\321\201\302\260\321\210\321\20
(note:
dt = r0 \321\2011\321\204).
124 Two-Dimensional Applications \302\246
Chapter 4

The excitation column entries consist of two components\342\200\224one from the exci-
excitation function/(r) and another from the boundary condition on c1\"\"\". Specifically
[again, these entries are not related to the b' in D.87)],

% = f f dxdy
\320\222\320\224\320\223 -Jp \\
Nf(r)(E'tan
\342\200\242
!)dt D.89)

where s2 is a segment on cmner belonging to the eth element. Substituting for/1\" and

E(an, a more explicit expression for b' is

= $ A
- \320\233 \\ Nf(r) *
[

-7*o NflrKx sin 0o - vcos\321\2040)

\302\246
?* D.901
f

which can be evaluated in closed form for each of the line segments and triangles.
The assembly of the elemental equations D.86) is carried out in the saint
manner as done for the TM waveguide mode analysis. The resulting global system
will be of the form

D\320\233

]|54[ ]|)w
where we used the superscripts \"interior\" and \"boundary\" to indicate the separation
of the node field column as done in D.65).
This system is similar in all respects to D.65) for the TM mode analysis
However, there is a major difference in that has
(\320\244) now been replaced with the
field itself on the boundary. Thus, the convenient decomposition to a pair of smaller
systems is no longer possible, nor is it needed. Since the ABC permitted the elimina-
no additional
of {\320\244},
elimination equations are required for the solution of (\320\257?\"\"}.Addition
of the two left hand matrices gives

[\320\231 |

where the entire matrix is sparse and the system can be solved using an iterative
solver (see Chapter 9).

4.4.3 ScatteredField Computation

Once or ?fM is computed from

Hfal a solution of D.92), the scattered field
outside the
computational domain and echowidth of the scatterer are two observa-
bles of [Link] echowidth is obtained from the far zone scattered field using the
relation

=
Jim 2*1--^-
(\302\246-\302\246\320\236\320\241
D.9.T,

The field outside the computational domain is obtained by application of the surface

equivalence principle or KirchhorTs integral equation. In either case, the pertinem

expression is (see Chapter 1, Section1.4)

Section 4.4 Two-Dimensional
\302\246 Scattering 125

?/*=>) = -I \302\246
(\302\253'V't/(r') G2/)(r. r') - \302\246
[\302\253'V Gw(j, r')] U(r')} dl' D.94)

where U is equal to \320\225\320\263

or H., depending on the excitation. The contour CK can be of

any shape or form and can be located at any distance from the scatterer. Also,
G^Cr,r') is the two-dimensional Green's function

. ,..u,-
..\342\200\236 -r'l) D.95)
4
where H^^ denotes the Hankel
zeroth-order function of the second kind. This
Green'sfunction can be
interpreted as the field generated by a line source at r'
since it satisfies the differential equation
klG1D= -8{r- r') VzG2l) + D.96)

The scattered field representation D.94) is identical to that given in Chapter 1.

Also, D.94) can be related to the radiated fields due to equivalent current sources.
For example, assuming U = E., from Maxwell's equations [see A.76)]

J = x H
\320\264 =
-j^p- n x [f x V?J D.97)

and since V?. = + ?(\320\264?./\320\250),

\320\277(\320\264\320\225./\320\222\320\277) it follows that x t
(\302\253 = f)

L D.98)
J'd4
Also, the equivalent magnetic current is given by

M = E x =
\321\217 ?? = lM, D.99)
and therefore D.94) can be rewritten as

?!\302\260V)
= i l-jk0Z0J.(r)G2D(r.,') + M,[n \302\246
V'C20(r. r')]} dl' D. 100)

Although C/c can be arbitrary, it is best to choose it so that the integrated fields
are most accurate. For purely metallic scatterers, the computed field and its deriva-
is
derivative most accurate near the metallic surface. Thus, it is appropriate to choose CK
to be near to or coincide with the metallic surface of the scatterer. For TM incidence,
when CK is coincident with the metallic surface. ?/(r)= 0 for \320\263
\320\261 and
\320\241\320\272, thus

\302\246 - D.101)
[n V'Ez(r')}H^(k0\\r r'\\)dl'

Likewise,for \320\242\320\225
incidence, h \302\246
VU(r)
- 0 on CK and D.94)reducesto
1-
\320\2575>)
= r'\\)]dl' D.102)
|
However, when dealing with coated metallic or purely dielectric scatterers, D.101)
and D.102) cannot be [Link] this case, the contour CA- is placed above the outer
surface of the dielectric (see Fig. 4.16)and the scattered must be computed from
D.100)or its dual.

To evaluate the scattered field using D.94) or D.100),we proceed by first

discretizing CK into Nx segments each less than 1/10 of wavelength or as dictated
126 Two-DimensionalApplications \302\246
Chapter 4

by the original FEM tessellation of the domain. For simplicity, let us choose Q to

be a circle of radius |r'| = rK. Then

D.1

with

^
\320\244(\320\263')
- 2rrK cos(<p
-
*'
D.104!
= \320\277
\320\244(\320\263)
\342\200\242
VU(t)

2, + r - 2rrK cos@
- <
lA2)(k0Jr2
_ \320\224\320\276
r r u(x'\\
V
4 AJo + ~
y/r2 4
D.1051

and in D.105) we used the simplifications2

\320\263-\320\263'

\320\257'2)(*\320\276|\320\263-\320\263'|)
= -\321\204')] D.106)
-\302\253\320\276 [\320\263\320\272-\320\263\321\201\320\276$(\321\204
|\320\263-\320\263'|

in which \320\257{2)denotes the first-order Hankel function of the second kind. Also, q
refers to the angle between r' and the .v-axis, as shown in Fig. 4.16.
For far zone computations (i.e.,r -*\302\246 oo), we can simplify the above integrands

by introducing the large argument approximations of the Hankel functions.

Specifically,

Figure 4.16 Illustration of the inicgraiioii

contour for computing
\320\241\321\206 the scattered fidi

2Notc: BGzp/ait = n \302\246

VG2* = ^ (n
\302\246
Wj2)
V\302\253) =
(*o \320\257)
^
\342\200\242
\320\224)
(\302\253 H^\\k0R), where VJ? = R -.
Seciion 4.4 Two-Dimensional
\302\246 Scattering 127

D.107)

and we approximate the radical as

77 TT
\342\200\224