FAST SOLVER FOR LARGE SCALE EDDY CURRENT NON-DESTRUCTIVE
EVALUATION PROBLEMS
By
Naiguang Lei

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Electrical Engineering – Doctor of Philosophy
2014

ABSTRACT
FAST SOLVER FOR LARGE SCALE EDDY CURRENT NON-DESTRUCTIVE
EVALUATION PROBLEMS
By
Naiguang Lei
Eddy current testing plays a very important role in non-destructive evaluations of conducting
test samples. Based on Faraday’s law, an alternating magnetic field source generates induced
currents, called eddy currents, in an electrically conducting test specimen. The eddy currents
generate induced magnetic fields that oppose the direction of the inducing magnetic field
in accordance with Lenz’s law. In the presence of discontinuities in material property
or defects in the test specimen, the induced eddy current paths are perturbed and the
associated magnetic fields can be detected by coils or magnetic field sensors, such as Hall
elements or magneto-resistance sensors. Due to the complexity of the test specimen and
the inspection environments, the availability of theoretical simulation models is extremely
valuable for studying the basic field/flaw interactions in order to obtain a fuller understanding
of non-destructive testing phenomena.
Theoretical models of the forward problem are also useful for training and validation of
automated defect detection systems. Theoretical models generate defect signatures that are
expensive to replicate experimentally. In general, modelling methods can be classified into
two categories: analytical and numerical. Although analytical approaches offer closed form
solution, it is generally not possible to obtain largely due to the complex sample and defect
geometries, especially in three-dimensional space. Numerical modelling has become popular
with advances in computer technology and computational methods. However, due to the
huge time consumption in the case of large scale problems, accelerations/fast solvers are

needed to enhance numerical models.
This dissertation describes a numerical simulation model for eddy current problems using
finite element analysis.

Validation of the accuracy of this model is demonstrated via

comparison with experimental measurements of steam generator tube wall defects. These
simulations generating two-dimension raster scan data typically takes one to two days on
a dedicated eight-core PC. A novel direct integral solver for eddy current problems and
GPU-based implementation is also investigated in this research to reduce the computational
time.

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . .
Chapter 1 Introduction . . . . . . . . . . . . .
1.1 Methods of Electromagnetic Non-destructive
1.1.1 Magnetic Flux Leakage Testing . . .
1.1.2 Eddy Current Testing . . . . . . . .
1.1.3 Microwave Non-Destructive Testing .
1.2 Principles of Eddy Current Testing . . . . .
1.3 Motivation and Organization . . . . . . . . .

Part I Finite Element
Non-Destructive Evaluation

Method

. . . . . . .
Evaluation
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .

in

.
.
.
.
.
.
.

.
.
.
.
.
.
.

Eddy

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

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

xi
1
1
1
2
3
3
6

Current
8

Chapter 2 Finite Element Method for Eddy Current
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . .
2.2 Problem Statement . . . . . . . . . . . . . . . . . . .
2.3 Formulations for Eddy Current Problems . . . . . . .
2.3.1 Formulations in Eddy Current Free Regions .
2.3.1.1 ψ Formulation . . . . . . . . . . . .
2.3.1.2 A Formulation . . . . . . . . . . . .
2.3.1.3 Ar Formulation . . . . . . . . . . . .
2.3.2 Formulations in Eddy Current Regions . . . .
2.3.2.1 T , ψ Formulation . . . . . . . . . . .
2.3.2.2 A∗ Formulation . . . . . . . . . . . .
2.3.2.3 A, V Formulation . . . . . . . . . . .
2.3.2.4 Ar , V Formulation . . . . . . . . . .
2.3.3 Coupling of Formulations . . . . . . . . . . . .
2.3.3.1 Direct Coupling . . . . . . . . . . . .
2.3.3.2 Coupling A to ψ . . . . . . . . . . .
2.3.3.3 Coupling Ar to ψ . . . . . . . . . . .
2.3.3.4 Coupling A to T , ψ . . . . . . . . . .
2.4 Governing Equations . . . . . . . . . . . . . . . . . .
2.4.1 Ar , V − Ar Formulation . . . . . . . . . . . .
2.4.2 System Equations . . . . . . . . . . . . . . . .
2.5 Modelling of Eddy Current Probes with Ferrite Core
iv

Simulation
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .

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

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

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

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

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

9
9
10
12
12
13
13
14
14
15
16
16
17
18
18
18
18
19
20
20
21
23

2.6

Parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

Chapter 3 Application to Steam Generator Tube Inspection . . . . . . . .
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Eddy Current Probes Simulation and Experimental Validations . . . . . . .
3.2.1 Rotating Probe Coil – Pancake Coil Probe . . . . . . . . . . . . . . .
3.2.2 Rotating Probe Coil – Plus Point Probe . . . . . . . . . . . . . . . .
3.2.3 Array Coils Probe - Array Probes . . . . . . . . . . . . . . . . . . . .
3.2.4 Quantitative Comparison between Experimental and Simulation Signals
3.2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Simulation Software Features . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4 New Probe Design and Simulation . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 The Rotating Field Eddy Current Probe Operational Principles and
Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25
25
26
27
28
30
33
34
34
35

Part II

42

Direct Integral Solver and Fast Algorithms

35
38
40
41

Chapter 4 Moment Method . . . . . . . . . . . . . . . .
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Volume Integral Equation Method . . . . . . .
4.1.1.1 Bounded Eddy Current Problems . .
4.1.1.2 Unbounded Eddy Current Problems
4.1.2 Surface Integral Equation Method . . . . . . .
4.2 Surface Integral Equations . . . . . . . . . . . . . . .
4.2.1 Electric Field Integral Equation . . . . . . . .
4.2.2 Magnetic Field Integral Equation . . . . . . .
4.3 High-order Basis Function . . . . . . . . . . . . . . .
4.3.1 RWG Basis Function . . . . . . . . . . . . . .
4.3.2 High-order RWG Basis Function . . . . . . . .
4.4 Moment Method Procedure . . . . . . . . . . . . . .
4.5 Singularity of Green’s Function . . . . . . . . . . . .
4.6 Low Frequency Breakdown of SIE . . . . . . . . . . .
4.7 Validation . . . . . . . . . . . . . . . . . . . . . . . .
4.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . .

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

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

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

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

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

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

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

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

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

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

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

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

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

43
43
43
43
45
53
54
55
57
58
58
59
60
61
63
65
70

Chapter 5 Fast Algorithms for Integral Equation
5.1 Matrix Compression Algorithms . . . . . . . . . .
5.1.1 Low-Rank Matrix . . . . . . . . . . . . . .
5.1.2 Adaptive Cross-Approximation . . . . . .
5.1.3 Randomized Algorithm . . . . . . . . . . .
5.1.3.1 Interpolative Decomposition . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

71
71
71
72
74
74

v

.
.
.
.
.
.

.
.
.
.
.
.

5.1.3.2

5.2

5.3

5.4

Converting Interpolative
Decomposition . . . . .
Adaptive Integral Method . . . . . . . .
5.2.1 FFT-based Method . . . . . . . .
5.2.2 Adaptive Integral Method . . . .
Fast Multipole Algorithm . . . . . . . .
5.3.1 Middle Frequency Regime . . . .
5.3.2 Low Frequency Regime . . . . . .
5.3.3 Multilevel Implementation . . . .
Conclusion . . . . . . . . . . . . . . . . .

Decomposition into
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .
. . . . . . . . . . .

Chapter 6 Direct Integral Solver . . . . . . .
6.1 Matrix Representation for Fast Algorithms
6.1.1 Single level Matrix Representation
6.1.2 Multilevel Matrix Representation .
6.2 Direct Integral Solver . . . . . . . . . . . .
6.3 Analysis of Computation Complexity . . .
6.3.1 Computational Scalability . . . . .
6.3.2 Matrix Compression . . . . . . . .
6.4 Conclusion . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

Singular Value
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .

75
76
76
77
79
80
81
84
85

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

86
86
86
87
88
89
90
91
91

Chapter 7 Summary and Future Work . . . . . . . . . . . . . . . . . . . . .
7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 Finite Element Method . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1.1 Refactoring the existing code . . . . . . . . . . . . . . . .
7.2.1.2 Optimization . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.2 Direct Integral Solver . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.2.1 GPU Implementation of the Matrix Compression Method
7.2.2.2 Optimization of Moment Method Galerkin Testing . . . .
7.2.2.3 Multi-body Problems . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.

93
93
94
94
94
95
95
95
96
96

APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
BIBLIOGRAPHY

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

vi

. . . . . . . . . . . . . . . . . 106

LIST OF TABLES

Table 3.1

Modelling geometry for pancake coil probe . . . . . . . . . . . . . .

27

Table 3.2

Modelling geometry for plus point probe . . . . . . . . . . . . . . .

28

Table 3.3

Modelling geometry for array probes . . . . . . . . . . . . . . . . . .

30

Table 3.4

Quantitative comparison of rotating probe coil . . . . . . . . . . . .

33

Table 3.5

Quantitative comparison of bobbin probe . . . . . . . . . . . . . . .

34

Table 3.6

Quantitative comparison of X-probe . . . . . . . . . . . . . . . . . .

34

Table 4.1

Parameters of test experiment No. 1 . . . . . . . . . . . . . . . . . .

66

Table 4.2

Parameters of test experiment No. 2 . . . . . . . . . . . . . . . . . .

67

Table 6.1

Computation time of matrix A . . . . . . . . . . . . . . . . . . . . .

90

Table 6.2

Matrix compression: memory usage and compression time . . . . . .

91

vii

LIST OF FIGURES

Figure 1.1

Magnetic flux leakage testing . . . . . . . . . . . . . . . . . . . . . .

2

Figure 1.2

Principles of eddy current testing . . . . . . . . . . . . . . . . . . .

5

Figure 2.1

An arbitrary geometry with closed surface . . . . . . . . . . . . . . .

11

Figure 2.2

Modelling method in case of ferrite core . . . . . . . . . . . . . . . .

24

Figure 3.1

Illusion of steam generator . . . . . . . . . . . . . . . . . . . . . . .

26

Figure 3.2

Impedance of pancake coil probe on a raster scan around the defect

27

Figure 3.3

Illustration of plus point probe . . . . . . . . . . . . . . . . . . . . .

28

Figure 3.4

Impedance of plus point probe on a raster scan around the defect
(2-D plot) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Impedance of plus point probe on a raster scan around the defect
(3-D plot) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

Figure 3.6

Illustration of array probes . . . . . . . . . . . . . . . . . . . . . . .

30

Figure 3.7

Real part of axial channel impedance of array probes . . . . . . . . .

31

Figure 3.8

Imaginary part of axial channel impedance of array probes . . . . .

31

Figure 3.9

Real part of circumferential channel impedance of array probes . . .

32

Figure 3.10

Imaginary part of circumferential channel impedance of array probes

32

Figure 3.11

Principle of rotating field . . . . . . . . . . . . . . . . . . . . . . . .

36

Figure 3.12

Rotating field eddy current probe structure . . . . . . . . . . . . . .

38

Figure 3.13

Signal of rotating field eddy current probe with bobbin pickup coil .

39

Figure 3.5

viii

Figure 3.14

Signal of rotating field eddy current probe with GMR array sensors

39

Figure 3.15

Signals of parametric studies . . . . . . . . . . . . . . . . . . . . . .

40

Figure 4.1

A conducting body excited by a line conductor carrying current . . .

44

Figure 4.2

Unbounded geometry examples . . . . . . . . . . . . . . . . . . . . .

47

Figure 4.3

An arbitrary geometry with closed surface . . . . . . . . . . . . . . .

53

Figure 4.4

Illusion of RWG basis function . . . . . . . . . . . . . . . . . . . . .

59

Figure 4.5

Vector plot of first-order RWG basis function . . . . . . . . . . . . .

59

Figure 4.6

Diagram for the measurement of ∆Z due to a surface breaking slot .

65

Figure 4.7

Variation of coil impedance change with scanning location above an
aluminium alloy plate . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Geometry of a rotating probe coil inside an inconel steam generator
tube with an axial through-wall notch . . . . . . . . . . . . . . . . .

67

Variation of coil impedance change with scanning location inside an
inconel steam generator tube . . . . . . . . . . . . . . . . . . . . . .

68

Figure 5.1

Projection of triangular RWG basis onto rectangular grids . . . . . .

78

Figure 5.2

Direction from rn to rm broken into three segments: d1 , D, d2 . . .

79

Figure 5.3

Illusions of groups . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

Figure 7.1

A two-body surface integral equation problem . . . . . . . . . . . .

96

Figure 1

Block diagram of simple operation of steam generator tube simulator
software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

Figure 2

Screen-shot showing various commands available in the simulation
software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

Figure 3

Analysis screen for the rotating probe coil . . . . . . . . . . . . . . . 103

Figure 4.8

Figure 4.9

ix

Figure 4

Display of the 1-D coil analysis screen (bobbin probe) . . . . . . . . 104

x

KEY TO SYMBOLS AND ABBREVIATIONS

EM

electromagnetic/electromagnetics

NDE

non-destructive evaluation

NDT

non-destructive testing

PDEs

partial differential equations

MFL

magnetic flux leakage

EC

eddy current

ECT

eddy current testing

PoD

possibility of detection

FE

finite element

FEM

finite element method

SG

steam generator

RPC

rotating probe coil

RoFEC

rotating field eddy current

GMR

giant magneto-resistance

MoM

moment method/method of moments

VIE

volume integral equation

SIE

surface integral equation

BEM

boundary element method

EFIE

electric field integral equation

MFIE

magnetic field integral equation

RWG

Rao-Wilton-Glisson

SVD

singular-value decomposition

CA

cross-approximation

ACA

adaptive cross-approximation
xi

ID

interpolative decomposition

DFT

discrete Fourier transform

FFT

fast Fourier transform

AIM

adaptive integral method

FMA

fast multipole algorithm

MLFMA

multilevel fast multipole algorithm

RRQR

rank-revealing QR

xii

Chapter 1
Introduction
Electromagnetic (EM) methods of non-destructive evaluation (NDE), also called
non-destructive testing (NDT), comprise a full spectrum of techniques ranging from magnetic
flux leakage (DC, static fields) to eddy current (quasi-static frequencies) and microwave (high
frequency). In each of these techniques, the underlying physics is fundamentally different as
the fields are described by different classes of partial differential equations (PDEs). In the
static regime, we deal with elliptic PDEs, whereas in eddy current problems the underlying
fields are described by parabolic PDEs and in microwave NDE, we have hyperbolic PDEs.
Consequently, this affects the different aspects of NDE including the modelling of forward
problem and solution to corresponding inverse problems.

1.1

Methods

of

Electromagnetic

Non-destructive

Evaluation
1.1.1

Magnetic Flux Leakage Testing

Magnetic flux leakage (MFL) is a magnetic method of NDT that is used to detect corrosion
and pitting in steel structures, such as pipelines and storage tanks. The basic principle is
to magnetize the steel part using a permanent or electromagnet. At areas where there is
corrosion or missing metal, the magnetic field “leaks” from the steel. In an MFL tool, a

1

magnetic field detector is used to detect the leakage field.
flux leakage

Figure 1.1: Magnetic flux leakage testing.

As an MFL tool navigates the pipeline a magnetic circuit is created between the pipe wall
and the tool. Brushes typically act as a transmitter of magnetic flux from the tool into
the pipe wall. High field MFL tools saturate the pipe wall so that the reluctance is very
high. Strategically placed tri-axial Hall effect sensor heads can accurately measure the
three-dimensional magnetic flux density vector.

1.1.2

Eddy Current Testing

Eddy current (EC) NDE uses EM induction to detect flaws in conductive materials.
An alternating current source generates changing magnetic field which interacts with a
conducting specimen to induce eddy currents (ECs) in it.

The induced ECs and the

associated magnetic fields opposes the direction of the inducing current and field. This affects
the net flux linking the coil and hence the impedance of the probe coil. Variations in the
phase and magnitude of these ECs can be monitored by measuring the terminal impedance
of the probe coil. Variations in the electrical conductivity or magnetic permeability of the
test object, or the presence of any flaws, will cause a change in ECs and a corresponding
change in the phase and amplitude of the measured impedance.
EC testing (ECT) can detect very small cracks on or near the surface of the material, the
2

surfaces need minimal preparation, and physically complex geometries can be investigated.
It is also useful for making electrical conductivity and coating thickness measurements. The
testing devices are portable, provide immediate feedback, and are non-contact methods,
allowing rapid inspection of the part.

1.1.3

Microwave Non-Destructive Testing

The term microwaves refer to alternating EM waves frequencies in the range from 300
MHz to 100 GHz. Since the penetration of microwaves in good conducting materials is
very small, microwave NDT techniques are mainly used for non-metallic materials. The
spatial resolution of these techniques depends on the wavelength of the EM wave. In the
microwave band of 3-100 GHz, wavelength varies from 100 mm to 3 mm. These techniques
have advantages over other NDE methods (such as radiography, ultrasonics and EC) and
offer good penetration in non-metallic materials, good resolution and non-contact nature
feature of the microwave sensor (antenna).

This dissertation focuses on EC methods of NDT at quasi-static frequencies, where
the penetration depth in conducting parts is of the order of millimeters or a few centimeters.

1.2

Principles of Eddy Current Testing

The basic principle underlying ECT can be illustrated with a simple arrangement shown
in Figure 1.2a. When a coil carrying an alternating current is brought in close proximity
to an electrically conducting test specimen, the alternating magnetic field causes induced
currents in the conducting test specimen in accordance with Faraday’s law of EM induction.

3

The induced currents are called ECs since they follow closed circulatory patterns that are
similar to eddies found in water bodies. The alternating ECs, in turn, establishes a field
whose direction is opposite to that of the original or primary field. Consequently, the net
flux linkages associated with the coil decrease. Since the inductance of a coil is defined as
the net flux linkages per ampere, the effective inductance of the coil decreases relative to
its value if it were to be suspended in air. The presence of ECs in the test specimen also
results in a resistive power loss. The effect of this power loss manifests in the form of a small
increase in the effective resistance of the coil. An exaggerated view of the changes in the
terminal characteristics of the coil is shown in Figure 1.2b, where the variation in resistance
and inductance is plotted in the impedance plane. When a flaw or inhomogeneity, whose
conductivity and/or permeability differ from that of the host specimen, is encountered, the
current distribution is altered. Consequently, the impedance of the coil is different relative
to the value observed for flawless regions, as shown in Figure 1.2b. Systems that are capable
of monitoring the changes in impedance of probe coil can, therefore, be used to detect flaws
in a specimen.
ECs exhibit a unique phenomenon known as the “skin effect” which causes the current
density at a particular depth to decrease with an increase in the frequency of excitation.
This result can be derived very simply from the solution of a plane wave propagating into a
conducting half plane. Skin depth (δ), also called standard depth of penetration, is defined
as the depth at which EC density has decreased to reciprocal of Euler-Mascheroni constant
of the surface value. The skin depth can be computed as follows:

1
δ=√
πf µσ

4

(1.1)

where f is the excitation frequency of the circuit, µ is the magnetic permeability of the
target material, and σ is the electrical conductivity of the target material. The skin depth
is often used as a guideline to select the excitation frequency for testing a given specimen.
The variations in coil impedance caused by discontinuities in the test specimen are often
very small in comparison with the quiescent value of the coil impedance. The detection and
measurement of the small changes are often made possible with the use of bridge circuits [1,2].

(a)

(b)

Figure 1.2: Principles of eddy current testing [1, 2]: (a) eddy current generation and flow
in a conducting specimen; (b) impedance-plane trajectory of a coil over a conducting test
specimen.

Factors that influence the EC field, and therefore coil impedance are:
• separation between the coil and specimen surface, called lift-off ;
5

• electrical conductivity of the specimen;
• magnetic permeability of the specimen;
• frequency of the AC;
• design of the EC probe coil;
• specimen geometry factors;
• discontinuities, such as cracks, corrosion, pitting.
The central problem in NDT is flaw detection and characterization. Successful detection
and characterization of flaws can benefit from a clear understanding of field/flaw interaction
and effects of various probe and operational parameters. Correlations between operational
parameters and signal features help in signal processing design procedures to compensate
for these effects and eliminate undesired responses. Correlations between defect parameters
and signal features can help in designing more accurate defect detection and characterization
algorithms.

1.3

Motivation and Organization

Theoretical models of the forward problem are necessary for studying the basic field/flaw
interactions in order to obtain a full understanding of NDT phenomena. Theoretical models
generate defect signatures that are expensive to replicate experimentally. The model provides
solutions to the forward problem under controlled variation of experimental parameters.
Such studies are important to use the technology effectively and maximize the possibility of
detection (PoD) of critical flaws.

6

Modelling methods can be classified into two categories: analytical and numerical. Although
analytical approaches offer closed form solutions, it is generally not possible to obtain
an analytical solution in the case of complex sample and defect geometries, especially in
three-dimensional (3-D) space. Numerical modelling has become popular with advances in
computer technology and computation methods. However, due to the huge time consumption
for large scale problems, accelerations/fast solvers are needed for obtaining practical solutions
using numerical models. This is the motivation for the work carried out in this dissertation.
This dissertation is organized into two parts: part I. finite element methods in EC NDE, part
II. direct integral solver and fast algorithms. In the first part, the finite element method is
introduced, and the application to steam generator inspection and new EC probe design are
presented. In the second part, several fast algorithms for microwave problems are introduced,
and the fast algorithms and a novel direct integral solver for large scale EM NDE problems
are investigated.

7

Part I
Finite Element Method in Eddy
Current Non-Destructive Evaluation

8

Chapter 2
Finite Element Method for Eddy
Current Simulation

2.1

Introduction

Among various numerical modelling methods, the finite element method (FEM) is one of
the most commonly used techniques. The reasons for this are:
• FEM has well-established and relatively simple formulations;
• results in sparse and banded matrices that can be efficiently stored and solved.
The fundamental principle underlying finite element (FE) analysis is the replacement of a
continuous domain by a number of sub-domains (elements) in which the unknown function is
represented by simple interpolation functions with unknown coefficients. A set of algebraic
equations are obtained by applying the Ritz or Galerkin method [3]. Finally, solution of the
boundary value problem is achieved by solving the system of equations.
Since it was first proposed in the 1940s, the FEM has been developed and applied extensively
to structural analysis and other engineering fields. The FEM was popularized in the area
of EM NDT by Lord etc. [4–7]. Lord and Palanisamy [8, 9] developed a 2-D axisymmetric
model to study EC phenomena in steam generator (SG) tube inspection for the nuclear power
industry. Over the past five decades, many researchers have developed FE models, such as the
9

one reported in [10], for the 3-D simulation of EC problems. There are commercial software
packages, such as Multiphysics

produced by COMSOL, available for 3-D FE simulation of

general EC problems.

2.2

Problem Statement

In frequency domain, the Faraday’s law and Amp`ere’s law are given by

∇ × E = −jω B

(2.1)

∇ × H = J inc + J + jω D

(2.2)

B = µH

(2.3)

D= E

(2.4)

J = σE

(2.5)

∇ · J + jωρ = 0

(2.6)

with the constitutive relations

and the continuity equation is given by

where
E – electric field intensity
D – electric flux density

10

H – magnetic field intensity
B – magnetic flux density
J inc – source current density
J – conduction current density
ρ – electric charge density
– permittivity
µ – permeability
σ – conductivity
In general, the EC problem can be represented in terms of a conducting domain Ωc and
surrounding non-conducting domain Ωn , as shown in Figure 2.1.

n
ˆ

Ωc
Ωn

Figure 2.1: An arbitrary geometry with closed surface.

For EC problems, by ignoring the displacement current, the quasi-static regime equations
are represented as:
in Ωn ,

∇ × H = J inc

11

(2.7)

in Ωc ,

∇ × H = σE

(2.8)

∇ × E = −jωµH

(2.9)

∇·J =0

(2.10)

From equation (2.7), we can represent the magnetic field due to the source current as,

∇ × Hs = J inc
˚
J inc (r ) × ∇
Hs (r) =

2.3

(2.11)
1
dV
|r − r |

(2.12)

Formulations for Eddy Current Problems

There are several choices of formulations in terms of scalar and vector potentials for EC
problems [11]. The purpose of this section is to review the commonly used formulations of
EC problems. Appropriate formulation is selected based on the material property of the
test specimen, the quantity of interest, the limitations of computational resources and so on.
These formulations in both EC free regions and EC regions are introduced as following.

2.3.1

Formulations in Eddy Current Free Regions

In the EC free region, the quasi-static magnetic field can be represented in two ways: 1. a
magnetic scalar potential ψ; 2. a magnetic vector potential A.

12

2.3.1.1

ψ Formulation

From equations (2.7) and (2.11), we have

H = Hs − ∇ψ

(2.13)

and using the magnetic Gauss’ law (∇ · B = 0), we obtain

∇ · (µ∇ψ) = ∇ · (µHs )

(2.14)

ψ = 0 or n
ˆ · (µ∇ψ) = n
ˆ · (µHs )

(2.15)

with the boundary conditions

2.3.1.2

A Formulation

Considering the magnetic Gauss’ law, we can represent magnetic flux density as

B =∇×A

(2.16)

and using equations (2.3) and (2.7), we obtain

1
∇ × ( ∇ × A) = ∇ × Hs
µ

(2.17)

1
n
ˆ × A = 0 or n
ˆ × ( ∇ × A) = n
ˆ × Hs
µ

(2.18)

with the boundary conditions

13

2.3.1.3

Ar Formulation

To represent the source current on the right-hand side only, we introduce the reduced
magnetic vector potential Ar . Considering the magnetic Gauss’ law, we can represent
magnetic flux density as

B = µ0 Hs + ∇ × Ar

(2.19)

and using equations (2.3) and (2.7), we obtain

µ
1
∇ × ( ∇ × Ar ) = ∇ × Hs − ∇ × ( 0 Hs )
µ
µ

(2.20)

with the boundary conditions

1
n
ˆ × Ar = 0 or n
ˆ × ( ∇ × Ar ) = 0
µ

2.3.2

(2.21)

Formulations in Eddy Current Regions

In the EC region, the EM field can be represented by two vector potentials: 1. a current
vector potential T ; 2. a magnetic vector potential A. A magnetic scalar potential ψ is
necessary to be employed with the the current vector potential representation; however, an
electric scalar potential V can be, but not have to be, introduced with the magnetic vector
potential representation.

14

2.3.2.1

T , ψ Formulation

From equation (2.8), we represent magnetic field as

H = Hs + T − ∇ψ

(2.22)

thus, electric field and the induced ECs are

E=

1
1
∇ × Hs + ∇ × T
σ
σ

J =∇×T

(2.23)
(2.24)

and using the magnetic Gauss’ law, equations (2.3), (2.5) and (2.8), we obtain

1
1
∇ × ( ∇ × T ) + jωµT − jωµ∇ψ = −∇ × ( ∇ × Hs ) − jωµHs
σ
σ

(2.25)

∇ · (µT − µ∇ψ) = −∇ · (µHs )

(2.26)

ψ = ψ0 = constant or n
ˆ · (µT − µ∇ψ) = −ˆ
n · (µHs )

(2.27)

1
1
n × ( ∇ × Hs )
n
ˆ × T = 0 or n
ˆ × ( ∇ × T ) = −ˆ
σ
σ

(2.28)

with the boundary conditions

15

2.3.2.2

A∗ Formulation

Considering the magnetic Gauss’ law and equation (2.9), we can represent magnetic flux
density and the electric field as

B = ∇ × A∗

(2.29)

E = −jω A∗

(2.30)

and using equations (2.3), (2.5) and (2.8), we obtain

1
∇ × ( ∇ × A∗ ) + jωσ A∗ = 0
µ

(2.31)

1
n
ˆ × A∗ = 0 or n
ˆ × ( ∇ × A∗ ) = n
ˆ × Hs
µ

(2.32)

with the boundary conditions

2.3.2.3

A, V Formulation

Similarly to A∗ formulation, considering the magnetic Gauss’ law and equation (2.9), we can
represent magnetic flux density and the electric field as

B =∇×A

(2.33)

E = −jω A − ∇V

(2.34)

16

and using equations (2.3), (2.5) and (2.8), we obtain

1
∇ × ( ∇ × A) + jωσ A + σ∇V = 0
µ

(2.35)

with the boundary conditions

2.3.2.4

1
n
ˆ × A = 0 or n
ˆ × ( ∇ × A) = n
ˆ × Hs
µ

(2.36)

V = V0 = constant or n
ˆ · (jωσ A + σ∇V ) = 0

(2.37)

Ar , V Formulation

Similarly to A, V formulation, considering the magnetic Gauss’ law and equation (2.9), we
can represent magnetic flux density and the electric field as

B = µ0 Hs + ∇ × Ar

(2.38)

E = −jω As − jω Ar − ∇V

(2.39)

and using equations (2.3), (2.5) and (2.8), with As shown in equation (2.54), we obtain

1
µ
∇ × ( ∇ × Ar ) + jωσ Ar + σ∇V = −jωσ As − ∇ × ( 0 Hs )
µ
µ

(2.40)

with the boundary conditions

n
ˆ × Ar = −ˆ
n × As

1
µ
or n
ˆ × ( ∇ × Ar ) = −ˆ
n × ( 0 Hs )
µ
µ

V = V0 = constant or n
ˆ · (jωσ Ar + σ∇V ) = −ˆ
n · (jωσ As )

17

(2.41)
(2.42)

2.3.3

Coupling of Formulations

2.3.3.1

Direct Coupling

The coupling of corresponding quantities at the interface between non-conductive domain
and conductive domain is very important to enforce the underling physics. For certain
formulations as discussed in sections 2.3.1 and 2.3.2, we can simply set each quantity equal
to each other directly at the interfaces, such as A∗ –A; A, V –A; Ar , V –Ar and T , ψ–ψ.

2.3.3.2

Coupling A to ψ

If we split the non-conductive domain Ωn into two sub-domains ΩA , Ωψ and represent each
in terms of A and ψ respectively, for ΩA , Ωψ ∈ Ωn , on the interface ΓAψ we have
1
ˆ × Hs − n
ˆ × ∇ψ
n
ˆ × ( ∇ × A) = n
µ
1
n
ˆ · ( ∇ × A) = n
ˆ · Hs − n
ˆ · ∇ψ
µ

2.3.3.3

(2.43)
(2.44)

Coupling Ar to ψ

As in section 2.3.3.2, if we split the non-conductive domain Ωn into to sub-domains ΩAr , Ωψ
and represent each in terms of Ar and ψ respectively, for ΩAr , Ωψ ∈ Ωn , on the interface
ΓAr ψ we have
1
n × ∇ψ
n
ˆ × ( ∇ × Ar ) = −ˆ
µ
1
n
ˆ · ( ∇ × Ar ) = −ˆ
n · ∇ψ
µ

18

(2.45)
(2.46)

2.3.3.4

Coupling A to T , ψ

If we split the conductive domain Ωc into two sub-domains ΩA , ΩT and represent ΩA in
terms of A and ΩT in terms of (T , ψ) respectively, for ΩA , ΩT ∈ Ωc , on the interface ΓAT
we have

1
ˆ × Hs + n
ˆ×T −n
ˆ × ∇ψ
n
ˆ × ( ∇ × A) = n
µ
1
n
ˆ · ( ∇ × A) = n
ˆ · Hs + n
ˆ·T −n
ˆ · ∇ψ
µ

(2.47)
(2.48)

As discussed above, due to the challenges of coupling different quantities on the interfaces,
the most commonly used formulations is A, V –A (or Ar , V –Ar ) and T , ψ–ψ formulation.
The major advantage of T , ψ–ψ formulation over A, V –A (or Ar , V –Ar ) formulation is that a
scalar quantity ψ is used in the non-conducting domain. Thus, the number of the unknowns
in T , ψ–ψ formulation is normally less than that of A, V –A (or Ar , V –Ar ) formulation.
From equation (2.24), the T , ψ–ψ formulation is straight forward for calculating the induced
ECs in conducting domain, however, from equation (2.33), the A, V –A formulation is more
efficient for calculating the magnetic flux density in both conducting and non-conducting
domain. Therefore, if the quantity of interest is induced ECs, the T , ψ–ψ formulation should
be used, and if the quantity of interest is magnetic flux density, the A, V –A formulation is
a better choice.

19

2.4
2.4.1

Governing Equations
Ar , V − Ar Formulation

Among those formulations as discussed in section 2.3, the A, V − A formulation has been
used widely in EC problems [12–15]:

∇ × ν∇ × A + jωσ A + σ∇V = Js
∇ · (jωσ A + σ∇V ) = 0

(2.49)
(2.50)

where ν is the reluctivity, σ the conductivity, ω the angular frequency, and Js the source
current density. A is the magnetic vector potential in the whole solution domain whereas the
electric scalar potential V is used only in conductors. The formulation requires generating
mesh for excitation coil together with the test sample, thus, the re-meshing will be needed
with this formulation when the excitation coil has relative motion to the sample. Re-meshing
should be avoided because: it is difficult and cumbersome; and it results in large computation
error which may be larger than the signal.
The Ar , V − Ar formulation has been proposed to solve the above problem [11, 16–19]. In
this formulation, the magnetic field intensity H is decomposed into Hs and Hr where Hs
is the field intensity in air due to excitation current and Hr is the field intensity due to
induced and/or magnetization currents. Correspondingly, the magnetic vector potential A is
decomposed into As and Ar , where As is the vector potential in air due to excitation current
and Ar is the vector potential due to induced and/or magnetization currents, as shown in

20

following equations:

∇ × ν∇ × Ar + jωσ Ar + σ∇V = ∇ × Hs − νr ∇ × Hs − jωσ As
∇ · (jωσ Ar + σ∇V ) = −∇ · jωσ As

(2.51)
(2.52)

where, νr is the reciprocal of relative permeability. From Biot-Savart’s law, we have,

Hs = ν0 ∇ × As
˚
Js (r )
µ0
dΩ
As =
4π
Ωs |r − r |

(2.53)
(2.54)

The excitation is described by Hs and As only, which avoids mesh generation for the coil.
This has many advantages [11, 16] including, eliminating the need for re-meshing different
coil positions in EC NDE simulation. As a consequence, the resultant system matrix remains
unchanged at different coil positions.

2.4.2

System Equations

By choosing proper shape functions for the quantities Ar and V , we obtain

V (x, y, z) =

Ni (x, y, z)Vi

(2.55)

Ni (x, y, z)Ar,i

(2.56)

i

Ar (x, y, z) =
i

21

Applying Galerkin method on equations (2.51) and (2.52), we have
˚
Ω

i

Ni · ∇ × ν∇ × Ar + jωσ Ar + σ∇V dΩ
˚
=
i

˚
i

Ω

Ω

Ni · (1 − νr )∇ × Hs − jωσ As dΩ

(2.57)

˚

Ni ∇ · (jωσ Ar + σ∇V ) dΩ =

Ω

i

Ni −∇ · jωσ As dΩ

(2.58)

As a result, for the i -th equation,
˚
Ω

˚
Ω

ν(∇ × Ni ) · (∇ × Ar ) + jωσ Ni · Ar + σ Ni · ∇V dΩ−
¨
ν(ˆ
n × Ni ) · (∇ × Ar ) dΓ
∂Ω
˚
=
Ni · (1 − νr )∇ × Hs − jωσ As dΩ

(2.59)

jωσ∇Ni · Ar + σ∇Ni · ∇V dΩ−
¨
jωσˆ
nNi · Ar + σˆ
nNi · ∇V dΓ
∂Ω
˚
Ni ∇ · jωσ As dΩ
=

(2.60)

Ω

Ω

Thus, we get a linear system equations as following:

[G]n×n [x]n = [b]n



GA



GV A

GAV



xA





bA

(2.61)



   
  =  
V
G
xV
bV

22

(2.62)

where,
˚
GA
ij

=
Ω

˚
GAV
ij

=
Ω

ν(∇ × Ni ) · (∇ × Nj ) + jωσ Ni · Nj dΩ−
¨
ν(ˆ
n × Ni ) · (∇ × Nj )dΓ

(2.63)

σ Ni · ∇Nj dΩ

(2.64)

∂Ω

˚

¨

GVij A

=
jωσ∇Ni · Nj dΩ −
jωσˆ
nNi · Nj dΓ
Ω
∂Ω
˚
¨
V
σ∇Ni · ∇Nj dΩ −
σˆ
nNi · ∇Nj dΓ
Gij =
Ω
∂Ω
˚
A
Ni · (1 − νr )∇ × Hs − jωσ As dΩ
bi =
Ω
˚
V
Ni ∇ · jωσ As dΩ
bi =

(2.65)
(2.66)
(2.67)
(2.68)

Ω

2.5

Modelling of Eddy Current Probes with Ferrite
Core

When a probe has a ferrite core, then the situation is more challenging. If only the coil
winding is separated from the FE mesh, then the ferrite core needs to be meshed with the
test sample. In this case, re-meshing is necessary for each coil/core position and the method
discussed in section 2.4 cannot be directly applied. A new modelling method/scheme was
developed to solve this problem [20]. The meshes of the test sample and the ferrite core
are generated separately in this method. The coil is however not meshed. Decomposing
the magnetic field into three parts: H = Hcoil + Hcore + Hsample , where Hcoil , Hcore ,
and Hsample are magnetic field intensities due to coil, core, and test sample, respectively, as
shown in Figure 2.2, this problem can be solved by following procedure:

23

• in the core domain, use Hcoil as excitation and compute Hcore ); in sample domain,
use it as excitation to compute Hsample ;
• if the change of coil impedance calculated from the total field (Hcoil + Hcore + Hsample )
convergences a certain value, exit the iteration;
• otherwise, in core domain, use (Hcoil + Hsample ) as excitation to compute Hcore , and
repeat the above step.
Core
Coil

Hcoil

Hcore

Hsample

Sample
(b)

(a)

Figure 2.2: Modelling method in case of ferrite core [20]: (a) configuration; (b) illustration
of the solving procedures.

2.6

Parallelization

Typical ECT problems involve multi-position scan. The total number of scanning positions
varies from a few tens (for a typical 1-D scan) to hundreds (for a typical 2-D scan) or even
thousands for a high resolution scan. By applying the Ar , V − Ar formulation and the
technique described in section 2.5, calculations in each scanning position is independent of
other scanning positions. Using this property, the calculation of multiple positions can be
conducted in parallel.

24

Chapter 3
Application to Steam Generator Tube
Inspection

3.1

Introduction

Heat exchanger tubes are used in a variety of industries, including power stations,
petrochemical plants, oil refineries, and air conditioning and refrigeration units for
transferring heat to the fluid circulating outside the tube. An SG is a typical heat exchanger
used in nuclear power plants, as shown in Figure 3.1 [21]. Steam generators (SGs) transfer
heat from the primary loop to the secondary loop feed water circulating on the outside of
the tubes to produce steam that drives the turbines.
It is critical that the radioactive primary coolant does not leak into the secondary side.
The SG tubes are continuously exposed to harsh environmental conditions including high
temperatures, pressures, and material interactions resulting in various types of degradation
mechanisms such as mechanical wear between tube and tube support plates, outer diameter
stress corrosion cracking, pitting, volumetric changes, primary water stress corrosion
cracking, and inter-granular attack.

These flaws can result in tube thinning and/or

development of multiple crack-like flaws, thereby increasing the risk of contaminating the
steam on the secondary side. Consequently the SG tubes in nuclear power plants are required
to be inspected periodically for degradation. EC inspection has proved to be a fast and
25

(b)

(a)

Figure 3.1: Illusion of steam generator [21]: (a) steam generator; (b) steam generator in
nuclear power plant.
effective way to detect and size most degradation mechanisms that occur in SGs.

3.2

Eddy

Current

Probes

Simulation

and

Experimental Validations
Theoretical models simulating EC inspection are very valuable in evaluating performance
of probes in critical situations and also understanding the effect of different variabilities on
the probe signal. However, a theoretical model is useful only if it is validated by comparing
with experimental measurements. The three dimensional FE model was used to simulate
EC SG tube inspection using a number of commercially available probes. The simulation
geometry was chosen to represent the calibration standard samples. The calibrated signals
were compared with experimental signals. Typical results from three types of probes, namely,
the pancake coil probe, plus point probe and array probes are presented as following.

26

3.2.1

Rotating Probe Coil – Pancake Coil Probe

An unshielded 0.115 inch diameter pancake coil probe with ferrite core was modelled and
inspection of a tube with axial notch using the probe was simulated [22]. The geometrical
parameters for the tube and notch are summarized in Table 3.1.

Geometry
Tube outer diameter
Tube inner diameter
Defect type
Defect depth
Defect length
Defect width

Tube with Defect
0.877”
0.775”
Axial notch
63% Outer diameter (OD)
0.501”
0.006”

Table 3.1: Modelling geometry for pancake coil probe

(a)

(b)

(c)

(d)

Figure 3.2: Impedance of pancake coil probe on a raster scan around the defect: (a) real part
of simulation signal; (b) imaginary part of simulation signal; (b) real part of experimental
signal; (d) imaginary part of experimental signal.

The simulations were carried out at 100 kHz frequency and the normalized signal from a
27

two dimensional scan was compared with experimental measurements. Comparison of the
model predicted signals and experimental signals (real and imaginary parts of probe coil
impedance) are shown in Figure 3.2.

3.2.2

Rotating Probe Coil – Plus Point Probe

Inspection of a tube with axial notch using a plus-point probe with ferrite core was simulated.
An illustration of plus-point probe is shown in Figure 3.3, and the geometrical parameters
of the tube and notch are summarized in Table 3.2.

Figure 3.3: Illustration of plus point probe.

Geometry
Tube outer diameter
Tube inner diameter
Defect type
Defect depth
Defect length
Defect width

Tube with Defect
0.877”
0.775”
Axial notch
100% Through wall (OD)
0.500”
0.005”

Table 3.2: Modelling geometry for plus point probe

The simulations were carried out at 300 kHz frequency and signals from a two dimensional

28

(a)

(b)

(c)

(d)

Figure 3.4: Impedance of plus point probe on a raster scan around the defect (2-D plot):
(a) real part of simulation signal; (b) imaginary part of simulation signal; (b) real part of
experimental signal; (d) imaginary part of experimental signal.

(a)

(b)

(c)

(d)

Figure 3.5: Impedance of plus point probe on a raster scan around the defect (3-D plot):
(a) real part of simulation signal; (b) imaginary part of simulation signal; (b) real part of
experimental signal; (d) imaginary part of experimental signal.
scan around the defect were calculated. Comparison of the model predicted signals and
experimental signals (real and imaginary parts of probe coil impedance) are shown in

29

Figure 3.4 (2-D plot) and Figure 3.5 (3-D plot).

3.2.3

Array Coils Probe - Array Probes

Inspection of a tube with axial notch using array probes, as shown in Figure 3.6, was
modelled. The geometrical parameters for the tube and defect are as shown in Table 3.3.

(a)

(b)

Figure 3.6: Illustration of array probes: (a) top view; (b) plan view.

Geometry
Tube outer diameter
Tube inner diameter
Defect type
Defect depth
Defect diameter

Tube with Defect
0.750”
0.665”
Through wall hole
100% Through wall (OD)
0.052”

Table 3.3: Modelling geometry for array probes

The signals from a scan around the defect were modelled at 250 kHz frequency. Comparison
of the model predicted signals (real and imaginary parts of probe coil impedance) with
experiment are shown in Figures 3.7 and 3.8 (axial channel) and Figures 3.9 and 3.10
(circumferential channel).
30

(b)

(a)

Figure 3.7: Real part of channel impedance of array probes: (a) simulation; (b) experimental.

(a)

(b)

Figure 3.8: Imaginary part of axial channel impedance of array probes: (a) simulation; (b)
experimental.

31

(b)

(a)

Figure 3.9: Real part of circumferential channel impedance of array probes: (a) simulation;
(b) experimental.

(a)

(b)

Figure 3.10: Imaginary part of circumferential channel impedance of array probes: (a)
simulation; (b) experimental.

32

3.2.4

Quantitative

Comparison

between

Experimental

and

Simulation Signals
For different probes, we define the quantitative metric as following:
• rotating probe coil (RPC) and X-probe:
for each channel,

MAX|exp − sim|
MAX|exp|

(3.1)

MAX|magnitudeexp − magnitudesim |
MAX|magnitudeexp |

(3.2)

MAX|phaseexp − phasesim |
MAX|phaseexp |

(3.3)

maximum error =

• bobbin probes:

magnitude error =
phase error =

The quantitative comparison of RPC between experimental (exp) and simulation (sim)
signals for a flaw at 300 kHz are shown in Table 3.4.
Flaw dimension

Probe type

Axial notch, 100% TW,
+Point
Length 0.38”, Width 0.005”
Pancake
Axial notch, 57% OD,
+Point
Length 0.38”, Width 0.005”
Pancake
Tube dimension: ID 0.647”, OD 0.745”

Maximum error(%)
Horizontal Vertical
6.31
3.13
3.96
10.50
17.22
16.97
6.07
9.61

Table 3.4: Quantitative comparison of rotating probe coil.

The quantitative comparison of a .610” bobbin probe between experimental (exp) and
simulation (sim) signals for a flaw at 270 kHz are shown in Table 3.5.

33

Flaw dimension Magnitude error(%)
100% TW hole
12.63
60% TW FBH
6.05
Tube dimension: ID 0.668”, OD 0.75”

Phase error(%)
11.67
0.58

Table 3.5: Quantitative comparison of bobbin probe.
The quantitative comparison of X-probe between experimental (exp) and simulation (sim)
signals for a flaw at 250 kHz are shown in Table 3.6.
Flaw dimension

Channel
Axial Vertical
100% TW hole
Circumferential Vertical
Tube dimension: ID 0.664”, OD 0.75”

Maximum error(%)
5.47
13.90

Table 3.6: Quantitative comparison of X-probe.

3.2.5

Conclusion

This section describes the development of a FE model for simulating EC inspection of SG
tubes using a variety of EC probes, flaw and tube geometries. Validation results show that
the model accuracy is within the bounds specified by industry and hence can be used to
generate signals from a variety of degradation mechanisms. In fact, it is important to note
that error reported is largely due to experimental measurement noise. The model can be
used to study effect of different operational parameters on the probe signal which in turn
can be used to develop PoD models.

3.3

Simulation Software Features

In addition to the FE solver, a software package with a graphical user interface that allows
easy use of the model for generating signals from SG inspection has been developed [23]. The

34

user can select a probe, defect size and defect location using a pull down menu. The basic
features of the graphic user interface of this software package are described in section 7.2.2.3.

3.4

New Probe Design and Simulation

SG tube inspection using EC techniques has evolved over the years from a simple bobbin coil,
to RPC and array probes, in an attempt to improve the speed and reliability of inspection.
As the probe design has evolved, the new sensors have provided more data, resulting in higher
detection accuracy, resolution, sensitivity, speed and robustness. The bobbin probe serves
high performance speed for full tube inspection but is insensitive to circumferentially oriented
crack, and the RPC probe offers the best spatial resolution but involves complex mechanical
control and slow inspection speed whereas the array probe is a compromise between speed
and resolution but need complicated hardware implementation.
This section presents a design for a rotating field eddy current (RoFEC) probe that offers the
sensitivity of an RPC probe without need for mechanical rotation of the probe assembly and
speed of the bobbin probe. The feasibility of the design is demonstrated using simulation
model. [24]

3.4.1

The Rotating Field Eddy Current Probe Operational
Principles and Design

The rotating field principle is a fundamental property of all electrical machines and consists of
three identical coils located on axes physically 120 degrees apart and supplied from a balanced
three-phase supply (i.e. with difference between the phases of exactly 120 degrees). In the
volume close to the multiphase probe the instantaneous magnetic flux densities, resulting
35

from each of the coils is illustrated in Figure 3.11.

i

ia

ib

Let ia (t), ib (t) and ic (t) represent the

Bb

ic
t

Ba

Bc
(b)

Bv

M

α
Bh
(c)

(a)

Figure 3.11: Principle of rotating field : (a) three-phase currents; (b) flux due to currents in
three-phase winding; (c) horizontal and vertical components of the resultant magnetic flux.

excitation currents in the three windings of the probe as expressed by equation (3.4).

ia (t) = I sin(ωt)
2π
)
3
4π
ic (t) = I sin(ωt +
)
3
ib (t) = I sin(ωt +

(3.4)

where, I =maximum value of the current due to one phase alone.
The magnetic flux density associated with the three coils can be represented by
equation (3.5), and are perpendicular to the plane of each coil. On a 2-D cross-sectional
plane of the tube and probe coils, the horizontal and vertical components of the resultant

36

magnetic flux density can be expressed by equation equation (3.5). As a result, the total
magnetic flux density has a constant amplitude and its phase will change in the same manner
as the excitation frequency as expressed in equation equation (3.6).

Ba (t) = φ sin(ωt)
2π
)
3
4π
Bc (t) = φ sin(ωt +
)
3
Bb (t) = φ sin(ωt +

(3.5)

where, φ =maximum value of the magnetic flux due to one phase alone.

3
Bh = Ba − (Bb + Bc ) cos(60 deg) = φ sin(ωt)
2
3
Bv = Bb cos(30 deg) − Bc cos(30 deg) = φ cos(ωt)
2
3
M 2 = Bh2 + Bv2 = ( φ)2
2

(3.6)

A 2-pole probe operating at 60 Hz generates a magnetic field that scans the pipe
circumference at 3600 rpm. Higher excitation frequencies result in proportionately higher
circumferential scan rates. The field generated by the probe is largely radial that result
in induced currents that flow circularly around the radial axis rotating around the tube at
synchronous speed effectively producing induced ECs that are multi-directional; in fact the
induced currents are identical to that produced by rotating pancake coils. The probe will
consequently be sensitive to cracks of all orientations in the tube wall. The response signal
due to the field/flaw interaction can be picked up using a simple bobbin coil for detecting a
discontinuity in tube wall ans shown in Figure 3.12.

37

(a)

(b)

Figure 3.12: Rotating field eddy current probe structure: (a) top view; (b) plan view.

3.4.2

Simulation Results

A typical free span region of SG tube geometry was modelled. The outer and inner diameters
of the tube are 22.3 mm and 19.7 mm respectively. A defect of width 3.5 mm, length 4 mm
and depth of 60% tube wall thickness was introduced on the inside of the tube wall as shown
in Figure 3.13a. The probe was excited by using a 320 kHz source and the voltage induced
by the rotating fields in a bobbin coil is calculated as the probe moves axially within the tube
as shown in Figure 3.13b. To validate the performance of the RoFEC probe, we also present
the results of a normal differential bobbin coils excited at 320 kHz as shown in Figure 3.13c.

The probe was excited by using a 320 kHz source and the z-component of the magnetic flux
density induced by the rotating fields in an array of giant magneto-resistance (GMR) sensors
is also calculated as shown in Figure 3.14.
Parametric studies were conducted for axial notches of different depth and locations.
Simulation results for ID and OD defects of various are presented in Figure 3.15a. The
correlation between defect geometry and the signal phase and amplitude is evident from

38

(a)

(b)

(c)

Figure 3.13: Signal of rotating field eddy current probe with bobbin pickup coil : (a) geometry;
(b) rotating field eddy current probe with absolute pickup coil; (c) differential bobbin coil.
-10

-10

0

0

-10

-10
50

150

250

350

50

150

(a)

250

350

(b)

-10

0

-10
50

150

250

350

(c)

Figure 3.14: Signal of rotating field eddy current probe with GMR array sensors: (a)
magnitude; (b) horizontal channel; (c) vertical channel.

39

these simulated signals. In order to test the robustness of performance of the probe, defect
signals from centred and off-centred probes were calculated. These results for a 60% ID
notch are presented in Figure 3.15b. It is seen that the phases of the three signals are the
same. Parametric variations of the signals with respect to increasing frequencies were also
conducted. The signals from a 60% ID defect at frequencies of 100, 320 and 600 kHz shown
in Figure 3.15c.
80%
60%
40%
20%

ID
ID
ID
ID

20% OD
40% OD
60% OD
80% OD
100% TW

(a)

+1 mm
−1 mm

100kHz
320kHz
600kHz

centered

(b)

(c)

Figure 3.15: Signals of parametric studies: (a) 20%, 40%, 60%, 80%, 100%, I.D. and O.D.
defects; (b) 60% I.D. defect with +1, 0, -1 mm offsets; (c) 60% I.D. defect at 100, 320, 600
kHz operating frequencies.

3.4.3

Conclusion

This section presents a use of the computational model in the design of a novel EC probe
based on rotating magnetic fields for the inspection of SG tubes in nuclear power plants. The
simulation results show the advantages of the RoFEC probe, which are high performance
speed and high resolution for full tube inspection. The probe design is based on basic
principles of electrical machines and uses a simple excitation system. The performance is
40

relatively robust with respect to probe wobble about the axis.

3.5

Conclusion

The FE modelling provides an inexpensive and fast method to simulate realistic field
conditions and defect geometries without having to produce flawed tube specimens. The
modelling signals can assist in signal interpretation, serve as training tools, generate training
signals, assist in PoD calculations with applying noise models, qualify probe technique,
demonstrate probe performance demonstration, aid in probe design. The experimentally
validated method is widely in use in the SG inspection industry by utility engineers,
inspection vendors, probe developers and NDE instructors.

41

Part II
Direct Integral Solver and Fast
Algorithms

42

Chapter 4
Moment Method

4.1

Introduction

Integral equation method, used extensively for solving PDEs, is often referred to as method
of moments (MoM) in EM.
The integral equation method is a numerical computational method of solving linear PDEs
which have been formulated as integral equations (i.e. in boundary integral form). It has
applications in many areas of engineering and science including fluid mechanics, acoustics,
EM and fracture mechanics.

4.1.1

Volume Integral Equation Method

Generally, there are two categories of volume integral equation (VIE) methods for EC
problems: 1. free space Green’s function based VIE method for bounded geometries; 2.
numerical Green’s function based VIE method for unbounded geometries.

4.1.1.1

Bounded Eddy Current Problems

For bounded geometries, due to the need for discretizing the whole 3-D geometry, the
computation cost is prohibitive. A typical problem is a conducting body excited by a line
source [25,26] as shown in Figure 4.1 Similar to the reduced potential formulation for FEM,
the vector potential A can be split into three parts: As due to excitation current, Ae due
43

Figure 4.1: A conducting body excited by a line conductor carrying current [25, 26].
to induced currents and Am due to induced magnetization. From Amp`ere’s and Faraday’s
laws of time varying quasi-static EM fields, we have [27, 28]:

A(r) = As (r) + Ae (r) + Am (r)
˚
µ0
Js (r )
As (r) =
dΩ
4π
Ωs |r − r |
˚
−µ0 σ
jω A(r ) + ∇ V (r )
dΩ
Ae (r) =
4π
|r − r |
Ωe
˚
µ0
1
Am (r) =
M (r ) × ∇
dΩ
4π
|r − r |
Ωm

(4.1)
(4.2)
(4.3)
(4.4)

where, V is electric scalar potential. Applying the Coulomb gauge (∇ · A = 0), we obtain:

−µ0 σ
4π

˚
[jω A(r ) + ∇ V (r )] · ∇
Ωe

44

1
dΩ = 0
|r − r |

(4.5)

Combining equation (4.1) and equation (4.5) and using Galerkin method, we form a complete
system as following:


  

A 0 0 D Ax  As,x 

  


  

 0 B 0 E  Ay  As,y 

  


  = 


  

 0 0 C F  Az  As,z 

  


  

G H I J
V
0

(4.6)

From equations (4.2) to (4.4), it is obvious that volume of the excitation current region Ωs ,
the induced current region Ωe and the induced magnetization region Ωm are required to be
finite. If the geometry of interest is unbounded, in other words the volume of the geometry of
interest is infinite, we will not be able to discretize such geometry. In this case, this approach
is not applicable.
Similar to the fast solvers in microwave problems, several fast schemes for bounded VIE
problems based on adaptive integral method, fast multi-pole method etc. have also been
reported [29–41].

4.1.1.2

Unbounded Eddy Current Problems

For unbounded geometries, it is not possible to discretize the infinite volume, therefore,
new Green’s functions satisfying the boundary conditions of different cases, such as, half
space conductor [42–45], layered slab conductor [46–53], wedge conductor [54,55] and layered
rod/bolt-hole conductor [56–59] need to be derived. This is one of the disadvantages of this
approach.

45

We can represent Faraday’s law and Amp`ere’s law in the form as:

∇ × E = −jωµor H − jω(µ − µor )H = −jωµor H − Mδ

(4.7)

∇ × H = jω or E + [jω( − or ) + (σ − σor )]E + Js = jω or E + Jδ + Js

(4.8)

where the subscript or is short for original.
↔

Then, assuming we have the suitable dyadic Green’s functions G that satisfies the problem
boundary conditions, we can obtain the fields due to the source current Js and equivalent
electric current Jδ and magnetic current Mδ in the following form:
˚

↔

E inc

= −jωµor
Js · GdΩ
Ωs
˚
↔
inc
H
=
(∇ × Js ) · GdΩ
Ωs
˚
↔
E1 (Jδ ) = −jωµor
Jδ · GdΩ

Ωδ
Ωδ

(4.11)

↔

(∇ × Jδ ) · GdΩ

˚
E2 (Mδ ) =

(4.10)

Ωδ

˚
H1 (Jδ ) =

(4.9)

(4.12)

↔

(∇ × Mδ ) · GdΩ
˚

H2 (Mδ ) = −jωµor

Ωδ

(4.13)

↔

Mδ · GdΩ

(4.14)

Finally combining the equations (4.9) to (4.14) in the form of:

E = E inc + E1 (Jδ ) + E2 (Mδ )

(4.15)

H = H inc + H1 (Jδ ) + H2 (Mδ )

(4.16)

46

and using Galerkin method, we obtain a complete system of equations.
The Green’s functions for slab conductor [47] and borehole [57] EC problems will be shown
as example in the following sections.

Ωδ

σ=0
σ = σ0
σ=0

(a)

σ=0
σ = σ0
σ=0

Ωδ

(b)

Figure 4.2: Unbounded geometry examples: (a) slab conductor; (b) borehole geometry.

Dyadic Green’s functions for slab conductor
For slab conductor (−c < z < 0) geometry as shown in Figure 4.2a, the relation between
electric dyadic Green’s function and magnetic vector potential dyadic Green’s function is:

↔

jω GE (r; r ) =

↔
1
∇ × ∇ × GA (r; r )
µ0 σ0

47

(4.17)

The magnetic vector potential dyadic Green’s function fitting the slab conductor geometry
is:

↔

↔

↔

↔

↔

↔

GA (r; r ) = I φ(r; r ) + I φ(r; r1 ) + I φ(r; r2 ) + I φ(r; r3 ) + I φ(r; r4 )
↔

↔

+∇ × zˆ∇ × zˆV1 (r; r ) + ∇ × zˆ∇ × zˆV2 (r; r ) + I W1 (r; r ) + I W2 (r; r )

(4.18)

↔

↔

where, r1 = r −2z zˆ, r2 = r −2(z +c)ˆ
z , r3 = r −2cˆ
z , r4 = r +2cˆ
z , I = xˆxˆ + yˆyˆ+ zˆzˆ, I =
xˆxˆ + yˆyˆ − zˆzˆ , and

↔

φ(r; r ) =

I−

∇∇
k02

e−jk|r−r |
4π|r − r |

(4.19)

The remaining terms are given by:

1
V1 (r; r ) = 2
4π

ˆ∞ ˆ∞ (e)
Γ − Γ(m) γ(z−z −2c)
e
+ e−γ(z−z +2c) eiu(x−x )+iv(y−y ) dudv
2
2κ γ

−∞ −∞

(4.20)
1
V2 (r; r ) = 2
4π

ˆ∞ ˆ∞
−∞ −∞

Λ(e) − Λ(m) γ(z+z )
e
+ e−γ(z+z +2c) eiu(x−x )+iv(y−y ) dudv
2
2κ γ
(4.21)

ˆ∞ (m)
Γ
eγ(z−z −2c) + e−γ(z−z +2c) J0 (κρ)κdκ
2γ

1
W1 (r; r ) =
2π

−∞
ˆ∞

−1
W2 (r; r ) =
2π

Λ(m) γ(z−z −2c)
e
+ e−γ(z+z +2c) J0 (κρ)κdκ
2γ

−∞

48

(4.22)

(4.23)

where,

Γ(e) (κ) =

(γ − κ)2
−1
(γ + κ)2 + (γ − κ)2 e−2γc

1
−1
1 − e−2γc
γ 2 − κ2
Λ(e) (κ) =
+1
(γ + κ)2 − (γ − κ)2 e−2γc
Γ(m) (κ) =

1
+1
Λ(m) (κ) = −2γc
e
−1

(4.24)
(4.25)
(4.26)
(4.27)

and κ2 = u2 + v 2 , γ 2 = κ2 − k 2 , ρ2 = (x − x )2 + (y − y )2 .

Dyadic Green’s functions for borehole geometry
For the borehole geometry as shown in Figure 4.2b, by using the transverse electric potential
Wa and the transverse magnetic potential Wb , for piecewise homogeneous regions of a
structure in cylindrical coordinate system, the magnetic vector potential can be derived as:

A = ∇ × W = ∇ × (ˆ
z Wa + zˆ × ∇Wb )

(4.28)

Therefore, the electric field can be represented as:

E = ∇ × (ˆ
z Wa ) + ∇ × ∇ × (ˆ
z Wb ) −

1
P
σ0

(4.29)

where σ0 is the host conductivity, P = (σ − σ0 )E is an electric dipole density which can
be consider as the equivalent source of the flaw field, and σ is the flaw conductivity. By

49

substituting equation (4.29) into the electric field equation:

∇ × ∇ × E − k 2 E = iωµP

(4.30)

where, k 2 = iωµσ0 we can show that the potentials satisfy:

  
zˆ · ∇ × P

Wa  
(∇2 + k 2 )   = 

1
Wb
2 zˆ · ∇ × ∇ × P

(4.31)

k

A set of scalar Green’s functions can be used as a solution of equation (4.31), which satisfies:




1 0
Gaa Gab 
(∇2 + k 2 )∇2t 
=
−

 δ(r − r )

0 1
Gba Gbb
2

(4.32)

2

where, ∇2t = ∂ 2 + ∂ 2 , Gij = −∇2t Uij , i, j = a, b.
∂x
∂y
By representing the electric field by using the dyadic Green’s function, we obtain,
˚
E(r) =

E inc (r) + iωµ0

↔

P (r) · G(r; r )dV

(4.33)

V

where the dyadic Green’s function is:

↔
1
G =(∇ × zˆ)∇ × (ˆ
z Uaa ) + 2 (∇ × zˆ)∇ × ∇ × (ˆ
z Uab )
k
1
+ (∇ × ∇ × zˆ)∇ × (ˆ
z Uba ) + 2 (∇ × ∇ × zˆ)∇ × ∇ × (ˆ
z Ubb )
k

50

(4.34)

Represent the scalar kernels in Fourier transform as:
ˆ∞
∞
1
˜ m (ρ; ρ , κ)eiκ(z−z ) dκ
G(r; r ) = 2
eim(φ−φ )
G
4π m=−∞
m=−∞

(4.35)

˜ m (ρ; ρ , κ) = −γ 2 U˜m (ρ; ρ , κ), γ 2 = κ2 +k 2 . The Fourier transformed kernels for an
where, G
equivalent source in the borehole conductive region can be represented as the summation of
an unbounded domain kernel and a term expressing outward propagation from the interface

˜ ij (ρ; ρ , κ) = G
˜ (0) (ρ; ρ , κ) + G
˜ (Γ) (ρ; ρ , κ)
G
ij

(4.36)

The unbounded domain term is

˜ (0) (ρ; ρ , κ) = Im (γρ< )Km (γρ> )
G

(4.37)

where, Im and Km are associated Bessel functions. The field propagating from the interface
is represented by

(Γ)
U˜m,ij (ρ; ρ , κ) = Γm,ij Im (γρ)Km (γρ )

51

(4.38)

where, Γm,ij are the boundary reflection coefficients as following:

Γm,aa =

Im (γb) −m2 µr κ2 + M (γb)[µr κ2 Λ(γb) − γ 2 Λ(|κ|b)]
·
Km (γb) −m2 µr κ2 + M (γb)[µr κ2 M (γb) − γ 2 Λ(|κ|b)]

(4.39)

Γm,ab =

−m2 µr κk 2
1
·
[Km (γb)]2 −m2 µr κ2 + M (γb)[µr κ2 M (γb) − γ 2 Λ(|κ|b)]

(4.40)

Γm,ba =

1
−m2 µr κ
·
[Km (γb)]2 −m2 µr κ2 + M (γb)[µr κ2 M (γb) − γ 2 Λ(|κ|b)]

(4.41)

Γm,bb =

Im (γb) −m2 µr κ2 + Λ(γb)[µr κ2 M (γb) − γ 2 Λ(|κ|b)]
·
Km (γb) −m2 µr κ2 + M (γb)[µr κ2 M (γb) − γ 2 Λ(|κ|b)]

(4.42)

xI (x)

xK (x)

where, Λ(x) = I m(x) , M (x) = K m(x) . By using a reciprocity theorem [60], with the
m
m
dipole density calculated, the impedance change of the excitation coil due to the flaw can be
determined from:
˚
I 2 ∆Z

P (r)E inc (r)dV

=−

(4.43)

V

From the two examples discussed above, we can clearly see that the Green’s functions need
to be derived for each geometrical case. Furthermore, as shown in the equations (4.20)
to (4.23) for slab geometry and equation (4.35) for borehole geometry, the Green’s functions
require numerical evaluation of infinite integrals, which is always very challenging. For more
irregular geometries, it’s even difficult to find the Green’s functions in either analytical or
numerical form. In such cases, this approach is not applicable.

52

4.1.2

Surface Integral Equation Method

The surface integral equation (SIE) method, also called boundary element method (BEM),
has been used in EC problems [61–78].
Snc

n
ˆ

Vc
Vn

Figure 4.3: An arbitrary geometry with closed surface.

A typical BEM problem in a bounded arbitrary geometry is shown in Figure 4.3, where, n
ˆ
is the outer normal direction, Vn is free space, Vc is the domain carrying EC and Snc is the
boundary between Vn and Vc .
Deriving from Faraday’s law and Amp`ere’s law, the Helmholtz equations for electric field
and magnetic field are:

∇ × ∇ × E − k 2 E = −jωµJ inc

(4.44)

∇ × ∇ × H − k 2 H = ∇ × J inc

(4.45)

where the square of the wave number in EC region and free space are k 2 = kc2 = ω 2 µ −
jωµσ ≈ −jωµσ and k 2 = kn2 = ω 2 µ0 0 , respectively. Applying the Huygens’ equivalent

53

principle on Helmholtz equations [79], we obtain:
"

↔

E=−

↔

[(ˆ
n × E) · (∇ × G) − jωµ(ˆ
n × H) · G]dS
"Snc

↔

(4.46)
↔

H=−

[(ˆ
n × H) · (∇ × G) − (−jω + σ)(ˆ
n × E) · G]dS
Snc
"
↔
↔
inc
E=E −
[(ˆ
n × E) · (∇ × G) − jωµ0 (ˆ
n × H) · G]dS
Snc
"
↔
↔
H=−
[(ˆ
n × H) · (∇ × G) + jω 0 (ˆ
n × E) · G]dS

(4.47)
(4.48)
(4.49)

Snc

where, equations (4.46) and (4.48) are called electric field integral equation (EFIE) for EC
region and free space respectively, equations (4.47) and (4.49) are called magnetic field
integral equation (MFIE) for EC region and free space respectively, with Green’s function
↔
↔
G = I − ∇∇
2

k

e−jkR
4πR

↔

↔

−jkR

in EC region and Green’s function G = I e 4πR in free space.

The incident electric and magnetic fields are represented as:
˚
E inc

↔

= −jωµ0
J inc · GdV
inc
V
˚
↔
inc
(∇ × J inc ) · GdV
H
=
V inc

(4.50)
(4.51)

The detailed discussion about adapting EFIE and MFIE to fit the EC problem and
corresponding numerical issues are discussed in the following sections.

4.2

Surface Integral Equations

As generally discussed in section 4.1.2, the SIE method, has been developed and used in EC
problems. In this section, we outline the governing SIEs used to formulate EC problems. The
derivation of such integral equations starting with Maxwell’s equations will not be presented

54

here. One can reference any textbook on advanced EM theory [79–86].

4.2.1

Electric Field Integral Equation

From Faraday’s law and Amp`ere’s law, the Helmholtz equation for electric field is:

∇ × ∇ × E − k 2 E = −jωµJ inc

(4.52)

with k 2 = ω 2 ˜µ, where ω is angular frequency, µ is the permeability and ˜ is the equivalent
σ in conductive region.
permittivity with ˜ = in non-conductive region and ˜ = − jω

By introducing the three-dimensional Green’s function G that satisfies the Helmholtz
equation,

∇2 G(r; r ) + k 2 G(r; r ) = −δ(r − r )

(4.53)

and the Sommerfeld radiation condition, we have a well-known solution:

G(r; r ) =

e−jk|r−r |
4π|r − r |

(4.54)

Next, by applying the scalar-vector Green’s theorem [87],
˚
[b∇ × ∇ × a + a∇2 b + (∇ · a)∇b]dV
V
"
=
[(ˆ
n · a)∇b + (ˆ
n × a) × ∇b + (ˆ
n × ∇ × a)b]dS
S

55

(4.55)

letting a = E and b = G, we obtain,
"
E inc (r)−

[ˆ
n · E(r )]∇G(r; r ) + [ˆ
n × E(r )] × ∇G(r; r ) + jωµ[ˆ
n × H(r )]G(r; r ) dS
S

=




E(r) r ∈ V



0

(4.56)

r∈
/V

where,
˚
E inc (r)

= −jωµ0

V inc

1
J inc (r )G(r; r ) + 2 ∇ · J inc (r ) ∇G(r; r ) dV
k

(4.57)

To write equation (4.56) in a compact form, we define two operators
"
L X(r) = −jk
"
K X(r) =

S

1
X(r )G(r; r ) + 2 ∇ · X(r ) ∇G(r; r ) dS
k

X(r ) × ∇G(r; r ) dS

(4.58)
(4.59)

S

and the equivalent surface currents, Js = n
ˆ × H and Ms = E × n
ˆ , we can represent
equation (4.56) as

E inc − ZL(Js ) + K(Ms ) =





E(r) r ∈ V



0

where Z is the equivalent surface impedance defined as Z =
surface EFIE.

56

(4.60)

r∈
/V

µ
˜.

Equation (4.60) is called

4.2.2

Magnetic Field Integral Equation

From Faraday’s law and Amp`ere’s law, the Helmholtz equation for magnetic field is:

∇ × ∇ × H − k 2 H = ∇ × J inc

(4.61)

with k 2 defined as same as in equation (4.52). Using the Green’s function defined in
equation (4.54) and applying the scalar-vector Green’s theorem defined in equation (4.55)
and letting a = H and b = G, we obtain,
"
H inc (r)−

[ˆ
n · H(r )]∇G(r; r ) + [ˆ
n × H(r )] × ∇G(r; r ) − jω [ˆ
n × E(r )]G(r; r ) dS
S

=




H(r) r ∈ V



0

(4.62)

r∈
/V

where,
˚
H inc (r)

=−

V inc

J inc (r ) × ∇G(r; r ) dV

(4.63)

Using the same operators defined in equations (4.58) and (4.59), and the equivalent surface
currents, Js = n
ˆ × H and Ms = E × n
ˆ , we can represent equation (4.62) as

H inc −





H(r) r ∈ V

1
L(Ms ) − K(Js ) =

Z


0

with the equivalent surface impedance Z =

µ
˜.

57

(4.64)

r∈
/V

Equation (4.64) is called surface MFIE.

4.3

High-order Basis Function

In general, there are two categories of basis function for integral equation method. The
first category is defined over the entire solution domain like the plane wave basis set etc.
The second category is defined over piecewise sub-domain of the solution domain. For
most EM problems, it’s very difficult to find a set of basis functions which can completely
approximate the EM quantities over the entire domain. Hence the second category of
piecewise approximation basis functions are used.

4.3.1

RWG Basis Function

As in FEM [3], the locally defined sub-domain basis function has been used very commonly.
In particular, the triangular rooftop function, also known as Rao-Wilton-Glisson (RWG)
basis function [88], is one of the most popular basis functions that is used extensively in EM
problems.
The RWG basis function is defined over two joined triangles at their common edge l as

Λ(r) =





 l + ρ+

r ∈ T+




 l − ρ−

r ∈ T−

2A
2A

(4.65)

where T ± are the joined triangles, A± denote the areas of those triangles, l is the length of
the common edge, and ρ± denote the vectors defined as shown in Figure 4.4a. The vector
plot of RWG basis function is shown in Figure 4.4b.
The most important feature of this basis function is that the normal component to the edge
is a constant and zero at other edges; therefore, it is called zeroth-order RWG basis function.

58

✯
✟
✟
✟
✟
✟
ρ−✟

T+

l
T−

ρ+
✟

✯
✟
✟
✟
✟
✟

(a)

(b)

Figure 4.4: Illusion of RWG basis function: (a) two joined triangles; (b) vector plot of the
zeroth-order RWG function.

4.3.2

High-order RWG Basis Function

As in FE analysis [3], the MoM using the RWG and the corresponding higher-order basis
functions needs a well-connected mesh. Recently, an interesting set of higher-order RWG
basis functions has been developed for both FEM and MoM. These basis functions are
defined over each patch using the Lagrange interpolation polynomials [89]. The Lagrange
interpolation points are chosen to be the same as the nodes of the well-developed Gaussian
quadrature points in triangle. As a result, evaluation of the integrals in the MoM is greatly
simplified. Numerical results exhibited higher-order convergence. As an example, the vector
plots of first-order RWG basis functions are shown in Figure 4.5.

(a)

(b)

(c)

Figure 4.5: Vector plot of first-order RWG basis function: (a) basis function associated with
the common edge; (b) the other basis function associated with the common edge; (c) basis
function associated with the triangle.

59

4.4

Moment Method Procedure

To solve the surface integral equations (4.60) and (4.64), we firstly discretize the entire
surface S into small triangular patches. Then we apply basis functions to expand the surface
current density as:
N

Js =

In Λn (r)

(4.66)

Kn Λn (r)

(4.67)

n=1
N

Ms =
n=1

where N is the number of unknowns.
Since the MoM procedures for both EFIE and MFIE are similar, the MFIE is taken as
an example here. For a homogeneous body, using the same basis functions as the testing
functions, and applying the Galerkin method to the MFIE as defined in equation (4.60) on
both sides of the surface, we obtain the matrix equations:
N

N

Amn In +
n=1
N

(4.68)

Dmn Kn = gm

(4.69)

n=1
N

Cmn In +
n=1

Bmn Kn = fm

n=1

60

in which, m = 1, 2, · · · , N , and
"
Amn = Z0
Λm · Λn + n
ˆ × K(Λn ) dS
S
"
Bmn =
Λm · n
ˆ × L(Λn ) dS
S
"
Cmn = Z
Λm · Λn + n
ˆ × Ki (Λn ) dS
S
"
Dmn =
Λm · n
ˆ × Li (Λn ) dS
S
"
fm =
Λm · n
ˆ × H inc dS

(4.70)
(4.71)
(4.72)
(4.73)
(4.74)

S

gm = 0

(4.75)

Thus, we can form a complete linear system for In and Kn that should be solved.

4.5

Singularity of Green’s Function

Both operators defined in equations (4.58) and (4.59) include Green’s function which is
defined in equation (4.54):

G=

e−jkR
e−jk|r−r |
=
4π|r − r |
4πR

(4.54 revisited)

In order to evaluate the matrices defined in equations (4.70) to (4.73) numerically, the
diagonal terms, when m = n, need to be evaluated at the singularity points when R =
0. Moreover, due to the behaviour of the Green’s function, when R → 0, the numerical
quadrature methods such as Gaussian quadrature etc. will suffer from low accuracy in these
evaluations [90–104].
One popular method to solve this problem is Duffy’s method [103]. Basically, Duffy’s method
61

consists of choosing two sets of non-overlapping Gauss-Legendre quadrature points in the
triangle to avoid numerical division by zero. However, as mentioned above, the accuracy
is always very bad as R → 0. Other methods involving transformation of the integrals to
cylinder coordinates [99] have not shown significant improvement either.
Another method has been presented in [104] to fully solve this problem on a plane triangle.
Expanding the Green’s function in Maclaurin series,

G=

1 1
(−jk)2
(−jk)3 2
e−jkR
=
+ (−jk) +
R+
R + ···
4πR
4π R
2
6

(4.76)

1 is the reason which causes the accuracy issue.
we can clearly see that the first term R
2
1 and third term (−jk) R lead to additional accuracy problems when
Additionally, the first R
2

evaluating the gradient of Green’s function. Thus, we need to remove the singular points
and evaluate the integral by reconstructing the Green’s function as

G=

G−

1
k2
1
k2
+
R +
−
R
4πR 8π
4πR 8π

(4.77)

1
k2
1 1
k2
+
R +
∇ −
∇R
4πR 8π
4π R 8π

(4.78)

As a consequence,

∇G = ∇ G −

The first term in equations (4.77) and (4.78) is very smooth and can be evaluated accurately
by using Gauss-Legendre quadrature. The second and third terms in those equations can be
analytically evaluated [104]. This approach has shown considerable promise.

62

4.6

Low Frequency Breakdown of SIE

Firstly, we introduce the regimes of EM physics:
• low frequency, kd

1;

• middle frequency, kd ∼ O(1);
• high Frequency, kd

1.

where, k = 2π
λ is wave-number and d is typical size of the domain of interest.
The EFIE suffers from a low-frequency breakdown [105–112] for perfect electrical conductors
(PEC) problems, and both EFIE and MFIE suffer from low-frequency breakdown for
penetrable body problems. This problem, however, is not specific to integral equations.
It is also germane to differential equation solutions of EM. In other words, a numerical EM
solver, formulated for mid-frequency, may not work well for low frequency.
We can understand the low-frequency breakdown problem by looking at the Maxwell’s
equations when ω → 0 [105, 112]. At zero frequency, they become

∇×E =0

(4.79)

∇×H =J

(4.80)
∇·J
ω→0 jω

∇ · D = ∇ · E = ρ = lim
∇ · B = ∇ · µH

(4.81)
(4.82)

From equation (4.80), the current density J that produces the magnetic field must be
divergence free. However, the current density J that produces a surface charge density
ρ in equation (4.81) cannot be divergence free, but in order to keep the right-hand side of
63

equation (4.81) bounded, the J ∼ O(ω), when ω → 0.
To solve this problem, we can decompose the current density into a solenoidal (divergence
free) part and a irrotational (curl free) part. Consequently, the respective basis functions
need to be decomposed into those parts, according to loop-star and loop-tree decomposition
[34–39,105,112]. However, this kind of decomposition is not robust and may not be resolved
properly for certain discretization.
As a solution [67–70, 73, 74], especially for EC problems, we have,
inside the conducting region,
"
Ei =

S

[Ms × ∇Gi + jωµi Js Gi ]dS

"
Hi =

S

[σ Ms Gi − Js × ∇Gi +

(4.83)

∇ · Ms
∇Gi ]dS
jωµi

(4.84)

and in free space,
"
Ho =

H inc

∇ · Ms
∇Go ]dS
jωµo

(4.85)

[Ms × ∇Gi + jωµi Js Gi ]dS

(4.86)

[Js × ∇Go −

+
S

Therefore, we obtain,
"
Ms = −ˆ
n×
∇ · Ms
=n
ˆ·
jωµi
Js

−n
ˆ × H inc

S

"
S

[σ Ms Gi − Js × ∇Gi +

"

=n
ˆ×

[Js × ∇Go −
S

∇ · Ms
∇Gi ]dS
jωµi

∇ · Ms
∇Go ]dS
jωµo

(4.87)
(4.88)

With proper choice of basis functions, we can easily discretize any kind of geometry and
form a linear system of equations that should be solved.

64

4.7

Validation

In this section, the validation of EC problems by applying the formulations discussed in
section 4.6 has been presented and discussed.
The following two problems have been considered: i) an air core pancake coil is scanned,
along the direction of the length of a rectangular notch in an aluminium alloy plate; ii) an
air core RPC is scanned, along the axial direction inside an inconel 600 SG tube with a
longitudinal through-wall notch.
In the first case, the coil lift-off is fixed and the changes in coil impedance, ∆Z, is measured
as a function of coil-center position. The parameters for this test experiment are listed in
Table 4.1. This is a benchmark problem reported in [113–117]. The comparison between the
experiment and simulation signals is shown in Figure 4.7.

Figure 4.6: Diagram for the measurement of ∆Z due to a surface breaking slot [115].

65

Coil
Inner radius (a2)
9.34 ± 0.05 mm
Outer radius (a1)
18.4 ± 0.05 mm
Length (b)
9.00 ± 0.20 mm
Number of turns (N)
408
Lift-off (l)
2.03 ± 0.05 mm
Test specimen
Conductivity (s)
3.06 ± 0.02 × 107 S/m
Thickness
12.22 ± 0.02 mm
Defect
Length (2c)
12.60 ± 0.02 mm
Depth (h)
5.00 ± 0.05 mm
Width (w)
0.28 ± 0.01 mm
Other parameters
Frequency
7000 Hz
Skin depth at 7000 Hz
1.09 mm
Isolated coil inductance 3.96 ± 0.10 mH
Table 4.1: Parameters of test experiment No. 1.

Figure 4.7: Variation of coil impedance change with scanning location above an aluminium
alloy plate.
In the second case, the coil lift-off is fixed and the changes in coil impedance, ∆Z, is measured
as a function of coil-center position. The parameters for this test experiment are listed in
Table 4.2. This is a World Federation of NDE Centers (WFNDEC) 2012 EC benchmark
66

problem [118]. The comparison between the experiment and simulation signals is shown in
Figure 4.9.

Figure 4.8: Geometry of a rotating probe coil inside an inconel steam generator tube with an
axial through-wall notch [118].

Coil
Inner radius (r2)
Outer radius (r1)
Height (x2-x1)
Number of turns (N)
Lift-off (λ)

1.529
3.918
1.044
305
1.235

± 0.004 mm
± 0.003 mm
± 0.005 mm
mm

Tube
16.64 ± 0.025 mm
18.99 ± 0.025 mm
Inconel 600

Inner diameter (a)
Outer diameter (b)
Material
Defect
Length
Depth
Width

12.20 mm
through-wall
85 µm
Other parameters
Frequency
25, 50, 100 kHz
Isolated DC coil inductance 465 µH
Isolated coil resistance
19.00 Ω
Table 4.2: Parameters of test experiment No. 2.

67

(a)

(b)

Figure 4.9: Variation of coil impedance change with scanning location inside an inconel steam
generator tube at: (a) 25 kHz; (b) 50 kHz; (c) 100 kHz.

68

Figure 4.9 (Cont’d)

(c)

As shown in Figures 4.7 and 4.9, the simulation signals are qualitatively matched with the
experiment signals. The mismatch between the simulation and experimental signals is mainly
caused by two reasons: mesh quality and Galerkin testing. For mesh quality, if a curvilinear
surface has been discretized by using planner triangular patches as we present in the case
two of a tubing geometry, inaccuracy of simulation signals is higher. For instance, simulation
signals in the first case of slab geometry compares more closely with experiment signals than
that in the second case, since the mesh of a plate geometry is more regular than that of a
tube geometry. For Galerkin testing, proper choice of the testing procedure, as discussed in
section 7.2.2.2, could considerably improve the accuracy of modelling geometries with thin
notches. This is an issue discussed in future work.

69

4.8

Conclusion

Using the higher-order RWG basis function, we use equations (4.86) to (4.88) to solve the
EC problem using the MoM as discussed in section 4.4. Due to the localized nature of
the EC fields, the near-field EM interaction plays a major role. Applying the singularity
removal techniques discussed in section 4.5, we have computed the EC field distributions
with adequate accuracy.

70

Chapter 5
Fast Algorithms for Integral Equation

5.1
5.1.1

Matrix Compression Algorithms
Low-Rank Matrix

Consider a group of basis functions and a group of testing functions, whose interactions
are represented as a sub-block in MoM matrices, Amn , Bmn , Cmn and Dmn , formulated
in section 4.4. When the two groups, the group of basis functions and the group of testing
functions, are near to each other, the EM interactions between any pair of testing and basis
functions will be sensitive to each other; on the contrary, when they are far-apart, the EM
interactions between any pair of testing and basis functions will be similar. In other words,
the corresponding sub-block of the matrices, which represents the far-apart groups, can be
re-described with fewer parameters. This kind of matrix is considered to be a low-rank
matrix, and the sub-blocks which represents near field interactions are full-rank matrices.
For instance, if an M-by-M sub-block [a] is low-rank, it can be represented as:

[a]M ×M = [u]M ×R [v]R×M + [e]M ×M

(5.1)

where, [e] is the error matrix and R is called rank of matrix [a].
To take advantage of the low-rank of the sub-blocks, the well-established singular-value

71

decomposition (SVD) method [119] can be used to derive a compressed form of the MoM
matrices. However, SVD method is computationally expensive. There are several methods
that can be applied to compress a low-rank matrix. One of the most popular methods,
adaptive cross-approximation (ACA), and another recently developed method, namely,
randomized algorithm, is introduced in the following sections.

5.1.2

Adaptive Cross-Approximation

Similar to the LU decomposition method with full pivoting, in cross-approximation (CA)
method, matrices [u] and [v] are constructed by minimizing the error matrix with given
tolerance. This method, however, needs a full knowledge of the matrix [a], and therefore it
is computationally expensive like the SVD method.
To minimize the computational requirements in CA method, the ACA method completes
the same task using the following procedure:
• first iteration
– randomly pick a row number I1 ∈ [1, M ],
– set v(1, :) = a(I1 , :),
– find column number J1 of the largest element in v(1, :),
a(:,J )

– set u(:, 1) = v(1,J1 ) ;
1
• second iteration
– find a row number I2 of the largest element in u(:, 1),
– set v(2, :) = a(I2 , :) − u(I2 , 1)v(1, :),
– find column number J2 of the largest element in v(2, :),
72

– set u(:, 2) =

a(:,J2 )−u(:,1)v(1,J2 )
;
v(2,J2 )

• k-th iteration
– find a row number Ik of the largest element in u(:, k − 1),
– set v(k, :) = a(Ik , :) −

k−1
i=1 u(Ik , i)v(i, :),

– find column number Jk of the largest element in v(k, :),
– set u(:, k) =

a(:,Jk )−

k−1 u(:,j)v(j,J )
k
j=1
.
v(k,Jk )

To terminate the iterations, we need to compute the Frobenius norm of the error matrix and
the original matrix at the k-th step. We can approximate both norms as:

||e|| ≈ ||e(k) || ≈ ||u(:, k)|| · ||v(k, :)||

(5.2)

||z||2 ≈ ||z (k) ||2 = ||z (k−1) ||2 + ||u(:, k)||2 + ||v(k, :)||2
k−1

|uT (:, j)u(k, :)| · |v(j, :)v T (:, k)|

+2

(5.3)

j=1

and stop the iteration once ||e|| < ε||z||.
We can see that the ACA method requires partial knowledge of the matrix [a] to find its
approximated compressed form. The entire process requires O(R2 M ) operations to obtain
the compressed matrix with a memory requirement of O(RM ), and once the approximation
is generated, O(RM ) operations will be needed to calculate each matrix-vector multiplication
[120, 121].

73

5.1.3

Randomized Algorithm

Recently, a fast randomized algorithm for the approximation of low-rank matrices has been
presented [122–124]. There are two variations of this algorithm: first, compress a matrix
by the randomized algorithm called interpolative decomposition (ID); second, use the ID to
obtain a much more efficient SVD method. By using this algorithm, we can construct a rank
k approximation [Z] from [A] at a cost proportional to O(m2 log(k) + 2l2 m).

5.1.3.1

Interpolative Decomposition

Find an m-by-k complex matrix [B] with k columns appropriately chosen from [A], and
construct a k-by-m matrix [P], such that

||P ||

4k(m − k) + 1

(5.4)

√
m kσk+1

(5.5)

and

||BP − A||

We pick an integer l close to k with k < l < m, and
• using a random number generator, form a l-by-m matrix [R] whose entries are
independent and identically distributed (i.i.d.) Gaussian random variables of zero
mean and unit variance, and compute the l-by-m product matrix: Y = RA;
• using the algorithm of [122], form a complex l-by-k matrix [Z] whose columns constitute
a subset of the columns of [Y ], and a complex k-by-n matrix [P ], such that some subset

74

of the columns of [P ] makes up the k-by-k identity matrix, and

4k(m − k) + 1ηk+1

||ZP − Y ||

(5.6)

where ηk+1 is the (k + 1)-th greatest singular value of [Y ];
• considering the fact that the columns of [Z] constitute a subset of the columns of
[Y ], for any j = 1, 2, · · · , k, an integer ij is existing such that the j-th column of
[Z] is the ij -th column of [Y ]. Form a real m-by-k matrix [B] from [A] that for any
j = 1, 2, · · · , k, the j-th column of [B] is the ij -th column of [A].

5.1.3.2

Converting

Interpolative

Decomposition

into

Singular

Value

Decomposition
Now we construct an SVD of matrix [A], such that

||U ΣV ∗ − A||

√
m kσk+1

(5.7)

where [U ] and [V ] are two column-wise orthonormal complex m-by-k matrices, [Σ] is a
non-negative real diagonal k-by-k matrix.
By using the ID, we can compute an approximated SVD of [A] as following:
• construct a lower triangular complex k-by-k matrix [L], and a complex m-by-k matrix
[Q] with orthonormal columns, and make those two matrices satisfy that P = LQ∗ ;
• compute the m-by-k product matrix C = BL;
• construct an SVD of [C], which is C = U ΣW ∗ , where [U ] is a column-wise orthonormal

75

complex m-by-k matrix, [Σ] is a real diagonal k-by-k matrix with non-negative entries,
and [W ] is a column-wise orthonormal complex k-by-k matrix;
• compute the m-by-k product matrix V = QW .
For the purpose of compressing the matrices, it will be sufficient to use ID only, and the
computational efficiency complexity of ID is O(m2 log(k) + lkm).

5.2
5.2.1

Adaptive Integral Method
FFT-based Method

This method is based on the well-known fast Fourier transform (FFT), where we calculate a
discrete convolution in O(N logN ) steps instead of a direct calculation in O(N 2 ) steps. For
certain geometries, it is possible to perform the calculation in equations (4.68) and (4.69)
using discrete convolutions. As a consequence, when we are using iterative methods, like
the stationary methods or Krylov subspace methods etc., to solve the linear system, we
can easily accelerate each matrix-vector multiplication by reducing O(N 2 ) double-floating
computations to O(N logN ) double-floating computations [125–141].
Once the geometry is meshed into uniform grids, the element values of MoM matrices,
Amn , Bmn , Cmn and Dmn , formulated in section 4.4 will only depend on the difference
between m and n. Such type of matrices is called Toeplitz matrix [142] and the matrix-vector
multiplication with this matrix can be calculated very efficiently using FFT. For example,
the multiplication Amn In can be written in the form of a cyclic convolution:
N

Amn In = A1n ⊗ In
n=1

76

(5.8)

This convolution can be computed by using the discrete Fourier transform (DFT):

p

[A]{I} = F −1 {F{An }

p

F{In }}

where F −1 and F denote inverse DFT and DFT respectively. Operator
p

(5.9)

is the Hadamard

p

product, {An } and {In } are expanded vectors defined as:

p

1,n

An =

p

In =





A

n = 1, 2, · · · , N
(5.10)




A1,2N −n n = N + 1, N + 2, · · · , 2N − 1




In n = 1, 2, · · · , N



0

(5.11)

n = N + 1, N + 2, · · · , 2N − 1

The DFTs can be easily calculated using FFT, and the requirements of memory and the
computation are proportional to O(N ) and O(N logN ) respectively. As discussed in [134]
for conducting plate and [141] for dielectric objects, this method shows great efficiency in
both aspects of computation time and memory requirement.

5.2.2

Adaptive Integral Method

For most problems, the solution domain is not easy to be meshed into uniform grids. To
remove this difficulty, the adaptive integral method (AIM) has been developed [29, 30, 143–
151]. The fundamental idea of AIM is firstly to project the sub-domain basis functions onto
uniform grids and then apply the FFT to the matrix-vector multiplication as described in
above section.
As shown in Figure 5.1, the projection of basis function ψm (x, y, z) can be calculated by the
77

Figure 5.1: of triangular RWG basis onto rectangular grids.
multipole moment approximation:
˚

(M +1)2

(xp − x0

)u (y

v
p − y0 ) (zp

− z0 )w Tmp

=
V

p=1

ψ(x, y, z)(x − x0 )u (y − y0 )v (z − z0 )w dV
(5.12)

where, 0 ≤ u, v, w ≤ (M +1)2 , r0 = (x0 , y0 , z0 ) is the reference point, and M is the projection
order. Examples are shown in Figure 5.1 with M = 3 and M = 4. With this definition in
equation (5.12), we can form a set of equations to solve for Tmp . Considering Amn as example
again, we can approximate it as:
(M +1)2 (M +1)2

Amn ≈

f ar
Amn

Tmp G(rp ; rq )Tnq

=
p=1

(5.13)

q=1

This projection only offers good accuracy when the basis and testing functions are far away
from each other. To correct the near fields relations, we define:

f ar

Anear
mn = Amn − Amn

78

(5.14)

where Anear
mn will be a sparse matrix when its small elements are neglected. As result, we
can calculate the matrix-vector multiplication as following:

[A]{I} = [Anear ]{I} + [T ]F −1 F{Gp }

[T ]T {I}

F

p

(5.15)

The number of computations of each matrix-vector multiplication is proportional to
O(N log N ) [144].

5.3

Fast Multipole Algorithm

In addition to the matrix compression techniques and AIM, another algorithm to speed up
the matrix-vector multiplication computation is fast multipole algorithm (FMA), which is
one of the most popular methods used in solving various EM problems [125, 152–173].
To use the FMA, we need to divide all the basis functions into groups, as shown in
Figure 5.2, where rp is the center of group Gp , rq is the center of group Gq , the distance
between rm belonging to Gp and rn belonging to Gq is broken in to three segments:
d1 = r q − r n , D = r p − r q , d2 = r m − r p .

rm − rn
n

m
d2
p

D
d1 q

Gp
Gq

Figure 5.2: Direction from rn to rm broken into three segments: d1 , D, d2 .
All groups are divided to two types with respective to each group: neighbouring groups and
far-field groups.
79

Meanwhile, we revisit the regimes of EM physics:
• low frequency, kd

1;

• middle frequency, kd ∼ O(1);
• high Frequency, kd

5.3.1

1.

Middle Frequency Regime

By using the addition theorem [125], the three-dimensional free space Green’s function can
be derived as:
e−jk|d+D|

∞

jk
(2)
ˆ
(−1)l (2l + 1)jl (kd)hl (kD)Pl (dˆ · D)
=−
4π
4π|d + D|
l=0

(5.16)

(2)

where d < D, jl (x) and hl (x) are spherical Bessel function and the second kind spherical
Hankel function respectively, and Pl (x) is Legendre polynomial.

From the elementary

identity [27]:
l
ˆ = j
jl (kd)Pl (dˆ · D)
4π

"

ˆ 2 kˆ
e−j k·d Pl (kˆ · D)d

(5.17)

we can represent the Green’s function as:
e−jk|rm −rn |
e−jk|d1 +d2 +D|
=
4π|rm − rn |
4π|d1 + d2 + D|
"
∞
jk
(2)
−j k·(d1 +d2 )
ˆ 2 kˆ
=−
e
(−1)l (2l + 1)hl (kD)Pl (kˆ · D)d
(4π)2

G(rm ; rn ) =

l=0

80

(5.18)

Thus,

1
G(rm ; rn ) ≈
jk

"

e−j k·(d1 +d2 ) α(k, D)d2 kˆ

(5.19)

where,

α(k, D) =

L
k 2
(2)
ˆ
(−1)l (2l + 1)hl (kD)Pl (kˆ · D)
4π

(5.20)

l=0

As a consequence, we can simplify the matrix-vector multiplication as:
"

N

Amn In ≈
n=1

Amn In +

Vmp (k, d1 ) ·

m∈Bp n∈Gq

Vqn (k, d2 )d2 kˆ

αpq (k, D)
m∈B
/ p

n∈Gq

(5.21)
where Bp denotes the neighbouring groups of Gp including Gp itself.
The FMA and its multilevel implementation, called multilevel FMA (MLFMA), offer
O(N 1.5 ) and O(N log N ) computational complexity respectively [174], and are shown to
achieve acceleration in mid-frequency EM scattering problems [152–156].
Additionally, the truncation number L of the series in equation (5.20) has been well studied by
1

analysing the truncation error [157–159], and a good choice of L is given by L = kD+6(kD) 3 .

5.3.2

Low Frequency Regime

As discussed above, in the low-frequency regime including EC problems, kd

1, the choice of

L is small which cause instability in the evaluation of the equation (5.19). As a consequence,
the representation shown in equation (5.18) is not suitable for low-frequency regime [160–

81

163].
Using the addition theory for spherical functions,

jl (kr)Ylm (θ, φ) =
∞

∞

l

4π(−i)l +l −l jl (kr )Yl m−m (θ , φ )

jl (kr )Yl m (θ , φ )
l =0

l =0 m =−l


· (−1)m



l
(2l + 1)(2l + 1)(2l + 1)  l l l   l


4π
0 0 0
−m m


l
m−m




(5.22)

(2)

hl (kr)Ylm (θ, φ) =
∞

l

l =0 m =−l

(2)
h (kr
l

∞

4π(−i)l +l −l jl (kr )Yl m−m (θ , φ )

)Yl m (θ , φ )
l =0


· (−1)m

f or

l
(2l + 1)(2l + 1)(2l + 1)  l l l   l


4π
0 0 0
−m m

r <r




l
m−m




(5.23)

82

we can derive the Green’s function as,

G(rm ; rn ) =

e−jk|d1 +d2 +D|
e−jk|rm −rn |
=
4π|rm − rn | 4π|d1 + d2 + D|
jk
=−
4π

∞

l

(2)
(−1)l (2l + 1)hl (kD)Y˜lm (θD , φD )

l=0 m=−l
l1
∞

·
l1 =0 m1 =−l1
l+l1

·
l2 =|l−l1 |



(−1)m (2l1 + 1)jl1 (kd1 )Y˜l∗ m (θd1 , φd1 )
1 1

il1 +l2 −l (2l2 + 1)jl2 (kd2 )Y˜l∗ m−m (θd2 , φd2 )
2
1




l
l l l   l
·

0 0 0
−m m

l
m−m

(5.24)




where, the spherical harmonics

Ylm (θ, φ) =

2l + 1 ˜
Y (θ, φ)
4π lm

(5.25)

Y˜lm (θ, φ) =

(l − m)! m
P (cosθ)e−imφ
(l + m)! l

(5.26)


l l
Plm (x) is given by associated Legendre polynomials and 
m m


l 
 is the Wigner 3-j
m

symbol [175].
With the expansion in equation (5.24), as in equation (5.21), we can represent the
matrix-vector multiplication as [160–167],
N

Amn In ≈
n=1

Amn In + Vmp (k, d1 ) ·
m∈Bp n∈Gq

αpq (k, D)
m∈B
/ p

83

Vqn (k, d2 )
n∈Gq

(5.27)

with,

jk
αpq (k, D) = −
4π

5.3.3

p

l

(2)
(−1)l (2l + 1)hl (kD)Y˜lm (θD , φD )

(5.28)

l=0 m=−l

Multilevel Implementation

To further enhance the advantages of using FMA, when the groups are really far away from
each other, we can aggregate a few far-field groups together into a larger group and transfer
the far-field EM interactions from the centres of the smaller groups to the center of the larger
group. As a result, we can represent the EM interactions using the large group instead of
the small groups to other far-away groups. For this purpose, a tree structure will be needed
as shown in Figure 5.3.

Figure 5.3: Illusions of groups – far-field groups for 3-level.

In the mid-frequency regime, the number of expansion terms L in equation (5.20) will vary
1

between different levels because of the varying distance D in the relation L = kD + 6(kD) 3 .
To solve this problem, the interpolation and anterpolation between levels have been discussed
84

[168–173]. For low frequency regime, we can maintain the number of expansion terms as
derived in equation (5.28) for different levels [160–163], thus we only need to regroup the
testing and basis functions.

5.4

Conclusion

Using the feature of FFT, AIM can be a powerful technique. However, for a 3D SIE problem,
AIM still needs a 3D uniform grid to take the advantage of FFT, which is a major drawback
of AIM. By taking the advantage of the low-rank matrix, several matrix compression methods
have been developed. As discussed in this chapter, both ACA and ID method are easy to be
computed in parallel. Meanwhile, the MLFMA can directly calculate the compressed form
of the MoM matrices and this method is also easy to be parallelized. To take advantage
of the compressed form of the MoM matrices using either matrix compression methods or
MLFMA, we can further construct a direct integral solver. Matrix representation of MLFMA
and construction of direct integral is discussed in the following chapter.

85

Chapter 6
Direct Integral Solver

6.1

Matrix Representation for Fast Algorithms

6.1.1

Single level Matrix Representation

We can split the MoM matrices, for example [A], into two parts, near-field interactions
[D], which is very sparse, and far-field interactions [Af ar ]. By using FMA as discussed in
section 5.3, [Af ar ] is a low-rank matrix which can be represented in compressed form as
following,

A ≈ D + Af ar = D + LSR

(6.1)

where, [L] is the multipole expansion of grouped sources, [R] is the multipole expansion
of grouped observations, and [S] is the kernel expansion as shown in equation (5.20) for
mid-frequency and equation (5.28) for low-frequency.
The accelerated matrix-vector multiplication can be clearly seen:

Ax = Dx + LSRx

86

(6.2)

For the linear system, Ax = b, furthermore, we can represent the compressed matrices in a
more compact matrix form,

D


R



   
L
 x  b 
   
   
 = 
−I 
 y  0
   
0
z
−I S

(6.3)

with z = Rx and y = Sz. With this form, we can use the sparsity of the matrix and
precondition them very efficiently for further saving computational resources.

6.1.2

Multilevel Matrix Representation

By using the same ideas as discussed above, for MLFMA, we can represent the MoM matrix
[A] as:

A ≈ D(1) + L(1) [D(2) + L(2) (· · · D(λ) + L(λ) SR(λ) · · · )R(2) ]R(1)

(6.4)

where, [L(1) ] is the multipole expansion of grouped sources, [R(1) ] is the multipole expansion
of grouped observations, [D(λ) ] is the near-field interactions of λ-th level, and [S] is the
kernel expansion of the sparsest level. In the mid-frequency regime, with λ > 1, [L(λ) ] and
[R(λ) ] are the interpolation and anterpolation matrices between (λ − 1)-th level and λ-th
level respectively, whereas in the low frequency regime, [L(λ) ] and [R(λ) ] are the regrouping
permutation matrices between (λ − 1)-th level and λ-th level.

87

Again, we can represent the compressed matrices in a more compact matrix form as following:






 

D(1) L(1)
  x  b


  


  
 R(1)
  y (1)  0
−I


  


  


  
(2) L(2)

  z (1)  0
−I
D


  


  



..  =  .. 
..
..
(2)



  
.
.
R

  .  .


  

  (λ−1)   
..
. D(λ) L(λ)

 z
 0


  


  


  
R(λ)
−I   y (λ)  0



  


  
−I
S
z (λ)
0

(6.5)

with z (1) = R(1) x and y (λ) = S (λ) z (λ) . As a result, the number of non-zero entries of the
matrix on the left-hand side of equation (6.5) is proportional to O(N log N ), whereas [A] is
an N -by-N full matrix.

6.2

Direct Integral Solver

Once we represent the MoM matrices in the form of equation (6.1), we can directly calculate
the inversion of [A] as following [176],

A−1 ≈ D−1 − D−1 L(RD−1 L)−1 RD−1
+ D−1 L(RD−1 L)−1 [(RD−1 L)−1 + S]−1 (RD−1 L)−1 RD−1

(6.6)

For the multi-level representation in equation (6.4), we can obtain [176],

A−1 ≈ D(1) + L(1) [D(2) + L(2) (· · · D(λ) + L(λ) SR(λ) · · · )R(2) ]R(1)

88

(6.7)

where,

D(λ) = [D(λ) ]−1 − [D(λ) ]−1 L(λ) (R(λ) [D(λ) ]−1 L(λ) )−1 R(λ) [D(λ) ]−1

(6.8)

L(λ) = [D(λ) ]−1 L(λ) (R(λ) [D(λ) ]−1 L(λ) )−1

(6.9)

R(λ) = (R(λ) [D(λ) ]−1 L(λ) )−1 R(λ) [D(λ) ]−1

(6.10)

S = [(R(λ) [D(λ) ]−1 L(λ) )−1 + S]−1

(6.11)

If the matrix [D] is not in the form of block diagonal, the computation of (RD−1 L)−1 will
be significantly expensive. Thus, using the matrix form as in equations (6.3) and (6.5) will
be a suitable choice. Otherwise, as studied in [176–180], the direct method is ideal for the
case of varying right-hand side.

6.3

Analysis of Computation Complexity

Typical EC problems always involve hundreds of simulations at each scan position of
the probe on test specimen as mentioned in section 2.5. In effect this implies that we
have hundreds of right-hand sides of the system of equations.

Therefore, there is a

significant advantage in deriving the direct inversion of matrix A as shown in equation (6.7).
Consequently, it is better to choose matrix compression methods over MLFMA for this case.
By using the matrix compression algorithms, the following procedure is needed:
1. form MoM matrix A from the test specimen;
2. calculate the multiple right-hand sides B = [b1 , b2 , · · · , bNscan ];
3. compress the MoM matrix and use equation (6.7) to form its direct inversion A−1 ;

89

4. obtain the solution by calculating X = [x1 , x2 , · · · , xNscan ] = A−1 B.

6.3.1

Computational Scalability

From section 4.4, we can clearly see that the first two steps is highly scalable, since the
calculation of each entry in MoM matrix A and right-hand sides matrix B are independent
of each other. Similarly, the last step is also highly scalable, because the calculation of each
entry in matrix X is obtained from matrix–matrix multiplication, which is also independent
of each other. Further in the GPU implementation, by using the feature of independent
computation of each entry, we can easily program the computation of each entry in matrices
A, B and X as kernel functions by CUDA. The computation time for the first step, namely
forming matrix A, is shown in Table 6.1.
Unknowns
2920
6890
15500

CPU (1 core)
122s
782s
3956s

CPU (10 cores)
12s
82s
405s

CPU (20 cores)
7s
41s
207s

GPU (2496 cores)
<1s
∼3s
∼9s

Table 6.1: Computation time of matrix A.

In this test, the CPU that is considered is Intel

Xeon

Processor E5-2670 v2 (25M Cache,

2.50 GHz), and the GPU that is considered is nVIDIA Tesla K20 (5GB, 706MHz).
From the timing results, we can see that the speed-up ratio is proportional to the number
cores under same core clock frequency. Although the core clock frequency of the GPU is
lower than that of CPU, due to larger amount of cores on GPU, we still obtain significant
speed-up.

90

6.3.2

Matrix Compression

As discussed in section 6.1.2, theoretically, the memory usage of the compressed form of
matrix A is proportional to O(N log N ). However, both the practical memory usage and the
compression time are proportional to O(N 1.8 ) as shown in Table 6.2.
Unknowns
2920
6890
15500
26030
59590

Compression
levels
2
3
4
5
5

Memory Usage
Original
Compressed
130.103 MiB 127.078 MiB
724.367 MiB 163.395 MiB
3.58 GiB
1.34269 GiB
10.0964 GiB 3.53157 GiB
52.0293 GiB 19.0287 GiB

Compression Time
RRQR
ID
6.69745s 6.65206s
27.5766s 29.6826s
393.642s 410.042s
1699.68s 1665.99s
15064.6s 14911.4s

Table 6.2: Matrix compression: memory usage and compression time.

In this test, the CPU that is considered is Intel

Xeon

Processor E5-2690 (20M Cache,

2.90 GHz).
In theory, the ID method as described in section 5.1.3.1 is much faster than the rank-revealing
QR (RRQR) decomposition. Several studies indicate that ID method is on an average 10∼11
times faster than RRQR [123, 181, 182]. However, the efficiency of ID method is based on
the apriori knowledge of the actual rank of the target matrix. Once a trial rank is smaller
than the actual rank of the target matrix, re-applying ID method with a larger trail rank is
necessary, which will significantly draw down the efficiency of ID method. Consequently, the
performance of ID method is about the same as RRQR decomposition as shown in Table 6.2.

6.4

Conclusion

Either matrix compression algorithms or MLFMA can be applied to form the direct integral
solver for large-scale EC problems. Although the MLFMA is more memory efficient than

91

matrix compression algorithms, it is impractical to resolve a direct inversion of the MoM
matrix, and the evaluations in MLFMA is also expensive. On the other hand, the matrix
compression algorithms is less memory efficient. However, it is sufficient to form a direct
inversion of the MoM matrix, which is a considerable advantage when we have multiple
right-hand sides such as that encountered in a scanning probe test situation.

92

Chapter 7
Summary and Future Work

7.1

Summary

The contributions of this thesis consists of two parts. In the first part, a FEM based
computational model has been developed for simulating EC inspection of SG tubes using
a variety of EC probes.

Validation results demonstrate that the model predicts the

experimental signals from a variety of degradation mechanisms very accurately.
This FEM model has also been applied to assist a novel design of EC probe based on rotating
magnetic fields for the inspection of SG tubes in nuclear power plants. Based on this new
design, we could achieve high performance speed and high resolution measurement for full
tube inspection.
In the second part, the MoM with higher-order basis function has been presented to solve the
EC problem. Because the near-field EM interaction plays a major role, singularity removal
techniques has been applied to compute the EC field distributions with adequate accuracy.
To accelerate the MoM, fast algorithms, namely FFT-based algarithm, AIM, ACA method,
ID method and MLFMA have been reviewed, and the ID method and MLFMA has been
studied. Due to low-frequency break down for EC problems, an alternative low-frequency
MLFMA has been presented. Furthermore, a newly proposed direct integral solver has been
introduced to solve both low-frequency and mid-frequency problems. The direct integral
solver has been developed using both MLFMA and RRQR/ID matrix compression methods
93

in this study.

7.2
7.2.1

Future Work
Finite Element Method

As shown in part I, we have a well-written FORTRAN implementation of the FEM solver.
It produces the correct results and is widely used in various applications in the NDE Lab.
However, the current solver has several major drawbacks:
1. it supports only certain types of elements;
2. it is single threaded;
3. it is not possible to launch it on distributed machines;
4. it is difficult to accelerate the code on the modern hardware.
In order to make the code more efficient, many improvements can be made.

7.2.1.1

Refactoring the existing code

Implementation of this FEM solver on GPU will require utilization of some external
optimized libraries. In theory, it is possible to link the FORTRAN codes with C, but for our
applications it is complicated due to implementation aspects. Hence it will be advantageous
to translate the code to C programming language.
The coding standards and techniques in C/C-style C++/C++ languages are different from
the ones used in FORTRAN. It is necessary to switch to the modular architecture to make
the code more clear and decrease the optimization time. In addition, it can also provide a
basis for further improvement using vector elements and higher-order scheme.
94

7.2.1.2

Optimization

In general, the FE solver takes two steps: form the stiffness matrix and then solve the sparse
system of equations.
To accelerate the computation of the stiffness matrix, colouring scheme should be considered.
The basic idea of the colouring scheme is to separate the elements into several groups. Within
each group, make sure no element connect to any other elements in this group. Thus, the
calculation of all the elements in a same group are independent of each other, and such
independent calculations can be easily conducted in parallel. In theory, this scheme is highly
parallel and especially good for the GPU implementation.
The acceleration of solving the sparse system of equations could be initially achieved by
implementing the latest fast sparse matrix solver, SuiteSparse, in the existing FE solver.
Based on this sparse matrix solver, further optimizations can be considered.

7.2.2

Direct Integral Solver

As discussed in part II, direct integral solver shows significant advantages in several aspects.
However, as a new algorithm, there are several issues that will need to be further studied.

7.2.2.1

GPU Implementation of the Matrix Compression Method

As shown in section 6.3, there are four major steps of direct integral solver.

GPU

implementation of the first two steps and the last step is discussed in section 6.3. The
GPU implementation of third step, matrix compression algorithm, will need to be further
studied.

95

7.2.2.2

Optimization of Moment Method Galerkin Testing

As discussed in [183], the accuracy of the SIE may be affected by MoM Galerkin testing.
Proper choice of the quadrature points could considerably improve the accuracy of modelling
certain geometries with sharp edges. Accurate evaluation of the SIE will need to be further
studied, especially for EC problems.

7.2.2.3

Multi-body Problems

If we have more than one object of interest, we need to construct a system of equations to
solve these problems. A two-body problem as shown in Figure 7.1 will be discussed as an
example in the following.

Figure 7.1: A two-body surface integral equation problem.

96

Similar to section 4.4, apply basis functions to expand the surface current density as:

(1)
Js

(1)

Ms

(2)

Js

(2)

Ms

N

(1) (1)

=

In Λn (r)
n=1
N

=

(7.1)

(1) (1)

Kn Λn (r)
n=1
N

=

(7.2)

(2) (2)

In Λn (r)
n=1
N

=

(7.3)

(2) (2)

Kn Λn (r)

(7.4)

n=1

for the first object,

(1)
Amn
(1)

Bmn
(1)

Cmn
(1)

Dmn
(1)

fm

"

(1)

(1)

(1)

= Z0
Λm · Λn + n
ˆ × K(Λn ) dS
S
"
(1)
(1)
Λm · n
ˆ × L(Λn ) dS
=
S
"
(1)
(1)
(1)
Λm · Λn + n
ˆ × Ki (Λn ) dS
=Z
S
"
(1)
(1)
Λm · n
ˆ × Li (Λn ) dS
=
S
"
(1)
Λm · n
ˆ × H inc dS
=

(7.5)
(7.6)
(7.7)
(7.8)
(7.9)

S

(1)

gm = 0

(7.10)

97

for the second object,

(2)
Amn
(2)

Bmn
(2)

Cmn
(2)

Dmn
(2)

fm

"

(2)

(2)

(2)

= Z0
Λm · Λn + n
ˆ × K(Λn ) dS
S
"
(2)
(2)
=
Λm · n
ˆ × L(Λn ) dS
S
"
(2)
(2)
(2)
=Z
Λm · Λn + n
ˆ × Ki (Λn ) dS
S
"
(2)
(2)
=
Λm · n
ˆ × Li (Λn ) dS
S
"
(2)
=
Λm · n
ˆ × H inc dS

(7.11)
(7.12)
(7.13)
(7.14)
(7.15)

S

(2)

gm = 0

(7.16)

between two objects,

(12)
Amn

"

(1)

(2)

(2)

= Z0
Λm · Λn + n
ˆ × K(Λn ) dS
S
"
(1)
(2)
(12)
Λm · n
ˆ × L(Λn ) dS
Bmn =
S
"
(21)
(2)
(1)
(1)
Amn = Z0
Λm · Λn + n
ˆ × K(Λn ) dS
S
"
(2)
(1)
(21)
Λm · n
ˆ × L(Λn ) dS
Bmn =

(7.17)
(7.18)
(7.19)
(7.20)

S

Thus, we can form,


A(1)

B (1)



 (1)
C
D(1)


 (21)
A
B (21)


0
0

A(12)
0
A(2)
C (2)



(1)

 f 

 

  (1)  





0  K   0 


=

  (2)   (2) 
(2)




B   I  f 


 

(2)
(2)
D
K
0

B (12)

98



I (1)



(7.21)

The application of the direct integral algorithm to this problem should be conducted. Further
more, if those objects are moving relatively to each other, such as a scanning coil with a
ferrite core in EC problems as described in section 2.5, the procedure of applying the direct
integral algorithm and corresponding GPU implementation will need to be studied.

99

APPENDIX

100

Graphic User Interface (GUI)
A procedural overview of the use of the software to predict SG inspection signals is
summarized in the following steps:
Select
Geometry

Select
Defect

Select
Probe

Run
Simulation

Generate
Signal

Figure 1: Block diagram of simple operation of steam generator tube simulator software.

The different steps are provided as commands in the tool bar as shown in Figure 2. User
can select from two predefined tube geometries: free span and support plate (as shown in
Figure 2). A rectangular defect of specific length, depth and width can be introduced in
the tube wall on the ID or OD along axial or circumferential direction. Probe coil and
material properties in different parts of the domain can also be selected from a pull-down
menu options on the commands tool bar as seen in Figure 21 .

1 From

Materials Evaluation, Vol. 69, No. 12. Reprinted with permission of the American
Society for Nondestructive Testing, Inc.
101

Figure 2: Screen-shot showing various commands available in the simulation software.

102

Figure 3: Analysis screen for the rotating probe coil.

103

Figure 4: Display of the 1-D coil analysis screen (bobbin probe).
104

When the user selects “Run”, the graphic user interface passes geometrical information to the
mesh generator. The latter generates a 3-D FE mesh of the tube (and support plate), core
and coil, forms probe scanning grid and saves mesh and grid data. The graphic user interface
passes material property values to the solver. After the simulation is done, the calculated
induced voltages at each probe position are stored, calibrated and displayed. Screen-shots
of the displays for bobbin and rotating coil probes are shown in Figures 3 and 42 .

2 From

Materials Evaluation, Vol. 69, No. 12. Reprinted with permission of the American
Society for Nondestructive Testing, Inc.
105

BIBLIOGRAPHY

106

BIBLIOGRAPHY
[1] S. S. Udpa and P. O. Moore, eds., Nondestructive Testing Handbook, vol. 5:
Electromagnetic Testing. American Society for Nondestructive Testing, 3rd ed., 2004.
[2] W. Lord and R. Palanisamy, “Development of theoretical models for NDT eddy
current phenomena,” in Symposium on Eddy Current Characterization of Materials
and Structures, (Gaithersburg, Maryland), NBS, September 1979.
[3] J.-M. Jin, The Finite Element Method in Electromagnetics. Wiley-IEEE Press, 2nd ed.,
May 27 2002.
[4] J. H. Hwang and W. Lord, “Finite element modeling of magnetic field/defect
interactions,” ASNT Journal of Testing and Evaluation, vol. 3, no. 1, 1975.
[5] J. H. Hwang and W. Lord, “Magnetic leakage field signatures of material
discontinuities,” in Proceedings of the 10th Symposium on NDE, (San Antonio, Texas),
1975.
[6] W. C. Yen, “Finite element characterization of residual leakage fields,” m.s. thesis,
Colorado State University, 1978.
[7] S. S. Udpa and W. Lord, “Finite element modeling of residual magnetic phenomenon,”
in International Magnetics Conference, (Boston), Apr. 1980.
[8] W. Lord and R. Palanisamy, “Magnetic probe inspection of steam generator tubing,”
Materials Evaluation, vol. 38, May 1980.
[9] R. Palanisamy, Finite Element Modeling of Eddy Current Nondestructive Testing
Phenomena. PhD thesis, Colorado State University, 1980.
[10] N. Ida, Three Dimensional Finite Element Modeling of Electromagnetic NDT
Phenomena. PhD thesis, Colorado State University, 1983.
[11] O. B´ır´o, “Edge element formulations of eddy current problems,” Computer Methods
in Applied Mechanics and Engineering, vol. 169, 1999.

107

[12] J. B. Manges and Z. J. Cendes, “A generalized tree-cotree gauge for magnetic field
computation,” IEEE Trans. On Magnetics, vol. 31, May 1995.
[13] K. Preis, O. B´ır´o, H. Reisinger, K. Papp, and I. Ticar, “Eddy current losses in large
air coils with layered stranded conductors,” IEEE Trans. On Magnetics, vol. 44, June
2008.
[14] R. Albanese and G. Rubinacci, “Magneto-static field computations in terms of
two-component vector potentials,” International Journal for Numerical Methods in
Engineering, vol. 29, 1990.
[15] O. B´ır´o, K. Preis, and K. R. Richter, “On the use of the magnetic vector potential
in the nodal and edge finite element analysis of 3D magneto-static problems,” IEEE
Trans. On Magnetics, vol. 32, May 1996.
[16] Z. Zeng, X. Liu, Y. Deng, and L. Udpa, “Reduced magnetic vector potential and electric
scalar potential formulation for eddy current modeling,” Przeglad Elektrotechniczny,
vol. 83, June 2007.
[17] O. B´ır´o and K. Preis, “An edge finite element eddy current formulation using a reduced
magnetic and a current vector potential,” IEEE Trans. On Magnetics, vol. 36, Sept.
2000.
[18] K. Preis, I. Bardi, O. Biro, C. Magele, W. Renhart, K. R. Richter, and G. Vrisk,
“Numerical analysis of 3D magneto-static fields,” IEEE Trans. On Magnetics, vol. 77,
Sept. 1991.
[19] O. B´ır´o, K. Preis, W. Renhart, K. Richter, and G. Vrisk, “Perpormance of
different vector potential formulations in solving multiply connected 3-D eddy current
problems,” IEEE Trans. On Magnetics, vol. 26, Mar. 1990.
[20] Z. Zeng, L. Udpa, and S. S. Udpa, “An efficient finite element method for modeling
ferrite-core eddy current probe,” International Journal of Applied Electromagnetics
and Mechanics, vol. 33, no. 1-2, 2010.
[21] M. Melton, “Materials degradation management - NEI 03-08 applications to
steam generator management program,” in Steam Generator Management Program:
Proceedings of the 26th Steam Generator NDE Workshop, (Big Sky, MT), Electric
Power Research Institute, July 2007.
[22] Z. Zeng, L. Udpa, and S. S. Udpa, “Finite-element model for simulation of ferrite-core
eddy-current probe,” IEEE Trans. on Magnetics, vol. 46, March 2010.
108

[23] N. Lei, C. Bardel, L. Udpa, S. S. Udpa, J. Benson, and R. Williams, “Simulation for
eddy current steam generator inspection,” Materials Evaluation, vol. 69, Dec. 2011.
[24] N. Lei, L. Udpa, S. Udpa, and Z. Zeng, “Rotating field eddy current (RoFEC)-probe
for steam generator inspection,” International Journal of Applied Electromagnetics and
Mechanics, vol. 33, pp. 1279–1285, Oct. 2010.
[25] O.-M. Kwon, M. V. K. Chari, S. J. Salon, and K. Sivasubramaniam, “Development
of integral equation solution for 3-D eddy current distribution in a conducting body,”
IEEE Trans. On Magnetics, vol. 39, Sept. 2003.
[26] O.-M. Kwon, M. V. K. Chari, and S. J. Salon, “Three-dimensional (3D) eddy current
time-domain integral equation with Coulomb gauge condition,” Journal of Applied
Physics, vol. 97, May 2005.
[27] J. A. Stratton, Electromagnetic Theory. New York: McGraw-Hill, 1941.
[28] B. H. McDonald and A. Wexler, Finite Elements in Electrical and Magnetic
Field Problems: Mutually constrained partial differential and integral equation field
formulations. New York: Wiley, 1980.
[29] G. Rubinacci, A. Tamburrino, S. Ventre, and F. Villone, “A fast algorithm for solving
3-D eddy current problems with integral formulations,” IEEE Trans. On Magnetics,
vol. 37, Sept. 2001.
[30] E. Cardelli, A. Faba, R. Specogna, A. Tamburrino, F. Trevisan, and S. Ventre,
“Analysis methodologies and experimental benchmarks for eddy current testing,” IEEE
Trans. On Magnetics, vol. 41, May 2005.
[31] G. Rubinacci, A. Tamburrino, S. Ventre, and F. Villone, “A fast 3-D multipole method
for eddy-current computation,” IEEE Trans. On Magnetics, vol. 40, Mar. 2004.
[32] G. Rubinacci, A. Tamburrino, and F. Villone, “Circuits/fields coupling and multiply
connected domains in integral formulations,” IEEE Trans. On Magnetics, vol. 38, Mar.
2002.
[33] G. Rubinacci, A. Tamburrino, and F. Villone, “A novel integral formulation for the
solution of Maxwell equations,” IEEE Trans. On Magnetics, vol. 39, May 2003.
[34] G. Rubinacci, A. Tamburrino, and S. Ventre, “An efficient numerical model for a
magnetic core eddy-current probe,” IEEE Trans. On Magnetics, vol. 44, June 2008.
109

[35] G. Rubinacci and A. Tamburrino, “A broadband volume integral formulation based on
edge-elements for full-wave analysis of lossy interconnects,” IEEE Trans. On Antennas
and Propagation, vol. 54, Oct. 2006.
[36] L. D. Negro, G. Miano, G. Rubinacci, A. Tamburrino, and S. Ventre, “A fast
computation method for the analysis of an array of metallic nanoparticles,” IEEE
Trans. On Magnetics, vol. 45, MARCH 2009.
[37] G. Rubinacci and A. Tamburrino, “Automatic treatment of multiply connected regions
in integral formulations,” IEEE Trans. On Magnetics, vol. 46, AUGUST 2010.
[38] G. Miano, G. Rubinacci, A. Tamburrino, and F. Villone, “Linearized fluid model for
plasmon oscillations in metallic nanoparticles,” IEEE Trans. On Magnetics, vol. 44,
June 2008.
[39] G. Miano, G. Rubinacci, and A. Tamburrino, “Numerical modeling for the analysis
of plasmon oscillations in metallic nanoparticles,” IEEE Trans. On Antennas and
Propagation, vol. 58, Sept. 2010.
[40] G. Rubinacci and A. Tamburrino, “Efficient modeling of the skin-effect in low frequency
regime,” in International Conference on Electromagnetics in Advanced Applications,
2007.
[41] G. Rubinacci, A. Tamburrino, S. Ventre, and F. Villone, “Fast computational methods
for large-scale eddy-current computation,” IEEE Trans. On Magnetics, vol. 38, Mar.
2002.
[42] T. Theodoulidis and J. Bowler, “Eddy-current interaction of a long coil with a slot in
a conductive plate,” IEEE Trans. On Magnetics, vol. 41, April 2005.
[43] J. R. Bowler, S. A. Jenkins, L. D. Sabbagh, and H. A. Sabbagh, “Eddy-current probe
impedance due to a volumetric flaw,” Journal of Applied Physics, vol. 70, Aug. 1991.
[44] J. R. Bowler and N. Harfield, “Evaluation of probe impedance due to thin-skin
eddy-current interaction with surface cracks,” IEEE Trans. On Magnetics, vol. 34,
March 1998.
[45] J. R. Bowler, “Time domain half-space dyadic Green’s functions for eddy-current
calculations,” Journal of Applied Physics, vol. 86, Dec. 1999.
[46] S. J. Norton and J. R. Bowler, “Theory of eddy current inversion,” Journal of Applied
Physics, vol. 73, Jan. 1991.
110

[47] Y. Yoshida and J. R. Bowler, “Vector potential integral formulation for eddy-current
probe response to cracks,” IEEE Trans. On Magnetics, vol. 36, March 2000.
[48] J. R. Bowler and F. Fu, “Time-domain dyadic Green’s function for an electric source
in a conductive plate,” IEEE Trans. On Magnetics, vol. 42, Nov. 2006.
[49] F. Fu and J. Bowler, “Transient eddy-current driver pickup probe response due to a
conductive plate,” IEEE Trans. On Magnetics, vol. 42, Aug. 2006.
[50] J. P´av´o and D. Lesselier, “Calculation of eddy current testing probe signal with global
approximation,” IEEE Trans. On Magnetics, vol. 42, April 2006.
[51] Y. L. Bihan, J. P´av´o, M. Bensetti, and C. Marchand, “Computational environment
for the fast calculation of ECT probe signal by field decomposition,” IEEE Trans. On
Magnetics, vol. 42, April 2006.
[52] B. Rejaei, “Mixed-potential volume integral-equation approach for circular spiral
inductors,” IEEE Trans. On Microwave Theory and Techniques, vol. 52, Aug. 2004.
[53] D. Chen, K. R. Shao, and J. D. Lavers, “Very fast numerical analysis of benchmark
models of eddy-current testing for steam generator tube,” IEEE Trans. On Magnetics,
vol. 38, Sept. 2002.
[54] T. Theodoulidis and J. R. Bowler, “Interaction of an eddy-current coil with a
right-angled conductive wedge,” IEEE Trans. On Magnetics, vol. 46, April 2010.
[55] T. P. Theodoulidis and J. R. Bowler, “Eddy current coil interaction with a right-angled
conductive wedge,” Proc. R. Soc. A, vol. 461, Oct. 2005.
[56] H. Sun, J. R. Bowler, and T. P. Theodoulidis, “Eddy currents induced in a finite length
layered rod by a coaxial coil,” IEEE Trans. On Magnetics, vol. 41, Sept. 2005.
[57] J. R. Bowler and T. Theodoulidis, “Boundary element calculation of eddy currents in
cylindrical structures containing cracks,” IEEE Trans. On Magnetics, vol. 45, March
2009.
[58] A. Skarlatos, G. Pichenot, D. Lesselier, M. Lambert, and B. Duchne, “Electromagnetic
modeling of a damaged ferromagnetic metal tube by a volume integral equation
formulation,” IEEE Trans. On Magnetics, vol. 44, May 2008.

111

[59] S. K. Burke and T. P. Theodoulidis, “Impedance of a horizontal coil in a borehole: a
model for eddy-current bolt-hole probes,” J. Phys. D: Appl. Phys., vol. 37, 2004.
[60] J. R. Bowler, “Eddy-current interaction with an ideal crack. I. the forward problem,”
J. Appl. Phys., vol. 75, June 1994.
[61] P. Alotto, M. Guarnieri, and F. Moro, “A boundary integral formulation on
unstructured dual grids for eddy-current analysis in thin shields,” IEEE Trans. On
Magnetics, vol. 43, April 2007.
[62] M. H. Pham and A. J. Peyton, “A model for the forward problem in magnetic induction
tomography using boundary integral equations,” IEEE Trans. On Magnetics, vol. 44,
Oct. 2008.
[63] I. R. Ciric and R. Curiac, “A single-source surface integral formulation for eddy-current
problems,” IEEE Trans. On Magnetics, vol. 40, March 2004.
[64] J. P´av´o, “Approximate methods for the calculation of the ECT signal of a crack in a
plate coated by conducting deposit,” IEEE Trans. On Magnetics, vol. 40, March 2004.
[65] M. H. Pham and A. J. Peyton, “Computation of 3-D sensitivity coefficients in magnetic
induction tomography using boundary integral equations and radial basis functions,”
IEEE Trans. On Magnetics, vol. 44, Oct. 2008.
[66] P. Beltrame and N. Burais, “Computing methods of hypersingular integral applied to
eddy-current testing,” IEEE Trans. On Magnetics, vol. 38, March 2002.
[67] K. Ishibashi, “Eddy current analysis by boundary element method utilizing impedance
boundary condition,” IEEE Trans. On Magnetics, vol. 31, May 1995.
[68] K. Ishibashi, “Eddy current analysis by BEM utilizing edge boundary condition,” IEEE
Trans. On Magnetics, vol. 12, May 1996.
[69] K. Ishibashi, “Eddy current analysis of a conductor with a groove by surface integral
equation with one unknown,” IEEE Trans. On Magnetics, vol. 41, May 2005.
[70] K. Ishibashi, “Eddy-current analysis of a conductor with a crack by line integral
equations with loop electric and magnetic currents as unknowns,” IEEE Trans. On
Magnetics, vol. 42, April 2006.

112

[71] J. Smajic, Z. Andjelic, and M. Bebendorf, “Fast BEM for eddy-current problems using
H-matrices and adaptive cross approximation,” IEEE Trans. On Magnetics, vol. 43,
April 2007.
[72] Y. Higuchi and M. Koizumi, “Integral equation method with surface impedance model
for 3D eddy current analysis in transformers,” IEEE Trans. On Magnetics, vol. 36,
July 2000.
[73] H. Fujita and K. Ishibashi, “Nonlinear eddy current analysis of thin steel plate by
boundary integral equations,” IEEE Trans. On Magnetics, vol. 44, June 2008.
[74] K. Ishibashi, “Numerical analysis of eddy current testing by integral equation method,”
IEEE Trans. On Magnetics, vol. 37, Sept. 2001.
[75] C.-H. Im and H.-K. Jung, “Numerical emulator for walk-through metal detectors using
3-D indirect boundary integral equation method,” IEEE Trans. On Instrumentation
and Measurement, vol. 54, June 2005.
[76] D. Ioan and M. Rebican, “Numerical model for eddy-current testing of ferromagnetic
steel parts,” IEEE Trans. On Magnetics, vol. 38, March 2002.
[77] Y. Fujishima and S. Wakao, “Surface charge analysis in eddy current problems,” IEEE
Trans. On Magnetics, vol. 39, May 2003.
[78] K.Satake, M.Tanaka, N.Shimizu, Y.Araki, and K.Morimoto, “Three-dimensional
analysis on eddy current testing for S/G tubes,” IEEE Trans. On Magnetics, vol. 28,
March 1992.
[79] R. F. Harrington, Time-Harmonic Electromagnetic Fields.
September 2001.

Wiley-IEEE Press,

[80] C. A. Balanis, Advanced Engineering Electromagnetics. Wiley, 2nd ed., January 2012.
[81] J. A. Kong, Electromagnetic Wave Theory. EMW Publishing, 2008.
[82] E. J. Rothwell and M. J. Cloud, Electromagnetics. CRC Press, 2001.
[83] J.-M. Jin, Theory and Computation of Electromagnetic Fields. Wiley-IEEE Press,
November 2010.

113

[84] L. Tsang, J. A. Kong, and K. H. Ding, Scattering of Electromagnetic Waves: Theories
and Applications. Wiley-Interscience, 2000.
[85] L. Tsang, J. A. Kong, K. H. Ding, and C. Ao, Scattering of Electromagnetic Waves:
Numerical Simulations. Wiley-Interscience, 2001.
[86] L. Tsang and J. A. Kong, Scattering of Electromagnetic Waves: Advanced Topics.
Wiley-Interscience, 2001.
[87] C.-T. Tai, Dyadic Green Functions in Electromagnetic Theory. IEEE Press, 2nd ed.,
January 1994.
[88] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by surfaces
of arbitrary shape,” IEEE Trans. Antennas Propag., vol. 30, no. 3, 1982.
[89] R. D. Graglia, D. R. Wilton, and A. F. Peterson, “Higher order interpolatory vector
bases for computational electromagnetics,” IEEE Transactions on Antennas and
Propagation, vol. 45, March 1997.
[90] I. Bogaert and D. D. Zutter, “High precision evaluation of the selfpatch integral
for linear basis functions on flat triangles,” IEEE Transactions on Antennas and
Propagation, vol. 58, May 2010.
[91] R. D. Graglia and G. Lombardi, “Machine precision evaluation of singular and nearly
singular potential integrals by use of Gauss quadrature formulas for rational functions,”
IEEE Transactions on Antennas and Propagation, vol. 56, April 2008.
[92] M. A. Khayat and D. R. Wilton, “Numerical evaluation of singular and near-singular
potential integrals,” IEEE Transactions on Antennas and Propagation, vol. 53, Oct.
2005.
[93] L. Rossi and P. J. Cullen, “On the fully numerical evaluation of the linear-shape
function times the 3-D Green’s function on a plane triangle,” IEEE Transactions on
Microwave Theory and Techniques, vol. 47, April 1999.
[94] R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D
Green’s function or its gradient on a plane triangle,” IEEE Transactions on Antennas
and Propagation, vol. 41, Oct. 1993.
[95] D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Buncak, and C. M.
Butler, “Potential integrals for uniform and linear source distributions on polygonal
114

and polyhedral domains,” IEEE Transactions on Antennas and Propagation, vol. 32,
March 1984.
[96] A. Tzoulis and T. F. Eibert, “Review of singular potential integrals for method of
moments solutions of surface integral equations,” Advances in Radio Science, vol. 2,
2004.
[97] S.-W. Lee, J. Boersma, C.-L. Law, and G. Deschamps, “Singularity in Green’s function
and its numerical evaluation,” IEEE Transactions on Antennas and Propagation,
vol. 28, March 1980.
[98] I. H¨anninen, M. Taskinen, and J. Sarvas, “Singularity subtraction integral formulae
for surface integral equations with RWG, rooftop and hybrid basis functions,” Progress
In Electromagnetics Research, vol. 63, 2006.
[99] W. Cai, Y. Yu, and X. C. Yuan, “Singularity treatment and high-order RWG basis
functions for integral equations of electromagnetic scattering,” Int. J. Numer. Meth.
Engng., vol. 53, 2002.
¨ ur Erg¨
[100] L. G¨
urel and Ozg¨
ul, “Singularity of the magnetic-field integral equation and its
extraction,” IEEE Antennas and Wireless Propagation Letters, vol. 4, 2005.
[101] R. E. Hodges and Y. Rahmat-Samii, “The evaluation of MFIE integrals with the use
of vector triangle basis functions,” Microwave and Optical Technology Letters, vol. 14,
Jan. 1997.
[102] S. Caorsi, D. Moreno, and F. Sidoti, “Theoretical and numerical treatment of surface
integrals involving the free-space Green’s function,” IEEE Transactions on Antennas
and Propagation, vol. 41, Sept. 1993.
[103] M. G. Duffy, “Quadrature over a pyramid or cube of integrands with a singularity at
a vertex,” SIAM Journal on Numerical Analysis, vol. 19, Dec. 1982.
[104] P. Yl¨a-Oijala and M. Taskinen, “Calculation of CFIE impedance matrix elements
with RWG and n×RWG functions,” IEEE Transactions on Antennas and Propagation,
vol. 51, Aug. 2003.
[105] W. C. Chew, M. S. Tong, and B. Hu, Integral Equation Methods for Electromagnetic
and Elastic Waves. Morgan and Claypool Publishers, 2009.
[106] J. R. Mautz and R. F. Harrington, “An E-field solution for a conducting surface small
or comparable to the wavelength,” IEEE Trans. Antennas Propag., vol. 32, Apr. 1984.
115

[107] E. Arvas, R. F. Harrington, and J. R. Mautz, “Radiation and scattering from
electrically small conducting bodies of arbitrary shape,” IEEE Trans. Antennas
Propag., vol. 34, Jan. 1986.
[108] W. Wu, A. W. Glisson, and D. Kajfez, “A comparison of two low-frequency
formulations for the electric field integral equation,” in Tenth Ann. Rev. Prog. Appl.
Comput. Electromag., vol. 2, (Monterey, California), Mar. 1994.
[109] M. Burton and S. Kashyap, “A study of a recent, moment-method algorithm that is
accurate to very low frequencies,” Appl. Comput. Electromagn. Soc. J., vol. 10, Nov.
1995.
[110] F. Yuan, “Analysis of power/ground noise and decoupling capacitors in printed circuit
board systems,” in IEEE Electromagn. Compat. Symp. Proc., (Austin, TX), Aug. 1997.
[111] G. Vecchi, “Loop-star decomposition of basis functions in the discretization of EFIE,”
IEEE Trans. Antennas Propag., vol. 47, Feb. 1999.
[112] J. S. Zhao and W. C. Chew, “Integral equation solution of maxwell’s equations from
zero frequency to microwave frequencies,” IEEE Trans. Antennas Propag., James R.
Wait Memorial Special Issue, vol. 48, Oct. 2000.
[113] J. R. Bowler, S. A. Jenkins, L. D. Sabbagh, and H. A. Sabbagh, “Eddy-current probe
impedance due to a volumetric flaw,” J. Applied Physics, vol. 70, pp. 1107–1114,
August 1991.
[114] “The ACES collection of canonical problems: Set 1,” pp. 3–8, The Applied
Computational Electromagnetics Society, Spring 1990.
[115] S. K. Burke, “A benchmark problem for computation of δZ in eddy-current
nondestructive evaluation (NDE),” J. Nondestructive Evaluation, vol. 7, pp. 35–41,
1988.
[116] D. M. McKirdy J. Nondestructive Evaluation, vol. 8, pp. 45–51, 1989.
[117] W. S. Dunbar J. Nondestructive Evaluation, vol. 7, pp. 43–53, 1988.
[118] H. Xie, Y. Ji, and J. R. Bowler, “Eddy current pancake coil measurements on a
longitudinal through-wall notch in an inconel tube,” ftp://cnde.iastate.edu/.

116

[119] G. H. Golub and C. F. V. Loan, Matrix Computations. Baltimore, MD: The John
Hopkins Univ. Press, 1996.
[120] S. Kurz, O. Rain, and S. Rjasanow, “The adaptive cross-approximation technique for
the 3-D boundary-element method,” IEEE Trans. Magn., vol. 38, no. 2, 2002.
[121] K. Zhao, M. N. Vouvakis, and J.-F. Lee, “The adaptive cross approximation algorithm
for accelerated method of moments computations of EMC problems,” IEEE Trans.
Electromagn. Compat., vol. 47, no. 4, 2005.
[122] H. Cheng, Z. Gimbutas, P. G. Martisson, and V. Rokhlin, “On the compression of low
rank matrices,” SIAM J. SCI. COMPUT., vol. 26, no. 4, 2005.
[123] E. Liberty, F. Woolfe, P.-G. Martinsson, V. Rokhlin, and M. Tygert, “Randomized
algorithms for the low-rank approximation of matrices,” PNAS, vol. 104, Dec. 2007.
[124] F. Woolfe, E. Liberty, V. Rokhlin, and M. Tygert, “A fast randomized algorithm for
the approximation of matrices,” Appl. Comput. Harmon. Anal., vol. 25, 2008.
[125] W. C. Chew, J. M. Jin, E. Michielssen, and J. M. Song, eds., Fast and Efficient
Algorithms in Computational Electromagnetics. Norwood, MA: Artech House, 2011.
[126] T. K. Sarkar, E. Arvas, and S. M. Rao, “Application of FFT and the conjugate gradient
method for the solution of electromagnetic radiation from electrically large and small
conducting bodies,” IEEE Trans. Antennas Propagat., vol. 34, no. 5, 1986.
[127] T. J. Peters and J. L. Volakis, “Application of a conjugate gradient FFT method to
scattering from thin planar material plates,” IEEE Trans. Antennas Propagat., vol. 36,
no. 4, 1988.
[128] C. Y. Shen, K. J. Glover, M. I. Sancer, and A. D. Varvatsis, “The discrete Fourier
transform method of solving differential-integral equations in scattering theory,” IEEE
Trans. Antennas Propagat., vol. 37, no. 8, 1989.
[129] K. Barkeshli and J. L. Volakis, “On the implementation of the conjugate gradient
Fourier transform method for scattering by planar plates,” IEEE Antennas Propagat.
Mag., vol. 32, no. 2, 1990.
[130] M. F. Catedra, J. G. Cuevas, and L. Nuno, “A scheme to analyze conducting plates of
resonant size using the conjugate gradient method and fast Fourier transform,” IEEE
Trans. Antennas Propagat., vol. 36, no. 12, 1988.
117

[131] A. P. M. Zwamborn and P. M. van den Berg, “A weak form of the conjugate gradient
FFT method for plate problems,” IEEE Trans. Antennas Propagat., vol. 39, no. 2,
1991.
[132] J. L. V. J. M. Jin, “A biconjugate gradient FFT solution for scattering by planar
plates,” Electromagn., vol. 12, 1992.
[133] Y. Zhuang, K. Wu, C. Wu, and J. Litva, “A combined full-wave CG-FFT method for
rigorous analysis of large microstrip antenna array,” IEEE Trans. Antennas Propagat.,
vol. 44, no. 1, 1996.
[134] C. F. Wang, F. Ling, and J. M. Jin, “A fast full-wave analysis of scattering and
radiation from large finite arrays of microstrip antennas,” IEEE Trans. Antennas
Propagat., vol. 46, no. 10, 1998.
[135] D. T. Borup and O. P. Gandhi, “Fast-Fourier transform method for calculation of SAR
distributions in finely discretized inhomogeneous models of biological bodies,” IEEE
Trans. Microwave Theory Tech., vol. 32, no. 2, 1984.
[136] M. F. Catedra, E. Gago, and L. Nuno, “A numerical scheme to obtain the RCS of
three-dimensional bodies of resonant size using the conjugate gradient method and the
fast Fourier transform,” IEEE Trans. Antennas Propagat., vol. 37, no. 5, 1989.
[137] A. P. M. Zwamborn, P. M. van den Berg, J. Mooibroek, and F. T. C. Koenis,
“Computation of three-dimensional electromagnetic fields distributions in a human
body using the weak form of the CG-FFT method,” ACES J., vol. 7, 1992.
[138] A. P. M. Zwamborn and P. M. van der Berg, “Computation of electromagnetic
fields inside strongly inhomogeneous objects by the weak-conjugate-gradient
fast-Fourier-transform method,” J. Opt. Soc. Am. A, vol. 11, 1994.
[139] H. Gan and W. C. Chew, “A discrete BCG-FFT algorithm for solving 3D
inhomogeneous scattering problems,” J. Electromagn. Waves Appl., vol. 9, 1995.
[140] J. M. Jin, J. Chen, H. Gan, W. C. Chew, R. L. Magin, and P. J. Dimbylow,
“Computation of electromagnetic fields for high-frequency magnetic resonance imaging
applications,” Phys. Med. Biol., vol. 41, 1996.
[141] C. F. Wang and J. M. Jin, “Simple and efficient computation of electromagnetic fields
in arbitrarily-shaped, inhomogeneous dielectric bodies using transpose-free QMR and
FFT,” IEEE Trans. Microwave Theory Tech., vol. 46, no. 5, 1998.
118

[142] R. H.-F. Chan and X.-Q. Jin, An Introduction to Iterative Toeplitz Solvers.
Fundamentals of Algorithms, Society for Industrial and Applied Mathematics,
September 2007.
[143] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “AIM: Adaptive integral method
for solving large-scale electromagnetic scattering and radiation problems,” Radio Sci.,
vol. 31, 1996.
[144] F. Ling, C. F. Wang, and J. M. Jin, “Application of adaptive integral method
to scattering and radiation analysis of arbitrarily shaped planar structures,” J.
Electromagn. Waves Appl., vol. 12, 1998.
[145] S. S. Bindiganavale, J. L. Volakis, and H. Anastassiu, “Scattering from structures
containing small features using the adaptive integral method (AIM),” IEEE Trans.
Antennas Propagat., vol. 46, no. 12, 1998.
[146] F. Ling, C. F. Wang, and J. M. Jin, “An efficient algorithm for analyzing large-scale
microstrip structures using adaptive integral method combined with discrete complex
image method,” IEEE Trans. Microwave Theory Tech., vol. 48, no. 5, 2000.
[147] C. F. Wang, F. Ling, J. M. Song, and J. M. Jin, “Adaptive integral solution of combined
field integral equation,” Microw. Opt. Tech. Lett., vol. 19, no. 5, 1998.
[148] J. R. Phillips and J. K. White, “A precorrected-FFT method for electrostatic analysis
of complicated 3D structures,” IEEE Trans. Comput. Aided Des., vol. 16, no. 10, 1997.
[149] X. C. Nie, L. W. Li, N. Yuan, T. S. Yeo, and Y. B. Gan, “Precorrected-FFT solution of
the volume integral equation for 3-D inhomogeneous dielectric objects,” IEEE Trans.
Antennas Propagat., vol. 53, no. 1, 2005.
[150] C. H. Chan, C. M. Lin, L. Tsang, and Y. F. Leung, “A sparse-matrix/canonical grid
method for analyzing microstrip structures,” IEICE Trans. Electron., vol. E80C, no. 11,
1997.
[151] S. M. Seo and J. F. Lee, “A fast IE-FFT algorithm for solving PEC scattering
problems,” IEEE Trans. Magn., vol. 41, no. 5, 2005.
[152] K. C. Donepudi, J. Song, J.-M. Jin, G. Kang, and W. C. Chew, “A novel
implementation of multilevel fast multipole algorithm for higher order Galerkin’s
method,” IEEE Transactions on Antennas and Propagation, vol. 48, Aug. 2000.

119

[153] J. L. Volakis and K. Sertel, Integral Equation Methods for Electromagnetics. SciTech
Publishing, Feb. 2012.
[154] H. Wall´en, S. J¨arvenp¨a¨a, P. Yl¨a-Oijala, and J. Sarvas, “Broadband M¨
uller-MLFMA
for electromagnetic scattering by dielectric objects,” IEEE Transactions on Antennas
and Propagation, vol. 55, May 2007.
[155] C. C. Lu and W. C. Chew, “Fast algorithm for solving hybrid integral equations,” in
IEE. Proceedings-H, vol. 140, Dec. 1993.
[156] E. Darve, “The fast multipole method: Numerical implementation,” Journal of
Computational Physics, vol. 160, 2000.
[157] E. Darve, “The fast multipole method I: Error analysis and asymptotic complexity,”
SIAM J. Numer. Anal., vol. 38, no. 1, 2000.
[158] M. L. Hastriter, S. Ohnuki, and W. C. Chew, “Error control of the translation operator
in 3D MLFMA,” Microwave and Optical Technology Letters, vol. 37, May 2003.
[159] S. Koc, J. Song, and W. C. Chew, “Error analysis for the numerical evaluation of
the diagonal forms of the scalar spherical addition theorem,” SIAM J. Numer. Anal.,
vol. 36, no. 3, 1999.
[160] L. J. Jiang and W. C. Chew, “A mixed-form fast multipole algorithm,” IEEE
Transactions on Antennas and Propagation, vol. 53, Dec. 2005.
[161] J. S. Zhao and W. C. Chew, “Applying LF-MLFMA to solve complex PEC structures,”
Microwave and Optical Technology Letters, vol. 28, Feb. 2001.
[162] J. S. Zhao and W. C. Chew, “MLFMA for solving integral equations of 2-D
electromagnetic problems from static to electrodynamic,” Microwave and Optical
Technology Letters, vol. 20, Mar. 1999.
[163] J. S. Zhao and W. C. Chew, “Three-dimensional multilevel fast multipole algorithm
from static to electrodynamic,” Microwave and Optical Technology Letters, vol. 26,
July 2000.
[164] E. Darve and P. Have, “Efficient fast multipole method for low-frequency scattering,”
Journal of Computational Physics, vol. 197, 2004.

120

[165] H. Cheng, W. Y. Crutchfield, Z. Gimbutas, L. F. Greengard, J. F. Ethridge, J. Huang,
V. Rokhlin, N. Yarvin, and J. Zhao, “A wideband fast multipole method for the
Helmholtz equation in three dimensions,” Journal of Computational Physics, vol. 216,
2006.
[166] J. S. Zhao and W. C. Chew, “Applying matrix rotation to the three-dimensional
low-frequency multilevel fast multipole algorithm,” Microwave and Optical Technology
Letters, vol. 26, July 2000.
[167] L. Greengard, J. Huang, V. Rokhlin, and S. Wandzura, “Accelerating fast multipole
methods for the Helmholtz equation at low frequencies,” IEEE Computational Science
and Engineering, vol. 5, no. 3, 1998.
¨ ur Erg¨
[168] Ozg¨
ul and L. G¨
urel, “Enhancing the accuracy of the interpolations and
anterpolations in MLFMA,” IEEE Antennas and Wireless Propagation Letters, vol. 5,
2006.
[169] I. Chowdhury and V. Jandhyala, “Integration and interpolation based on fast
spherical transforms for the multilevel fast multipole method,” Microwave and Optical
Technology Letters, vol. 48, Oct. 2006.
[170] J. Song and W. C. Chew, “Interpolation of translation matrix in MLFMA,” Microwave
and Optical Technology Letters, vol. 30, July 2001.
¨ ur Erg¨
[171] Ozg¨
ul and L. G¨
urel, “On the lagrange interpolation in multilevel fast multipole
algorithm,” IEEE Antennas and Propagation Society International Symposium, 2006.
¨ ur Erg¨
[172] Ozg¨
ul and L. G¨
urel, “Optimal interpolation of translation operator in multilevel
fast multipole algorithm,” IEEE Transactions on Antennas and Propagation, vol. 54,
Dec. 2006.
[173] J. Sarvas, “Performing interpolation and anterpolation entirely by fast Fourier
transform in the 3-D multilevel fast multipole algorithm,” SIAM J. Numer. Anal.,
vol. 41, no. 6, 2003.
[174] S. Aluru, “Greengard’s N-body algorithm is not order N,” SIAM Journal on Scientific
Computing, vol. 17, May 1996.
[175] J. J. Sakurai and J. J. Napolitano, Modern Quantum Mechanics. Addison-Wesley,
2nd ed., July 2010.

121

[176] A. Gillman, P. M. Young, and P.-G. Martinsson, “A direct solver with O(N) complexity
for integral equations on one-dimensional domains,” in Frontiers of Mathematics in
China, vol. 7, April 2012.
[177] K. L. Ho and L. Greengard, “A fast direct solver for structured linear systems by
recursive skeletonization,” SIAM J. SCI. COMPUT., vol. 34, no. 5, 2012.
[178] L. Greengard, D. Gueyffier, P. Martinsson, and V. Rokhlin, “Fast direct solvers for
integral equations in complex three-dimensional domains,” Acta Numerica, vol. 18,
May 2009.
[179] Y. Chen, “A fast, direct algorithm for the Lippmann-Schwinger integral equation in
two dimensions,” Advances in Computational Mathematicsm, vol. 16, 2002.
[180] P. Martinsson and V. Rokhlin, “A fast direct solver for boundary integral equations in
two dimensions,” Journal of Computational Physics, vol. 205, 2005.
[181] J.-G. Wei, Z. Peng, and J.-F. Lee, “A fast direct matrix solver for surface integral
equation methods for electromagnetic wave scattering from non-penetrable targets,”
RADIO SCIENCE, vol. 47, 2012.
[182] X.-M. Pan, J.-G. Wei, Z. Peng, and X.-Q. Sheng, “A fast algorithm for multiscale
electromagnetic problems using interpolative decomposition and multilevel fast
multipole algorithm,” RADIO SCIENCE, vol. 47, 2012.
¨ ur Erg¨
[183] Ozg¨
ul and L. G¨
urel, “Improved testing of the magnetic-field integral equation,”
IEEE Antennas and Wireless Propagation Letters, vol. 15, October 2005.

122