s; a: 33.9.
E . AfiAxXLux
1| ‘q‘o 1. 31 ‘ ‘1 Y.‘
.. 6. r .3. .rn...
. .O)‘
,z.

‘ 7 . .
t. x. (it
.i ; .1.

. : ~\- in.)

«strunmtufirtih

.m. irhwhflw

45.".—
mmh .

3.3%....»th r.
L. :07. Q ‘10:}.
zmufi.phr fl ...
1.3. $552
14 : 2.,

 

 

:\
‘5‘ .2.
.« I

駧 gagxé

.fi‘ 5.. I. ‘I

 

11.31

.1005]

 

‘ LIBRARY
Michigan State
University

 

 

 

This is to certify that the
dissertation entitled

APPLICATION OF THE WAVEFIELD TRANSFORM TO
NONDESTRUCTIVE EVALUATION

presented by

YONG TIAN

has been accepted towards fulfillment
of the requirements for the

Doctoral degree in Electrical Engineering

 

 

MIL MN 4% Cw dart/(bu:

Major ProfessBr’s Signature

Nov-1 II; 4120!

Date

MSU is an Affirmative Action/Equal Opportunity Institution

PLACE IN RETURN Box to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

APPLICATIONS OF THE WAVEFIELD TRANSFORM TO
NONDESTRUCTIVE EVALUATION
By

Yong Tian

A DISSERTATION

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

DOCTOR OF PHILOSOPHY
Department of Electrical and Computer Engineering

2005

ABSTRACT

APPLICATIONS OF THE WAVEFIELD TRANSFORM TO
NONDESTRUCTIVE EVALUATION
By

Yong Tian

Eddy current nondestructive techniques offer many attractive benefits such as
reduced inspection time, low cost and reproducibility. Nevertheless, they are not used
in many industrial applications, primarily due to the difficulty associated with the
lack of simple and physically meaningful interpretation techniques. In contrast, wave
propagation phenomena based non-destructive evaluation (N DE) techniques employ
a host of physical intuitive concepts, among which a prominent one is time—of—flight
(TOF), i.e. the time duration between the excitation pulse and its response. TOF not
only provides position and other information relating to the flaw in an “explicit” way,
but also enables the use of mapping algorithms based on wave propagation. There is
thus a clear need to study means in which techniques for analyzing wave propagation
based NDE data can be applied to eddy current testing (ECT) data. This research
aims at inverting ECT data from the perspective of wave propagation phenomena
and presenting the inversion results in the same format as that obtained from wave
propagation based testings, thereby facilitating possible future data fusion processes.

Towards this goal, it is necessary to obtain a comprehensive understanding of the
distinct physical characteristics, incompatible test data formats, and different math-

ematical tools required for analysis and data processing. To this end, we employ a

wavefield transform, also called the Q-transforrn, a mapping relating wave fields to
diffusive fields, to retrieve TOF information from ECT data. In this research, the
TOF information is extracted directly from ECT data. Thus, we overcome the insta-
bility associated with numerical inversion of the Q-transform, which is widely adopted
in traditional Q-transforrn based TOF extraction methods. In order to demonstrate
the effectiveness of the proposed Q-transform approach, we estimate the TOP for
a host of canonical examples, in both the time and frequency domain. In addition
to demonstrating the merit of these models through numerical simulations, an ex-
perimental set-up was built for validating the concept. The measured data shows
excellent agreement with theoretical predictions. The experimental data was also
used to estimate the source position successfully.

The Q-transform based TOF extraction methods presented in this work demon-
strates the potential for retrieving the TOP information from test objects with com-
plex geometries. The extracted TOF data possesses the same format as that mea-
sured from wave propagation based NDE techniques. This makes it possible to fuse
these measurements together for improving inspection accuracy and reliability. In a
broader context, the successful development of the Q-transform approach may inspire
and encourage future research on methods for addressing ECT conductivity imaging

problems.

To my parents, sister and wife.

iv

ACKNOWLEDGMENTS

Foremost, I would like to express my gratitude to my major advisor, Dr. Satish
Udpa, for his consistent support in the past few years. He offered not only his keen
insight and guidance in my research area, but also his generous encouragement. I
have no doubt that I wouldn’t reach this point without his help. I would like to
thank Dr. Antonello Tamburrino, who also has been my major advisor, for his great
contribution to this study. His guidance and inspiration have been a key source of
motivation. I would like to thank all my committee members, Dr. Richard C. York,
Dr. Lalita Udpa, Dr. Shanker Balasubramaniam, Dr. Leo Kempel and Dr. Edward
Rothwell for taking the time to serve on my committee and for providing invaluable
comments and suggestions to enhance the quality of this dissertation.

I am grateful to my colleagues. Mr. Naveen Nair contributed to both my theoret-
ical and experimental work in the final stage of this work, Mr. Robert Clifford, Mr.
Michael Chan and Mr. Brian Wright provided significant support for my experimental

work.

TABLE OF CONTENTS

LIST OF TABLES ix
LIST OF FIGURES xiii
1 INTRODUCTION 1
1.1 Nondestructive Evaluation ........................ 1
1.2 Eddy Current Testing and Time—of-Flight ................ 3
1.3 Preliminaries ............................... 5
1.4 Literature Review ............................. 8
1.5 Organization of This W'ork ........................ 12

2 FUNDAMENTALS OF ELECTROMAGNETIC THEORY AND THE F I-

NITE ELEMENT METHOD 14
2.1 Basics of Electromagnetic Theory .................... 14
2.1.1 Maxwell’s Equations ....................... 15
2.1.2 Constitutive Relations ...................... 17
2.1.3 Boundary Conditions ....................... 17
2.1.4 Time-Harmonic Fields ...................... 18
2.1.5 Quasi-Static Approximation ................... 20
2.1.6 Potential Functions ........................ 22
2.2 Potential Formulations for Eddy Current Problems .......... 23
2.2.1 A — V Formulation ........................ 25
2.2.2 A Formulation .......................... 29

vi

2.3 Finite Element Method .......................... 30

2.3.1 Discretization ........................... 31
2.3.2 Galerkin’s Method ........................ 34
2.3.3 Solution Methods ......................... 37
THE Q-TRANSFORM 39
3.1 Derivations of the Scalar Q-transform .................. 39
3.2 Mathematical Properties of the Q-transform .............. 43
3.3 Numerical Implementation of the Scalar Q-transform ......... 47
3.3.1 Numerical Approximation of the Q-transform ......... 47
3.3.2 Truncation of the Infinite Integral ................ 50
3.4 Derivations of the Vector Q-transform .................. 52

TIME-OF -F LIGHT EXTRACTION FOR TRANSIENT DIFFUSIVE FIELDS 58
4.1 A Scalar Diffusion Equation Problem .................. 58

4.1.1 Time-of-F light and Peak Value of the Eddy Current Measurement 59

4.1.2 Analytical Solutions of the Scalar Problem ........... 62
4.1.3 Design of the Excitation Signal ................. 66
4.2 Numerical Simulation: A Scalar Problem ................ 68
4.2.1 Numerical Implementation .................... 68
4.2.2 Numerical Results ........................ 71
4.3 A Vector Diffusion Equation Problem .................. 76

TIME-OF-F LIGHT EXTRACTION FOR HARMONIC DIFFUSIVE FIELDS 82

5.1 Time-of—Flight in Harmonic Propagative Fields ............. 82

vii

5.2 A Cylindrical Configuration .......................
5.2.1 Solution of the Radiation Problem ...............
5.2.2 Time-of—Flight and the Solution of the Radiation Problem

5.3 Interface Removal and Time-of—Flight Estimation ...........
5.3.1 Interface Removal Using Linear Filtering ............
5.3.2 Time-of—Flight Estimation in the Frequency Domain ......

5.3.3 Numerical Simulations ......................

6 EXPERIMENTAL VERIFICATION S
6.1 Interface Removal for a Planar Configuration ..............
6.2 Experimental Set-up ...........................

6.3 Experiment Results ............................
7 CONCLUSIONS

A Q-TRANSFORM PAIRS AND THE Q-TRANSFORM
A.1 Q-transform Pairs .............................

A.2 The Q-transform .............................

BIBLIOGRAPHY

viii

84

92

94

94

96

101

101

103

106

118

121

121

122

125

3.1

LIST OF TABLES

Maximum support for numerical evaluation of the Q-transform . . . .

ix

52

1.1

1.2

1.3

2.1

2.2

3.1

3.2

4.1

4.2

4.3

4.4

4.5

4.6

4.7

LIST OF FIGURES

A typical N DE system .......................... 2
Data fusion approach using the Q-transform .............. 9
Alternative data fusion approach using the Q-transform ........ 11
A typical solution domain of an eddy current problem ........ 23
A triangular element ........................... 32
Weighting function of the Q-transform in linear scale. ........ 57
Weighting function of the Q—transform in logarithmic scale ...... 57
A small anomaly in the vicinity of a point source ........... 59

The relation between time-of-fiight go and peak position tmax of the
eddy current measurement ....................... 61
Time domain and fictitious time domain excitation signals as given by
Eq. (4.23) and (4.24), respectively (q,- equals zero in all cases). . . . 66
Left: A schematic of the geometry for the numerical simulation; Center:
The finite element mesh generated by FEMLAB®; Right: Expanded
view of the mesh in the region containing the anomaly. ........ 70
The plot of "Us (0, t) (solid) together with the plot of Q {us (0, t)} (*).
The predicted peak position is at t = 0.009 s ............. 72
The plot of vs (0, t) (*) together with the plot of the approximate re-
sponse (solid). case n = 1 ........................ 73
The plot of v, (0, t) (*) together with the plot of the approximate re-

sponse (solid). case n = 2 ........................ 74

4.8

4.9

4.10

4.11

5.1

5.2

5.3

5.4

5.5

Shifting effect of truncation on excitation signal, q, = 0.15. Both
ideal response 1),; (0, t) (solid) and shifted response 2); (0,t) (dashed)
are normalized. ..............................
Truncation error. Left: Relative error with respect to peak value, 72. = 1
(solid line) and n = 2 (dashed line); Right: Relative error with respect
to peak position, 72 = 1 (solid line) and n = 2 (dashed line) .....
Error in the predicted position (*) and amplitude (o) of the peak as
function of 7 for n = 1 (dashed line) and n = 2 (solid line), respectively.
Case X = —0.25. .............................
Error in the predicted position (*) and amplitude (o) of the peak as

function of 7 for n = 1 (dashed line) and n = 2 (solid line), respectively.

Two time-of—flight terms included in the Fourier coefficients .....
Geometry for the reference problem. Left: 3D View; Right: 2D Cross-
section at z = 2:0 .............................
Geometry for the Radiation Problem ..................
Finite element mesh in a neighborhood of the cylinder. Additional
artificial boundaries have been enforced to adjust the mesh density

V(r,6n; jw), Upper: real part; Lower: imaginary part. (Solid line:
analytical solution given by Eq. (5.57); ’ X’: High frequency approxi-

mation given by Eq. (5.58); ’0': Numerical solution constructed with

xi

77

84

85

88

99

6.1

6.2

6.3

6.4

6.5

6.6

6.7

6.8

6.9

6.10

6.11

6.12

6.13

6.14

6.15

6.16

6.17

6.18

6.19

Geometry of the planar configuration for interface removal. Left: in

presence of the plate; Right: inside an unbounded homogeneous medium

101

Experimental Set-up ........................... 103
Real part of Byp (GMR sensor) ..................... 108
Imaginary part of Byp (GMR sensor) ................ .. . 108
Phase of Byp (GMR sensor) ....................... 109
Magnitude of Byp (GMR sensor) .................... 109
Complex plane plot of Byp (GMR sensor) ............... 110
Real part of Byp (Inductive coil) .................... 110
Imaginary part of BW (Inductive coil) ................. 111
Phase of 8,”, (Inductive coil) ...................... 111
Magnitude of BW (Inductive coil) .................... 112
Complex plane plot of BW (Inductive coil) ............... 112
Real part of Byh (x, yh) (filtered from GMR data) ........... 113
Imaginary part of Byh (2:, yh) (filtered from GMR data) ........ 113
Complex plane plot of By}, (1:, yh) (filtered from GMR data) ..... 114
Real part of By}, (T, yh)(filtered from inductive coil data) ....... 114
Imaginary part of By}, (1:, yh) (filtered from inductive coil data) . . . 115

Complex plane plot of Byh (51:, yh) (filtered from inductive coil data) . 115
Error function used for recovering the source location (—a, a), unit: m,

(contour plot, ‘GMR I’) ......................... 116

xii

6.20 Error function used for recovering the source location (—a, a), unit: m,
(mesh plot, ‘GMR I’) .......................... 116
6.21 Error function used for recovering the source location (—a, a), unit: m,
(contour plot, ‘Coil I’) .......................... 117
6.22 Error function used for recovering the source location (—a, a), unit: m,

(mesh plot, ‘Coil I’) ........................... 117

xiii

CHAPTER 1

INTRODUCTION

1.1 Nondestructive Evaluation

NDE techniques are used very widely in a number of industrial areas to control
product quality, identify component failures, monitor manufacturing processes, etc.
A key aspect of NDE lies in the fact that inspections are performed in a manner
that does not affect the future usefulness of the object or material. As such, NDE
provides balance between quality control and cost effectiveness. Some of the major
applications can be found in aerospace industry, gas transmission pipeline industry
and nuclear power plants, where the detection of defects prior to device failure often
translates to enormous savings in terms of both economical costs and human lives.

A NDE technique usually introduces some form of physical energy, such as X-
rays, elastic waves, or electromagnetic fields, into the test specimen. The nature of
the interaction between the energy and the test specimen is a function of several
variables, including the energy type, material properties, defects and inhomogeneities
in the material and so on. The interaction is sampled through a transducer and
the response of the transducer is analyzed to characterize the material properties
of the specimen. As shown in Figure 1.1, a typical NDE system employs an input
transducer with an excitation source to expose the test specimen to the physical field.
The response of the test object is then sampled by a measurement transducer and
passed on to a computing unit for signal/ image processing, recognition, and fusion.

Common NDT techniques include radiography testing (RT), ultrasonic testing (UT),

 

Excitation
Source

 

 

 

 

 

Signal/Image Display
Recognition {Storage

 

 

 

 

 

 

 

 

 

U Input Transducer

Specimen v
\-._a-/

 

 

Measurement D

Transducer Signal/Image

Processing

 

 

 

 

 

Figure 1.1: A typical N DE system

electromagnetic testing (ET) and acoustic emission testing (AB) [1].

Since a single technique may not be capable of providing all the necessary NDE
information, the idea of judiciously combining measurements from multiple NDE
techniques, or from the same kind of sensors with different parameters has been
proposed as a means for garnering additional information [2]. We can benefit from
the enlarged pool of complementary information regarding the test specimen, thus
offering great improvement in the detectability and reliability of NDE inspections. For
example, in ECT testing, fields from eddy current probes are limited to the surface of
a conductive specimen. For a given set of material parameters, the depth of energy
penetration is inversely proportional to the square root of the excitation frequency
used. We employ low excitation frequencies for detecting deep flaws. Shallow flaws
can be detected with greater sensitivity by using a high excitation frequency. If

shallow as well as deep flaws are to be detected, it may be necessary to employ a wide

range of frequencies. The complexity and cost of such multi-mode NDE inspections

can be substantially reduced with advances in automated data fusion tools [3].

1.2 Eddy Current Testing and Time-of-Flight

ECT is an important NDE technique employing low—frequency electromagnetic
fields where the magnetoquasi-static approximation holds, i.e., the displacement cur-
rent is negligible. Therefore, the governing Maxwell’s equations are in the form of
vector diffusion equations.

ECT applications find widespread adoption in industry. ECT provides nearly
instantaneous measurements with great reproducibility for metallic materials. More—
over, modern ECT offers a low-cost method for high-speed and large—scale inspections,
even in extreme environmental situations, such as under high pressure and in high
temperature conditions.

The basic idea underlying ECT is to use a time varying magnetic field (excita-
tion source) to induce electrical currents in the conductive test sample. The induced
currents, known as eddy currents, generate a time varying magnetic field which in-
teracts with the original field. According to Lenz’s law, the magnetic field generated
by eddy currents always opposes a change in the net magnetic flux. The intensity of
the eddy current reaches its largest value at the surface of the specimen and decays
exponentially with the depth below the surface, giving rise to the so—called skin effect.
This effect constrains the usage of ECT to surface or subsurface inspections. Eddy
current techniques are sensitive to magnetic permeability and geometric properties of

the specimen.

The excitation sources (typically coils) associated with conventional ECT systems
are time harmonic in nature. The measurement employs sensors which are usually
either coils or magnetic field transducers. In swept-frequency eddy current (SFEC)
testing, a number of different frequencies are used to excite the probing coil. Since the
diffusion of the eddy currents into metals is governed by the skin effect, the responses
at different frequencies reflect the internal structure at different depths. Therefore,
SFEC can be used to obtain the “depth profile” of the specimen. Another prominent
form of ECT is pulsed eddy current (PEC) testing, which is a time-domain method
that uses a step-function as the excitation signal. Compared with other methods, the
PEC method offers the potential advantages of greater penetration and also serves as

a ready means of multi—frequency measurement.

When wave propagation NDE methods are used, the term time-of-flight (TOF)
refers to the time interval between the transmission of an excitation pulse and the
arrival of its echo. In ultrasonic testing, for example, the TOP is used as a basis for
locating the scatterer. Based on the knowledge of TOF, it is possible to determine the
location of the defects and generate reflectivity images about the subsurface struc-
ture of the specimen. The TOF measuring techniques also offer intuitive physical
interpretation of the measured data and are simple to implement.

Unfortunately, TOF measurements are not. physically meaningful in the case of
methods that generate a diffusive field, such as ECT [4], since no definition of pulse
travel~time and no concepts of wave propagation exist therein. Alternatively, the size

and location information of the defects are often obtained through complex decision-

making algorithms (such as so—called model-based metht'xls) for ECT measurements.
For the same reason, well-developed ultrasonic imaging methods cannot be employed
for analyzing ECT data. This suggests that one stands to benefit by exploring the
possibility of extracting information such as TOF that are traditionally associated

with wave propagation based NDE methods from diffusion based (such as eddy cur-

rent) NDE methods.

1.3 Preliminaries

This work aims at developing time-of-flight extraction methods from eddy current
measurements. The fundamental mathematical tool exploited in this research is the
wavefield transform, which is an integral operator that is capable of transforming
solutions of wave equation problems to solutions of diffusion equation problems. The
wavefield transform [5] is also called the Q-transform [6], or diffusion-to-wave trans-
form [7]. For the sake of simplicity, we choose to call it the Q-transform throughout
this work.

Suppose I} (x, t) is the solution of the diffusion equation [8, 9]:

V21) (x,t) — c(x)6tv (x,t) = F (x, t) in 9, fort Z 0 (1.1)
a (x) v (x, t) + B (x) 3"?) (x, t) = G (x,t) on 89, for t 2 0 (1.2)
v (x,0) = h (x) in R (1.3)

and u (x, q) is the solution of the wave equation:

V221 (X. (1) - C(X) 3qu (MI) = f (X. a) in 9. for q 2 0 (1-4)

a (X) H (x, 61) + B (X) 8w (x, q) = 9 (x. q) on (99, for q 2 0 (1-5)
u (x,0) = 0 in e (1.6)

aqu (x,0) = h (x) in e, (1.7)

where x is a spatial variable, I) is the domain of RN (N 2 1) with a smooth boundary
852, c (x) is a real and nonnegative function depending on the material properties, 8,,
denotes the normal derivative on ('99, and oz and B are two known scalar functions.

The Q-transform operator, denoted by Q, is defined as [6, 8, 9, 10]:

2

Q: u(x,q) —+ v(x,t) E him/0 00qexp (—%) u(x,q)dq. (1.8)

 

One can show that

v (x,t) = qu 09(1)} (0 (19)

provided that

F(Xat) = Q{f (x, (1)} (t) and G(x,t) = Q{9(X,q)} (t)-

It is evident from (1.8) that the Q-transform serves a bridge between the wave phe-
nomenon u (x, q) and the diffusion phenomenon v (x, t). Taking Fourier transform

on both sides of Eq. (1.8) [11], we obtained the Q-transform in the frequency (w)

domain:

fQ:u(x,q)—>V(w)Ev/O+ooexp (—\/j—wfi)u(x,q)dq (1.10)

With the help of Q-transform, we will demonstrate that the time-of—fiight can be
calculated from diffusive fields. In another words, measurements of the diffusive field
established by the eddy current probe can be interpreted from the perspective of wave
propagation phenomena.

The major benefits we stand to accrue include:
0 Interpretation of ECT data in a clear and intuitive manner.

0 Ability to estimate TOF from diffusive field measurements that correspond to
those that are typically obtained from wave propagation based NDE techniques,

for data fusion.

0 Ability to apply a wide variety of image reconstruction techniques that are

typically applied to wave propagation measurements.

In this work, the basic strategy of extracting the time-of-flight from transient ECT

measurements is described:

0 With respect to every diffusive field 2) (x, t) of interest, a fictitious wave field

it (x, q) is established analytically using the Q-transform.

0 Although the time-of-flight term appears explicitly in the fictitious wave field,
we wish to associate the TOF with certain distinct features of the time do—
main response. To this end, manipulations of the form of excitation signal are

carefully investigated.

o Parametric forms of diffusive field responses containing the TOF information,
are derived by transforming the fictitious wave response using Q-transform,
which, unlike inverse Q-transform, is not ill-posed. Extraction methods are

then developed to recover the TOF that is treated as a parameter.

In case of harmonic fields, material discontinuities could significantly increase the
complexity of the TOF retrieval process. To overcome this difficulty, a two—step

procedure is followed:

0 A linear filter is constructed to transform the field in presence of material dis-
continuities to a field excited in an infinite homogeneous background medium.

This filter is typically independent of the nature of the excitation source.

0 Based on knowledge of the field inside a homogeneous background, which is
much easier to handle with than fields resulting from the presence of mater-
ial discontinuities, we describe a technique for retriving T OF information. The
TOF information can be fused or combined with equivalent estimates from wave
propagation based NDE methods such as ultrasoninc of microwave nondestruc-
tive testing techniques to improve the TOF estimate. This is an example of a

phenomenological approach to data fusion.

1.4 Literature Review

NDE techniques are governed by different partial differential equations. It is
important to recognize the mathematical connections among these partial differential

equations [10]. In particular, the linkage between hyperbolic and parabolic equations,

UT Data

 

 

IL

Data Fusion >

Fused Data
:> Inverse fi

Q-transform

 

 

 

 

 

 

 

 

 

 

ECT Data

Figure 1.2: Data fusion approach using the Q-transform

i.e. wave propagation and diffusion equations, is of special interest here. Bragg and
Dettman [8, 9, 12, 13] were among the earliest to build these connections through
the Laplace transform. Other mathematical works addressing this issue were pursued
by Reznitskaya [14], Lavrent’ev [15] and Romanov [10]. Specifically, they developed
an integral transformation, that are called Q-transform later [6], that are capable of
mapping solutions of hyperbolic problems into solutions of parabolic problems.
Some authors, encouraged by developments in ultrasonic imaging techniques, such
as holographic methods and diffraction tomography methods, have attempted to ex-
tend these imaging methods to eddy current testing [17, 18, 19, 20, 21]. A systematic
way to treat diffusive measurements by using wave propagation based imaging meth-
ods has been made possible through the use of the Q-transform. The underlying
idea is to invert the Q-Transform (this requires the solution of a first kind Fredholm
integral equation) to cast eddy current data into the framework of a wave propa-
gation problem. The philosophy of this approach is illustrated in Figure 1.2. Such

attempts were made for the purpose of analyzing petroleum reservoir exploration

data [22] and magnetotelluric soundings over a stratified earth [23, 24]. The two last
references do not make use of the Q-Transform directly but refer to the underlying
connection between wave propagation and diffusion equations. The primary obstacle
for these methods arises from the ill—posedness of the inversion of the Q-Tfansform
(the Q-transform is a compact operator). This ill-posed problem has been addressed
[7, 25, 26] along with a discussion on regularization algorithms.

The Q-transform was extended to vector diffusion problems by Lee et a1. [6],
where the term Q-transform was coined. Details about how to implement the Q-
transform for electromagnetic imaging in conductive media were elaborated by Ger-
shenson [27, 28]. Experimental work was carried out to validate the time domain
(diffusive fields) to fictitious wave domain (fictitious wave fields) transform and in-
vestigate the constraints on bandwidth and signal-to-noise ratio [29].

Among recent developments, a nonlinear tomographic inversion scheme was ap-
plied to reconstruct the two-dimensional spatial distribution of the conductivity start-
ing from the TOF information associated to the fictitious wave propagation data,
using the inverse Q-Transform on eddy current measurements [5].

In [26], Gibert and others approximate the frequency domain eddy current re-
sponse through the Q-Transform. Specifically, they postulate that the response for
the fictitious field consists mainly of weighted and delayed replica of the driving (fic-
titious) waveform. Each term in this sequence can be regarded as a reflection, or
multiple reflection, from certain discontinuities. Then, they obtain the approximate
frequency domain eddy current response by applying the Q-transform to the ficti-

tious wave propagation response. The values of TOF’s and weights involved in the

10

 

UT Data

:[> Q—transform fl

 

 

 

 

 

 

 

 

Data Fusion >

jf Fused Data

 

 

 

ECT Data

Figure 1.3: Alternative data fusion approach using the Q-transform

approximated response can be estimated numerically from ECT measurements. Levy
et al. [30] indicate that this method works well if the desired response is a spike train
signal, which is often the case for wave reflection in layered structures.

An alternative is to represent the fictitious wave field either with a ray series or
its leading term [31]. The determination of the corresponding coefficients and TOF
can be very difficult. When only the leading term is used, the complexity of this
parameter estimation problem is reduced but multiple reflections cannot be taken
into account.

The Q-transform has also been used in another way to transform ultrasonic mea-
surements into equivalent ECT measurements. This enables one to transform ultra-
sonic measurements into the ECT format and, therefore, to fuse ultrasonic and ECT
measurements from the ECT standpoint as shown in Figure 1.3. A relevant problem
that has to be addressed to successfully fuse the data is the so—called registration
problem, i.e. a time-shift affecting the ultrasonic data may significantly affect its

Q-transform [16].

11

Several papers summarizing the work done to date have been published. TOF ex—
traction for time-domain fields was proposed in [32, 33]. The corresponding numerical
analyses of scalar and vector cases can be found in [34] and [35]. The treatment of
interfaces has been developed in [36] and its numerical validation has been presented

in [37].

1.5 Organization of This Work

This dissertation is organized into seven chapters.

Chapter 1, introduces the concept of nondestructive evaluation. The physics un-
derlying ECT techniques and diffusion phenomena is discussed. After a brief descrip-
tion of the time-of-flight and Q-transform concepts, the basic strategy proposed for
solving the ECT inverse problem (TOF extraction) using mathematical connections
between diffusion and wave propagation phenomena are presented. A short literature
review is also included in the chapter.

Chapter 2, summarizes the basics of electromagnetic fields. The eddy current
formulation based on magnetic vector potential A and electric scalar potential V is
presented. The fundamentals of the finite element method are also described.

In Chapter 3, the definition of the Q-transform is introduced through detailed
derivations for both the scalar and vector cases. Important mathematical proper-
ties are presented and numerical implementations for computing the transform are
discussed.

In Chapter 4, the constraints associated with (fictitious) wave front recovery in

the time domain are introduced. These constraints are then taken into consideration

12

for designing appropriate excitation waveforn'is for eddy current testing. The TOF
information is extracted afterwards by measuring special features of the received
signals, such as the position of the peak value. This strategy is demonstrated through
an example involving localization of a small anomaly embedded in an unbounded
homogeneous medium. The scalar diffusive case is investigated first for its theoretical
simplicity and significance. The investigation strategy is then extended to the vector
diffusive case. The approach is validated using data from numerical models simulating
some simple test geometries. The results obtained reveal good agreement with the
theoretical predictions.

In Chapter 5, the problem of material interfaces is treated. A cylindrical structure
is presented as an illustrative example. To remove the effect of multi-reflection on the
material interfaces, a stable linear filtering method is developed. The time-of-flight is
then extracted by solving a one-dimensional single—parameter minimization problem.

In Chapter 6, the idea of interface removal is adapted to a planar configuration.
The corresponding experimental set-up is built and tested. The experimental mea-
surements obtained agree with numerical as well as analytical predictions very well.
The measured data is then used successfully to recover the position of point sources.

In Chapter 7, a brief summary is presented.

13

CHAPTER 2

FUNDAMENTALS OF ELECTROMAGNETIC THEORY

AND THE FINITE ELEMENT METHOD

This chapter reviews the basics of electromagnetic theory, including Maxwell’s
equations as well as a short discussion covering constitutive relations and boundary
conditions. Some of the popular potential formulations for eddy current problems
are then discussed. The last section focuses on a short introduction to one of the
most important numerical techniques, i.e. finite element method (FEM). More de-
tailed discussion on electromagnetic theory, finite element method and computational

electromagnetics can be found in references such as [38, 39].

2.1 Basics of Electromagnetic Theory

The problem of electromagnetic analysis on a macroscopic level can be addressed
by solving Maxwell’s equations, subject to certain boundary and interface conditions.
Various scalar and vector potentials can be introduced to reformulate the Maxwell’s

equation into preferred forms depending on specific applications.

14

2.1.] M axwell ’3 Equations

For time-varying electromagnetic (EM) fields, the differential form of Maxwell’s

equations is given by

VXE(x,t) = flag—fl (2.1)
VxH(x,t) = Jam-625i) (2.2)
V-D(x,t) = p(x,t) (2.3)
V-B(x,t) = o, (2.4)

where E (volts/ meter) is the electric field intensity, H (amperes / meter) is the mag-
netic field intensity, D (coulombs/meterz) is the electric flux density, B (Tesla or
webers/meterz) is the magnetic flux density, J (amperes/meter2) is the electric cur-
rent density, p (coulombs/meter3) is the electric charge density and x = x (x, y, z) is
the spatial position. Here, Eq. (2.1) and (2.2) are referred to as Faraday’s law and
Maxwell-Ampere’s law, respectively. At the same time, Eq. (2.3) and (2.4) are two
different forms of Gauss’s law, i.e. the electric and magnetic form, respectively.

Another fundamental equation, the equation of continuity, can be derived from

Eqs. (2.2) and (2.3), which is

 

= 0. (2.5)

It is indeed a mathematical form of the law of the conservation of charges, which
means that the net flow of electric current out of a small volume equals the time rate

of decrease in electric charges.

15

Among the above five fundamental Eqs. (2.1)-(2.5), only three of them are inde—
pendent. Eq. (2.1)-(2.2) combined with either Eq. (2.3) or (2.5) and proper initial
conditions can be used to form an independent system, while others can be derived
from this system.

Although we will utilize the differential form of Maxwell’s equations in this re-
search, the equivalent integral form of Eqs. (2.1)-(2.4) is worth presenting here for
sake of their significance in both analytical and numerical treatments. The equations

are

jfE(x,t)-d1 = —/S§§(—§—t’it—)-ds (2.6)

C
. _ x 8D (x, t) .

CH(x, t) d1 —— [S[J( ,t) + —————at dS (2.7)
éD (x,t) - d3 = /Vp(x,t) dv (2.8)
fB (x, t) . d8 = 0. (2.9)

S

Equation (2.6), Faraday’s integral law, states that the circulation of the electric field
E around a contour C is determined by the time rate of change of the magnetic flux
linking the surface S enclosed by contour C. Equation (2.7), Ampere’s integral law
requires that the circulation of the magnetic field intensity H around a closed contour
C is equal to the net current passing through the surface S spanning the contour plus
the time rate of change of the net displacement flux density dD/O‘t through the
surface. Equation (2.8) and (2.9) are called Gauss’ integral law of electric flux and
magnetic flux, respectively. The first one indicates that the net electric flux crossing

a surface S that encloses a volume V is equal to the charge contained in this volume.

16

The second one states that the total magnetic flux passing through a closed surface

is zero, i.e. the magnetic flux is conservative.

2. 1. 2 Constitutive Relations

Assume the medium is isotropic, linear and nondispersive, the constitutive rela-
tions describing the material properties of the medium are usually combined with

Maxwell’s equations. They are

D (x, t) = e (x) E (x, t) (2.10)
B (x, t) = a (x) H (x, t) (2.11)
J (x, t) = a (x) E (x, t) , (2.12)

where e ( f arads / meter) is the electric permittivity of a medium, ,u (henrys/ meter)
is the magnetic permeability and o (siemens/ meter) is the electrical conductivity.
In a more general case, i.e. the medium may be anisotropic and nonlinear, these

property parameters are tensors while their values depend on the field values.

2.1.3 Boundary Conditions

Maxwell’s equations should be satisfied everywhere in space. However, certain
boundary conditions are necessary to determine a unique solution for a specific prob—
lem.

At interfaces between two media, the following conditions can be derived from the

17

integral form of Maxwell’s equations:

fix(E1—E2) = o (2.13)
fi-(Di—D2) = 93 (2-14)
fix(H1—H2) —_— J3 (2.15)
fi-(Bl—Bg) = 0, (2.16)

where fl is defined as the unit vector normal to the interface and pointing from
medium 2 to medium 1. This implies that only the tangential component of the
electric field intensity E and the normal component of the magnetic flux density B are
always continuous on inter-material boundaries. The discontinuity of the tangential
component of H and the normal component of D are characterized by the surface
current density J 3 and the surface charge density p3. Meanwhile, only two of these
four boundary conditions are independent. One of Eqs. (2.13) and (2.16), together
with one of Eqs. (2.14) and (2.15), form a set of two independent conditions.

In addition, the interface condition for the current density is

1:1.(31 — J2) = —3p3/Bt. (2.17)

2.1.4 Time-Harmonic Fields

In many engineering applications involving only linear materials, it is sufficient to
deal with only the steady-state solutions of the EM fields when produced by sinusoidal

currents [40, 41]. Such fields are said to be sinusoidally time-varying or time harmonic,

18

i.e. they vary at a sinusoidal frequency w. An arbitrary sinusoidal field F (x, t) can

be expressed as

F (x, t) = Re [Fp (x;jw) exp (jwt)] , (2.18)

where Fp (x; jw) is the phasor form of F (x, t) and is in general a complex value, Re []
indicates “taking the real part of” the quantity within brackets, and w is the angular
frequency (in rad/s) of the sinusoidal excitation.

Using the phasor representations of the EM quantities and replacing the time

derivatives o/at by jw, the Maxwell’s equations, in the harmonic case, become

V X Ep (X001) = —jWBp (X;jw) (2-19)
V X Hp (X;jw) = Jp (XU'W) +jWDp (X;jw) (220)
V : Dp (mica) = Pp (X;jw) (2-21)
V - Bp (x;jw) = 0. (2.22)

Also, interface conditions (2.13)-(2.16) are converted to:

a x (Elp — Egp) = o (2.23)
a (D1,, — Dgp) = p5,, (2.24)
a x (Hlp — ng) = Jsp (2.25)
a (131,, — 32p) = o. (2.26)

In the above equations, the time-space dependent Maxwell’s equations are reduced

to space dependent equations and the differential operations are reduced to alge-

19

braic operations. Thus, time—harmonic analysis could simplify the solution process
significantly.

Furthermore, considering the fact that w is one element of a whole frequency
spectrum, we can represent a nonsinusoidal field in terms of its time-harmonic com—

ponents, i.e.

F (x,t) = /00 Fp (x,w) ejmdw. (2.27)

--00
Therefore, if a time-harmonic field is known for any frequency w, the corresponding
nonsinusoidal field can be obtained through Fourier analysis. In later sections, we

will drop the subscript p when no misunderstanding could occur.

2.1.5 Quasi-Static Approximation

For eddy current NDE, the so—called magneto—quasistatic (MQS) approximation
holds. Thus the general electromagnetic wave propagation phenomena represented
by the original Maxwell’s equations can be simplified.

As can be seen from Eq. (2.2), there exist two kinds of currents: conduction
current and displacement current. The conduction current .1 is proportional to the
electric field intensity, as stated by Ohm’s law (2.12). The displacement current J d is

defined as the time varying rate of electric flux density, i.e.

Jd (x, t) E 219%2 (2.28)

Under the assumption that the time varying rate is low enough such that a >> we,

the displacement current can be neglected and the EM field can be obtained by

20

considering stationary current at every instant. This fact implies that the equation

of continuity can be rewritten as

V - J (x, t) = 0 (2.29)

and the time derivative of the electric displacement 8D (x, t) /8t can be disregarded

in Maxwell-Ampere’s law to yield

V x H (x,t) = J (x, t). (2.30)

By neglecting the displacement current, the Maxwell’s equations and corresponding
constitutive relations can be divided into two sub—systems of equations, i.e. one

representing an electrostatic field system

V - D (x,t) = p(x, t) (2.31)

D (x, t) = e (x) E (x, t) (2.32)

and the other a magnetodynamic field system [42]

V x H (x,t) = J (x, t) (2.33)
V - B (x, t) = 0 (2.34)
V x E (x, t) = __8§é7xfl (2.35)

21

Eflmflzpkflflmfl can

J (x, t) = a (x) E (x, t). (2.37)

2.1.6 Potential Functions

It is often helpful to formulate the problems using potential functions. The fol-
lowing two vector identities are frequently applied to define the potential functions.
They are

VxVV=0 (2%)

and

V-VxA=0. (2.39)

These equalities state that the gradient of an arbitrary scalar function is irrotational
and the curl of an arbitrary vector function is solenoidal, provided these functions
are sufficiently differentiable. Two useful potential functions are the magnetic vector

potential A and the electric scalar potential V. They are given by the equalities

B = V x A (2.40)
8A
E — ——a—t — VV, (2.41)

which are direct consequences of the magnetic form of Gauss’s law and Faraday’s law,

respectively.

22

 

Figure 2.1: A typical solution domain of an eddy current problem

2.2 Potential Formulations for Eddy Current Problems

For eddy current problems, the quasi—static approximation is valid and potential
formulations are often introduced to simplify the governing Maxwell’s equations and
minimize computational costs. For a general three-dimensional eddy current problem
shown in Figure 2.1, the solution domain Q can be divided into a conducting region
{21 and a nonconducting region {22, which may include the excitation current source
J 3. On the outside boundary SH, the tangential component of the magnetic field
is zero while on the other part of the boundary, S 3, the normal component of the
magnetic flux density is zero. Also, we have 312 as interface between the conducting
region and the nonconducting region.

The Maxwell’s equations and corresponding boundary conditions for eddy current

23

problems can be written as:

VtzoE

X at
VXHZJS
V-B=0
B-fi=0
foi=0
fi-(Bl—B2)=0

le(H1—H2)_—"0

in $21
in $21
in 91
III 92
III 92
on 53
on SH
on 312

on 512.

(2.42)
(2.43)
(2.44)
(2.45)
(2.46)
(2.47)
(2.48)

(2.49)

where fl is the outer normal on the corresponding surface, and the subscript 1 and 2

refer to quantities in the region 91 and fig, respectively.

Providing boundaries SB and S H are simply connected, the uniqueness of B and

E as the solutions of these equations can be proved [43]. There are a number of

potential based methods aimed at reformulating the basic Maxwell’s equations. Each

of them is associated with certain advantages and disadvantages [44]. In this section,

Eqs. (2.42)-(2.50) will be reformulated in terms of magnetic vector potential A and

electric scalar potential V.

24

2.2.1 A — V Formulation

In A— V formulation, the magnetic vector potential A and electric scalar potential
V are used to represent the electromagnetic field in the conducting region S21, while
in the nonconducting region 92, only A is used. This method is also referred as
A — V — A formulation [45].

Invoking Eq. (2.40) and (2.41), Eqs. (2.42) - (2.50) can be expressed as

1 '9
VX—VXA+0£:IA:+0VV=0 infll (2.51)
a dt
V x (1V x A) 2 JS in 02 (2.52)
a
fi-VXA=0 onSB (2.53)
I
fiVxAxfi=0 onSH (2.54)
fi.(VxA1—VXA2)=O on512 (2.55)
1 I
(—V x A1— —V x A2) x n = 0 on 812. (2.56)
#1 #2

Notice that Eq. (2.43), (2.44) and (2.46) are automatic satisfied. Eq. (2.55) is
satisfied when A is continuous across the interface 812 [45].

Equation (2.40) specifies the curl of the magnetic vector potential as the magnetic
flux density. Helmholtz theorem states that a vector field is uniquely determined when
both its curl and divergence are specified. Therefore, the divergence of the magnetic
vector potential must be specified to fix the additional degrees of freedom associated
with A (called gauging). This value may be specified freely without affecting the

physical problem. In many circumstances, the divergence free condition is imposed,

25

which is

V - A = 0. (2.57)

1

This is the well-known Coulomb gauge . If we adopt the Coulomb gauge, A will

have a unique solution in a closed region providing the boundary conditions

fiXA=0 on SB (2.58)

fl - A = 0 OI] SH (2.59)

are satisfied [47]. Consequently, VV is determined uniquely [43]. Eq. (2.53) can be
derived from Eq. (2.58).
In order to incorporate the Coulomb gauge, we may add a term -—V (1 / u) V - A

on the left-side of both Eq. (2.51) and (2.52), which leads to [43, 48]

1 1 '
VX;VXA—V;V~A+U%?+UVV=O infll (2.60)

I I
V X (;V X A) — V—V ' A = J3 in 92. (2.61)
,u

Now, the equation of continuity (2.29), which is automatically satisfied in Eq. (2.51),

no longer holds in (2.60). Hence, it must be enforced explicitly as

V - (—o%?— — UVV) = 0 in (21. (2.62)

 

lOther gauges could be used in eddy current problems too, such as the Lorentz gauge [46]

V-A= —poV.

26

Furthermore, we can assume that the media are linear, isotropic and homogeneous.

Accordingly, the magnetic permeability is a scalar. And, by using the vector identity

V><V><A=V(V-A)—V2A (2.63)

E’

and the Coulomb gauge, we may rewrite Eqs. (2.51) and (2.52) as:

711w A + 0%? + oVV = 0 in 12, (2-64)

—V2A = J5 in 92, (2.65)

where V2 E V - V is the vector Laplacian operator.
The normal component of the current density must be constrained to be zero on

the interface 512 since eddy currents are restricted to the conducting region only, i.e.

a 0%? — W) . f1 = 0. (2.66)

In addition, the following constraints are added to boundary conditions on SH and

S 3 to satisfy the Coulomb gauge [45]:

1
12,-VA: O OI] 53 (2.67)

(%V x A) x f1 = 0 on SH. (2.68)

In summary, the A — V formulation of eddy current computation can be expressed

27

{'IS

1 8A
——V2A+U,— +0VV = 0

a tit
8A
V- —— VV 2
0(8t+ ) 0
V2A=J3
fi-A=0
1 -
(—VxA)xn=0
u
fixA=0
1
—V-A=0
u
A1=A2
1
—VA1=-1—V A2
#1 #2
1 1 -
(—VXA1——VXA2)XHZO
#1 #2
a(—%—VV) fi=0
at

in 91
in 91
in 92
OH SH
01] SH
on 8;;
on 33
on 512
on 812

on 312

on 512.

(2.69)
(2.70)
(2.71)
(2.72)
(2.73)
(2.74)
(2.75)
(2.76)
(2.77)

(2.78)

(2.79)

With A — V formulation, the uniqueness of A and V is guaranteed, which leads to

stable numerical solutions for any quasi-static field problem. The handling of internal

interface conditions is trivial and applicable to multiply connected conductors. At

the same time, the current sources can be incorporated easily. The major drawback

lies in the fact that the number of unknowns is high, i.e. four unknowns in conducting

region and three unknowns in nonconducting region at each node.

28

2. 2. 2 A Formulation

In order to eliminate the scalar electric potential V in A — V formulation, we

em lov the following Gauge transformation
. O O C)

A=A*+V\II
8‘11

V=V*———.
at

If ‘11 satisfies the Laplace equation

v51; = 0,

Eq. (2.51) and (2.62) hold the same form for A" and V‘, which are

 

 

V x 1V x A” = —08,A — aVV"
u dt
At
V - (_068t — UVV‘) = 0.

Suppose the function \I’ is defined as

t
\Ilz—f th,

(2.80)

(2.81)

(2.82)

(2.83)

(2.84)

(2.85)

then V“ equals zero. By dropping the =1: sign, Eqs. (2.69)-(2.71) can be combined to

yield
8A

1
V —VA —-—=J.
X” +06% ,3

29

(2.86)

Equation (2.86) holds in both the conducting domain {21 and nonconducting domain
(22, and is referred to as the A formulation.

We must point out that the electric scalar potential V can be eliminated only if
the conductivity is constant everywhere [48]. When regions of different conductivi-
ties exist, the introduction of the electric scalar potential is necessary to ensure the

continuity of current at the interface of different regions.

2.3 Finite Element Method

Analytical solutions for problems governed by the Maxwell’s equations are not
available for most engineering applications. Hence, the implementation of numerical
methods referred to as computational electromagnetics (CEM) dominate electromag-
netic field analysis. Some of popular numerical methods include finite difference (FD)
method, boundary element method (BEM), method of moments (MOM) and finite
element method (FEM). Although each method has its own advantages and disadvan-
tages, FEM is considered more powerful and flexible in handling problems involving

complex geometries and inhomogeneous media [50]. FEM is also regarded as the most

Widely applied numerical simulation method in various engineering applications and

is not limited to electromagnetic field analysis. Moreover, FEM attracts researchers
in all fields and a large number of well-developed software packages are currently
available. The computer codes developed for one particular discipline can be easily
adapted to another discipline. In this section, the basic principles underlining the
finite element method are briefly presented.

The basic strategy of FEM is to discretize the solution domain into a finite number

30

of interconnected subdomains, called elements. The shape and the size of the elements
can be selected arbitrarily, which provide superior capabilities in fitting different ob—
jects and interfaces. In every element, the solution to be sought is approximated with
low-order polynomial functions. For nodal—based FEM methods, this approximation
is calculated on the nodes of the elements and is sufficient to approximate the solu—
tion on the whole element by interpolating with the basis functions, also called the
shape functions. We use the weighted residuals approach or the variational principle
to obtain a system of linear equations that are related to the governing differen-
tial equations and the applications of proper boundary conditions. Once this linear
system of equations, which is typically a sparse, banded, symmetric and positive
definite matrix equation, is solved for the interpolation coefficients in each element,
the solution is uniquely determined throughout the solution domain. The procedures
can be used for solving time-dependent, linear or nonlinear and two—dimensional or

three-dimensional problems [51].

2. 3. 1 Discretization

As mentioned before, the first step in FEM is to divide the solution domain into
small elements. A continuous solution domain is thus replaced with a number of sub-
domains in which the unknown is represented by simple interpolation functions with
unknown coefficients. Thus, the solution of the entire system is approximated by a
finite number of unknown coefficients to be solved. Typical elements used are trian—
gles and quadrilaterals in two-dimensional case, and tetrahedrons and hexahedrons in

three-dimensional case. The element type and the number of elements must be care-

31

lxlsyl.

lx3,y3l 2

Ix; ,y2 I

vN

 

Figure 2.2: A triangular element

fully selected to meet specific requirements on numerical accuracy and computational

resources.

Typically, the solution u of the governing differential equation is approximated

in the local element with a set of linearly independent polynomial functions , i.e.

shape functions. Let us consider a two-dimensional triangular element with vertex

coordinates (221,311), (ac-2,312) and (11:3,y3), which are designated as node 1, 2, and 3,

respectively (Figure 2.2). Also, the solutions on the nodes 1, 2 and 3 are ul, u2 and

u3, respectively. The approximated solution a (as, y) within the element can then be

expressed as [51, 52, 53]

U1

U(I,y) : (21013330 ¢2($7y) ¢3(xay) U2

 

U3

2 u1¢1($,y)+ ”2452 (55,11) + U3Q23 (33’ y) ’

32

I

 

(2.87)

where 61 (2:, y), (15-2 (2:, y) and (15;; (:17, y) are first order polynomial that satisfy:

952' (131,91) 2 515- (2-88)
Therefore,
O(r')-A1— (+5 +1)
.I “3y _ A _2A a] 1x Cly
A 1
452 (Isl/I = 13 = 53(612 + 52113 + C23!)
, A 1
(0303.31) = KS- = 2—5 (as +b3$+63yl

where A is the area of the whole element, A1, A2 and A3 are areas of the subordinate

triangles as shown in Figure 2.2, respectively. Thus, we have

 

 

 

 

1 171 311 1 :1: y
_ 1 _ 1
A—E 15152 y2 ’Ai“§ 1333' 313'
1 2:3 313 1 33k Elk
and
at = $191: - my), bi = yj — yk, Ci = 33k - 351' (289)

where (i, j, k) are cyclical permutations of (1,2,3).

For example, if the node coordinates are given by

(flail/1) : (010)? (17:27:92) : (130), ($32y3):(191)1

33

the basis functions become

[61(r.y)=1-x. 452(x.y)=x-y, ¢3(I,y)=y) (2-90)
and _ ,
U1

fl($79)=[1—x x—y y[' ug - (2‘91)
1%.

 

 

In general, the unknown solution in an element e can be written as

116m) = Z 8: (4.1/>24: = {5. (my)? {us}. (2.92)
[:21

where n is the number of nodes in the element, u]: is the value of u at node k of
the element and Sf: is the interpolation function for node k, which is also known as
nodal function or basis function. The superscript T denotes transpose operation. In

common, the order N of the element e is defined as the highest order of the functions

3; ($.31)-

2. 3. 2 Galerkin’s Method

The next step is to formulate the system of equations, which can be done using
either the weighted residuals approach or the variational principle. Here, we will
adopt one of the weighted residual methods, i.e. the Galerkin’s method, due to its
simplicity.

Let us consider a general differential equation problem defined in a domain Q

34

enclosed by surface 80 = (991 U 0522 with 8521 ft 8522 = 0, which is

Cu 2 f in Q (2.93)
u = g on 891 (2.94)
')
{—1—1 + au 2: h on 892, (2.95)
(in.

where [I is a differential operator up to the second order, f is the source term, u
is the unknown field and ft is the normal vector. Equation (2.94) is known as the
Dirichlet boundary condition, or essential boundary condition, which specifies the
value of u on surface 091. Equation (2.95) is called the Cauchy boundary condition,
or mixed boundary condition, which is a combination of u and its derivative. If a = 0,
Eq. (2.95) yields the Neumann boundary condition (also named as natural boundary
condition)
an

7—: = h on 0522, (2-96)
On

which depends on the derivative of u only.

For a given approximate solution it, the residual R is defined as
R = £21 — f. (2.97)

We wish to estimate an optimal a by letting the integral of the so—called “weighted

residual” to be zero, i.e.

R,=/S,Rdn=0, l=1,...,N (2.98)
9

35

where Q denotes the whole solution domain and S, is the l-th nodal function.
Following Jin’s formulation [39], we consider a single element, i.e. the eth element.
Then, the weighted residual

H; (2.90)

£0
II
7:542

and
Rf :/ Sf (Cue -— f)dQ 2': 1,2,3, ...,n, (2.100)
98

where n is the number of nodes in the element. By incorporating Eq. (2.92) into Eq.

(2.100), we obtain

35:] s:£{se}Tde{ue}—/ fodQ i=1,2,3,...,n, (2.101)
96 9e

or in its matrix form as

{Re} = [K8] {if} - {be}, (2.102)
where
{Re}=[ng, 3...,ij (2.103)
K5,: Qesrcswe (2.104)
5;: nefodQ. (2.105)

Taking into account Eqs. (2.98) and (2.102)-(2.105), we have

[7?] [a] = [6] , (2.106)

36

where [Tf] is a sparse matrix.

Boundary conditions specified in Eq. (2.94)-(2.95) must be incorporated to guar-
antee the uniqueness of the solution. In finite element formulations, the natural
boundary conditions are automatically satisfied while the implementation of the es-

sential boundary conditions must be imposed explicitly through some modifications

to matrix [17] and {5} in Eq. (2.106).

2. 3.3 Solution Methods

The matrix Eq. (2.106) derived from the Galerkin’s method often involves a very
large coefficient matrix [_], called the stiffness matrix. [7] is not necessarily a
symmetric matrix since the operator .C is not required to be self-adjoint. Therefore,
methods for solving such a linear equation system must be chosen carefully. Typical
solution methods can be classified as either direct solvers or iterative solvers.

Direct solution methods, such as Gaussian elimination and its symmetric variant
called Cholesky decomposition (triangular decomposition), lead to accurate solutions
of Eq. (2.106) providing the round off errors are negligible. However, storage re-
quirements can be excessive and deter the use of such methods when a large stiffness
matrix is involved. The problem could be minimized by choosing a node number-
ing scheme that results in a highly sparse and handed stiffness matrix, i.e. a great
number of the elements of the matrix are zeros and most of the nonzero elements are
distributed around the diagonal of the stiffness matrix. Under these conditions, it is
possible to use methods [49] that discard all zero elements outside a diagonal band.

Thus the storage requirement is significantly reduced.

37

Iterative solvers including the conjugate-gradient method, over—relaxation method,
and so on, are quite efficient with regard to memory usage since all nonzero elements
are discarded and the structure of the stiffness matrix does not affect memory re-
quirements. The major concern associated with iterative solvers is in regard to their
convergence properties. In principle, the rate of convergence is closely related to
the condition number of the stiffness matrix, i.e. the ratio of the highest eigenvalue
and the smallest one. Ideally, this ratio should be as close to unity as possible. In
order to achieve faster convergence, a preconditioner is usually applied before the
final solution step to lower the condition number of the stiffness matrix. Common
preconditioners include incomplete LU, diagonal scaling, symmetric successive over—
relaxation (SSOR) and so on [39]. The iterative solvers can also be accelerated using
the multigrid method and easily implemented in a parallel computing environment.

Some other key factors relevant to the success of FEM implementation include re-
liable mesh generation, handling open—boundary domains, etc. In addition, a balance

between the available computational resources and the accuracy must be achieved.

38

CHAPTER 3

THE Q-TRANSFORM

In this chapter, the Q—transform will be presented in a scalar form in the first
section and extended to vector form in the last section. Some of the more impor-
tant mathematical properties of the Q-transform are delineated in the second section
while numerical implementations of the scalar Q-transform are discussed in the third

section.

3.1 Derivations of the Scalar Q-transform

In this section, we focus our attention on the Q-transform in scalar form only. It
is worth noting that familiarity with the scalar Q-transform is essential for using the
Q-transform for analyzing electromagnetic fields governed by vector equations such as
Maxwell’s equations. The validation of this statement will be provided in subsequent
sections.

Consider the following two initial boundary value problems [10]:

V22) (x, t) - c (x) 8m (x, t) = F (x, t) in Q, for t Z 0 (3.1)
a (x) v (x, t) + 5 (x) an?) (x, t) = C (x, t) on 89, for t 2 0 (3.2)
v (x,0) = h (x) in $2 (3.3)

39

 

and

V2u (x,q) — C(X) aqqu (x, q) = f (x, q) in e, for q 2 0 (3.4)

a (x) u (x, q) + 6 (x) an (x, q) = g (x, q) on on, for t 2 0 (3.5)
u (x, 0) = 0 in e (3.6)

oqu (x, 0) = h (x) inQ, (3.7)

where x denotes the spatial variable, 52 E ERN (N _>_ 1) with a smooth boundary 09,
c (x) is a real and nonnegative function, 8,, denotes the normal derivative on the 89
and the two scalar functions a, 6 are known.

To derive the relation between a solution of diffusion equation problem (3.1)-
(3.3) and a solution of wave equation problem (3.4)-(3.7), we follow the method
presented by Romanov [10]. Taking the Laplace transform1 of Eqs. (3.1)-(3.3) and
Eqs. (3.4)-(3.7) (transforming t to s and from q to p, respectively) and denoting

Laplace transformed variables with an over-hat, we have

A

F (x, s) (3.8)

V2?) (x, s) — c (x) [30 (x, s) — v (x, 0)]

a(x)6(x,s)+6(x)a,0(x,s) = C(x,s) (3.9)

 

1The Laplace transform and the inverse Laplace transform are defined as

00 t 1 7 +joo t
F(s) E/ f(t)e—3 dt and f(t) E —/ . F(s)e3 ds,
0 2m 7 _ 300

respectively, where 'y is a vertical contour in the complex plane chosen so that all singularities of the
F (s) are to the left of the contour.

40

and

V211 (Xe?) -C(X) [19221 (Km) —PU(Xa0) -6qu 090)] = f (x,t)) (3-10)

0(X)fl(x,p)+/3(X)01.21(X.p) = 510920)- (3-11)

Invoking initial conditions given by Eqs. (3.3) and (3.6)-(3.7), we obtain

V26 (x, s) — C(X) so (x, s) + c (x) h (x) = F (x, s) (3.12)

a(x)u(x,s)+fi(x)8nu(x,s) = C(x,s) (3.13)

and

A

V211 (X110) - C 0019211 (x,t)) + C (X) h (X) f (X. P) (314)

0(X)fi(x,p)+fl(X)8nfi(xap) = @0910)- (315)

If we let 3 = p2 and if D is the difference between i) (x,p2) and u (x,p), we can

subtract Eq. (3.12) from Eq. (3.14) to obtain

v21“) (x, p) — c (x) p2D (x, p) = 0, (3.16)

providing that source terms F (x,p2) and f (x, p) satisfy

F(x,p2) = f(x,p). (3.17)

41

Similarly, by subtracting Eq. (3.13) from Eq. (3.15), we have

6. (x) D (x, p) + 6 (x) 6,1“) (x, p) = 0, (3.18)
given that boundary conditions C (x,p2) and Q (x, p) satisfy
C (x,p2) = Q (x,p). (3.19)
Observing Eq. (3.16), it is clear that D (x, p) must equal to zero, or

6 (x, p2) = 6 (x, p). (3.20)

Equation (3.20) can be further written in an integral form as

oo oo

[0 v (x, t) exp (—p2t)dt = [0 u (x,q) exp (—pq)dq. (3.21)

Then, if we take the inverse Laplace transform of both sides of Eq. (3.21), with

s = p2, we obtain the Q-transform relation from u (x, q) to v (x, t), i.e.

 

1 °° q2
’0 (xi t) = ZWA qexp(-E)u (xlq)(1qi (3'22)

which is a one—to—one transform. Comparing Eq. (3.17) and Eq. (3.19) with Eq.
(3.20), it is easy to see that the same relation, i.e. Eq. (3.22), exists for the source
and boundary conditions in the q and t domains.

In addition, the inversion of Eq. (3.22), called the inverse Q—transform, can be

42

 

achieved by taking the inverse Laplace transform of Eq. (3.21) [7], which gives

1 ”7+ioo oo ‘
u (x, q) = a] / v (x, q) exp (—p2t) dt exp (pq)dp, (3.23)
l ,7 0

—ioo

where 7): is defined in the same manner as the contour associated with the Laplace
transform defined earlier.

From now on, we will frequently use the following simplified notation to denote

Eqs. (3.22) and (3.23), which are

‘00) = Q{U(q)}(t) and ”(1(0): Q"{v(t)}(q),

respectively (we drop the spatial variable x). Some important properties of the Q-

transform are presented in Appendix A.

3.2 Mathematical Properties of the Q-transform

The mathematical properties of the scalar Q-transform have been derived by Tam-
burrino [32] and others [6, 16]. We begin by rewriting the definition of the scalar

Q-transform adopted from Eq. (3.22) here, which is

2

W) = Q{U(Q)}(t) = 2‘27, AwleM-Z—tw ((1)617.

 

Some of the more important properties of the scalar Q—transform are:

43

Invariance of the Dirac function
Q{5(€1)}(t) = 5 (t)
Linearity

Q{au1(q)+ bu2 (q)} (t) = aQ {at (0)} (t) + bQ {Hz (61)} (t),

where a and b are arbitrary constants.

Translation

Q{u(q— 40)} (t) = {fig—77m (—j—Z)} * 20(4)} (0,

where :1: stands for the convolution operation.

Scaling

Q{U(a(1)}(t) = aQ {u (9)} (0%) ,
where a is a positive scaling factor.

Q—transform of the first derivative

Qtu' (4)} 0) = $98124 (4)} (t) — e {W } (t) _ ME.

2v at3

q

44

(3.24)

(3.25)

(3.26)

(3.27)

(3.28)

0 Time derivative

dQ {u ((1)) (t) _ i 3

(It. — 41:2 Qiq2'u(‘1)l(t)“ 5Q {u (0)} (t)

 

o Laplace transform
5 {Q {u (0)} (0} (S2) = 5 {u (0)} (8)
o Fourier transform
+00
45 {e {u (a) (0} (w) = f u (q) exp (—q¢eZ)dq
0 Relationship with the Fourier sine transform

We define the Fourier sine transform of function u (q) as:

Us (w) E (E [Omit (q) sin (wQ)dq,

we have

0 Relationship with the Fourier cosine transform

45

(3.20)

(3.30)

(3.31)

(3.32)

(3.33)

(3.34)

We define the Fourier cosine transform of a function u (q) as:

+00
UC (w) E fl/O u (q) cos (wq)dq, (3.35)

we have

Qtu (q)} (0 = (E [”0 V1777 [1 — 2e2t<1>(1;§;—p2t)] awe». (3.36)

where (D () is a confluent hypergeometric function [54].

Q-transform of functions with finite support in the q domain

If supp{u (q)} g [q0,q1], where supp{-} refers to the support of a function,

we have:

Q{u (0} = 71,—; [u (q) exp (4:3)]: + 3% jg: u' (4) exp (fig-)3

Q-transform of functions with finite support in the to domain

If supp{U3 (w)} E [w0,w1], we have:

 

1 2 “’1— _1_ WI '—weX) —w2 w
Q{u(t)}= m [Us(w)exp(—w t)]w;+ mfg}; US( ).1( t)d.
(3.38)

In Appendix A, some of useful Q-transform pairs are provided for quick reference.

46

Also given, is a short list of mathematical properties of a related transform, i.e. the

C-transform .

3.3 Numerical Implementation of the Scalar Q-transform

The scalar Q—transform defined in the previous sections reveals the unique map—
ping relation between solutions of wave equation problems and those of diffusion
equation problems. By evaluating Eq. (3.22) directly, we can convert solutions of the
scalar wave equation problems to solutions of the scalar diffusion equation problems.
Unfortunately, the inverse of this process, i.e. the conversion of solutions of the scalar
diffusion equation problems to those of scalar wave equation problems suffers from
numerical instability as pointed out by Gibert [26] and Ross [7] This instability
originates from the mathematical nature of Eq. (3.22), which is a Fredholm integral
of the first kind [55]. We also notice that Eq. (3.22) is an infinite integral from zero
to infinity, which requires some sort of truncation to be imposed during numerical
evaluation. This issue, as well as methods of numerically approximating Eq. (3.22),

will be discussed in this section.

3.3.1 Numerical Approximation of the Q-transform

The following numerical schemes have been proposed in [32] to evaluate the Q-
transform. The numerical approximation of a function u (q) could be expressed in

the form of
N

u (9) = 2: f (at) 3,, (q), (3.39)

‘:1

47

where u)c is a given sequence with k = 1,2, . ..,N, 0 S q1< £12 < < qN, and N is
an integer. f (uk) and sk (q) are unknown functions of uk and q is to be found. The

Q-transform of u (q) is given by

N

Qiu (C7)} (I) = 21(01):) Q{3k (Gll (U- (3-40)

11:21

A discussion on two methods of approximation follows.

Let u (q) be a piecewise constant. Then,

u((1) = U1IH(q - 9k) — H (q — Qk+1)la (3-41)

1132

where one more point up.“ = 0 for qk+1 > qk is added for convenience and H (q) is
the so—called Heaviside step function. Taking the Q-transform of both sides of Eq.

(3.41), we obtain

Q {u (q)} (t) = ZuxQ {H (q - (1k) - H (q — Qk+l)} (t) ,

which yields

Q {U (9)} (t) = Zak [Fe ('5) - Fk+1(t)lv (342)
where (See section 3.2)
F. (t) a Q{H (q — 41.)}(t) = 3%” (~33). (3.43)

48

Rearranging terms, Eq. (3.42) can be rewritten as

N+I

Q {u (q)} (t) = 2(1),. — 6.4) F. (t), (3.44)

kzl

where we also set u0 = 0. Obviously, Eq. (3.44) is more efficient than Eq. (3.42) in

terms of computational costs.

Using a similar strategy, employing a piece-wise linear approximation of u (q):

2

 

u (q) = [up + (u... — up) ——q " qt )[H (q — q.) — H (q — 41.1)]
k=1 Clk+1—(Ik
”‘1u(—)—u(— )
= “I q q" " q “*1 [Hews—Hewett. (3.45)
k=1 qk+l—Qk

we have,

N—I

95330) = Z—“fl—Qttq—qthtq—q.)—H(q-q...)1}.(3.46>

k=1 Qk+1 - (1k

Further manipulations of Eq. (3.46) lead to cumbersome semi-analytical expressions
imposing an extremely heavy computational burden.

Instead, if the support of u (q) is [q1,qN], one can show that (See Eq. (3.37))

1 2 2
Q{u (4)} (t) = __ U1 eXI) _Q_1_ — UN BXP —q—N + 1(1), (3.47)
at 4t 4t
where,
1 N“ 9k+1 (I2
I (t) E —- ujc/ exp (-——-)dq, (3.48)
7ft k=1 qk 4t

where u}, is the first order approximation of the derivative of u(q), i.e.,

I uk+l _ uk
u. = —————. 3.49)
A Qk+1 — (It (

By introducing the complementary error function, which is defined as:

+00

erfc (:r) E % exp (—y2)dy, (3.50)

the integral included in Eq. (3.48) can be evaluated and I (t) can then be shown as:

231;} [..f.( 7,) (3%)]. (3...)

Again, with the rearrangements of terms,

N
f(t) )_=_
k:

2(6 —u3,_ 1) )erfc (21%), (3.52)

1

El“
H.

subject to the condition that us 2 u’N = 0.
The piece-wise linear approximation presented in Eq. (3.52) results in similar com-
putational costs as the piece-wise constant approximation, i.e. Eq. (3.44). However,

the calculation accuracy improves in general.

3.3.2 Truncation of the Infinite Integral

In order to numerically evaluate the Q-transform, which is an integral from zero
to infinity in the q domain as specified before, we must truncate the integral in ways

that offer sufficient accuracy and, hopefully, involve reduced computational costs at

50

the same time. This mission can also be interpreted as an optimization task in order
to determine the support of the Q—transform.
As reported in [32], the Q-transform can be represented in forms that differ from

Eq. (3.22), such as

 

v(t)=Q{u(q)}(t)= 31—7/ °°.,(q,,(,¢2-,)dq, (3.53)

where

10(2) 2 fix exp (£23) . (3.54)

As shown in Figure 3.1 and Figure 3.2, the function 10(2) reaches its maximum at
a: = 1 and then decays very fast as :1: becomes larger than one. This fact reveals that,
when the integral in Eq. (3.53) is evaluated at a certain time to, the most relevant
part of the q domain signal u (q) is constrained within a region q E [q1, ([2] ‘centered’
at q = J23}, or equivalently, a: = 1.

Since q = mm, the lower and upper bounds of the region that contributes most

to the value of Q {21 (q)} (t), i.e. [q], q-z], can be determined by solving equation
w (as) = 2:0, (3.55)

where :60 is a predetermined threshold based on the behavior of function u (q) that is
applied and the expected level of numerical accuracy. The smaller solution 2:1 of the
Eq. (3.55) can be obtained analytically with the help of the so-called Lambert’s W
function [56] while the bigger solution 272 can be reached numerically. Since the :61

is often very small, we can set it to zero while some typical values of 2:2 are listed in

51

Table 3.1.

Table 3.1: Maximum support for numerical evaluation of the Q—transform

 

2:0 1E—4 1E—5 IE6 IE7

 

 

 

 

 

 

321.2 6.9 7.6 8.1 8.7

 

 

In our experience, the piece-wise linear approximation appears to be more efficient
and reliable in most applications. Although the piece-wise constant approximation
is equally efficient, it is less accurate. The higher order approximation could be
very complex in implementation and is not necessary. The interval for truncating
the improper integral Eq. (3.53) in a numerical evaluation can be determined by
solving Eq. (3.55) numerically. Beyond these considerations, the size of interval of

the evaluation points is also critical for a successful numerical evaluation.

3.4 Derivations of the Vector Q-transform

The Q—transform was extended from scalar form to vector form by Lee in low-
frequency EM field analysis [6]. The underlining strategy is similar to that of the
scalar case. However, some special considerations had to be taken into account. For
the sake of thoroughness, this derivation will be sketched here.

If we consider only linear, isotropic and charge free media, the Maxwell’s equations
(Eqs. (2.1)-(2.4)) together with the constitutive equations (Eqs. (2.10)-(2.12)) can

be written as:

V x E (x,t) = —u (x) Eli—(Sill (3.56)
VXH(x,t)=o(x)E(x,t)+e(x)9m—é:(’—t)-+Js(x,t). (3.57)

52

Assuming that u is constant, Eqs. (3.56)-(3.57) implies

.2

0 c) 8
V x V x E (x, t) + no (x) 5E (x,t) + as (x) 55E (x, t) = —ua—tJ3 (x, t). (3.58)

In eddy current applications, the EM field is considered magnetoquasistatic (MQS)
and the displacement current term in Eq. (3.57) can be neglected. Then, by substi-
tuting Eq. (3.56) into Eq. (3.57), we obtain the diffusive governing equation for the

electric field, which is
V x V x E (x, t) + ya (x) g—tE (x, t) = —u(%Js (x, t). (3.59)
The corresponding initial and boundary conditions are given by
E(x,0) 2 0, By 2 E(xb,t) t > 0,
where F is the boundary of the solution domain at x = xb. On the other side, for a

lossless problem, we have

2

V x V x E (x,t) + #5 (x) %E (x, t) = —ug—tJS (x, t). (3.60)

Therefore, referring to the diffusive electric field shown in Eq. (3.59), we can
construct a fictitious wave field in the q domain by replacing the first-order time

derivative with a second-order derivative with respect to q , i.e.

82 (92 ,
V x V X U (x,q) + #0 (x) (3—q2U (x, q) = —;1,5(?2—j3(x,q), (3.61)

53

together with

('3
U(X,0)=3-qU(X,q)|q=o=0, UF=U(Xb»CI) q>0,

where js is the corresponding external source for the fictitious field. It is clear that,
by comparing Eq. (3.61) with Eq. (3.60), the fictitious field U (x,q) behaves as a
propagating wave with a velocity of (ua)_1/ 2 while the independent variable q has
the dimension of square root of time.

Again, we take the Laplace transform to both Eqs. (3.59) and (3.61) from t to 3

and q to p, respectively. We obtain

V X V x E (x, s) + ya (x) 5E (x, s) = —usjs (x, 3) (3.62)
Ely = E (xb, s) — 7r/2 < arg (s) < +7r/2 (3.63)
and
v x v x e (M?) + #0 (We (m) = wfi. (m), (3.64)
filp = f] (xb,p) — 7r/2 < arg (p) < +7r/2. (3.65)

By observing the similarities between Eq. (3.62) and (3.64), and selecting s = p2 and
limiting arg (p) to the domain (—7r/4, +7r/4), we may subtract Eq. (3.62) from Eq.
(3.64) to obtain:

v x V x f‘ (w) + w (x) 2221? (w) = o (3.66)

filp = 0, —7r/4 < arg (p) < +7r/4, (3.67)

54

where E (x,p) is the difference between E (x,p2) and f} (x,p), while conditions for

source and boundary terms

3 (x, 192) =3‘(x, p) (3.68)

A

E (Xb,p2) = fl (xb,p) (3.69)

are satisfied. Then, by multiplying both sides of Eq. (3.66) with the complex con-
jugate of E (x, p) and integrating over the region V enclosed in boundary P, we can

obtain (with some additional manipulations)

/ IV x f‘ (w) I2 + #0 (x) p2lf‘ (w) 1% = o. (3.70)
v
which implies

1509122) = 13 (Km), (3-71)
or

E (x, s) = f] (x, ,5) . (3.72)

The integral relation between the electric field E (x, t) and its fictitious wave domain

counterpart U (x, q) is then presented, by taking the inverse Laplace transform to

both sides of Eq. (3.72), as we did before:

 

E(x,t) = 2%]; ooqexp(—%)U(x,q)dq. (3.73)

This vector relation is identical as Eq. (3.22) for scalar case. A similar transforming

55

relation can be derived from Eq. (3.68) and (3.69) for the corresponding source and

boundary conditions.

56

 

I

0.9

 

 

 

0.8— . , .i. , .4
O.6~ ‘

’>'<‘

‘5 0.5” ................ .4
0.4" .4
0.3* ‘
o.2~ . . ‘

O l l l
0 1 2 3 4 5
x

 

 

 

 

 

 

10‘15 —
—20 9

1O 1 l l l l
1o“° 10‘8 o 6 o" 10'2 10°

1 _ Iog(x) ‘

Figure 3.2: Weighting function of the Q-transform in logarithmic scale

57

CHAPTER 4

TIME-OF-FLIGHT EXTRACTION FOR TRANSIENT

DIFFUSIVE FIELDS

In this chapter, the main issues concerning the extraction of time-of-flight informa—
tion from measurements of transient diffusive fields are discussed with reference to a
scalar configuration. We will show that time-of—flight information associated with po-
sition of a “small” scatterer can be extracted providing appropriate excitation signals
are applied. In addition, numerical simulations are presented along with theoretical
analysis. We will also introduce a vector diffusive field configuration briefly. The ap-
proaches discussed in section 4.1 and 4.3 have been proposed in [32, 33]. Numerical

validation and analysis has been developed in [34, 35].

4.1 A Scalar Diffusion Equation Problem

As we pointed out in the previous chapters, the time-of—flight concept has no
physical meaning in the case of diffusion phenomena. However, with help of Q—
transform, a fictitious wave field could be established with respect to the diffusive
field of interest. The time-of—flight information could be estimated in the fictitious
wave domain and converted back to obtain its equivalent in the diffusion domain. In
general, this process of estimating the time-of—flight could be very complicated. In
this section, we will introduce a simple reference scalar diffusive field problem that
shows how the position of a “small” scatterer embedded in a homogeneous medium

with a point source can be estimated.

58

 

 

 

__ anomaly
A— X D
d
4\ - r
Point Source

Figure 4.1: A small anomaly in the vicinity of a point source

Referring to Figure 4.1, assume that the host medium is a conductor of conductiv-
ity 0'0, permeability M0 and that the displacement current can be neglected. Also, we
assume that the anomaly is confined within the ball BR (x0) & {x E W] [X — Xol S R},
where R << Ixol. A point source is located at the origin of the coordinate system and
the measured quantity is the reaction field (evaluated at the origin of the coordinates
system) due to the anomaly. The unknown quantity that is sought in this canonical
problem is the distance from the scatterer to the origin, i.e. d£ |x0| shown in Figure
4.1. In order to show the main ideas underlying the proposed approach, we consider
a scalar diffusion equation first. We further indicate that the position of the anomaly
is arbitrary. The anomaly is on the z-axis, as shown in Figure 4.1. However, it could

be anywhere else.

4.1.1 Time-of-Flz'ght and Peak Value of the Eddy Current Measure-
ment

To make time—of-flight information “visible” in eddy current measurements, our

strategy is to build a link between the time-of-fiight for the fictitious field and certain

59

characteristic properties of the measured ECT signal in the time domain.
We notice that the measurement from a. wave receiver is zero before a certain time
qo due to the propagation delay of the excitation signal. The Q-transform of such a

fictitious wave domain signal u (q) can be expressed as (See Eq. (3.37) in section 3.2):

 

em = qu (q)} (e) = u (:5? exp (:5) + % jjwu' <q>exp (—§}:)de, (4.1)

where u’ (q) is a unknown function. If the first term of the r.h.s of the Eq. (4.1) is

dominant, i.e.
>>

[0 00 u’ (q) exp (—%>dql, (4.2)

v (t) can then be written in a more explicit form as

 

 

q2
u (q?) exp (”4%)

 

 

and obtain its maximum

+
lvlmax = L (qo )l V —2- (4.4)
90 7T8

at tnm. = q8 / 2. This relation between the time-of-flight go and peak position tmax of
the time domain measurement is illustrated in Figure 4.2.
It has been shown [32] that, if u (q) is finite and different from zero only for q > qo,

u (q) is differentiable for q > qo, lu’ (q)| is bounded by a constant M, and
M << Iu (qJH/qo, (4.5)

60

 

 

 

 

 

 

° ‘10 ‘I

Figure 4.2: The relation between time-of—flight go and peak position tmax of the eddy
current measurement

then Eq. (4.2) can be satisfied [33]. Similarly, in a neighborhood of tmax, when u (q)

contains a Dirac pulse such as

U(q) = (MW-(1;) +h(q-qf)e (4-6)

where a is a scaling constant, [2 (q) vanishes for q < 0 and Ill (q)| is bounded by a

constant M1, and

MI << a/qo. (4.7)

Eq. (4.2) can also be satisfied. The corresponding peak value and peak position are

lvl _ la] 1 6 3
max _ qu 7T 8

and tmax = q3 / 6, respectively. Therefore, under proper conditions (4.5) or (4.7), the
position of the peak of Iv (t)| is, simply, proportional to the square of time-of-flight

(10- If additive noise n (t) is present, the peak is observable providing its maximum

61

value is much bigger than noise level |n (t)|, i.e.

|u(qa‘)l ‘2 .
(10 \/7r:e->>|n(t)| (4-8)

01'

lal l (9)3 >> |n(l)|. (49)

fire

It is worth noting that Eq. (4.8) not only sets a limit on the acceptable noise level
but also imposes an upper limit on the retrievable time-of-flight go for a given u (q; )
and noise level.

The necessary conditions (4.5) and (4.7) are stated in terms of the fictitious wave
u (q). The only way to impose those conditions is to design the source of the fictitious
wave problem properly and the excitation signal to be used during the eddy current

test.

4.1.2 Analytical Solutions of the Scalar Problem

The mathematical description of the eddy current problem shown in Figure 4.1 is

in the form of a diffusion equation, which is

V221 (x,t) - s (x) 8tv(x,t) = G(x,t) in 323, for t 2 0

p (x,0) = 0 in ER3 (410)
Illim u(x,t)=0 for tZO

where s (x) is a function of material properties and G (x,t) = 6 (x) g (t) is a point

source, located at the origin of the coordinate system. The fictitious time domain

62

counterpart of Eq. (4.10) can then be written as

V2u (x, q) — s (x) dqqu (x,q) = F (x, t) in 5R3, for q 2 0

u (x, 0) = 8,121 (x,0) = 0 in 923 (4.11)
lllim u(x,q)=0 for qZO

where F (x, q) = (5 (x) f (q) is a point source in fictitious time domain. Let a (x) be the
conductivity in the region of interest, X (x) £0 (x) /00 ~— 1 be the contrast function,
and let .9 (x) £ [1 + X (x)] /c3, where c0 = 1/,/,uooo is the velocity of the fictitious

wave. With the imposition of the additional constraint, as discussed in Chapter 3,

N) = Q HUN ('5), (4.12)

v (x, t) and u (x, q) constitute a Q-transform pair, i.e.

v (x,t) = Q {11094)} 00- (4-13)

It is straightforward to obtain an analytical solution of the wave equation problem
given by Eq. (4.11) while Eq. (4.13) can then be applied to obtain the analytical
solution of the diffusion equation problem given by Eq. (4.10).

If we define uo (x, q) as incident wave field (which is solution of Eq. (4.11) if the
contrast function X (x) is zero everywhere, i.e. anomaly doesn’t exist), then we have

for the scattered field a, (x, q)

us (x,q) = u (er1) - Up (x, (Il- (4-14)

63

It can be shown that uS (x, q) is the solution of the wave equation problem

V2us (x, q) — s (x) oquu, (x, q) = #0qu (x, q) in 333, for q 2 0
'0
u (x, 0) = dqu (x, 0) = 0 in ER3
lim u (x, q) = 0 for q 2 0

IXI-‘OO

(4.15)
and the solution of the corresponding diffusion equation problem, i.e. the measured
signal in the time domain, is v, (x, t) = Q {us (x, q)} (t).

Eq. (4.15) can be further simplified by introducing the first-order Born approxi-
mation, which says that the unknown field inside the anomaly can be replaced with
the incident field at the same location providing that the scattered field is “weak”.
In addition, we can replace the source term in the r.h.s of Eq. (4.15) with a point

source located at the center of the anomaly. This yields

1
Vzus (x, q) — gaqqus (x, q) = 5 (X - Xe) f3 (<1) in 5R", for q 2 0
0
u (x, 0) = aqu (x, O) = 0 in 993 (4-16)
Illim u(x,q)=0 for qZO

 

where
A K
fs (q) =gaqquo (Xe. q) . (4.17)
0
Ki/ x(x)dx, (4.18)
anxol
and
— x c.
110(XOJI) = _f(q I OH 0)- (4-19)
47rlxol

64

Now, the scattered field u, (x, q) can be treated as a propagating wave launched by a

point source, which should be in the same form as Eq. (4.19), i.e.

fs ((1 — IX - Xol /00)
47r|x — x0] .

 

us (Xoeq) = — (4-20)

Subsequently, by substituting Eq. (4.17) and (4.19) into Eq. (4.20), we can express

the scattered field as:

1 K

8 , = a _ _ — . 4.21
u(xq) 4,|X_XO,034,,IXOI qqf(q Ix Xol/Co Ion/eo) < )

 

To reduce the complexity of the above expression, we set our measuring point at

x = O. This allows us to rewrite the scattered field as:

K
us ) = _ Qaqq — X0 . ,
(0 q) ————(47r |x0| CO) f (q 2| l/CO) (4 22)

The following approach of designing the excitation signal consists of developing a
‘simple’ analytical model for the measured quantity in the fictitious time domain and
thus imposing either constraint (4.5) or (4.7) on f (q). Finally, 9 (t) can be obtained
as the Q-transform of f (q). Notice that the design of the source terms f (q) and g (t)

depends on the geometry under consideration.

65

f(q}

 

 

 

 

 

 

 

 

 

n=1
0 q
:36)
1 n=
n=1
0 t 0 t

Figure 4.3: Time domain and fictitious time domain excitation signals as given by
Eq. (4.23) and (4.24), respectively (q,- equals zero in all cases).

4.1.3 Design of the Excitation Signal

The design of the excitation signal in the fictitious time domain was discussed in

[32]. A candidate signal is:

f ((1) = (q - Ch)" H (q - (11), (4-23)

where q,- is a nonnegative constant. It can be shown that Eq. (4.7) is automati-
cally satisfied for n =1 whereas Eq. (4.5) is automatically satisfied for n =2. The

corresponding time domain excitation signals are

g<t> = omen} (t) = erfc ((1,/m) n =1
Alfie“) (‘Z—i) — 2Q1erfc (zq—C/E) n = 2 .

Examples of both time domain and fictitious time domain excitation signals are shown

(4.24)

in Figure 4.3.

66

Using these excitation signals, the scattered field in time domain can be obtained
by substituting Eq. (4.23) into Eq. (4.22) and applying the Q-transform to the result.

This yields

Bq exp —q2 4t /\/47rt3 n=1
p, (0,1): f ( f/ ) , (4.25)

2Bexp (—q?/4t) Affi n = 2

where BiK/(47rc0 Ixol)2 and qfiq, + 2 |x0|/CO. Finally, the positions of the peak
values of the scattered signal as (0, t) can be shown as

2
q 6 n = 1
tmax = f/ . (4.26)

qfi/2 n=2

The peak values are

m : \/54/(7re3)/q% n = 1
B J‘s/(_xa/ef n = 2

(4.27)

and the solutions |x0| of the inverse problem are

|x0|= (”Gt’“ win/2 ”:1 . (4.28)

(\/2tmax - q,) c0/2 n = 2

We want to indicate that other choices of excitation signals may exist. Also,
it is possible to associate other properties of the time domain measurements with

time-of—fiight.

67

4.2 Numerical Simulation: A Scalar Problem

This section presents a finite element method based numerical simulation to show
the effectiveness of the proposed method. Most of the results presented herein have

been published in [34] except for part of section 4.2.2.

4.2.1 Numerical Implementation

The numerical solutions for diffusion equation problems given by equation (4.10)
have been obtained. In addition, the numerical solutions for the corresponding wave
equation problem given by Eq. (4.11) have also been obtained for the purpose of
cross validation. All numerical calculations were carried out using a commercial
package called FEMLAB®. In order to further simplify the problem and reduce the
computational burden, which is always a concern when finite element techniques are
employed, several important simplifications have been made.

Instead of solving Eqs. (4.10)-(4.11) directly for total fields u (x, q) and v (x, t), we
will solve Eq. (4.16) and its time domain equivalent for as (x, q) and U, (x, t). This
strategy can ease the process of introducing source excitation and, more critically,
avoid numerical difficulties that arise when subtracting two quantities with small
differences (as in the case of Eq. (4.14)).

Notice that the direction from the origin to the small anomaly is not relevant,
we can reduce the original three-dimensional problem to a two-dimensional axisym-

metric formulation. The governing equations for u, (x, q) and u, (x, t) in cylindrical

68

coordinates (r, z) are

021), (x,t) 16hr, (x, t.) + 82'128(x,t) 1 81), (x,t) X (x)
,2

(9,2 ,. 0,. 8 - — 237 “ ’gQ {dqqu (Xoqul (429)

2 _ ‘2 2
W + 1w + 0 Us (X’Q) 1 M X (x)(9qqu (x0, q), (4.30)

(9T2 r 07' (92.2 cf, 8q2 c3

respectively. In the above equations, uo (x, q) is determined analytically using Eq.
(4.19).
Additional simplification is achieved by setting the delay q,- of excitation signal to

zero, which reduces the q domain excitation signal to:

f ((1) = WV ((1), n = 1,2 (4-31)
and the time domain excitation signal to:

1 n = 1
g (t) = . (4.32)

4 t/7r n: 2
Other specifications include: the anomaly is a ball of radius R = 0.1 mm and
the distance between its center and the origin of the coordinate system is |x0]=10
mm. The conductivity 00 of the surrounding material is 3.54 X 107S/m (copper) and

the relative permeability is 1, which corresponds to a fictitious wave velocity of 0.15

m/fi. The anomaly has a conductivity of 2.655 X 107 S, which implies a contrast

69

   

       

v '1 v v v- 7 7
same: sex
Z 'A‘L‘avflg‘. .Is A
l V "17‘ n" ‘-

  
  
 
 

 

7
&‘

m 1" ‘5
é'tifiife'é'fiwmé
“ :3.“.£;,Iiva§§\'flg‘

   

"é
E
v.-
I
C
r- ,4
01
4%

   

 

    
   

  
   
     

 

    

    

    
 
  
 
 
  
 

   
   
     
    

av

4" 4
l ‘JI ‘ ‘vo
6w P‘E’q’fi ‘fl (Hatfiflmuuv $5.
07 ”‘0 Iv wwwmumva 1
53" 945530”: “i'y‘g’rmmmv 3““
AV ‘ mam ‘ixv‘¢"“'-'«'4mw
“A" 0 “_ 0 Quip r 4}
.1. 709‘“ a 35“};th ”4'4 fia‘il‘ay
(A. I r ‘ 1‘) Te‘h'mfiaih
a r 1 (D N ‘x'rm‘evx'mf‘
4fl 5.. ‘5 "v '3‘ :5‘6'
b a: v E 5W"- an: 4
.‘1 7 5"45"
‘ ‘1’, «'4
£59.34"
in ‘
3i

b‘ ’u
'5'
IA

A

4%“ r

1 1"
u
154

A ‘

75%

A
47"

.\
{v}

    
  

V
v

 

  

o :4

7
Z

i

' .A

449:3
4:2

.4

a

   

  
  

gt!
V

32'
D A"
17.3 t

  
 
  
 
 
  

3, .
k

‘t‘i

YA

    
       
   
  
 
   
    

VA?

 

 

 

4‘
4 7‘
pg

 
  

      
  

     
 
   

 
 
       

 
    

A

v i 4
4 9. 4 V" 9":

(A v e a“ 5 " . .. *‘9 51""
'5; 2.2-24:13:: ”5%:an

a n‘ m- :3 15,4 51:73:. assume... 0,“

b ‘V )5?» a, 4 00001; Sgiglu'o'fifigm- 3b

‘ .1900!“ <' wiqu “you.“ "we? ..A‘
i A V, q, n 0 ‘ p “’u“’ fln'h” 1.15 fig:
we 1 «easing: atxxe'lexzewes’sa‘aai
E's '0..po b . 'W 'e" 'e' ,

\V
sex
p
'Afé

1""
'

$9 v
47"
mug

V“
“w
an
7A8

 

A

 

 

 

 

 

Figure 4.4: Left: A schematic of the geometry for the numerical simulation; Center:
The finite element mesh generated by FEMLAB®; Right: Expanded view of the
mesh in the region containing the anomaly.

X = —0.25 in 8;; (x0).

As shown in Figure 4.4 (left), the small dark area represents the anomaly (not in
scale), which is considered the source region in a scattering problem formulation and
its strength is specified by the r.h.s. of Eq. (4.15). A zero Neumann boundary con—
dition is imposed on the left border of the solution domain to simulate the symmetry
axis and a zero Dirichlet boundary condition is imposed on the other three borders to
simulate the open boundaries of an infinite domain. The basic requirement for using
a zero Dirichlet boundary condition is that the size of the solution domain must be
large enough so that the boundary reflections are negligible. We also have an addi-

tional requirement that the duration of the excitation signal must be long enough so

70

that the peak of the scattered field (in time domain) is detectable. However, if the
size of the solution domain is too large, the finite element mesh could be too coarse
to obtain good numerical results considering the limitation of available computer re—
sources. Thus, attention rnust be paid concerning the choice of the solution domain
so that all these requirements can be satisfied. A trial and error approach, although
cumbersome and time-consuming, is often best for optimizing the mesh design. Our
final mesh includes 6517 nodes and 12778 triangular elements as shown in Figure 4.4
(Middle and Right). In addition, we employ a denser mesh and higher order elements
(quadratic elements) in the vicinity of the anomaly (shadowed area) to enhance the
accuracy of the solution. Using a PC equipped with a 2.2 GHz Pentium®IV CPU
and 1 GB physical memory, the major computational burden is the huge memory

consumption when the time-stepping algorithm is used.

4.2. 2 Numerical Results

A comparison between the time domain scattered field as (0, t) and the fictitious
time domain response us (0, t) was carried out for the n = 2 case to highlight the
mapping property of the Q-transform. Figure 4.5 shows agreement between the (nu-
merical) solution of the diffusion problem (Eq. (4.29)) and the numerically computed
Q-transform of the (numerical) solution of the wave propagation problem (Eq. (4.30)).
Moreover, the peak of the time domain response (the observable quantity vs (0, t)) is
very close (within 2.3%) to the theoretical value of 0.0098 predicted by Eq. (4.26).
This confirms that lxol can be estimated using Eq. (4.28) once the peak position has

been determined from the measured waveform v8 (0, t). The numerical solutions of

71

 

v (0,t)/B (s-1 ’2)
M L?) k 01 0') NJ

—b

 

D 1 1 1 L 1 1 1 1

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
1(3)

 

 

Figure 4.5: The plot of us (0, t) (solid) together with the plot of Q {us (0, t)} (*). The
predicted peak position is at t = 0.009 s

diffusion equation problem are presented together with the corresponding analytical
solutions (see Eq. (4.25)) for cases n = 1 and n = 2 in Figure 4.6 - 4.7, respectively.

Two major sources of error are of concern. One is caused by the finite duration
of the excitation signal while the other results from the Born approximation.

We know that the duration 7’ of an excitation signal cannot be infinite long in
numerical simulation. Hence, we need to determine the minimum value (Tmin) of r
such that the error caused by truncation of the signal is negligible. Since the Q-
transform is a linear operation, the scattered field u; (0, t) resulting from a truncated

excitation signal is (for n = 1) (see Eq. (4.25) in section 4.1.3)

2 2 2

, qf <1 (1+A)q; (1+4) 4;

t B = —— — ———_ _—
v.(0, )/ exp( 4) 2 7t, exp 4,

— t
_ q! (1} (2A '1’ A2) q; . .
_ Wexp (71?) [1— (1+ A)exp (— 4t )] 1(4-33)

[\3

fi

(a.
0.3

 

72

 

 

 

 

04 ‘

 

0 0.006 0.01 0.015 0.02
t (s)

Figure 4.6: The plot of v, (0, t) (*) together withthe plot of the approximate response
(solid). case n = 1

where )1 is defined as the truncation constant:
Air/qr (4.34)

It can be show that, at the point t = tmax = qfi/6, the ratio of of, (0,t,,m) and

v3 (0, tmax) is

3A2 GA
—+—) . (4.35)

1"” C1: v;(01tmax)/vs(01tmax)=1_(1+ A) 8X13 (_ 2

Accordingly, we can show that the relative error (1 with respect to the peak value

is less than 1% if the truncation constant A Z 1.14, which yields rm," = 1.14qf.

73

 

(h

k
I

 

{30,078 (5'1’2)

[\J
L

—L
l

 

 

 

 

0 0.02 0.04 0.06 0.08 0.1
t (8)

Figure 4.7: The plot of us (0, t) (*) together with the plot of the approximate response
(solid). case n = 2

Similarly, for n = 2, we have

, _ 2 fi __2_ _(1+/\)2q}
v3(0,t)/B — mexp(—4 t) mexp< _—4t )

= exp 42—) (ex- Warn]

X2+2A
2

 

 

1— C2 = u; (0, tmax) /us (0, tmax) = 1 — exp (— (4.37)

and the relative error (2 with respect to the peak value is less than 1% if A Z 2.20, or
Tmm = 2.20qf. The effect of the truncation of the excitation signal is shown in Figure
4.8 and the relative errors with respect to the peak value are plotted in Figure 4.9
(Left). The error with respect to peak position has been calculated numerically and

presented in Figure 4.9 (Right). We notice that the relative errors with respect to

74

Figure 4.8: Shifting effect of truncation on excitation signal, q, = 0.15. Both ideal
response vs (0, t) (solid) and shifted response 12; (0, t) (dashed) are normalized.

.0
cn

Normalized geak Value

Figure 4.9: Truncation error. Left: Relative error with respect to peak value, 72 = 1
(solid line) and n = 2 (dashed line); Right: Relative error with respect to peak

.0
a.

.0
M

in

Normalized responses

 

 

 

 

'5.
0.8 - ”x s
| '\.~
I "s,
0.5 "\
l 'x.
I ."s,
I ‘U
0.4 l ‘3“ 4
0.2 . 1 “i
I
0 1 l g 1 1
0.5 1 1 .5 2 2.5
tl(qf2/6)

 

A
I

 

 

 

 

Normalized Peak Position

 

 

 

 

 

 

position, n = 1 (solid line) and n = 2 (dashed line)

75

 

both the peak value and peak position decrease with increasing A.

The Born approximation was invoked to obtain Eq. (4.16) and the related in—
version formula (4.28). We notice that, if q,- = 0, the relative error caused by Born
approximation depends only on two dimensionless parameters: y-é—R/d and the con-
trast X. In Figure 4.10 and Figure 4.11, we show numerically calculated relative errors
with respect to peak value and peak positions using dimensionless coordinates x’ and
t’ defined by x’=x/ |x0| and t’ =t /T, where Ti- (Ixol/Col2, respectively. The value of
contrast used in Figure 4.10 is X = —0.25 while that of Figure 4.11 is X = —0.99.
These values of contrast have been chosen because they are representative of the
range of interest, which is the interval from —1 (perfectly insulating anomaly) to 0
(no anomaly). Our results show that the error in the estimate given by Eq. (4.28) is
very small for 7 S 0.1 and the relative error grows faster than a linear rate for larger
values of 7. Also, the relative errors on both peak value and peak position are larger

for higher contrast.

4.3 A Vector Diffusion Equation Problem

Most of ECT problems are governed by the vector diffusion equations. The treat-
ment for vector diffusion problem has been developed in [32]. Numerical validation
results and additional analysis can be found in [35].

Let us consider the reference problem solved in sections 4.1-4.2. Keeping the same

material properties, the electromagnetic field generated, under quasi-magnetostatic

76

 

20 x e .

 

 

 

 

Error (%)

 

 

 

 

0.1 0.15 0.2 0.25 0.3
Y

Figure 4.10: Error in the predicted position (*) and amplitude (o) of the peak as
function of 7 for n = 1 (dashed line) and n = 2 (solid line), respectively. Case
X = -0.25.

20

 

15,.........,...:. . , q ......... ._

 

 

Error (%)

 

 

 

 

0.1 0.15 0.2 0.25 0.3
7

Figure 4.11: Error in the predicted position (*) and amplitude (o) of the peak as
function of ’7 for n = 1 (dashed line) and n = 2 (solid line), respectively. Case
x = —0.99.

77

approximation, in the form of the Maxwell’s equations, is given by:

V X E = —lt531H in Q
VXH=J0+0E inf)
H(x,t-—-0)=0 inQ

E(|x| ——>oo,t)=0

(4.38)
(4.39)
(4.40)

(4.41)

where the domain {2 2 ER3 is enclosed by surface 89 and fl is the outward normal on

surface 89. Additionally, the boundary condition (4.41) is equivalent to the radiation

condition at infinity. Accordingly, a fictitious wave problem

Vxe=—u8qh inf)
VXh=jo+ane 1119
h(x,q=0)=0 inf?

e(x,q=0)=0 inf?

can be established to fulfill

E(X7t) = Q{<9qe(x,(1)}(t)
H(Xet) = Q{h(X.Q)}(t)

78

(4.42)
(4.43)
(4.44)
(4.45)

(4.46)

(4.47)

(4.48)

if the source terms satisfy

Jo (X7 0 = Q Up (X) (1)} (0- (4-49)

For our fictitious time domain problem (Eq. (4.42)-(4.46)), we assume for a point
source, a magnetic dipole m (q), placed at the origin of the coordinate system and
oriented along the z—direction. Again, we choose our measurable quantity for this
problem as the scattered field evaluated at the origin. The governing equations for

the scattered fields are given by:

V X e, = —u8qh3 in Q (4.50)
v x h8 = j, + aooqe, in a (4.51)
h, (x, q = 0) = 0 in a (4.52)
e, (x, q = 0) = 0 in Q (4.53)
e, (le ——> oo,q) = 0 (4.54)
Where
js = ooXaqe (4.55)

is the source term. By using the linear Born approximation for a small anomaly, the
source j 3 can be treated as an electric dipole located at x0. Note that e in Eqs. (4.55)
is replaced by e0 and

_usinfi m” (q — d/co) + m’ (q - 61/60) i
47r cod (12 it)

 

(4.56)

 

80 (X0, q) :

79

is the solution of Eq. (4.42)-(4.46) in case of o = 00 and L, is a direction vector
pointed in the circumferential direction with respect to the direction of the magnetic
dipole in spherical coordinates (r, 6, (15) [38]. We can denote the approximate source

as

is = UoKaqeo (X010) - (4.57)

Then, we have the scattered field

lsine' Egg—dc Ego—ea .
hs (01g): 477 (Cod/O)+ ([2/ ) 1W? (4°58)

 

 

 

where l is the distance between two electric charge and a new spherical coordinate
system (7", o’fl’) is established with its origin located at x0. Then, by substituting

Eq. (4.57) into Eq. (4.58), we have

I!” (

_HoooKlsinflsinO’ m q-qf) +2mm(q-CI/l +fl(q_—_qL)_ i ,
(4W)2 (Codlz Cod3 ‘14 ¢ .
(4.59)

ha (0, q) =

 

 

where qf = 2d/co.
Now, we need to determine a suitable choice for m (q) so that the last two terms
in the bracket of Eq. (4.59) can be neglected. As proposed in [32], by choosing 777. (q)

proportional to q", Eq. (4.59) can be reduced to

_ lay 11000 Kl sin 6 sin 9’
(40001702

 

h. (o, q) = H (q — we.) . (4.60)

80

Applying the Q-transform to the above equation, we obtain

 

H, (0, t) =

_i¢,)ioUoKlsin08in 6’ 1t 'p( (I?) , (4.61)

2 — ex
(4c0d7r) 7F

whose maximum value (2d / CO) is located at t = q; / 2.

In this chapter, an inversion method based on time-of-flight for diffusive phenom-
enon has been shown to be effective for detecting anomalies in conductive materials.
The design of the excitation waveform is critical in order to estimate the time-of—flight
properly. Numerical modeling has made it possible to investigate the limits of validity

of the inversion method.

81

CHAPTER 5

TIME-OF-FLIGHT EXTRACTION FOR HARMONIC

DIFFUSIVE FIELDS

Q-transform based tiIne—of—flight extraction techniques can be implemented not
only for transient diffusive fields, but also for harmonic cases. In this chapter, in
contrast to the inversion methods discussed in chapter 4, the time-of-fiight informa-
tion is estimated from nonlocal measurements first. The TOF extraction scheme is
based on a novel approach that involves the removal of the material interface. Apart
from 5.2.1, the methods presented in this chapter have been proposed in [36] and

numerically validated in [37].

5.1 Time-of-Flight in Harmonic Propagative Fields

In the time domain, time-of-flight for the fictitious field appears as a delay in the
measured field at a given field position (local quantity). In the frequency domain,
the time—of-flight appears in terms of the factor e‘j’”, where k = w/c is wave number
at frequency w, c is the wave velocity and 1" represents distance. Alternatively, the
time—of—flight information can also be extracted from the Fourier coefficients (nonlocal
measurements) [36]. To introduce this concept, we begin by investigating, in a full
electrodynamics sense, the canonical problem of radiation from an infinite current-
carrying line embedded in an unbounded homogeneous medium.

Let the current-carrying line be parallel to the z axis and z invariant. We know

82

that the electric field radiated by the line is given by

#00110

Em...) = — H52) (km. (5.1)

where i0 = ioiz is the strength of the electric current, H62) () is the Hankel function
of second kind of zero order and R is the distance from a point x to the current-
carrying line [57]. Since this field is independent of z, the field is given by E (x; jw) =
izEz (r, 6; jw), where the subscript denotes the 2 component of the field. Also, we

denote the projection of current-carrying line on the my plane as x0 (r0, 00). Then, the

electric field can be expanded in a Fourier series, which is (the subscript is dropped)

+oo

E(r,0;jw) = 2 E, (r; jw)ej"9 (5.2)
. w i . _.n

En(r;Jw) = —&::—Hl(:|)(k7‘) [20J|n|(kr0)e , 90], (5.3)

where Jn (-) is the nth order Bessel function of the first kind [58].
Using the asymptotic approximations of H112) (hr) and Jn (kro) for large arguments

(kr >> 1, kro >> 1), we have [54, 60]

H?) (:r) = gexp [—j (a: — n7r/2 — 7r/4)] + O (lb/a?) (5.4)
Jn(a:) = \/;-§-:COS(IE—n7T/2—7T/4)+O(I/\/;§). (5.5)

83

Circle of Evaluation

o;\ -
r‘L’O 4!}, ¢0)—————>

 

Figure 5.1: Two time-of—fiight terms included in the Fourier coefficients

Then, the coefficient E, (r; jw) can be expressed as:

#0010

= _ [1062.0 { -—jk(T—To) _1 In] - —jk(r+ro)} ~jn00 5 6
_27r rro e + ( ) 3e 6 . ( . )

As shown in Eq. (5.6), the time—of-flight e‘j’" found in local measurements is split
into two separate terms. One term, e‘j"(”’°), represents the time-of-flight from the
source position x0 (r0, 60) to the closest point of the evaluation circle while the other
term, e“jk("+r°), represents the time-of-flight from the source to the most distant point

on the evaluation circle (Figure 5.1).

5.2 A Cylindrical Configuration

In order to illustrate our time—of—flight extraction strategy for harmonic eddy cur-

rent fields, we take the cylindrical structure shown in Figure 5.2 as our reference

84

10

 

t-

 

 

V

 

Figure 5.2: Geometry for the reference problem. Left: 3D View; Right: 2D Cross-
section at z = 20

problem throughout this chapter. Such cylindrical geometries have been investigated
extensively in a variety of disciplines [59, 61].

In this problem, a long conducting cylinder is placed in an unbounded homoge-
neous non-conducting space. A long current-carrying wire of current intensity [0 is
embedded inside the cylinder. The axes of both the cylinder and the wire are parallel
with the z-axis.

We assume that all materials involved are linear, isotropic, non-dispersive and

nonmagnetic. Since the displacement current is negligible, the electric field E satisfies

V X V X E (X;J'w) + jwrtan (X;jw) = -jwquo (X;jw), (57)

where x = xi, + yiy is the position of a point.
Obviously, at the central part of the cylinder, the excited diffusive fields have

no 2 dependency. Thus, we can simplify this three-dimensional problem into a two-

85

dimensional problem on a cross-section of the cylinder, such as 2 = 0. At this
cross-section, the projection of the current-carrying wire is located at (330,110). The
conductivity of the cylinder is 00 = 3.57 X 107S/m (pure aluminum).

Now, with Jo (x;jw) = .102 (ngw)iz and E (xgjw) = E2 (x;jw)iz, Eq. (5.7) can

be reduced to a scalar equation in the form of

V2192 (x; W) -J'w1t00(X) E2 (X;J'w) = jwquOZ (X;jw), (58)

where J02 (x;jw) 2 [Oz (w) 6 (I - 170) 6 (y — yo) and 10(w) is the strength of the im-
posed current. Furthermore, we can drop the z subscript for simplicity and define

the contrast function as X (x) E o (x) /00 — 1 to obtain:

V213093141) - J'wuoao [1 + X (X)l E (X;J'w) = jwuolo (w) <5 (05 - 930) <5 (y - .40)- (5-9)

As stated in Chapter 2, the boundary conditions on the outer surface of the cylinder

are:

a x E (b‘,6;jw) = n x E (0+,6;jw) (5.10)

fiXH(b‘,6;jw) = fixH(b+,0;jw), (5.11)

where fl is the outward normal on that surface, superscripts ’—' and ’+’ denote the

medium outside and inside the cylinder, respectively. In the cylindrical coordinate

86

system, the above boundary conditions can be reduced to

SE (b’, 6,301) /(’)7‘ 2 (3E (b+,9;jw) /8r (5.12)

E (b‘,6;jw) = E (6+,6;ja). (5.13)

To establish a corresponding fictitious wave equation problem in the frequency do-
main, let c?) denote the fictitious frequency domain dependency, we have E (r, 19; jw) =

Q {E (r, 0; 363)}. The Q-transform relation in the frequency domain is

E(r,6; jw) = E (p, a, 71—5) , (5.14)

which follows from the Laplace transform relation given by Eq. (3.72) in section 3.4.

5. 2.1 Solution of the Radiation Problem

In order to understand the relation between fields inside and outside the cylinder,
we will solve this radiation problem analytically in the fictitious frequency domain
using the method of separation of variables.

The model of this radiation problem is shown in Figure 5.3, where the conductive
cylinder (r _<_ b) and the free space outside the cylinder (r > b) is denoted as domain
I and II, respectively. The conductivity 01 is equal to zero and 00 = 3.57 X 1078 / m.
To further simplify the problem, the point source is placed on the 113-axis at a distance

(relative to the origin) of r0 = V133 + yg without loss of generality.

87

‘5.
,0-
'.
*1

 

“V

 

 

Figure 5.3: Geometry for the Radiation Problem

The governing equations can be written as:

V217] (x;j&7) + ESE (x;j{IJ) = jfiuoj, (JD) for 0 S r < b
V2E(x;j&7) = 0 for b < r
68 (016351) /ap = of? (5+, 6, ya) /8r
E (mama) = E (Mega)

lim E(r,0;jc:2) = 0,

r—v+oo

where .7, (3'6?) E To (3'07) 6 (:1: — $0) 6 (y — yo).

(5.15)
(5.16)
(5.17)
(5.18)

(5.19)

The general solutions of Eqs. (5.15)—(5.19) at an arbitrary evaluation point x =

(r, 0) are of the form [58]:

~

~ ~ 4
E (130,315) =

71:1

88

1‘25 (2) ~_ , +°° , . ~
H0 (kor)+nZflC1,ncos(n9)Jn(kor) er

+00
2 Cg,” cos (n6) r‘" r > b

(5.20)

where r’ 2: |x — x0], r 2 IX], C1,, and C2,), are coefficients of the solutions in domain

I and II, respectively.

By incorporating the interface conditions (Eqs. (5.17)-(5.18)) on the circle r = b

and using the identity [58],

H<2)(k0 7J2) Zn 71))H,(,2) (for) .I,, (for) COS (W),

71220

with

The system of equations for Cl,” and C2,, is

0.577“ (n)

J,, (1206) C1,, — 54102,, = é‘——'————H,<,‘2> (1'05) 1,0205)

4
n b'"C2 n 110(31,T(n)

jn (gob) 01,11 + a , = 4

The solutions of this system are

 

"4(7)?" ~
C1,: “—°————°" (k0 (a)H,(,2_)1 kob) /J,,_1 1405 n—0 1 2
H_0 "J __
02,”: T2773)“ 00 (koa) )_/J,, 1(kob) ,n—1,2,....

Hf?) (hob) Jn (hoa) .

(5.21)

(5.22)

(5.23)

(5.24)

(5.25)

(5.26)

Embedding Eqs. (5.25)—(5.26) into Eq. (5.20), the complete solution can be shown

89

to be:

 

 

’ ~~, , ~ +oo .1, Ea. H53) 7.1 b) J,, '12 r
”0:11 —H62) (kor’) + gTOz) ( 0 ) Jn_11 (2;) ( 0 ) cos (n6)
5(7‘,9;jfi)=< ”Sb

. 7.. +°° 5 " Jn ((500)
_Jfiv MO/Uo ”2:; T (n) (;) m cos (n6)

,r>b

 

(5.27)

Identities used to manipulate the Bessel functions are [54]

where Zn (2) can be any of the Bessel functions J,, (2), Y" (2), H5,” (2) and H62) (2).

The analytical solution given by Eq. (5.27) can be obtained in another way. Let
us consider a generalized case of a nonzero conductivity 0) in the region outside the
conductive cylinder. In this situation, the general solution in domain II satisfies the

Helmholtz equation instead of the Laplace equation since the radiated field exists in

90

 

that domain. The general solution then has the form:

 

7'1: , ~ +00 ~
.6“) H52) (kor') + 2 C1,” cos (71.6) J" (kor) r S b
E (no-.132) = e... "=0 , (528)
Z ngncos (n6) H32) (7.5.x) r > b
11:]

where l0] 2 52,/11001. Imposing the same interface conditions (Eqs. (5.17)-(5.18)) as

previously, the equations for solving the coefficients are

~

.1, (12,5) 01,, — 57:02,, = Wm,” (E05) .1, (an) (5.29)
J}, (1205) 01,, + 565-"sz = W171?) (1505) J,, (2305). (5.30)

The coefficients are obtained as

= quDJn (6900) T (n) E 01H?) (gob) H6231 (Elb) — 230115.231 (Rob) H722) (Rb);

 

 

1'" 4 EIJ, (7605) H533, (76.5) — EOJ,_. (205) Hi?) (15.5)”31)
C _ jHOJJJn (600) 710‘) f: 1 {5 32
2’" — 27") £11, (£05) H53), (7615) — EOJ.._, (E05) H5?) (12.5) ' )

If conductivity 5, goes to zero, the combination of the coefficients Cm, and Cg", with
Eq. (5.28) is shown to be identical with Eq. (5.27). The result makes use of following

asymptotic relations:

 

x_. . 2 n F .
H5?) (:13) —9 ] (—) (n) ,n > 0 (5.33)
:1: 7r
x_. , 2
H53) (:16) ——9 J (E) In — ,n = 0 (5.34)
77 :1:

91

and the formulas

 

 

‘EIC‘

) ,n > 0 (5.35)

”32) ’12,) ~. 2 2 2 2 b 0
T(~—)‘ ALP [j (—) In ~—] / [j (—) ln :—] :1: <—-) ,n = 0. (5.36)
115, l (1.5) W kl'r 7r kpb 7”

5.2.2 Time-of-F light and the Solution of the Radiation Problem

: m
:A 3A
N ‘3:
AA
0‘ 3‘
0‘ ‘3
:l
l
O
{—1
A
SWIM
i
V
3
fl
>12?
:2
I—_.J
\
[Ks]
/"\
SWIM
0"
V
3
"'1
>1;
L—J
II

In order to understand how the time-of-flight information is embedded in the solu-
tion of the radiation problem, we modify Eq. (5.27) slightly, to introduce the source
at position x = (r0,60) and replace cos(-) function with an exponential function.

Thus, the solution of the radiation problem can be rewritten in the form of:

 

of. ~ +°° ~
J‘O 'Hg") (kor’) + Z 13:, (r;j172)ejno ,r _<_ b
130:0; )0) = ,0, ~ "=-°° , (5.37)
2 En (rm?) 6an ,r 2 b

n=—oo,n;é0
where the angular Fourier coefficients E, (r; jib) and Ef, (r; jib) are

~ (2) ~
E: (M5) = Lia-J [le14 (6070) [No] Jlnl (for) ::|:11((;::)) (5-38)

Inl
~ .~ . CO b ~ ~ —'11 1 ~
En(r;]w) = fits—:5 G) [I,-J.,,. (horo) e 1 90] 11 I 1 (1.05). (5.39)

 

 

If the evaluation circle is coincident with the surface of the cylinder, i.e. r = b,

92

the coefficients can be expressed as:

 

 

~ ~ ,~1~, ~ ~ _
E (0,310) = "0: D,(jw)Fn(jw)e‘J”“0 (5.40)
~ In!

~ .~ #00011 b .~ —'0

Enbi' : — _ n" jno, °
( .792) J 27w (7,) G (lee (541)

where

0,023.3) = .1],l (horo)Hl(:|)_1 (205) (5.42)
11.00) = J...(ie'ob)/J..._.(Eob) (5.43)
(1,047) = J,,l (220m) /.1,,._1 (75,5). (5.44)

As proposed by Tamburrino [36], we can apply Eqs. (5.4) and (5.5) to Eqs. (5.42)-

(5.44) to show

Notice that Eq.

2 1 [e—J'Eow-TO) + (_1)lnl je-J'Eo(b+ro)] e—jnoo (5.45)

filibov 0T0
_ . b 1 +6—2j(Foro—(n|n/2—pr/4
_ JV 7'0 1 + e—2j(horo—|n|w/2+n/4

)
)
= _jfil+6—2j(koro—|n|7r/2—1r/4:. (547)

r01+ e—2j(koro—|n|1r/2+7r/4

~

e-Jko(b—ro) (5.46)

 

 

(5.45) is analogous to Eq. (5.6) that includes explicit time-of-

flight terms and F, (jib), G, (ij) are periodic distributions of the frequency Eb that

can be expanded in term of a Fourier series. Furthermore, it was shown that such an

expansion reveals the existence of the time-of-fiight term of type (b — TO + 2mb) /co

and (b + r0 + 2mb) /co in F, (ij), where m is a nonnegative integer. These kinds of

93

time—of—fiight terms represents the multiple reflections of the fictitious waves on the

interface r = b.

5.3 Interface Removal and Time-of-Flight Estimation

The time-of-flight terms identified in the solution of the radiation problem (Eqs.
(5.37)-(5.39)) are complex due to the multiple reflections from the boundary of the
conductive cylinder. In contrast, the fictitious field I7 (x; jab) radiated from a point

source embedded in an unbounded homogeneous medium given by,

‘7 (X; 1'97): -#0:1

 

H"2 (E, [x — xol) (5.48)

is related in a simple way to time-of-flight |x — x0] /c0, which is in turn related to
the distance between the evaluation point x = x (r, 6) and the source location x0 =
x (ro,60). Therefore, we need a procedure for removing the effect of the material
discontinuity and to cast the TOF extraction issue in terms of TOF estimation from

a “simple” field given by Eq. 5.48 (see [36]).

5.3.] Interface Removal Using Linear Filtering

To find the linkage between the fictitious field I7 (x; job) and E (x; job), we represent

~

V (x; jib) in a Fourier series

(r, 6 ,jw) 2+: V,( (r ,job) )6)“, (5.49)

n=—OO

94

~ ~

V, (r;jcb) = —E%J—JH(2) (kor) [1:41],” (Fore) e‘jnoo] . (5.50)

I"!

Let us write Eq. (5.39) in form of

 

 

~ ~ ,1 l l I 1 ~ ~ .
En (r;jw) = —jé::§ (7;) ~ [1,4,] (1mm) (”90] , r 2 b, (5.51)
J,,,_. (1405)

where an additional term for case of n = 0 is added to original equation for conve-
nience of mathematical manipulations. Therefore, for r 2 b

17.. (me) = 540%” (%)” Hfj,’ (top) .11..-. (4.4)] E. (p.16)

2 L—jlbo—gé (%)]n] Hlifl) (7601“) Jlnl—l (Rob)

~

En (r3197), (55?)

 

 

where r’ 2 b is the radius of the evaluation circle. Thus, by substituting Eq. (5.52)

into Eq. (5.49) and invoking the following expression

~ ~ 1 27f~ ~ . I
E,(r’;jw) = — / E(r',6';jw)e""0d6', (5.53)
27f 0
we have
~ .~ 3605 2“ +00 7" ln] (2) ~ ~ '(9-0’)
We...) ._. __4_0 .51.. (.5) H, (1,31,, (page
E (r', 9’; jib) d6’. (5.54)

The transform relation given above is a linear filtering process that maps the mea-
surement E (r’ , 6’; jib) taken on a circle of radius r’ onto the field V (r, 6; job) evaluated

on a circle of radius r. The linear filter involved is stable for r’ > r [36]. By further

95

investigating the linearity of this transform, we can expand the field radiated from
a point source to that of a line or a surface source. Therefore, this linear filtering
method for interface removing could also be applied to more realistic defects of finite

dimensions [37].

5. 3.2 Time-of—F light Estimation in the Frequency Domain

The corresponding linear filter for solutions of diffusion equation problems can be

obtained using the Q—transform. If E = Q {E} and V = Q {V}, we have
E (..,j...) = E (x, «172), (5.55)

which is equivalent to replacing (.b in solutions of fictitious wave equations with ,/w / j

in solutions of diffusion equations. Thus, we obtain

”if (%’)'"'x8 (2)714 (filel

n=—oo

 

- 21r
V (139mb) = ____b,/fi.3 ]

460 0
-E (r', 6';jw) d6’. (5.56)

Also, from Eq. (5.48), we have

,. _ flex/071109)) (2) glx—xol
V(r,6,jw)-—— 4J3 H0 (\[J' Co > (5.57)

 

and in case of 100 [x - x0] >> 1,

4 . . .
. - g _#0 vij, (J01) Co _\/jw|x—x0|
V (r, 6,302) _. 4 , / 71 Ix _ X0] exp ( Co . (5.58)

96

 

 

In practical situations, we have noisy measurements Ema, = E + AE, where AB is
the noise term. The associated noisy field Vmeas = V + AV can be obtained with Eq.
(5.56). The estimation of the time-of—fiight term qTOF = [X — x0] /c0 can be achieved

by solving the one parameter minimization problem:

NIH/511' (2) ( lw ) x a: -
— . 1! “Sq _ Vmeas T 76 1.7“) 1 559)

where (r‘, 6*) is a given observation point and q" is the estimate of qTOF.

argmin
qTOF : q'

 

 

The estimation scheme can be summarized as follows:
1. Measure the data Ema, (r’, 6’;jw) for a fixed r’ and 0 g 6’ < 27r;

2. Convert Ema, (r', 6’; jw) into Vmeas (rk, 6),; jw), where (rk, 6k) is the k-th obser-
vation point. Repeat this step for k = 1,. .. ,NT, where NT is the number of

observation points;

3. Estimate q}: by solving the minimization problem (Eq. (5.59)) for (r‘,6*) =

(rk,6),), Repeat for k = 1,. . .,NT;

4. Calculate the distance d)c = cog; from the observation point (rk,6k) to locate

the source

It is worth pointing out that the evaluation of filter coefficient through Eq. (5.56)
can be expedited using the Fast Fourier transform (FFT) if the sampling points
(r', 6’) are uniformly distributed on the circle r = r’, which means 6’k : 27rk/NT for
k = 0,. . . ,NT — 1. On the other hand, the one parameter minimization problem is

formulated under the assumption that the source strength I,- (jw) is known. Other-

97

wise, we need to consider the value of I, (jw) as an additional unknown in Eq. (5.59)

to make it a two parameter minimization problem [37].

5.3.3 Numerical Simulations

To validate the inversion method based on the interface removing technique, we
present a numerical example with reference to Figure 5.3. In this simulation, we
choose the radius of the cylinder to be b = 30mm. The position of the point source
is (r0, 60) = (10mm, 00) and conductivity of the cylinder is 00 = 3.57 X 1078/ m. The
radius of the circle where E (r’, 6’; jw) is measured is r’ = 32mm and the radius of
the circle where V (r, 6; jw) is calculated is r 2 35mm. The excitation frequency is
75Hz which provides a reasonable skin depth of 9.8mm for possible detection of the
electric field and also satisfies the high frequency requirements of Eq. (5.58). Also,
the excitation strength is set to be I, (jw) = 1 for simplicity.

The simulation of the two—dimensional electric field E (r, 6; jw) has been conducted
using the commercial finite element package FEMLAB®. A circular solution domain
of radius 1000mm is combined with a Dirichlet boundary condition on its outer border
to simulate an unbounded space. To achieve best accuracy, we employ a high mesh
density around the boundary of the conductive cylinder and the source by adding
artificial boundaries (the automatic mesh generator produces a mesh that is dense
around all boundaries). The final mesh (see Figure 5.4) employs 24497 nodes and
48800 triangular elements of the Lagrange Quadratic type. The mesh size is as high as

we can possibly handle on our Pentium®IV PC equipped with 1GB physical memory.

98

 

' A
:51! 21:!»

, 1 . u
vyfifi-me‘é

0';
1:41 1 1‘
‘ ‘ ‘ u "v x-"vffir" my §5 23'3““

7.
633541“ '4‘4:.':'$.nflvl::.’ ‘ '
i . 1:0"..v‘y 7r u, '575‘ ‘
4 (1'qu 'w‘y‘, '10.}. “"0,“ fi

“.1 ram; 2.13.1. :‘xmr
' ‘ ‘l' " mam W." V
Jug: V n. ‘
.

. ’5
an; 1 sex.
zav‘p. .,

"en: 7

«V :' A

>.‘
»

v a :-
;;..;v 25$
”A a

EVAVAVAV
A
AV

«'0 1 hit * 5
0 0 '
‘ 0216.14
15 .

I
v
n

'

.n»“ 4‘
, a
I"

A
.vA
'fi
:-
‘ fi
«‘23
t‘ '4
’8'4'
uv- &

9
1

A ‘V
Anva‘é‘4»
‘ AVAVAVJ a.
A o
'4' >0
.745
3‘
'3 ,
55v

1
LVA')‘

4 5‘ .
’4..."
_.

c a
n”

‘

. '0

I“

«'1! ,

it'- t
w
0.; '1

V

"V1.7

‘1‘

.-

My‘u‘

AV

4:7 7
A

«Y‘g 7 k

4'

.
‘V

.5!
1L.

‘
5
v .
A"
‘ v"
..‘u
VA‘
VA"
' fie" ‘
o 4
.5
AV
AVAVA‘V

..,J,
“'5.
a"; .
Vth,
v»

S r
b 3 '13 ='
«.9 .
‘~V>
e:
i b.
‘95!
Ar
'1 1
’4'
V‘v
47

Q h I 1

' 5 '1‘” .
on" '
’ ‘yp '1'1'1'

' 1
5‘5
- E’s
V}
C
"3.
v
‘ v
'5'.
'-
v 5‘
.4 .
1: .v
,VAV
'e'e
33:“:
..

I

'46‘1‘ ‘
4

.
(
‘w ;'a ‘
,. a-
3%
Av

a'
v
'4' (D
V) 5!:
i, ‘t
l- ‘L
‘5:
”0.-
V 1‘
‘ 9c:

I
. 5
Z V
' V

I"
Y t "as?
V iégfivAvg-géggsz; .
6 ”1:1
$52,133,145
2510442345

1‘”...

 

 

 

Figure 5.4: Finite element mesh in a neighborhood of the cylinder. Additional artifi-
cial boundaries have been enforced to adjust the mesh density

To utilize the Fast Fourier Transform, we sample the circle r = r’ uniformly
at 64 points to get E (r’,6’ = 27rk/64; jw). The numerically calculated V (r,6; jw)
is compared with the analytical value given by Eq. (5.57) and its high frequency
approximation (5.58) is shown in Figure 5.5. The results agree very well.

If only one observation point (NT = 1) is used to determine the source location,
the estimation error in the distance between the point source and the observation
point (r, 64) is 0.16 percent in the noise free case and 1.36 percent for the 10 percent

additive noise case [37].

99

 

 

 

 

 

 

 

 

I I Ix:Xx I I I
1....X
d1
1
o
1.
o
O
O
. . 1
13¢...l L l l l ..hAK
n 32 64
I I 1‘. I f I
P d
r -1
q
q
:00 o. .0..E
. O
Q 0
. .- 1
*- .... 0..
1 1 1 L H 1

Figure 5.5: V (130"; jw), Upper: real part; Lower: imaginary part. (Solid line: ana-
lytical solution given by Eq. (5.57); ’ x’: High frequency approximation given by Eq.
(5.58); ’0’: Numerical solution constructed with Eq. (5.56).

100

 

CHAPTER 6

EXPERIMENTAL VERIFICATIONS

The interface removal algorithm discussed in Chapter 5 can be easily adapted to
a planar configuration. In such a configuration, line-type current sources are placed
under an infinite conducting plate. The interfaces between the plate and the back-
ground medium can be removed with a linear filtering process. Then, a localization
algorithm is developed to recover the position of the sources. The advantage of the

planar configuration lies in the fact that it is easier to be implemented in a laboratory.

6.1 Interface Removal for a Planar Configuration

 

 

 

y}: Zyp
Al)" 3 a My
_ 00
01 0 ________ l _________ yzyh
------------------ y=yp
00 Z ‘1 00 2“ ;x
0 0
.10 o—Io 1h .10 ._I0
01-0 .(——>I
a 0.0

 

 

Figure 6.1: Geometry of the planar configuration for interface removal. Left: in
presence of the plate; Right: inside an unbounded homogeneous medium

Figure 6.1 shows cross-section of the plane configuration, where an infinite con-

101

ducting plate is placed on the X Z plane. Sinusoidal current sources are located on the
line y = —h and extended along the z direction to infinity. The thickness of the plate
is d. The conductivity of the plate and the background medium is 00 = 2.59 x 107
(S/ m) (Aluminum alloy 3003, temper H14) and 01 = 0 (air), respectively.

In this configuration, the magnetic field has no 2: dependency, allowing us to
represent B = 8,, (x, y)’i; + By (:22, y)iy. In addition, the derivative of the field with
respect to the z direction will be zero. Let E(kx,y) = F; {B(:c,y)}, where F1.
denotes the spatial Fourier transform with respect to :17, k; is a spatial Fourier domain
wavenumber and an overbar indicates the variable is in the spatial Fourier domain.
Then, Em, indicates the y component of magnetic flux density field excited in presence
of the plate while Eyh indicates the field resulting from the same source excitations
but with homogeneous background media of 00. The linear relation between Eyh and

Em, is given by [62]:

ml

yp (km: yp) = 4exp [J'Ao (yr. + ’0] exp [- lkxl (yp - h - d)l
yh (km: yh) (1 - J' lkxl />\o)2 exp (J'Aod) + (1 +J' Ikxl /)\0)2 exp (-J'>\od)’

 

(6.1)

mil

where A3 E —jwu000 — kg, yh and yp is the position of the line of measurement for the
homogeneous media problem and for the planar problem, respectively. The interface
removal filter (6.1) is independent of the form of source excitations, which allows us to
use arbitrary source excitations located on the y = —h plane. The stability condition

of the filter is yh Z yp [62].

102

 

 

Signal

Generator

 

Techron
Current

 

l

 

Impedance
Matching

 

:1) Multimeter

 

 

 

 

 

 

 

 

Amplifier

 

 

 

 

 

 

 

 

/

 

\Q)

 

 

 

Lock-in ’
Amplifier SensorU / r
/ / Y

Aluminum Plate

[
L/ '2
Wire

 

 

Motion
C ontrollcr

 

6:9

 

 

 

\

 

 

Figure 6.2: Experimental Set-up

6.2 Experimental Set-up

The current source used in this experiment is a two-wire system. These two wires
are parallel to each other and carry alternating current in opposite directions. A
304.8mm x 304.8 mm (12 inch by 12 inch) aluminum plate of thickness 2.54 mm (0.1
inch) is placed 4.96 mm above the plane of the two-wire system. The length of the
parallel segment of the Wires is 635 mm (25 inch). The wires are twisted together
rest of the length and connected in series with a impedance matching circuit (0.59
plus 0.5 mH) to a Techron 8604 current amplifier (GE). A motion controller is used
to move a magnetic field sensor (oriented to detect the y component of the magnetic

field) along a line of measurement above the plate and perpendicular to the wires.

103

The outputs of the magnetic field sensor is measured using a lock-in amplifier.

This configuration can be described with a two dimensional model. In such a
model, the parallel wires acts as two point sources. A two-wire system was adopted
instead of a simpler system consisting only one infinite wire along the z direction for
practical reasons. A two-wire system is less sensitive to the environmental noises.
Also, it generates a stronger magnetic field, which is important for maintaining a
high signal to noise ratio.

The magnetic fields generated above the plate is in the range of zero to 0.15
Gauss (magnitude) when the excitating frequency lies between 100 and 6000 Hz. Two
magnetic field sensors are used: an inductive coil and a GMR (Giant magnetoresistive)
magnetic field sensor (NVE corporation, AAH-002-02 model).

The inductive coil used in this experiment includes 106 turns of wires. The diam—
eter and thickness of the coil is 3 mm and 1 mm, respectively. The relation between

its output voltage and the applied magnetic field can be expressed as:

Vea-
Hem, = kc f’ x 104, (6.2)

 

where ch'z is the output voltage in Volts and Head is the applied magnetic field in
Gauss (both in rms values), f is the excitation frequency in Hz and he 2 554.48 is a
calibration constant determined experimentally using a well-calibrated inductive coil.
The GMR sensor is sealed inside a standard SOIC8 package of size 3.91 by 4.9 mm.
The dimensions of its active area is about 100nm by 200nm and is located in the
middle of the package. The sensor generates a 10 mV output at 0.25 Gauss with a

4.96 V DC supply. It offers 96% linearity from 10% to 70% of full scale with excel-

104

 

lent directivity characteristic. The GMR sensor is also characterized with omnipolar
outputs, which means it provides positive outputs for both directional positive and
negative applied fields. In order to obtain bipolar output, the GMR sensor was biased
with a small permenant magnet during the experiment. Unfortunately, the offset DC
fields generated inside the sensor could not be controlled very well due to memory

properties of the sensor. The field is determined using

HGA/IR : anGMRa (63)

where VGMR is the output voltage in mV and Home is the applied magnetic field in
Gauss and kn = 1 / 42.73 is the calibration constant.

Due to their small dimensions, both the GMR sensor and the inductive coil can
be approximated as point sensors. The inductive coil is quite simple to use and
is characterized by high linearity. However, it is costly and less sensitive for low
frequency fields. The GMR sensor is inexpensive and provides high sensitivity without
significant frequency dependency. Both sensors were employed and the repeatability
of their performance has been tested.

Several aspects of this experimental set—up have been investigated to reduce possi-
ble error sources. (1) The transverse dimension of the conducting plate is sufficiently
large so that the effect of the finite dimensions of the plate is negligible; (2) The
output signals of both magnetic field sensor were measured with a lock-in amplifier in
the differential mode to exploit its high common mode noise rejection capability (up
to 100 dB); (3) The environmental noise was reduced; (4) The current fluctuations

were reduced to a minimum by using constant current mode of the Techron current

105

 

amplifier. The magnitude of the current was recorded and used to account for the

effect of current fluctuations; (5) Efforts were made to reduce positioning errors.

6.3 Experiment Results

we used 100 Hz and 6 kHz as the excitation frequency in case of the GMR sensor
and the inductive coil, respectively. The lift-off is 3.81 mm (0.15 inch) for the GMR
sensor and 1.8 mm for the inductive coil. The magnitude / phase pair or real / imaginary
components of the field By was recorded at each measuring point. To set a reference
for recorded phase information, the field at center between two wires (:1: = O) was
measured experimentally and the phase at this point is set to 7r.

In an effort to check accuracy of the experimental data, a two-dimensional nu-
merical model was developed using a commercial finite element analysis package
F EMLAB®. The agreement between experimental data and results of numerical
simulations is very good.

In case of the GMR sensor, three data sets were obtained. The first data set, ‘GMR
I’, recorded measurements in the form of magnitude/ phase pairs while other two sets,
i.e. ‘GMR II’ and ‘GMR III’, recorded measurements in the form of real/ imaginary
component pairs. Data sets ‘GMR I’ and ‘GMR II’ were obtained using the same
magnetic bias level while ‘GMR III’ use a different bias level. The numerical simula-
tion results and experimental data are shown in Figure 6.3 through 6.7. The ‘GMR
III’ data shows better match with the numerical results suggesting that the bias level
is better suited.

In case of the inductive coil measurements, two data sets were obtained, called

106

r3

 

|
Wm ..- -. an...

‘Coil I’ and ‘Coil II’, respectively. Both data sets were recorded in form of magni-
tude/phase pairs and shown in Figure 6.8 through 6.12. We notice that mismatch
in phase which arises when the magnitude of the fields are lower than 0.002 Gauss
(which corresponds to 20 MV output of the inductive coil).

The experimental results obtained can be used as inputs to the linear filtering
process expressed by equation (6.1). By selecting yh 2 yp, the fields of homogeneous
media By), (Lg/,1) are obtained and compared with analytical solutions (Figures 6.13
through 6.18). Based on knowledge of the excitation fields in homogeneous media, it is
much easier to solve for the locations of the current sources. Figure 6.19 through 6.22
show the error functions in the localization of the current source for recovering the
source location (:ta, —h). Assume h is known, the relative positioning error is 3.28%
using the GMR sensor data and 0.13% using the inductive coil data, respectively.

The expected value of a is 0.01582 m.

107

 

 

0.1 . .

.0
o
01

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

 

   

GMR I

O GMR ll

0 GMR Ill
numerical

 

 

 

 

 

 

 

—015

“—0.1 0.05

—0.05

Figure 6.3: Real part of Byp (GMR sensor)

 

 

 

 

 

0.1

 

 

 

 

 

0.025 . , .
3 g - GMR I
(,1. o GMR u
0.02 ’1-"9 .............. ' . D GMR 'II ‘
v 3‘ O O .
f . numerical
g ..
0 II C"
v 001 ........................ | ,- ..... I ................ 1' ...... -
'_.‘ J r- 3 '
a
E 0' ',
E 0005 .................. J. ................. ........ 0 . . . .
_ ' c 3 ’4 .
. ;
IV 1 :
C. I Q .
.O ........... f‘ , . . A, .......... l'r‘a .............. ..
o 16:; Z i “:7.
01' . . f3 . “.33...
....III. : - . .il .....
0.05

-0.005

—O.1

-0.05

Figure 6.4: Imaginary part of Byp (GMR sensor)

108

0.1

3.5

 

2.5 r

phase, (radian)

 

     
 

I

 

° GMRI
o GMR II
D GMR Ill

 

numerical

l

 

 

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

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

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

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

.........

 

 

0.16

Figure 6.5: Phase of BW (GMR sensor)

 

 

0.14

l

.o
—L
N

.0
—L

magnitude, (Gauss)
p o
O O
O) 0)

.o
o
A

 

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

[C ‘J
1 (1
1: t
...... q, .

:1
000‘

 

GMR I
GMR ll
GMR Ill

 

 

 

numerical q

 

 

 

Figure 6.6: Magnitude of BW (GMR sensor)

109

 

 

 

- GMRI 5. 3 i
o GMRII _______ 0. ‘. ,_
D GMRIII : ~
numerical

 

  
   

 

 

 

 

 

 

    

—0.1 5 —0.1 —0.05 0 0.05 0.1
Re[Byp], (Gauss)

 

Figure 6.7: Complex plane plot of By? (GMR sensor)

 

 

0.01 T I I

 

 

 

 

:4. 3 2;. . Coil I
0 0 w “ O CO“ II
. I~ . \ (I o
-_ : . numerical
0005. ..... L ........ . .
I I! u I n "C
a v'
g 0 0
(U 0 .............................
22. 0
In ‘.
>
0
£_0.005 _. .......................... I. ....... a
(D .a
(I . j, .,
0‘? 0"“!

_0'o1__,. ., ,. ............. ........ ., .. .. .-

 

 

 

-0015 5 3
-01 —0.05 0 0.05 0.1

Figure 6.8: Real part of BMD (Inductive coil)

110

lm[B ], (Gauss)

phase, (radian)

  
   
 
   
  
 

 

I

0 Coil l
O Coil ll ‘
numerical _

 

 

 

 

 

YP

 

 

 

 

 

 

 

_1 ..................
_2 - g
_3 ................... g
_5 . I. .
—0.1 —0.05 O 0.05 0.1
x, (m
Figure 6.9: Imaginary part of BW (Inductive coil)
4 ‘. I ‘
‘ : 0 Coil l
O Coil ll -
numerical

 

 

 

 

 

 

 

 

 

 

—0.1 —-0.05 0 0.05 0.1

Figure 6.10: Phase of Byp (Inductive coil)

111

magnitude, (Gauss)

1 : 0 Coill
O Coil ll
numerical -
O ‘ r
-0.1 —0.05 0 0.05 0.1

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6.11: Magnitude of Byp (Inductive coil)

 

 

  
 

 

 

 

 

 

 

 

—3
10
5 x T T r
p . Coil l
4 O Coil ll
3 --- numerical
A 2 .............. I .....................................................
U)
8
m 1 ... ......................................
(D
‘-: 0 P. ........................
”a
a _1 .. .......................................................
E
_ _2 _,
_3 _ ........... .012) . .. 0O
_4 ............. é. & ................................. , _,
05015 0 01 0 i 5 0 0 005 0 01
' ' Re[%0yp], (Gauss) ' '

Figure 6.12: Complex plane plot of By? (Inductive coil)

112

0.015

0.01

O
O
O
01

yh

Fie[B ], (Gauss)

—0.01 --

—0.01 5

 

 

 

 

 

 

 

 

 

T I I f I
' 0 GMR I
O GMR III
P" ............................. ana'ytical q
’1“ ’0) - .
1") 0‘“ I
l' ’1‘
n . . . " . {g “ ...... _,
v 0'. i5 ‘
i \ . 0 .“
i . '1‘. 0 0 . ‘
l \ g" ‘0‘: . . . . .‘. ................. ‘ .
"H. k I i 9‘ 'au,
_. . .3 -
i E u f
i i \v I
'. I 0 I
I. ..... ., ...~....“0 ................ a
: 3 s I 2
s 2 ‘i s 2 2
l l I l l
—0.1 —0.05 0 0.05 0.1 0.15

Figure 6.13: Real part of By), (3:, yh) (filtered from GMR data)

0.05
0.04
0.03
0.02
0.01

th

0.01

lm[B ], (Gauss)

—0.02
-0.03
—0.04
—0.05

Figure

 

 

 

 

 

 

 

 

 

6.14:

113

I I
g . GMR I
............................... ............. O GMRIII 1
_. ............ ..... 3 ..... ., ...... analyt'ca' .
'0 . ‘5‘ .
t‘ ‘ \
... ........ fl 0. . . .Z. . . . '. ‘5 ~—I
(- .
.. t, . u ....... CD ‘\ ........................ .
..t
) \
'lo ‘4‘": . :1 "a, ................... \ "‘ “‘ """

' 0 . m- .
...................... .
........ .. 1'. .1 fl .
t ..... f ..... E ......... I. , a
_ i'k' ......... ............ ......... _

: 2 3 2
; g g 1 ;
—O.1 -0.05 0 0.05 0.1 0.15

Imaginary part of By), (:13, yh) (filtered from GMR data)

 

0.04

 

lm[Byh], (Gauss)

 

 

 

 

 

 

 

003 ,.. . m ..... m .1. ........ ....................... ..
.-. . . ‘ o . : "a .

0.02-9.6. . ................. '0'. ........ .
. . . o

001""~ .LO .................................. .‘:’.r. _

. I . ‘ O. I

o o,‘

0- ....." g9»... ” ‘J 1" ......... .

O < "

_o_o1-.-. .., ................. .
_0.02 .. . GMR ' ..... . .............. , ........ Q. .0 .......... ,
o GMR III g 3 . § 0
_0.03 -. analytical .................. ...: ..... .: ...... 9 .......... .

. I . I '
_0.04_...........;...6.....Q.’......<; ..... ' ................ _,
—0.05 i i

—0.01 5

—0.01

-0.005 0
Re[Byh], (Gauss)

0.01

Figure 6.15: Complex plane plot of By), (3, yh) (filtered from GMR data)

 

 

 

 

 

 

 

 

 

x 10'6
Co ; S5
4 ... . .C. 0 . . .,
A .'
g 2 ...................... ........ 3‘ ...... .... ~-1
I q, . \
((3 .l . 0 C 0
v 0 “ ID
a .. "‘1 ............... V . . . A.
mg 0‘. " " 0 0 D II ‘0'
'5‘ 2 u i u '0
a: 4. ..
" : 0 _
g 0 COII l
_4 ....................... “A ....... i ........ . .. O Coil ll 4
. . u
o 5 analytlcal
_ : <53
—0.05 0.05
x, m)

Figure 6.16: Real part of By), (cc, yh)(filtered from inductive coil data)

114

 

lm[Byh], (Gauss)

lm[Byh], (Gauss)

 

 

 

 

 

 

 

 

 

x 10—
6 I
. ' 0: 0 Coil l
v. 0 O CO" H
4 _. .m 0 .......... 0 g.. analytical ..
0
0
.4 q "
2 ... . . . . . . . . . . . ....... ..,
u
0 0
.. 0 I
a 4
.t. ’0 Jo ..... a. 1".“ ...... 1. . . . n
0 ’ 0
Q!
II .
_2 _ .............................................. «0 . _.
‘1
O 1|
D
C
_4 -. ..... . u . ........ .. . . . _.
0' " 0
G ©
: O
30.05 0.05
x, m)

 

 

 

 

 

 

 

 

 

—6
10
6x ' Y I ' '
: . Coill
f 0 Coilll
': analytical}
f i f o i
i i i be
_4 i d
; 200 :0 a. 8'5
: i ' C? . 05
—66 5. é g é ti is 8
_ — _ Re , Gauss _
[ Yb] ( ) X10 6

Figure 6.18: Complex plane plot of By), (:13, yh) (filtered from inductive coil data)

115

 

      

 

 

 

(\” /a

0 0.005 0.01 0.015 0.02 0.025 0.03
a. (m)

 

 

 

Figure 6.19: Error function used for recovering the source location (—a, a), unit: m,
(contour plot, ‘GMR I’)

-5

 

 

 

 

 

 

 

 

 

 

 

 

x 10 .
., . :
5 E
) ,._ .
' . . , ~~~k~
C 5 ~:~:~:~.
o I I I I
:- ‘::::555
O ‘ ' I'M:
c 4 “I
a
5 l
t 3
LIJ
2

 

  

' 0.03
0.02

 
 

—0.02 0.01

"a! (m) a! (m)

Figure 6.20: Error function used for recovering the source location (—a, a), unit: m,
(mesh plot, ‘GMR I’)

116

 

a. (m)

Figure 6.21: Error function used for recovering the source location (—a, a), unit: m,
(contour plot, ‘Coil I’)

  
     
    
  
     
    
    
     

 

  
   

S\\\\\ H‘ \\ I
3“\\l“\\\\‘“ . . ”7"”

  

 
   

  

Error function
A
A
(II

 
 
  
 
 

   

 

0“ .
— 11” fill/ll
4-4 ‘ i33"33 ~ ‘ ”3'
:: {3, o
_ =‘13'3'3‘3°3°3:3:3:3:;2;;=’ 1"",51"
4'35 \\\‘:‘ :zwflozdl/I/ gig?
“\t‘: ::::::::;:II II// .
4.3 ’ r" "5”
4.25

0.03

   

—0.01 0.02

—0.02 0.01

-a, (m) a, (m)

Figure 6.22: Error function used for recovering the source location (—a, a), unit: m,
(mesh plot, ‘Coil I’)

117

CHAPTER 7

CONCLUSIONS

Q-transform based time—of—fiight extraction methods are discussed in this disserta-
tion. The feasibility of the methods were demonstrated through the use of canonical
inverse problems in eddy current testing both in the time and frequency domain.
The basic methodology used involves construction of a fictitious wave field from
ECT measurements using the Q—transform and estimating the time-of-fiight from
the constructed fictitious wave field. The time-of-fiight estimate can then be used for
source/ defect localization and potentially be fused with time-of-flight measurements
obtained from the wave field directly. Numerical simulations based on the finite el-
ement method were implemented for validating the concept. Experimental results

validating the concepts were also obtained.

a In the time domain, time-of—flight associated with (fictitious) wave front can
be determined using carefully designed excitation signals. The approach was
illustrated by studying the scattering from a small anomaly embedded in an
infinite homogeneous medium. The position of the anomaly was estimated
using measurements of diffusion domain responses with both scalar and vector

formulations.

o In the harmonic domain, the existence of the time-of-fiight information in non-
local measurements of the Fourier coefficients was proved. A linear filtering
method was developed to remove the material discontinuities for an long con-

ductive cylinder embedded in an unbounded homogeneous space and the po-

118

sition of a long cylindrical hole inside the conductor was estimated through

measured time-of—flight information.

o In the experimental work, two current carrying wires are placed under a con-
ducting plate. Measurements of the magnetic flux density field (component
perpendicular to the plate) was obtained using a GMR sensor as well as an
inductive coil. These experimental measurements agree with theoretical pre-
dictions very well. The effectiveness of the interface removal algorithm is then
demonstrated using experimental measurements. The experimental data was

also employed to recover the location of the source currents.

This work summarizes the progresses made using Q-transform based methods in
recent years. The approach presennted herein, we believe, is the first phenomeno-
logical approach to extracting time-of—flight information from diffusion based NDE
techniques.

This research has presented a useful paradigm for exploring the time—of-fiight ex-
traction approach for potential NDE applications. Several canonical problems have
been systematically studied, and more importantly, basic methodologies and strate-
gies for solving these problems have been established and summarized. Capitalizing
on this attractive paradigm, it will be beneficial to explore new applications of prac-
tical interest. This study prompts future work to identify and study other NDE

problems. Possible future work may include:

0 Identifying TOF information from 3D complex structures based on real indus-

trial applications;

0 Localizing the positions of sources/defects, or mapping material parameters,

119

 

through advanced imaging technique developed for wave propagation applica-

tions;

0 Fusing retrieved TOF information with measurements obtained from wave prop-

agation based NDE techniques.

120

APPENDIX A

Q-TRANSFORM PAIRS AND THE Q-TRANSFORM

A.1 Q-transform Pairs

Here we give a concise list of some useful Q—transforrns pairs. Detailed derivations

of Eqs. (A.1)-(A.7) were presented in [6, 7, 16, 32].

QfiHWlUk=5fil

Q {sin (wq)} (t) = wexp (——w2t)

Q {cos <qu (t) = g [1 — 2WD (1. g —w2t)]

q2n+1 . _ tn
Q{(2n+1)!}(t) _ E

Q {ii [sin (wq) cosh (wq) — cos (wq) sinh (wq)]} (t) = sin (2w2t)

 

Q {i5 [cos (wq) sinh (wq) + sin (wq) cosh (wq)]} (t) = cos (2w2t)

Q{H<q — cm} (t) = —1—-exp (i)

 

121

(A.10)

Q {m - (10)? H (q — (10)} (t) = —2q0erfc (in) + 41¢; exp (%) (A.11)

3 - 2
Q {q3H (q - q0)} (t) = 6terfc (%) + ill-+73%! exp <—gfl) (A.1?)

Q {(q — 00)3 H (q - (10)} (t) = (6t + 3613) erfc (5%) - 6(10\/;eXP(-§§) (A-13)

Note that go is a nonnegative constant in Eqs. (A?) through Eq. (A.13) and

erfc () is the complementary error function [54].

A.2 The Q-transform

Q-transform is also an invertible mapping between solutions of parabolic equations
and those of hyperbolic equations, which is slightly different with Q-transform [14, 10].

Consider partial differential equations (PDEs) of the parabolic type

Btu (x,t) = Ev (x, t) + f (x, t) (A.14)
v (x, 0) = d) (x) (A.15)
an?) (x,t) = ’t/J (x, t) on S (A.16)

and of the hyperbolic type

aqqu(x,q) = £u(x,q)+g(x,q) (A.1?)

u(x,0) = q>(x) (A.18)
a,u(x,0) = 0 (A.19)
8nu(x,q) = go(x,q)onS, (A20)

122

respectively. Here L is a uniformly elliptical operator with continuous coefficients
depending only on the variable x, n is a conformal to S. g(x,q) and (p (x, q) are
smooth functions of their arguments and grow not faster that C eaq at q —-> oo. Assume

that conditions

f(x,t) = «#7? [0°09 (x,Q)exp (£15) dq (A21)

and

((x,t) = —1— 0.80 (x,q) exp (“2) dq (A22)

are satisfied. The solutions of Eqs. (A.14)-(A.16) and Eqs. (A.17)-(A.20) are related

with the Q-transform, which is defined as

at) = em (1)} (t) s —1——/°°u (q) exp ("3) dq. (A23)

We notice that function

a (t, q) = _exp (-3) (A24)

is the solution of the equation of heat conductivity 6', = qu and

2(t) = mum} (t) = Zita {121?} (t). (A25)

As an alternative to the Q-transform, Q-transform could be useful for data fusion
tasks once its mathematical properties are fully explored.

Corresponding to formulas presented in section 3.2 and A.1 for the Q-transform,
important mathematical properties of the Q-transform are derived and listed below for

reference. Similarities between Q-transform and Q-transform can be easily identified

123

(see Eqs. (A.2)-(A.3) and (A.32)-(A.33) ).

 

 

de{ud(tq)} (t) _ fimqaw _ genie)»
ewe» = Jig) + Q(u(q)}
é{sin (MM Z 3%q) (1; ; _wzt) : _—Q{COS (mm + _%

Q {cos (wq)} = exp (—w2t)

Q{q"} = Q{ gm}

 

n+1
~ 273 (277’)! n
Q{q }=T)!_t n20,n€Z
é{q2n+1}=_n!_(4t)n+1 n>0n€Z
2./7Tt — ’

124

(A26)

(A27)

(A.28)

(A29)

(A30)

(A31)

(A32)

(A33)

(A34)

(A35)

(A.36)

BIBLIOGRAPHY

[1] Cartz, L.. Nondestructive Testing, ASM International, 1995.
[2] Gros, X. E, NDT Data Fusion, Arnold, 1997.

[3] Chan, S. C., Udpa, L. and Udpa, 8., “Design of an NDE Integrated Data
Acquisition System (NIDAS)”, in Review of Progress in Quantitative Nonde-
structive Evaluation, Edited by Thompson, D. O. and Chimenti, D. E., New
York: Plenum Press, vol. 22, pp.585-592, 2002.

[4] Udpa, L. and Lord, W., “Diffusion, Waves, Phase and Eddy Current Imaging”,
in Review of Progress in Quantitative Nondestructive Evaluation, Edited by
Thompson, D. O. and Chimenti, D. E., New York: Plenum Press, vol. 4,
pp.499-506, 1985.

[5] Lee, K. H. and Xie, G., “A New Approach to Imaging with Low Frequency
Electromagnetic Fields”, Geophysics, vol. 58, pp.786—796, 1993.

[6] Lee, K. H., Liu, C. and Morrison, H. F., “A New Approach to Modeling the
Electromagnetic Response of Conductive Media”, Geophysics, vol. 54, pp.1180-
1192, 1989.

[7] Ross, 8., Lusk, M., and Lord, W., “Application of a Diffusion-to-Wave "Hans-
formation for Inverting Eddy Current Nondestructive Evaluation Data”, IEEE
Trans. on Magnetics, vol. 32, pp.535-546, 1996.

[8] Bragg, L. R. and Dettman, J. W., “Related Partial Differential Equations and
Their Applications”, SIAM J. Appl. Math., vol. 16, pp.459—467, 1968.

[9] Bragg, L. R. and Dettman, J. W., “Related Problems in Partial Differential
Equations”, Bulletin of America Society of Mathematics, vol. 74, pp.375—378,
1968.

[10] Romanov, V. G., Inverse problems of mathematical physics, translated by Yuz—
ina, L. Ya., Utrecht, The Netherlands: VNU Science Press BV, 1987.

[11] Lee, K. H., Xie, G., Habashy, T. M. and Torres-Verdin, C., “Wave Field Trans-
form of Electromagnetic Fields”, presented in 64th Ann. Internet. Mtg, Soc.
Expl. Geophys, Expanded Abstracts, pp633-635, 1994.

[12] Bragg, L. R. and Dettman, J. W., “An Operator Calculus for Related Partial
Differential Equations”, J. Math. Analysis and Applic, vol. 22, pp.261-271,
1968.

[13] Bragg, L. R. and Dettman, J. W., “A Class of Related Dirichlet and Initial
Value Problems”, Proc. Amer. Math. Soc, vol. 21, pp.50—56, 1969.

125

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

Reznitskaya, K. G., “The Connection Between Solutions of the Cauchy Problem
for Equation of Different Types and Inverse Problems”, Mat. Problemy Geofiz.
Vyp., vol. 5, Part 1, pp.55—62, 1974. (in Russian)

Lavrent’ev, M. M., Romanov, V. G. and Shishatskii, S. P., Ill-Posed Problems
of Mathematical Physics and Analysis, America Mathematical Society, Provi-
dence, RI, 1986. (translated)

Sun, K., Udpa, S., Udpa, L., Xue, T., and Lord, W., “Registration Issue in
the Fusion of Eddy Current and Ultrasound NDE Data Using Q-transform”,
in Review of Progress in Quantitative Nondestructive Evaluation, Edited by
Thompson, D. O. and Chimenti, D. E., New York: Plenum Press, vol. 15,
pp.813—820, 1996.

Hildebrand, B. P. and Fitzpatrick, G. L., “Inversion of Eddy Current Data
Using Holographic Principles”, in Review of Progress in Quantitative Nonde-
structive Evaluation, Edited by Thompson, D. O. and Chimenti, D. E., New
York: Plenum Press, vol. 4, pp.507-515, 1984.

Bohbot, R. de. 0., Lesselier, D. and Duchene, B., “A Diffraction Tomograo-
hic Algorithm for Eddy Current Imaging from Anomalous Fields as Fictitious
Imaginary Frequencies”, Inverse Problems, vol. 10, pp109-127, 1994.

Devaney, A. J ., “Linearised Inverse Scattering in Attenuating Media”, Inverse
Problem, vol. 3, pp389-397, 1987.

Zorgati, R., Duchene, B., Lesselier D. and Pons, F., “Eddy Current Testing of
Anomalies in Conductive Materials, Part I: Qualitative Imaging via Diffraction
Tomography Techniques”, IEEE Transactions on Magnetics, vol. 27, pp.4416—
4437,1991.

Zorgati, R., Lesselier D., Duchene, B. and Pons, F., “Eddy Current Testing of
Anomalies in Conductive Materials, Part II: Qualitative Imaging via Determin-

istic and Scochastic Inversion Techniques”, IEEE Transactions on Magnetics,
vol. 28, pp.1850—1862, 1992.

Pierce, A., “Wave Methods for an Inverse Problem in Diffusion”, Inverse Prob-
lem, vol. 2, pp.205—217, 1986.

Kunetz, G., “Processing and Interpretation of Magnetotelluric Soundings”,
Geophysics, vol. 37, pp.1005-1021, 1972.

Levy, 8., Oldenburg, D., and Wang, J ., “Subsurface Imaging Using Magnetotel-
luric Data”, Geophysics, vol. 53, pp.104—117, 1988.

McWhirter, J. G. and Pike, E. R., “On the Numerical Inversion of the Laplace
Transform and Similar Fredholm Integral Equations of the First Kind”, Journal
of Physics A, vol. 11, pp.1729—1745, 1978.

126

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

Gibert, D., Tournerie, B. and Virieux, J., “High-resolution Electromagnetic
Imaging of the Conductive Earth Interior”, Inverse Problems, vol. 10, pp.341-
351, 1994.

Gershenson, M. “Simple Interpretation of Time Domain Electromagnetic
Sounding Using Similarities Between Wave and Diffusion Propagation”, .93 SE G
Annual Convention Extended Abstract, pp.1342-1345, 1993, 882.34.

Gershenson, M., “Simple Interpretation of Time Domain Electromagnetic
Sounding Using Similarities Between Wave and Diffusion Propagation”, Geo-
physics, vol. 62, pp.763-774, 1997.

Das, K., Becker, A. and Lee, K. H., “Experimental Validation of the Wave-
field Transform of Electromagnetic Fields”, Geophysical Prospecting, vol. 50,
pp.441-451, 2002.

Levy, S. and Fullagar, P. K., “Reconstruction of a Sparse Spike Train from a
Portion of its Spectrum and Application to High-Resolution Deconvolution”,
Geophysics, vol. 46, pp.1235—1243, 1981.

Lee, T. J., Suh, J. H., Kim, H. J., Song, Y. and Lee, K. H., “Electromag-
netic Traveltime Tomography Using an Approximate Wavefield Transform”,
Geophysics, vol. 67, pp.68-76, 2002.

Tamburrino, A., Udpa, S. 8., “On the Q-transform for Solving Inverse Problems
for the Diffusion Equation”, Internal Report, Department of Electrical and
Computer Engineering, Michigan State University, 2001.

Tamburrino, A. and Udpa, S. 8., “Solution of Inverse Problems for Scalar Par-
abolic Equations Using a Hyperbolic to Parabolic Transformation: Time-of-
Flight Analysis”, submitted for publish.

Tian, Y., Tamburrino, A., Udpa, S. S. and Udpa, L., “Time-of—Flight Mea-
surements from Eddy Current Tests”, in Review of Progress in Quantitative
Nondestructive Evaluation, Edited by Thompson, D. O. and Chimenti, D. E.,
New York: Plenum Press, vol. 22, pp.593—600, 2002.

Tamburrino, A., Fresa, R., Udpa, S. S. and Tian, Y., “Three—Dimensional Defect
Localization from Time-of—F light / Eddy Current Testing Data”, IEEE Trans. on
Magnetics, vol. 40, pp.1148-1151, March 2004.

Tamburrino, A., “TOF in the Presence of Material Interfaces: A Scalar Axisym-
metric Geometry”, Internal Report, Department of Electrical and Computer
Engineering, Michigan State University, 2003.

Tamburrino, A., Tian, Y. and Udpa, S. S., “Eddy Current Testing of Conduc-
tive Materials Using Time-of-F light Measurements”, in E’NDE, Electromag-
netic Non-Destructive Evaluation (VIII), T. Sollier et al. (Eds), pp. 104-111,
108 Press, 2004.

127

[38]

[39]

[40]

[41]

[42]

[43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

Haus, Hermann A, Melcher, James R., Electromagnetic Fields and Energy, New
Jersey: Prentice Hall, 1989.

Jin, Jianming, The Finite Element Method in Electromagnetics, New York:

Wiley, 2002.

Demarest, K. R., Engineering Electromagnetics, New Jersey: Prentice Hall,
1998.

Guru, B. S. and Hiziroglu, H. R., Electromagnetic Field Theory Fundamentals,
Boston: PWS Publishing, 1998.

Ida, N. and Bastos, J. P. A., Electromagnetics and Calculation of Fields, New
York: Springer, 1997.

Csendes, Z. J ., Weiss, J. and Hoole, S. R. H., “Alternative Vector Potential For-
mulations of 3—D Magnetostatic Field Problems”, IEEE Trans. on Magnetics,
vol. MAG-18, pp.367—372, 1982.

Biro, O., Preis, K. and Richter, K. R., “Various FEM Formulations for the
Calculations of Transient 3D Eddy Currents in Nonlinear Media”, IEEE Trans.
on Magnetics, vol. 31, pp.1307-1312, 1995.

Biro, O. and Preis, K., “On the Use of the Magnetic Vector Potential in the
Finite Element Analysis of 3-D Eddy Currents”, IEEE Trans. on Magnetics,
vol. 25, pp.3145-3159, 1989.

Bryant, C. F., Emson, C. R. I. and Trowbridge, C. W., “A Comparison of
Lorentz Gauge Formulations in Eddy Current Computations”, IEEE Trans. on
Magnetics, vol. 26, pp.430—433, 1990.

Morisue, T., “A New Formulation of the Magnetic Vector Potential Method in
3-D Multiply Connected Regions”, IEEE Trans. on Magnetics, vol. 24, pp.110-
113, 1988.

Morisue, T., “Magnetic Vector Potential and Electric Scalar Potential in Three-
Dimensional Eddy Current Problems”, IEEE Trans. on Magnetics, vol. MAG-
18, pp.531-535, 1982.

George, A. and Liu, J., Computer Solution of Large Sparse Definite Systems,
New Jersey: Prentice Hall, 1981.

Kraus, J. D. and F leisch, D. A., Electromagnetics with Applications, New York:
McGraw-Hill, 1999.

Zhou, Pei—bai, Numerical Analysis of Electromagnetic Fields, Berlin, Germany:
Springer-Verlag, 1993.

128

[52] Li, Y., Edge Based Finite Element Simulation 0f Eddy Current Phenomenon
and its Application to 3D Defect Characterization, Ph’D thesis, Iowa State
University, 2002.

[53] Wait, R. and Mitchell, A. R., Finite Element Analysis and Applications, New
York: John Wiley & Sons, 1985.

[54] Abramowitz, M. and Stegun, I. A. (Eds), Handbook of Mathematical Func-
tions with Formulas, Graphs, and Mathematical Tables, 9‘“ printing. New York:
Dover, 1972.

[55] Tikhonov, A. N ., Solutions of Ill-Posed Problems, New York: Wiley, 1977.

[56] Corless, R. M., Gonnet, G. H., Harc, D. E. G., Jeffrey, D. J. and Knuth, D. E.,
“On the Lambert W Function”, Advances in Computational Mathematics, vol.
5, pp.329—359, 1996.

[57] Richmond, J. H., “Scattering by a Dielectric Cylinder of Arbitrary Cross Section
Shape”, IEEE Trans. on Antennas and Propagation, vol. 13, pp.334-341, 1965.

[58] Harrington, R. F., Time-Harmonic Electromagnetic Fields, New Jersey: IEEE
Press, 2001.

[59] Kishk, A. A., Parrikar, R. P. and Elsherbeni, A. Z., “Electromagnetic Scattering
from an Eccentric Multilayered Circular Cylinder”, IEEE Trans. on Antennas
and Propagation, vol. 40, pp.295—303, 1992.

[60] Schwart, L., Mathematics for the Physical Sciences, Paris: Hermann, 1966.

[61] Parrikar, R. P., Kishk, A. A. and Elsherbeni, A. Z., “Scattering from an
Impedance Cylinder Embedded in a Nonconcentric Dielectric Cylinder”, IEE
Proceedings-H, vol. 138, pp169—172, 1991.

[62] Tamburrino, A. and Nair, N., “Interface removal algorithm for a planar prob-
lem”, Private communication, 2005.

129