r: a
. ,tfinnuu V
1 3;. .vai
pan unfrrhka .

3. ‘ $1.21.

“an. ,

Lu. 4. 1.. u. .6
$3.31.)? «n ‘2.
'N‘ Pigl- ‘5‘ ‘h n
. = .. x. as. it.
. ”33.9.... in s . I. m
.1135... r... n: .i tilt. 3
« .
fimafl.‘
.21! vs: 3?: in”?!
.A‘ a z. t I
.6 ‘9
3.5 5.43.! In l. ‘.
0":- i‘l“ I

.. .
._. 2&9

.
. .
lit 2 i

. ,
55.1.1.1 20.1.5333: I I
((-5 I?‘ u

, s”
pl “} .i-m-IW‘YIU

tilt.
’05:}...
9.019.453»..-

.ifi.

.‘J...’.~}.
54“}... 2:...
111,051....
I!

:3... f.v.v
:4 ti
’il.

 

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the
dissertation entitled

FORWARD AND INVERSE MODELING USING MESHLESS

PhD.

METHOD FOR NDE APPLICATION

presented by

XIN LIU

has been accepted towards fulfillment
of the requirements for the

degree in Electrical Enweering

 

 

041$: 4am!

 

Major Professor’s Signature

’25? owéufl

Date

MSU is an Affirmative Action/Equal Opportunity Employer

 

 

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

oifiltzl 59 2%

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

908 K:IProiIAoc&Pres/CIRCIDatoDm. hdd

 

A_

FORWARD AND INVERSE MODELING USING MESHLESS METHOD
FOR NDE APPLICATION
By

Xin Liu

A DISSERTATION

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

DOCTOR OF PHILOSOPHY
Electrical Engineering

2009

ABSTRACT

FORWARD AND INVERSE MODELING USING MESHLESS METHOD
FOR NDE APPLICATION

By
Xin Liu

Nondestructive Evaluation (NDE) methods are used extensively in industry in
inspecting and maintaining product quality. As new NDE sensors and systems are
developed for new applications, the availability of a theoretical mOdel that can
simulate the system performance is extremely important for understanding the
underlying physics and optimizing the system parameters. For most real world
problems, where analytical solutions are not available, various computational
methods have been developed. Among these, Finite Element Methods (FEM) are well
known and used extensively for computing static and quasi-static electromagnetic
fields associated with NDE applications. FEM offers a robust model, and the
simulation results using fast solvers are stable and accurate. One of its disadvantages
is the need for a computational mesh. Most of the test geometries in NDE are
complex and high dimensional, and mesh construction in FE models is an extremely
labor intensive and time consuming process. Furthermore, in electromagnetic
problems that involve geometrical discontinuities and propagating tight cracks, the
use of an underlying mesh creates difficulties in the treatment of discontinuities and

decreases the accuracy of simulation results.

A major contribution in this dissertation is the development of Element-Free
Galerkin (EFG) model, for NDE applications, which belongs to the newly developed
class of Meshless Methods that do not involve the use of a mesh for discretizing the
domain. Several improvements to the basic EFG model are proposed for
electromagnetic field computations. One, two, and three-dimensional (1D, 2D, and
3D) models for Poisson and diffusion equations, describing static or low frequency
quasi-static problems, are presented as examples to illustrate the validity of the EFG
method. The model is also applied to real world problems for electromagnetic field
calculations in aircraft skin inspection and the simulation results clearly demonstrate
the feasibility and advantage of the method.

The implementation of EFG method for solving the inverse problem is also
presented. The objective of the inverse problem is to estimate parameters
characterizing defect profile, including location, size, and shape, based on the NDE
measured signals. New formulations for model based inversion using EFG method as
forward model is developed in this dissertation. 1D parametric and 2D non-parametric
inversion techniques using State Space search and gradient-based search method are
presented along with preliminary results on simulated data. The major advantage of
BF G technique over FEM in inverse problem solutions is the elimination of
re-meshing in each iteration, which saves computational time while maintaining
accuracy of solutions.

The challenge of high dimensionality in 3D inverse problems is also addressed
by extending the iterative State Space search approach described for 2D problems.

Preliminary results validating these modifications are presented.

To my husband, daughter, and parents

iv

ACKNOWLEDGEMENTS

I would like to express my sincere appreciation to Dr. Lalita Udpa for her
continuous guidance and encouragement throughout my Ph.D study. I am grateful for
her intellectual, financial support, inspiration, invaluable advice, and all the trust that
helps me to get into and survive in the world of engineering. Dr. Lalita Udpa teaches
me everything, from the basic English writing to the attitude and passion in research.
Dr. Lalita Udpa is a lifelong role model in my professionalism.

I would like to express my sincere gratitude to my committee members: Dr. Satish
Udpa, Dr. Fang Z. Peng, and Dr. Neil Wright for taking the time to serve on my
committee and for providing invaluable comments and suggestions to enhance the
quality of the dissertation.

I would like to express my gratitude to all my collaborators for the helpful
discussion, guidance, providing data and samples. I thank each person related to this
work directly or indirectly for the successful completion of the dissertation.

I would particularly like to thank my husband and colleague, Dr. Yiming Deng for
his support, encouragement, and love. I would like to thank my beloved parents, all

my friends and colleagues in NDEL. Thank you!

f

TABLE OF CONTENTS

 

 

LIST OF TABLES -- - - -- - -- - - ...... - - - _- - - - . viii
LIST OF FIGURES ..................................................................................... ix
CHAPTER 1. INTRODUCTION ............................................ - - 1
1.1 Introduction to Nondestructive Testing ............................................................... 1
1.2 Eddy Current Testing ........................................................................................... 2
1.3 Forward and Inverse Problems ............................................................................ 5
1.4 Need for Modeling ............................................................................................... 7
1.5 Introduction to Meshless Methods .................................................................... 10
1.6 Research Objectives .......................................................................................... 12

CHAPTER 2. BASIC ELECTROMAGNETICS FOR EDDY

 

 

CURRENT MODELING 15
2.1 Introduction ....................................................................................................... 15
2.2 Maxwell’s Equation ........................................................................................... 15
2.3 Time-Harmonic Fields and Potential Functions ................................................ 16

2.3.1 Magnetic Vector Potential ............................................................................. 17
2.3.2 Electric Scalar Potential ................................................................................. 17
2.4 Weak Form and Galerkin Method ..................................................................... 18

CHAPTER 3. ELEMENT FREE GALERKIN METHOD ................. 21
3.1 Introduction ....................................................................................................... 21
3.2 Moving Least Square Approximation ............................................................... 22
3.3 Boundary Conditions ......................................................................................... 24
3.4 Weight Function ................................................................................................ 28
3.5 Discontinuities Approximation .......................................................................... 32
3.6 The Orthogonal Basis ........................................................................................ 33

CHAPTER 4. FORWARD MODEL VALIDATION - .......... 38
4.1 Introduction ....................................................................................................... 38
4.2 1D Problem ........................................................................................................ 39
4.3 2D Problem ........................................................................................................ 42

CHAPTER 5. FORWARD MODEL IMPLEMENTATION .............. 47
5.1 Introduction ....................................................................................................... 47

5.2 Mathematical Formulation ................................................................................ 47

5.3 Magneto Optic Imaging Application ................................................................. 50
5.3.1 Introduction to M01 in Airframe Inspection Application ............................. 51
5.3.2 Modeling Geometry ....................................................................................... 54
5.3.3 Numerical Results I — Fastener hole without crack ....................................... 57
5.3.4 Numerical Results 11 - Fastener hole with 5 mm second layer tight crack....60
5.3.5 Conclusion ..................................................................................................... 64

CHAPTER 6. INVERSE PROBLEM USING EFG METHOD .......... 65

6.1 Introduction ....................................................................................................... 65

6.2 1D Problems — Parametric Approach ................................................................ 67
6.2.1 Problem Statement and Inversion Scheme .................................................... 68
6.2.2 Inversion Results and Discussion .................................................................. 70
6.2.3 Conclusion ..................................................................................................... 74

6.3 Inversion based on‘State Space Search Techniques for 2D Problem ................ 74
6.3.1 Problem Representation ................................................................................. 75
6.3.2 Tree Structure ................................................................................................ 77
6.3.3 The Search Procedure .................................................................................... 78
6.3.4 Initial Results ................................................................................................. 82
6.3.5 Choice of Cost Function ................................................................................ 86
6.3.6 Conclusion ..................................................................................................... 90

6.4 Inversion based on Gradient Search using Adjoint Equation for 2D Problem..91
6.4.1 The Inversion Scheme in 2D ......................................................................... 92
6.4.2 Results ........................................................................................................... 95
6.4.3 Effect of Noise Level ................................................................................... 101
6.4.4 Conclusion ................................................................................................... 108

CHAPTER 7. 3D INVERSE PROBLEM USING EFG METHOD .. 109

7.1 Introduction ..................................................................................................... 109
7.2 Problem Statement ........................................................................................... 110
7.3 3D State Space Representation ........................................................................ 1 11
7.4 Search Procedure ............................................................................................. 112
7.5 The Detailed Search Procedure ....................................................................... 114
7.6 Initial Results ................................................................................................... 1 17
7.7 Conclusion ....................................................................................................... 125
CHAPTER 8. CONCLUSIONS AND FUTURE WORK .................. 126
8.1 Original Contributions to the Field .................................................................. 126
8.2 Future Work ..................................................................................................... 127
REFERENCES .................................................................... 129

LIST OF TABLES

Table 1 Boundary conditions ............................................................................................. 56
Table 2 Summary of simulation results ............................................................................. 57
Table 3 Summary of simulation results ............................................................................. 61

viii

LIST OF FIGURES

Figure 1.1 A generic NDE system ....................................................................................... 2
Figure 1.2 System of eddy current testing ........................................................................... 3
Figure 1.3 The system approach to illustrate the forward and inverse problem ................. 6
Figure 3.1 EFG coupled with finite elements along essential boundaries ................... 27
Figure 3.2 Domain of influence (a) circular; (b) square .................................................... 29
Figure 3.3 The weight functions over 1D domain of influence ......................................... 31
Figure 3.4 The asymmetry of the domain of influence ..................................................... 32

Figure 3.5 Domain of influence of node adjacent to material discontinuity (rectangular
domain) ...................................................................................................................... 33

Figure 4.1 The comparison of exact and EFG solution with 1“, 2nd order polynomial
basis .......................................................................................... 40

Figure 4.2 The comparison of exact and EFG solution with 3rd order polynomial basis .41

Figure 4.3 The comparison of exact solution with various order orthogonal basis ........... 42
Figure 4.4 Modeling geometry for 2D problem ................................................................ 43
Figure 4.5 Induced eddy current plot using EFG method by polynomial basis function (a)
st nd

1 order (b) 2 order ................................................................................................ 44
Figure 3.6 Induced eddy current plot using EF G method by orthogonal basis function (a)

211 order (b) comparison with Figure 4.5 (a) ............................................................ 45
Figure 4.7 The convergence rate using polynomial basis and orthogonal basis ............... 46
Figure 5.1 (a) The schematic of the M01 system (b) an experimental image ............... 52
Figure 5.2 (a) The M01 308 system (b) using of M01 system .......................................... 53

Figure 5.3 M01 simulation geometry (both for FEM and EFG). (a) 3-D view, (b) 2-D top
view and (c) side view of the multilayer plate with a 6mm diameter fastener hole ..54

Figure 5.4 2D cross-section view of (a) FEM mesh and (b) EFG discretization of the
solution domain ......................................................................................................... 55

 

Figure 5.5 (a) 3D view, (b) 2D top view and (c) 2D side view of the normal component of
the magnetic flux density with FEM and (d), (e), (f) with EF G ................................ 58

Figure 5.6 Binary results obtained by thresholding simulation results in Figure 5.5. (a)
FEM (b) EFG and (0) experimental image ................................................................ 60

Figure 5.7 (a) 3D view, (b) 2D top view, and (c) 2D side view of the normal component
of the magnetic flux density with FEM and (d), (e), (f) with EF G ............................ 61

Figure 5.8 Binary results obtained by thresholding simulation results in Figure 5.7. (a)
FEM (b) EFG ............................................................................................................. 64

Figure 6.1 Schematic representation of the iterative approach of the inverse problem. . ..67
Figure 6.2 Definition of the defect parameters for 1D crack ............................................. 68
Figure 6.3 Cost function versus (a) r, (h) z and (c) 6 ........................................................ 71

Figure 6.4 Initial guess (dotted line), the true defect (solid line) and the prediction
(dashed line). (a) crack perpendicular to the current (b) crack at an angle of 55

degrees with the current direction .............................................................................. 73
Figure 6.5 The EFG model Geometry (a) 3D View (b) cross-section view ...................... 75
Figure 6.6 Examples of States (Defects) ........................................................................... 76
Figure 6.7 the tree structure ............................................................................................... 78
Figure 6.8 (a) Input signal (b) initial guess of top surface defect (0) initial guess of bottom

surface defect ............................................................................................................. 79
Figure 6.9 Candidate defect shapes used in the search in layer 2 ..................................... 81

Figure 6.10 (a) True defect profile and (b) initial guess of a top surface defect (c) initial
guess of a bottom surface defect (d) reconstructed results for top surface defect ..... 83

Figure 6.11 (a) True defect profile and (b) initial guess of a top surface defect (c) initial
guess of a bottom surface defect (d) reconstructed result for bottom surface defect 85

Figure 6.12 True defect profile (a) top surface (b) bottom surface ................................... 88

Figure 6.13 Reconstructed top surface defect profiles (a) with old cost function (b) with
new cost function ....................................................................................................... 89

Figure 6.14 Reconstructed bottom surface defect profiles (a) with old cost function (b)
with new cost function ............................................................................................... 90

Figure 6.15 (a) The desired defect profile. (b) the initial guess ........................................ 97

Figure 6.16 The updating status of the defect profile ........................................................ 99
Figure 6.17 True defect ................................................................................................... 102

Figure 6.18 The signal corrupted with (a) 5.0%. (b) 10.0%. (c) 15.0%. (d) 20.0%. (e)
25.0%. (t) 30.0% random noise ............................................................................... 103

Figure 6.19 (a) to (f) The reconstructed defect profile (dashed lines) when the signal is
corrupted with (a) 5.0%. (b) 10.0%. (c) 15.0%. (d) 20.0%. (e) 25.0%. (f) 30.0%
random noise ............................................................................................................ 106

Figure 7. l The model geometry (a) 3D view, (b) 2D top view and (c) cross-section

view ............................................................................................ 1 10
Figure 7.2 3x3 state space in each layer with three sub-state space ................................ 113
Figure 7.3 The candidate state space for sub-state in Figure 7.2 (a) ............................... 113
Figure 7.4 Plot of estimated defect region using real, imaginary and absolute data, ...... 115
Figure 7.5 The ideal signal ofBz (a) real part (b) imaginary part ................................... 117
Figure 7.6 Initial guess of a top surface defect (a) cross-section view (b) top view ....... 118

Figure 7.7 Initial guess of a bottom surface defect (a) cross-section view (b) top view.1 19

Figure 7.8 The reconstructed defect profile from top to bottom layer ............................ 120
Figure 7.9 The cost function vs. node ............................................................................. 122
Figure 7.10 The reconstructed defect profile from bottom to top layer .......................... 123
Figure 7.11 The cost firnction vs. node ............................ 124

xi

CHAPTER 1. INTRODUCTION

1.1 Introduction to Nondestructive Testing

Nondestructive Testing (NDT) or Nondestructive Evaluation (NDE) methods are
widely used by a number of industries to control and maintain product quality,
prevent catastrophic failure and improve reliability. It can be defined as the
assessment of structural integrity of a material or component without any physical
damage it [1]. NDT/NDE plays a crucial role in assuring safety and reliability in the
operation of aircraft, gas pipeline, bridge, nuclear power plant etc.

These inspection procedures enable us to locate a flaw or defect and also provide
information about defect such as size, shape, and orientation. Furthermore, NDE can
also be used to characterize material properties. Some of the commonly used NDT
techniques include ultrasonic (UT) [2], radiographic [3], electromagnetic (ET), and
thermographic [4] etc.

The implementation of NDT includes the major components as illustrated in Figure
1.1 [5]. An input transducer is used to inject energy (excitation source) into the test
specimen. The nature of the interaction between the energy and the test specimen is a
function of several variables, including the type of the energy, material properties,
defects, inhomogeneities and so on. The receiver transducer will capture the response
of the material-energy interaction and the resulting signal is processed, recognized
and analyzed to determine the existence and the properties of a defect in the test

specimen.

f

A number of energy sources including electromagnetic, ultrasonic and
thermographic have been extensively used in a variety of applications. Commonly
used electromagnetic techniques for NDT include magnetic flux leakage [6] and eddy
current methods [7]. Eddy current methods have been widely used to detect cracks,
seams, pits and other surface and subsurface flaws in conducting test specimen. This
dissertation mainly focuses on electromagnetic methods which are based on principles

of eddy current testing described in the next section.

Source

  

Signal/Image
Recognition

 
    
 
 
   

Input
nsducer

    

Results ,

 

Measurement 1%;
Transducer ''''

  
 

 

Processing

  

Figure 1.1 A generic NDE system

1.2 Eddy Current Testing

Eddy current testing (ECT) is an electromagnetic technique and based on Michael

Faraday’s laws of electromagnetic induction discovered in 1831. A conventional eddy

current testing system is shown in Figure 1.2 [8]. It includes a coil or a set of coils

excited by a time varying current. The time varying magnetic field generated induces
eddy currents inside the conducting specimen (Faraday’s Law). The eddy current also
generates a magnetic field (Ampere’s Law) which opposes the field produced by the

source current (Lenz’s Law).

Source

 

(Q)

Conducting plate

 

 

 

Figure 1.2 System of eddy current testing

The existence of defects in the test specimen changes the eddy current distribution
which in turn alters the net magnetic flux linking the coil. Consequently the presence
of a defect is detected as a change of the probe coil impedance [8].

Eddy current techniques are widely used in the inspection of aircraft structures
composed of aluminum layers for detecting hidden cracks and corrosion. A major
limitation of eddy current inspection technique is it can detect only surface or near
surface defects because of the skin effect, which is related to the depth of penetration
of fields. It can be shown that the eddy current density is maximum on the surface and
decays with the depth inside the conducting plate. The amplitude of fields at a depth

x could be expressed as

x

Ax = Aoe 5
(1.1)

where A0 is the amplitude of field on the surface of an infinite conducting plate.

5 = /_1_
(for!

And 6 is the skin depth given by [9]

(1.2)
0' : the conductivity of the plate.

,u: the permeability of the plate.

f : the frequency of the excitation.

The eddy currents are affected by the test specimen and experimental operating
conditions. The test specimen parameters that vary with electrical conductivity,
magnetic permeability, thickness, shape etc. The experimental operation conditions
include the design of the excitation coil, operating frequency, lift-off and so on.
Several variations of the basic ECT have been developed in recent years. The newly
developed inspection methods based on ECT are Magneto Optic Imaging (M01) and
Giant Magneto-Resistance (GMR) sensor which essentially are sensors that can
measure the fields associated with induced currents. Details of these two techniques

are provided in the next chapter.

1.3 Forward and Inverse Problems

In NDE, forward and inverse problems are two essential areas of research. The
forward problem involves simulating the system performance given an excitation
source and predicting the measured signal on a flawed test specimen. In contrast, the
inverse problem involves estimating the defect parameters based on the information
in the measured signal [10], [11], [12].

A systems approach to illustrate both the forward and inverse problems is shown in
Figure 1.3 [13]. In Figure 1.3 (a), the input and the system transfer function are
known, and the output is to be determined, which is forward problem. In Figure 1.3
(b), the system and output are known, and the input is to be determined, which is the
problem of signal inversion. One approach is to assume a set of possible input
(defects) as a priori, and from the measured output, the true input signal (defect) is
estimated using a model based approach or other techniques such as maximum
posteriori probability estimation. In Figure 1.3 (c), the input and output are known,
but the system is to be determined, which is also an inverse problem also known as

system identification problem.

 

INPUT SYSTEM OUTPUT

 

 

 

 

(KNOWN) (KNOWN) (To BE
COMPUTED)
INPUT SYSTEM OUTPUT
(TO BE (KNOWN) (KNOWN)
COMPUTED)
INPUT SYSTEM OUTPUT
(KNOWN) (TO BE COMPUTED) (KNOWN)

 

Figure 1.3 The system approach to illustrate the forward and inverse problem

The mathematical term of well-posed problem starts from a definition given by

Hadamard [14] by the following three properties:

(a) Existence: there exists a globally-defined solution for all reasonable data.

(b) Uniqueness: the solution is unique.

(c) Continuity (stability): which means the solution depends continuously on the
given input data.

A problem which is not well-posed is defined to be ill-posed.

Generally, the well-posed forward problems will lead to the availability of a stable
accurate solution. While, the inverse problems are often ill-posed, which means the
measured signals lacks the continuous and reliable dependence on the defect profiles,
non-uniqueness is a major issue to be concerned, which means a measured signal can
correspond to more that one defect profiles. These challenges lead to the development

of various constrained techniques to find the optimal inverse solution from the set of

possible solutions. Both the details of the forward problem and inverse problem will

be described in the following chapters.

1.4 Need for Modeling

Research in both forward and inverse problem depends strongly on the
availability of a good forward model that can simulate the system performance.
Firstly, the model is useful in visualization the material-energy interaction which is
important for us to understand the underlying physical principles. Secondly, by
providing quantitative values of field distribution, a numerical model can be used to
adjust and improve the system parameters such as operating frequency, excitation
source, sensor, probe parameters and so on. Thirdly, when the measured signal is not
available or not adequate, the output of the forward model can be substituted for
training data in developing signal inversion schemes.

The main uses of computational models are summarized below:

1. Solve forward problem

2. Understand the underlying physics of material-energy interaction
3. Visualize fields, such as induced currents, reflected waves

4. Optimize the system design and operational parameters

5. Testbed for generating defect signatures

6. Model based probability of detection (POD) estimation

7. Model based inverse problem

All electromagnetic phenomena, including ECT, are best described using
Maxwell’s equations. The nature of NDE applications leads to problems with high
dimensionality, nonlinear material properties and complicated solution domain and
boundary. Often, these problems involve dynamic geometries and the solution may be
a function of time and position. These complications have led to development of
various types of models, which include analytical and numerical modeling.

Analytical modeling derived from Maxell’s equations played an important role
before the development of computers. Today, the computer aided analytical modeling
is also valuable for validation based on some assumptions for simplification. The
drawback is that analytical approaches can only handle simple domain geometries and
hence have limited applications.

In contrast, numerical modeling can simulate problems with complex geometries,
high dimensionality, nonlinear, anisotropic and inhomogeneous material properties.
Numerical models, by virtue of the massive time and memory resources of computers
are far more powerful than the analytical models.

Recently, the development and implementation of the numerical modeling have
become a major research focus in NDE [15]. Several numerical modeling methods
have been developed for solving different governing equations, such as Finite
Difference Method (F DM), Finite Element Method (FEM), Boundary Element
Methods (BEM), Boundary Integral Methods (BIM) and Volume Integral Methods
(VIM). Each method has its own advantages and limitations and is appropriate for

different kinds of numerical problems [7], [16], [17], [18].

FDM is the simplest numerical modeling method to solve partial differential
equations (PDE). The advantage is its relatively ease and simple form for replacing
the partial derivatives by appropriate difference formulation. It is widely applied for
direct current, quasi-static, transient fields, and linear problems. The limitation is that
FDM has poor convergence property for irregular geometries due to its regular
discretization. Also, it is rather difficult to model distributed parameters such as the
current densities, conductivities and permeability [1 8].

FEM evolved in the late 1950s as a numerical technique in structural analysis and
quickly developed as a major numerical modeling method in various engineering
fields to solve PDE [19], [20], [21]. Because of its accuracy and efficiency in
modeling, it was applied to electromagnetic NDE in late 19703, especially for direct
current, low frequency fields, and permanent magnet problems [7]. Compared with
FDM, FEM has many advantages including its ease to impose the essential boundary
conditions and model complex geometries. FEM can handle higher order
approximation and lead to faster convergence and better accuracy. FEM has already
been developed for various 2D and 3D eddy current NDE problems [8]. The
drawbacks of FEM include large computer resources needed, especially for nonlinear
and time-dependent problem. Also, it is not well suited for open region problems [16].

Differential equations are commonly used for describing low frequency
electromagnetic field problems, such as electrical machines and eddy current testing.
However, in the area of antennas and electromagnetic wave propagation, integral
equations [22], [25] are more commonly used. BEM and VIM methods are based on

the integral equations [16], [l 7].

BEM solves Maxwell’s equations by using the given boundary conditions to
discretize the surface in integral equations, rather than calculating all the values in
solution domain. In the post-processing, the integral equation is used again to
calculate the required physical quantity at any point in the solution domain. Hence
BEM is more efficient than other models with regard to computer resources and
ensures good accuracy of solutions. It can be used to solve the fields in linear
homogeneous problems when the Green’s function is available. The drawback is that
its global stiffness matrix is a full matrix in contrast to sparse banded matrix
encountered in FEM, and hence requires more computational time [1 6].

In the case of inhomogeneous and nonlinear problems, VIM is generally introduced
to discretize the volumetric domain before solution. In VIM, the field is determined at
a point by summing the effects of the sources at all points, convolved with the
Green’s frmction. The advantage is it only necessary to construct a mesh over the test
sample and solve for the currents in the test sample. But this method also require the
Green’s function for the problem to be available and its global stiffness matrix is a
full matrix instead of banded matrix in FEM, which requires more computational time

for solver [17].

1.5 Introduction to Meshless Methods

The traditional FEM is a well established technique that has been successfully
applied to ECT modeling. The fundamental idea in FEM is to approximate a

continuous function over the entire solution domain by piecewise continuous

10

f

approximations, usually polynomials, over a set of sub-domains called finite elements
(FE). The interconnecting structure of the elements via nodes is called a finite element
mesh. The reliance of FEM on a mesh leads to some characteristic disadvantages.
Generation of a mesh for a complex 3D geometry is generally difficult and time
consuming [26]. For example, in electromagnetic problems that involve geometrical
discontinuities and propagating tight cracks, the use of a dense underlying mesh leads
to a high dimensional matrix equation and large computational error.

In recent years, an alternate numerical approach, known as meshless method [2 6]
has been developed that eliminates some of the disadvantages of FEM. In meshless
methods, the unknown function is approximated entirely in terms “local” functions
defined on a set of nodes. In this approach, elements and the usual relationship
between nodes and elements are not necessary to construct a discrete set of equations,
which obviates the need for generating a mesh in the solution domain. Consequently
meshless methods are especially useful in the problems with discontinuities or
propagating cracks.

The initial idea of meshless methods dates back to the smooth particle
hydrodynamics (SPH) method for modeling astrophysical phenomena [27] (Gingold
and Monaghan, 1977). The research into meshless methods has become very active
after the publication of the Diffuse Element Method [28] by Nayroles, Touzot &
Villon (1992). Several so-called meshless methods, including: Element Free Galerkin
(EFG) by Belytschko, Lu, and Gu (1994); Reproducing Kernel Particle Method
(EKPM) [29] by Liu, Chen, Uras, and Chang (1996); the Partition of Unity Finite

Element Method (PUFEM) [30] by Babuska and Melenk (1997); hp-cloud method

11

[31] by Duarte and Oden 1996; Natural Element Method (NEM) [32] by Sukumar,
Moran, and Belytschko (1998); Meshless Galerkin methods using Radial Basis
Functions (RBF) [33] by Wendland (1999) have also been reported in the literature.
The major differences in meshless methods come from the techniques used for
interpolation.

The principal attractive features of meshless methods for NDT are the possibility

of:

(i) Working with a cloud of points that describes the underlying structure
exactly instead of relying on a tessellation of the domain in forward
problems

(ii) Lack of need for remeshing in solving inverse problems or probe scanning

Element-Free Galerkin Method (EFG) is one kind of meshless methods, which

employs moving least square approximation with Galerkin formulation. The EFG
method is now widely applied to problems in fracture mechanics and

electromagnetics [26], [33].

1.6 Research Objectives

The major objective of this dissertation is the development of an efficient
numerical model based on the EFG method for solving the forward problem and a
scheme for model-based solution of the inverse problem. Although the EFG method
is shown to be convenient and robust in electromagnetic field calculations, the

theoretical basis is not completed developed. Some improvements to the basic EFG

12

method are proposed and discussed for enhancing accuracy of solutions to forward
problems. Our preliminary 1D, 2D and 3D results of the improved EFG formulations
are shown to be promising. These improvements include:
1. Selection of weight function and the regularity of the influence radius of the
nodes [41]
2. Selection of orthogonal basis functions [42]
3. Effect of different formulations for applying boundary conditions and the
formulation of the penalty function [43]
4. Treatments of discontinuities in the interface with visibility criterion [34]
Using the forward EFG model, new approaches for solving inverse problems are
also developed in this dissertation. The major advantage of using EFG model in
inverse problems is that re-meshing at each iteration is eliminated resulting in reduced
computational time and increased accuracy of solutions.
The inversion schemes implemented in this work are:
1. 1D parametric inversion of tight cracks
2. 2D nonparametric inversion based on state space search
3. 2D nonparametric inversion based on gradient search
4. 3D nonparametric inversion based on modified state space search
The dissertation consists of the following chapters. Chapter 2 briefly introduces the
necessary electromagnetic theoretical background. Chapter 3 describes the
formulation and improvements of the forward EFG method in details. Chapter 4
shows results of 1D and 2D model validations of the EFG method in electromagnetic

modeling problems. Chapter 5 presents the 3D implementations of the EFG method in

13

simulating M01 and GMR inspections. Chapter 6 includes the introduction of inverse
problems and the preliminary numerical results of the application of EFG method for
1D and 2D defect reconstruction. Chapter 7 presents the modifications for the 3D
inversion using state space search technique and preliminary results are included for

validating the approach. Chapter 8 covers the conclusions and future work.

CHAPTER 2. BASIC ELECTROMAGNETICS FOR EDDY

CURRENT MODELING

2.1 Introduction

A brief introduction to the basic electromagnetic background for the eddy current
NDE is presented in this chapter. Eddy current NDE phenomenon is governed by

Maxwell’s equations.

2.2 Maxwell’s Equation

All electromagnetic phenomena, including ECT, are best described using

Maxwell’s equations. The differential form of these equations is expressed as follows:

\7><E=-§E
at
(2.1)
VxH=J+Q
at
(2.2)
V-D=p
(2.3)
V-B=0

15

(2.4)

along with the continuity equation derived from equation (2.2):

V-J+a—p=0
at

(2.5)

where the variables are:
E: electric field intensity (volt/meter)

H: magnetic field intensity (ampere/meter)
D: electric flux density (coulumn/meter3)
B: magnetic flux density (tesla)

J: electric current density (ampere/meter3)

p: electric charge density (coulomb/meter3)

2.3 Time-Harmonic Fields and Potential Functions

When the excitation energy (source) is varying sinusoidally with time and we are
interested in the steady state response, the solution can be considered in the frequency
domain and the Maxwell’s equations (2.1) and (2.2) can be rewritten in terms of the
excitation frequency:

V x E = — j (OB
(2.6)

VxH=J+ij

l6

(2.7)
The governing partial differential equation is usually solved in terms of potential
functions, which could be the magnetic vector potential and/or the electric scalar

potential introduced below:

2. 3. 1 Magnetic Vector Potential

We define magnetic vector potential A as the following identity:
V - V x A E 0
(2.8)

From equation (2.4) and (2.8), B can be expressed as the curl of magnetic vector

potential A.
B=VXA
(2.9)
2. 3.2 Electric Scalar Potential
The electric scalar potential V is defined by the identity:
VxVV=O
(2.10)
And from equation (2.1) and (2.9), we get:
V x (E + 91) = 0
at
(2.11)

When combining equation (2.10), we obtain:

17

Ez—i’i—VV
at

(2.12)
By introducing the definition of the magnetic vector potential and electric scalar
potential, the governing equations, can be fully rewritten using an A-V formulation

and solved using numerical modeling methods.

2.4 Weak Form and Galerkin Method

The electromagnetic problems discussed in this dissertation are boundary-value
problems described by a governing differential equation in a solution domain Q:
Lu = f in Q
(2.13)

together with the homogeneous boundary conditions on the boundary S that encloses

the domain:
2:- = O on 81
(2.14)
it = 0 on 82
(2.15)

where L is a differential operator, f is the excitation source or force function, u is

the unknown quantity, and S =S1 USZ is the boundary of domain Q. In

18

electromagnetics, the governing differential equations can range from simple Poission
equations to diffusion and scalar/vector wave equations.

Ideally, one desired to obtain an exact solution for the boundary-value problems,
while in most practical engineering applications, only approximate solution is

tractable. Generally, the approximate solution is defined by the expansion
[1 n
u (X) = ZVJ'CJ'
1' =1

(2.16)

where vj is the shape function defined over the entire domain, cj is the

coefficients to be determined, and n is the total number of discretization nodes. In
I: n

FEM and EFG, u is written as u = Zoiu, , where (I), is shape function and u,-
i=1

h

is the nodal value. Substituting u for u in equation (2.13) will result in a nonzero

residual of the form
r = Luh — f at 0 in Q

h directly, the weak form seeks a solution that satisfies the

Instead of solving for u
equation in an average sense, which is called a weak form. Therefore, the optimal

h

approximation for u is the one that reduces the residual r to the least value at all

points in domain Q. So a weighted residual form is considered as
Rl '-'-" JQWirdQ = O

(2.17)

where R1 defines the weighted residual integrals and Wt is the weighting function.

In Galerkin method, the weighting functions are chosen to be the same as the shape
functions used in the approximation expansion. This usually leads to the most

accurate solution. Particularly, in FEM and EFG, the weighting function is w,- = (D,- .

So equation (2. 17) can be written as

n
R,- =[Q(q>,- 22(1),... —cr>,- -f)dn=o

i=1
(2.18)

Equation (2.18) can be expressed by the linear algebraic equation
[Gllul = [F]

(2.19)

with Gij=jn(<r>,.-Loj)dn and F,=[Q<I>,--fd§2. [G] is called the stiffness

matrix, [F] is the load (or source) vector.

For each local integral, the EFG method involves higher order terms than

conventional FEM, which thereby yields more accurate solution.

20

CHAPTER 3. ELEMENT FREE GALERKIN METHOD

3.1 Introduction

In Element free Galerkin method, a set of nodes is used to construct the discrete
system of equations approximating the solution. However, in order to implement the
Galerkin procedure, it is necessary to compute the integrals over the solution domain;
this is done by defining the support of the basis functions using either a set of
quadrature points or a background mesh.

Consider u(x) is a continuous function to be approximated. The theoretical

foundation of the approximation in the EFG method is the moving least squares
(MLS) approximation, which relies on three components: (i) a weight function, (ii) a
polynomial basis, and (iii) a set of position-dependent coefficients. The weight
function is nonzero only over a small subdomain of a node, which is defined as the
domain of influence. The domain of influence plays an important role in the
performance of EFG method in terms of accuracy of the solution.

In the following sections, we will present the basic formulation of EFG method
using MLS approximation. The modifications developed for efficient implementation
of boundary conditions and interface conditions, selection of weight function and use

of higher order orthogonal basis firnctions are discussed in detail.

21

3.2 Moving Least Square Approximation

The MLS is a method to obtain a differentiable approximate function in the domain

using known function values at discrete points through a weighted least squares fit.

In MLS approximation, the interpolant uh (x) is given by [34]

u” (x) = ij(x)a)(x) = pT(x)a(x)
j=0

(3.1)

where m +1 is the number of terms in the basis function, p j (x) are monomial

basis functions, and a ,- (x) are coefficients that depend on position x.

When x and .‘x' are not superposed, Lancaster and Salkauskas (1981) defined

the local approximation by:
h m r
u (xx) = 2p 1(E)a ,-(x) =1» (Edam
.=0

(3.2)

where x is the approximation point, i is a particular node. For example, in two
dimensions, uh (x) can be expressed in terms of either a linear or a quadratic basis
as:
uh (x, y) = a0 (x, y) + all (x, y)x + 02 (x, y) y (linear basis)
(3.3)

h
u (x. y) = ao(x, y) + 01(x,y)x + 02(x. y)y +

2 2 (quadratic basis)
03(x, y)x + a4(x. y)x)2 + 05(x. y)y

22

(3 .4)

The coefficients aJ- (x) are determined by minimizing the weighted error between
local approximation and the nodal values u -, i.e., by minimizing the following

quadratic form:
" h 2
J = zMx-xjxu (x,xJ-)—uj)
J'=1

= ZMx-xj>12p.(xj>a.(x)—u .12.
j=l 1-0

(3.5)

Here w(x — x j) is a weight function with compact support, n is the number of

nodes in the neighborhood of x where the weight function does not vanish. In

matrix notation, Equation (3.5) can be rewritten as

J = (Pa — u)T W(x)(Pa — u)

(3.6)
where uT = (u1,u2,- --,u,,) are the unknowns
P = [pi(xj)]mxn ,
(3.7)
W(x) = diag[W(x " x1),W(x _ x2),“',W(x _ xn)]°
(3.8)
The minimization of J with respect to a(x) leads to
B(x)a(x) — C(x)u = 0 ,
(3.9)

23

where

B(x) = PTW(x)P,

(3.10)
C(x) -.—. PTW(x).
(3.1 1)
It follows that
a(x) = B"1(x)C(x)u.
(3.12)

Substituting (3.12) into (3.2), and letting f = x , the MLS approximation can be

written as

u”(x) = Zoj(x)uj

j=1
(3.13)
where the shape functions (Dj are given by
m
<I>,-(x) = gnaw—Race».- = [FE-‘0,-
(3.14)

3.3 Boundary Conditions

In a numerical solution, it is necessary to choose a solution domain with a

surrounding boundary to simulate a particular problem. Therefore, it is essential for

24

the numerical model to consider the physical processes in the boundary region. The
boundary conditions are very important for the existence and uniqueness of the
numerical solutions. Different boundary conditions may result in very different
simulation results. Improper sets of boundary conditions may introduce erroneous
influences on the simulation system, hence imposing the proper set of boundary
conditions is very important.

It should be noted that the shape function derived in equation (3.14) does not

satisfy the Kronecker delta criterion at nodal locations: i.e. (I) j(xk) at (SJ-k; therefore,

uh (x j) at uj , which means when the MLS approximation is evaluated at a particular

point in the space, the value obtained is not the actual value of the field being
approximated, although they are related. This makes it difficult to impose essential
boundary conditions in EFG. In order to address this, special techniques have been
applied such as Lagrange multipliers [34], modified variational principles [44], or
coupling with standard finite elements [45] at the boundary. Note that the
construction of the shape function is theoretically identical in both two and three
dimensions. Although all these techniques are workable, each technique has its own
advantages and limitations.

Generally, Dirichlet type boundary conditions can be imposed through the use of
Lagrange multipliers. Furthermore, it is the most accurate method for imposing
Dirichlet boundary conditions and useful in relatively small problems, such as 1D or
smaller 2D geometries, where the cost to solve the problem is not big. But it adds to

the complexity in that the resulting stiffness matrix will have an increased size and

25

will no longer be positive definite or banded, which may lead to slower convergence.
The proof of this statement is shown in the following.
With Lagrange multipliers method for boundary condition, the original linear

algebraic equation (2.19) is changed to
G K a _ F
K T O ,1 _ q

with Gij=IQ(CDi-L(Dj)d(2, Kij=—[r<p,.-der, u is the unknown, ,1 is the

(3.15)

Lagrange multiplier, F,- =1 ody- - fdQ, qk z—irNk -uodl" , and Nk is the shape

function associated with the Lagrange multipliers. While the original stiffness matrix
G is sparse, positive definite and banded, the new sub-matrix K is a full matrix
associated with the boundary elements. So the overall stiffness matrix is no longer
positive definite or banded, resulting in an increased number of unknowns and an
awkward structure for the iterative solver.

Modified variational principle method can result in a banded equation; it can also
reduce the number of unknowns and therefore reduce the computational costs, but it
sometimes yields unstable solutions as shown in [47].

In the coupled FEM-EFG [48] approach, the elements are placed around the
boundary of the domain as shown in Figure 3.1 [48]. The essential boundary
conditions are applied to FEM nodes. In the implementation of this method we use
finite elements as background mesh to perform the quadrature of the weak form. In

this approach the FEM nodes adjacent to the essential boundaries are already

26

 

available and easy to apply boundary condition. The limitation of the coupled
FEM-EFG approach is that the modified shape function has a discontinuity in the

derivative at the interface, which affects the accuracy.

EF G region

    

Boundary region

 

L I I T Essential boundary

Figure 3.1 EFG coupled with finite elements along essential boundaries

 

An alternate method based on the constrained variational principle is proposed in
this dissertation, which is the penalty function method originally by Suzhen Liu and
Qingxin Yang [43]. Consider the typical linear boundary value problem where the
linear algebraic equation using EFG method is obtained in the equation (2. 19)

[G llu] = [F l
(3.16)
Each element in the final discrete equations with penalty boundary condition is

obtained as:
Gil- =[Q(<1>,. -L(Dj)d§2+a[rd>i(bjd1‘
(3.17)
F,- =[thi-fd9+a[ruo<1>idl‘.

(3.18)

27

where F is the associated boundary region, a is the penalty parameter, and uo is

the boundary value. With the penalty function method, the stiffiress matrix is positive
definite and banded which translates into a reduced matrix size of the problem.

The selection of penalty parameter has a great effect on the accuracy of solutions.
If a is too small, the boundary condition term is negligible and the boundary
condition is not imposed effectively; if a is too large, diagonal elements
corresponding to boundary nodes and some off diagonal elements associated with the
boundary nodes become very large making the stiffness matrix ill-conditioned and the

convergence slower. Although the penalty methods are very simple to implement, a

much be chosen appropriately. a is determined to lie typically in the range fi'om 105

to 1020. Detailed research was conducted on some test problems in the following

chapter.
3.4 Weight Function

The weight function plays an important role in the performance of EFG method in
terms of the accuracy of solutions, complexity of computation (coding) and rate of
convergence.

The selection of the weight function within the domain of influence is not
restrictive, but it must conform to the following conditions:

1. The weight function should be positive valued;

2. The weight function should be relatively large for the nodes x j close to J? ,

28

and small for more distant nodes. In other words, it should be decrease in

magnitude with the distance from f to x j.

3. The weight function is continuous together with its derivatives up to the
desired degree.

In the numerical modeling, the solution domain is fully covered by the domains of

influence of all nodes. The shape of the domain of influence is arbitrary; a circular or

rectangular domain is typically used according to the domain geometry. Figure 3.2

shows two examples of the domain of influences.

 

(a) (b)

Figure 3.2 Domain of influence (a) circular; (b) square

In a 2D domain of influence, the weight Motion is derived and expressed at any
given point as [40]
w(x-xi) = w. -wy = wen-way).
(3.19)
where w(ri) , for i = x, y , are typically either Gaussians, exponentials or cubic

splines. These functions take the following form, respectively

29

2
_( -1)
w(r) -{1 e r

0

l

Z—4r2 + 4r3
3

w(r)=(§-—4r+4r2——§r
o

 

4 3

for r S 1 (Gaussian),
for r > 1
for r .<_ 1 (Exponential),
for r > 1
for r S -1-
2

for; < r S 1 (Cubic splines),

forr>l

(3.20)

(3.21)

(3.22)

where ,8 = 0.5 yields best convergence for the exponential weight, and r;- , for

i = x, y , are calculated as

rJr =|x—xj|/dmx , ry =ly—yjl/dmy ,

(3.23)

(3.24)

Figure 3.3 shows one example of weight function over a 1D domain of influence.

30

V.: V

l 2 X 3 4

 

Figure 3.3 The weight functions over 1D domain of influence

In the above equations, rimx is a scaling factor, (cx,cy) is the difference
between node x j and its nearest neighbor. It is important that dmax is carefirlly

chosen. If dmax is too small, then too few nodes are in the local domain of influence,
and the matrix B becomes singular; if dmax is too large, then too many nodes are
present in the local domain of influence, and the local character of approximation will

be affected, the accuracy will decrease and the calculation time will increase. In

general, the optimum choice of aimax is 1.4 to 3.5.

The continuity of the shape function is tightly related to the weight function. When

the weight function w(ri) and its (k — l) derivatives are continuous, the shape

function (DJ-(x) and its (k — 1) derivatives are also continuous. In 3D, the

definition of weight function is a natural extension of that presented for 2D solution
domain.

In the previous work, uniformly distributed background nodes were used. Further
the symmetric weight function makes the domain of influence symmetric as well.
When modeling large solution domain, non—uniformly distributed nodes are used for

efficiency, so the domain of influence is asymmetric. Hence the asymmetric weight

31

function should be modified and implemented accordingly so that it is continuous.

For example, when the cubic splines weight function is applied, it is itself continuous
. t d . . . . .

as well as Its 1S and 2H order denvatrves, even when 1t 18 asymmetric. In the case of

Gaussian and Exponential weight function, they are both higher order weight

functions and hence continuous as well. Figure 3.4 shows the asymmetry of the

domain of influence.

 

 

 

 

 

Figure 3.4 The asymmetry of the domain of influence

3.5 Discontinuities Approximation

In NDE applications the test geometry typically involves multiply connected
regions with each region characterized by its own material properties, such as
conductivity and permeability. This results in a discontinuity of the normal
component of current density at the interface boundaries which in turn implies that
the derivatives of the shape flmction or the shape function itself will be discontinuous

at the interface. Since continuity of shape functions is inherited from continuity of the

32

weight function, it is necessary to introduce discontinuity into the weight function.
One technique used to realize this is by using the visibility criterion [40].

The visibility criterion is illustrated in Figure 3.5, for the rectangular domain. If
there is no material discontinuity, the domain of influence is the total area of the
square. However, in the presence of a discontinuity, the domain of influence of node
x j shrinks to the area covered by dashed horizontal line, and the weight function
vanishes outside of that area. This procedure directly results in the discontinuity of

weight function, which in turn introduces the discontinuity of shape function and its

derivatives.

 

Figure 3.5 Domain of influence of node adjacent to material discontinuity
(rectangular domain)

3.6 The Orthogonal Basis

From equation (3.14), we can see that the basis function is one factor in the overall

shape function and hence affects the accuracy of solution. In theory, increasing the

33

order of basis functions improves the accuracy of solution. But in the EFG method it
is often seen that higher order (m 2 4) polynomial basis fimctions suffers from a lack
of convergence procedure because the associated shape function matrix becomes
ill-conditioned. The 1D proof of this statement is presented below.

In the EFG method, the shape function defined in equation (3.14) involves the
inversion of matrix B , which is defined in (3.10). Here W is the matrix of weight

function and P is the basis function coefficient vector.

 

B(x) = PTW(x)P

(3.25)

So B can be equivalently expressed as:

n T
B(x) = Zw(x —x,-)[1,x,-, ...... ,xlm] [l,x,-, ...... ,xfn].
i=1

(3.26)

Each element in matrix B is written in the form of an inner product:

" I
(p, (x), pk (x)) = Zw(x — xi)x,- x? , where l,k = 0,1, ...... ,m
i =1

(3.27)

So

<P1(x).P1(x)> (P1(x),Pm(x)>
B(x) = . . . . . . . .

(pm(x),pi(x)) (pm(x).pm(x))

34

(3.28)

Rewrite the inner product in integral form, we get
n I k
(p1(x),pt(x)) = In 2 W(x - xi)x x,- dx
i =1

(3.29)

When L2 (a,b) space is unit of L2(—1,1) space, which means the solution
domain is [-l, 1], then matrix B is Hilbert matrix, which is ill-conditioned and

difficult to invert numerically.

i 1 1/2 l/(m+l)-

1/2 1/3 1/(m+2)
B(x)=H=

 

 

_1/(m+1) 1/(m+2) 1/(2m+1)]
(3.30)

For the general L2 (a,b) space, the structure of matrix B is similar to Hilbert

matrix, which is ill-conditioned and hard to invert. This makes calculation of shape
function difficult. This is the reason that the higher order (m 2 4) polynomial basis
fImctions lead to lack of convergence in EFG procedure.

Choosing an orthogonal basis [42] can solve this problem. The definition of

orthogonal basis is shown below:
Vector p(x)=[p1(x),..., p j(x),..., pm(x)]T is an orthogonal basis, if the inner

product satisfies

35

N o 1 k
(Pl(x)apk(x))=ZMx—xi)p,(xi)pk(xi)={ ( it )

 

 

 

 

 

H C, at O (l at k)
(3.31)
Therefore, its shape matrix becomes
“(pilx),pi(x>) ‘
B(x) =
_ (R... (x). p... (x))_
(3.32)
The recursive relationship of orthogonal basis is given by,
P107) =1
P2 (x) = (x _ a2)Pl(x)
Pk+l (x) = (x _ ak+l)Pk (x) - flkPk—l (x)
(3.33)
where k = 1,2,...,m. The parameters a , ,6 are Obtained as,
. _ (wk(x),pt(x))
ak+l _
i <pk(x)spk(x)>
flk = (pk(x):pk(x))
[ (Pk—1(x),pk—1(x)>
(3.34)

where k =0,1,...,m—1. Substituting equations (3.32)-(3.34) into (3.13), the shape

functions are then expressed as:

36

"’-W(x xit)p,-(x)

 

 

 

 

 

 

 

(Di( )= ( )
x ,2—1 <Pj(x >p,-(x>) p x
(3.35)
The partial differential form of the shape function is
d<D,-(x) _ imdfflpflxiWfl-x)
2. j-) (pjtx) p,<x))
iwix- x.)p,-(x.-) dp’”
+
)h (mum-(x)?
0’ (x). -(x)
iwtx—xtipjtxtipjoi 00’ hp] >
j=1 (ir2,-(x),p,-(x))2
(3.36)
And the partial differential of inner product form is
d<Pj(x):Pj(x)>_ N Zdw(x— xi)
.1. ZMx— —x.->[p (xt>]’= .1. —Ln (xiii
i=1
(3.37)

These equations (3.3l)-(3.37) indicate that the computations in EFG method

eliminate the need for inversion of the shape function matrix when the orthogonal

bases are used. This in turn makes the implementation of EFG method simpler [48],

increases the accuracy of the solution by allowing higher order bases and also

converges in fewer numbers of iterations.

37

CHAPTER 4. FORWARD MODEL VALIDATION

4.1 Introduction

In Chapter 3, the basic formulation of EFG method and modifications to improve
its performance have been proposed. The improved formulation will be validated by
static and quasi-static electromagnetic field computation via comparison against
analytical models and the conventional EFG. Results of performance of the forward
model on both 1D and 2D problems are presented in this chapter.

The implementation of static and low fiequency electromagnetic fields is the basis
of eddy current inspection of conducting test specimen. To simulate this underlying
physics of material-energy interaction, the fields are typically expressed in terms of
potentials and the governing equations which are Poisson or diffusion equations are
solved together with appropriate boundary and interface conditions. In the following,
we will analyze the application of the improved EFG formulation to the numerical
solution of these equations.

In general, consider a domain of interest denoted by (2 that is bounded

by S] USZ = 80 , where S] is the Dirichlet boundary and S2 is the Neumann

boundary.

38

4.2 1D Problem

A 1D electrostatic Poisson problem for the potential function u is described as

follows:

V2u=—x (Q:0<x<l)
u(0)=0 u'(l)=0

(4.1)

To obtain the discrete system equations, it is necessary to use a weak form of the
equilibrium equation and boundary conditions. Galerkin method is used for the choice
of test function and Lagrange multipliers are applied for incorporating boundary
conditions. The final discrete function can be obtained by substituting the trial

functions u(x) and test function into the weak form, yielding the system of linear

algebraic equations.

3

The exact solution is u(x) = éx—%x for comparison with the numerical

solutions.

Figure 4.1 shows the exact solution and the numerical solution using EFG method

with 1St and 2nd order polynomial basis. Although the solution with 2nd order

polynomial basis convergences, it starts to show some instability as explained in
Chapter 3, that the higher order polynomial basis produces numerical oscillations of

the solution.

39

0.35

 

 

-ll— exact
0-3 ' —1— 1st order
—9— 2nd order

  

 

 

 

0.25 -

0.2 -

0.15-

0.1 -

0.05 -

 

 

   

 

l I l J

0.2 0.4 0.6 0.8 1
X

Figure 4.1 The comparison of exact and EFG solution with 1“, 2nd order
polynomial basis

Figire 4.2 shows the exact solution as well as the solution using EFG method with

d . . . .
3r order polynomial basrs. The numerical solutlon no longer converges. It means the

oscillations become more severe so that the numerical solution is invalid.

4O

 

 

0.35

 

—II- exact
0'3 ' —-l-— 3rd order

0.25

 

 

 

 

0.2 -

z 0.15 -
0.1 -
0.05 -

 

‘-
qr-

 

 

 

'0.05 r r r 1

Figure 4.2 The comparison of exact and EFG solution with 3rd order polynomial
basis

In Figure 4.3, the numerical solution using EFG method with orthogonal basis
function is presented. The order of the basis function varies from 1 to 4 and the
numerical solutions are seen to be stable and convergent, which demonstrates the

effectiveness of higher order orthogonal basis functions.

41

 

 

0.35 . . t .
—*— exact arr/AA

0'3 ' —+- 1 st order . #1:“ q

_ —9-2nd order
0'25 HER-— 3rd order
4111 order

 
 
   

 

 

 

 

 

 

Figure 4.3 The comparison of exact solution with various order orthogonal basis

4.3 2D Problem

The 2D eddy current modeling geometry is shown in Figure 4.4. Assume an

infinite-size 2D conducting plate with a tight (zero width) crack along y direction

(length = 1.0 cm) placed as shown. A time-varying harmonic current source is applied

using a thin foil in the x direction.

42

 

 

 

 

 

-i .0:5 0 0.5 1

Figure 4.4 Modeling geometry for 2D problem

The induced current will be strictly along - x direction in the absence of crack.
The induced currents will bend around the comers in the presenceof the crack. This
geometry is typically encountered in 2D eddy current problems. The plot of the

induced eddy current predicted by the model is shown in Figure 4.5.

Figure 4.5 (a) shows the induced eddy current with the crack obtained using 1St

order polynomial basis function and (b) shows the image obtained using 2nd order

polynomial basis function. Due to the complexity of the problem, the EFG model

usrng 2n order polynorrnal basrs function rs no longer convergent and solutlon lS

incorrect.

43

 

 

 

 

 

 

 

 

 

 

 

Q—Q—P
v-a-r-
o_5._._..,,,/\\ ._._.
> o.
$.5P‘hhk\\///’~_*—4
l—‘—‘~Mk.—Q—fW'—.——-o—
.1"
bhhkkhvaflkh
-1 -0.5 0 0.5 1
X
(a)
1..
0.5“-
> 0/
-0.5*
.1" _ ‘ '1
-‘I -05 0 05 1
X

Figure 4.5 Induced eddy current plot using EFG method by polynomial basis
function (a) 1”t order (b) 2mi order

Figure 4.6 (a) shows the induced current plot using EFG method with 2nd order
orthogonal basis and (b) shows the difference compared with Figure 4.5 (a). It is clear

that the EFG model with 2nd order orthogonal basis function still converges and only

. . . . t
rrunor differences are observed at the boundaries compared wrth the 18 order

 

 

 

 

 

 

 

 

 

 

polynomial basis function.

l—T-Iv—o-P-r-l‘w—‘a-v-fi‘fi—hfi-Ihh;
1;,E,,E_M“hahi
._.._._.‘,.‘_,/‘-u-..xh ._
0.5~—e—-e.—»//\\--e—e—J
...._...,., I xxx”..._
a or‘TTT ' ' ‘ -.._._.__
).—.—-—-—~ - .2...—._._
L_._.~\\ \ 1”“...—
.0.5~—-—--\\//.er~—~—-
__._...kx\/fr._.._._
1p-hhhhhhrrf-po_
TLTTT‘T‘TTTrTrTe-TTTTT‘IT
.1 -0.5 0 0.5 1

x

(a)
1, —« 1 ~33 ..v- w.
.. , , ‘ - r \ \ I .‘
' ’ " "‘ "‘ " "\
0_5_- . . . . . , , - , _ x-
’ ~ \ /
> 0 II‘C‘
_o_......,......5\
‘1’, . . . 1 , t , . . . .4
-1 -05 0 05 1

x

(b)

Figure 4.6 Induced eddy current plot using EFG method by orthogonal basis
function (a) 2"‘1 order (b) comparison with Figure 4.5 (a)

45

Figure 4.7 shows the iteration error in each step using lSt order polynomial basis

and 2nd order orthogonal basis. From the following figure, we can conclude that the

EF G model developed in this dissertation with orthogonal basis converges faster than

the one with conventional polynomial basis.

 

 

 

 

 

10“ ............. . .............................. .,
I?33553533321535323255335355553255553225315+ Poly-Basis E:
::-..‘::::::2::::::::::::::::::::::::::::3: -°- Who Basis r‘

10.1 n::::: £::: :s::::: :: 3%": as: :- .:::a: is: as: .::: .::: "shear: .::::2: :1

In

E 2

u, 10‘

C

2

g 10:3
10"-
10‘5

 

 

 

 

Figure 4.7 The convergence rate using polynomial basis and orthogonal basis

For the 1D problem in section 4.2, the Lagrange Multipliers is used to impose the
essential boundary condition, because of the small size of the solution domain. In the
above 2D eddy current problem, the penalty method is used to impose the essential
boundary condition. The size of the discrete matrix is reduced while the accuracy is

retained.

46

CHAPTER 5. FORWARD MODEL IMPLEMENTATION

5.1 Introduction

This chapter presents the application of the EF G method to electromagnetic field
computations in a 3D NDE test geometry. In quasi-static frequency range, the
governing equations can be simplified to a diffusion equation in terms of magnetic
vector potential and electric scalar potential. The numerical results are validated via

comparison against experimental data and results obtained using conventional FEM.
5.2 Mathematical Formulation

Let ()1 and 02 be partitions of the solution domain, where Q] denotes the
conducting region and 02 represents the surrounding free space. Let A be the
magnetic vector potential in both 91 and 02 , and V be the electric scalar
potential in (21. ,u and 0' are the permeability and conductivity of the media
respectively. Derived from Maxwell’s equations, the governing diffusion equation for

conventional eddy current problem in Cartesian coordinates can be written as [50]

inVxA+jcoaA+0VV=Js in Q
,u

(5.1)

V-(jcocrA+oVV)=0 in O,

47

I '1

l-

(5.2)

where J S is the assigned source current density.

To obtain the numerical solutions of equation (5.1) and (5.2), we first expand the
potentials (magnetic vector potential and electric scalar potential) as a linear
combination of shape functions associated with the nodes, we get:

A = 2(ijojd, +ij<I>jiiy + Azjcp 1d,)
,-

(5.3)
,- ,-

(5.4)

where ij , A Azj are the three components of the magnetic vector potential at

yj’

A

node j; Vj is the scalar potential; (DJ- is the shape firnction, and 8,, fiy, az

are the Cartesian unit vectors.

Homogeneous Dirichlet boundary conditions are applied at external boundaries and
current continuity conditioan = 01—ij —VV)«r‘r = 0 is applied on interfaces to
obtain a unique solution. Because EF G shape function does not satisfy the Kronecker
delta property, essential boundary condition can not be imposed directly as is done in
FEM. As discussed earlier, the penalty method is used for this purpose.

To obtain a system of linear algebraic equations, the Galerkin method is used.
Substituting the above expansions (5 .3) and (5.4) into the governing equations (5 .1)

and (5.2), we get:

48

3N
j=l

N 3N N
+Enact?)-V<Dde]VJ-+a2[[rk¢i-¢de]Aj+aZ[[rk(Di-<l>de]Vj
j=l j=1 j=1

3N 3N’ -N
:12] 0J5-¢de+ajZ=l[Irkuo-¢de]+aJZ=l[JrkuooCDJ-dr]

(5.5)

3N N
Zl[Inja)oV(Di .OdeMj + ZIUQOVQ‘. .chjdguyj z 0
J: }=

(5.6)

In (5.5), a is the penalty parameter. Generally, the penalty parameter should be

carefully chosen so that 0.92 equals 104 or 105 times greater than S at the

boundary nodes in the equation (5.5).

Equations (5.5) and (5.6) can be written in matrix form as

GA = Q
(5.7)

where the matrix G is a complex and sparse matrix; A is the vector of unknowns;
four unknowns are at each node, including the electric scalar potential and the three
components of the magnetic vector potential; Q is the load vector incorporating the
current source.

Equation (5.7) is solved using (TF) QMR [53]. Only non—zero entries of G need
to be stored in the iterative solution procedure. Compared with Lagrange multipliers,

the matrix obtained using penalty function is of a reduced size and positive definite,

49

which makes the convergence faster.
The solution of (5 .7) is then used to compute various measurable physical
quantities such as flux density B , and induced eddy current density in the area of

interest.

5.3 Magneto Optic Imaging Application

Magneto-optic imaging (MOI) is a relatively new sensor application of bubble
memory technology to NDT and produces easy-to-interpret, real-time analog images
[49]. In this technique, an induction foil is used induce eddy currents in the specimen
and a magneto-optic sensor is used to detect the magnetic flux density associated with
the induced currents. However, since the M01 data is binary, the information
provided is qualitative in nature. A theoretical model that can simulate the MOI
system can provide a more quantitative value of the magnetic fields and is hence
useful for the optimization of the magneto-optic (MO) sensor and system [51], [52].
A linear FE model utilizing the A — v formulation was initially developed for this
purpose. However the inclusion of infinitesimal air gaps between layers and the tight
crack makes mesh generation very cumbersome and decreases the accuracy of FEM

solutions. A more accurate EFG model is implemented and applied to MOI inspection

geometry.

50

5. 3.1 Introduction to M0! in Airframe Inspection Application

In aircraft industry, for the airframes to stay in service much longer than the
designed life, thousands of fasteners and bonded joints on the fuselage need to be
inspected. In order to cover the large inspection area, fast, accurate and cost effective
NDT inspection methods are clearly needed. Some of the challenges encountered in
ECT are:

(1) Detection of corrosions or cracks in the multi-layer structures

(2) Detection of cracks under the fastener (CUP)

(3) Detection of surface and subsurface defects close to edges

Conventional eddy current inspection method is time consuming due to the small
probe size and large inspection area. Furthermore, it requires well-trained operators
for data interpretation. Development of new techniques for rapid and accurate
inspection is of considerable interest to the aircraft industry. MOI provides analog
images that are fast and easy to interpret.

Briefly, in MO inspection, an induction foil is used to induce eddy currents in the
conducting sample. The magnetic fields associated with the induced eddy currents are
detected using a MO sensor. Polarized light passed through the MO sensor undergoes
a rotation of the plane of polarization in the presence of local magnetic fields. When
the reflected light is viewed through an analyzer, an analog image of local
magnetization is seen on a monitor. Figure 5 .1 (a) shows a schematic of the M01
device and Figure 5.1 (b) shows an experimental image. However the MO image is

binary because of the physical properties of the garnet film sensor and does not

51

provide a quantitative value of the measured signal. A simulation model for M0
inspection of multilayer aircraft geometries for predicting the continuous valued
magnetic flux density that produces the observed experimental MO images has been

developed.

Light source Q
I

Polarizer

 

Bias coil

Sensor 4 Induction foil

Test sam 1e

1

 

(a)

 

(b)
Figure 5.1 (a) The schematic of the MOI system (b) an experimental image

52

Figure 5.2 (a) shows a M01 308 system and (b) shows inspection being performed

using MOI system with ahead mounted display.

   
   
 
      

  

Personal
Video System

Low Frequency
Eddy-current
Attachment

 

 

(b)
Figure 5.2 (a) The M01 308 system (b) using of M01 system

53

5.3.2 Modeling Geometry

A typical 3D MO inspection geometry used in both FEM and the EFG methods is

shown in Figure 5.3.

 

 

 

 

 

 

 

 

 

 

 

 

(a) (b)

-J
h=1mm
I I d=3mm

| I

.....___..l L...—
, ,/ I I
all’ gap ’

 

 

 

 

(C)

Figure 5.3 M01 simulation geometry (both for FEM and EFG). (a) 3-D view, (b) 2-D
top view and (c) side view of the multilayer plate with a 6mm diameter fastener hole

The geometry consists of two layers of aluminum plates, of 3 mm total thickness
and of infinite width & length. An infinitesimally small air gap exists between the two

layers. An infinite inductive foil of 1 mm thickness, carrying a linear sinusoidal

excitation current density of amplitude 108 A/m 2 and frequency 3 kHz is placed

54

above the conducting multilayer plate at a distance of 1mm. The solution region is of
dimensions [—12,12]mm x[—12,12]mm x [—8, 8] mm . A fastener hole of 6 mm
diameter is introduced extending through the two layers and a second layer tight crack
of 5 mm radial length is also included. The number of nodes used for discretization
using the EFG method is 18 x18 x 12 whereas that used in FEM is 28 x 28 x 24 and

the cross-section is illustrated in Figure 5.4.

Figure 5.4 2D cross-section view of (a) FEM mesh and (b) EFG discretization of the
solution domain

0.015
0.01

0.005

-0.005

-0.01

 

.0015 -0.015 -0.01 -0.005 0 0.005 0.01 0.015

(a)

55

Figure 5.4 continued

 

0.015 T : = ? 3 : ? : 3 1T7 : 3 ? c 3 § : 3 4T

0.01

l
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O

I

0,005 i" O O O O O O O O O O O O O O O O O O Q —1

-0.005

I
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O

l

-0.01

I
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O
O

l

 

 

t)- o
0
I

-0.015

 

I o

_o

O
Co- 0
C
80- e
01

O
bo- e
..L

O
31)- o
01

15 20111

oo- o

-0. 05

(b)

The boundaries of the solution region are chosen in air and the imposed boundary

conditions are summarized in Table 1.

Table 1 Boundary conditions

 

Upper boundary Ax = Ay = A, = 0, V = 0
Lower boundary Ax = Ay = A2 = 0, V = 0
Left boundary Ax = 0, V ._. 0, Aya A2 are free
Right boundary Ax = 0, V = 0, Ay, A2 are free
Front boundary Ax = 0, V = 0, Ay, A2 are free

Rear boundary Ax = 0, V = 0, Ay, A2 are free

 

 

 

 

 

 

 

 

 

56

 

5. 3.3 Numerical Results I — Fastener hole without crack

Simulation results obtained using FEM and EFG for a defect free fastener hole are
shown in Figure 5.5. The normal component of the magnetic flux density 32 is plotted
at the MO sensor layer placed above the induction foil. A comparison of the
computational time and mesh parameters are summarized in Table 2. Figure 5.6

shows a comparison of model predictions after thresholding and the experimental MO

 

 

 

 

 

 

 

 

image.
Table 2 Summary of simulation results
FEM EFG
Peak value of magnetic 34.223 35.927
flux density (Gauss)
Computation time (hours) 3 0.5
Number of nodes 28x28x24 = 18816 18x18x12 = 3888
Discretization complexity complex (non-uniform simple (uniformly
elements) distributed nodes)

 

57

 

Figure 5.5 (a) 3D view, (b) 2D top view and (c) 2D side view of the normal
component of the magnetic flux density with FEM and (d), (e), (f) with EFG

 

 

 

 

 

 

 

 

(b) (C)

58

Figure 5.5 continued

NCO-h

 

0.01

 

 

 

 

 

 

(e) (f)

59

 

 

 

 

 

 

 

(C)

Figure 5.6 Binary results obtained by thresholding simulation results in Figure 5.5.
(a) FEM (b) EF G and (c) experimental image

5.3.4 Numerical Resulm II - Fastener hole with 5 mm second layer tight crack

Simulation results obtained using FEM and EFG for a fastener hole with 5 mm
radial crack are shown in Figure 5.7. The normal component of the magnetic flux
density is plotted at the MO sensor layer placed above the induction foil. A
comparison of the computational time and mesh parameters are summarized in Table

3. Figure 5.8 shows a comparison of the model predictions afier thresholding.

Table 3 Summary of simulation results

 

 

 

 

 

 

 

 

FEM EFG
Peak value of magnetic 35.467 37.779
flux density (Gauss)
Computation time (hours) 3 0.5
Number ofnodes 28x28x24=18816 18xl8x12=3888
Discretization complexity complex (non-uniform simple (uniformly
elements ) distributed nodes)

 

Figure 5.7 (a) 3D view, (b) 2D top view, and (c) 2D side view of the normal
component of the magnetic flux density with FEM and (d), (e), (f) with EFG

 

 

61

 

Figure 5.7 continued

 

 

 

 

 

 

(d)

62

Figure 5.7 continued

 

(6)

 

 

 

63

(f)

 

 

 

   

 

 

 

 

 

 

(a) (b)

Figure 5.8 Binary results obtained by thresholding simulation results in Figure 5.7.
(a) FEM (b) EFG

5.3.5 Conclusion

The EFG method, gives us an efficient fi'amework for modeling airframe geometry
in its full complexity. A benefit of the EFG method is that it is inherently a high order
scheme, which makes it more accurate than the conventional linear FEM. Another
advantage of the EFG lies in its simplicity of domain discretization because it relies
on only a cloud of nodes and does not require an underlying tessellation to describe
the 3D domain. The infinitesimal air gaps between layers and the tight cracks are
modeled by a layer of nodes alone. Initial results indicate that MO images predicted
by the EFG model are much closer to the experimental images than those of the FE

model with the same number of unknowns.

CHAPTER 6. INVERSE PROBLEM USING EFG

METHOD

6.1 Introduction

Inverse problem in NDE involves the estimation of true defect profile on the basis
of the information contained in a measured NDE probe signal. As mentioned in
chapter 1, inverse problems are generally ill-posed, and the full analytical solutions to
these problems are seldom possible. Therefore, practical solutions employing a
variety of constrained search techniques are used to find the optimal solution from the
set of possible solutions. These techniques range fiom simple calibration procedure to
pattern recognition algorithms such as neural networks. These approaches include

direct and iterative procedures.

In the following, several well-known methods to solve inverse problem are

reviewed:

Calibration methods [54] use calibration curves (or tables), which are obtained by
generating a set of signal parameters either from experiments or numerical models for
a variety of defect parameters to estimate (fit) a curve. The curve fitting technique can
use polynomial approximation or other functional approximation method. These
curves are then used to predict defect parameters for a given measured signal

parameter.

65

Pattern recognition algorithms using neural networks involve mapping the
measured signal directly to the defect profile [55]-[59]. In this case, signal inversion
is treated as a function approximation problem and the mapping from the measured
signal to the defect profile is learned by a neural network or other signal processing

algorithm using training data.

Model-based solution to inverse problems [60], [61] employs a forward model that
can predict the measured signal for a given defect profile in an iterative framework.
The iteration starts with an initial estimation of the defect profile represented by the
defect parameters, and forward problem is solved to determine the corresponding
signal. The cost function defining the error between the measured and predicted
signals is minimized iteratively by updating the defect profile. When the error is
below a pre-defined threshold, the iteration ends and the updated defect parameters

represent the desired solution.

Iterative methods generate accurate reconstruction solutions, however due to the
high-dimensional (2D or 3D) numerical models in each iteration; they are
computationally intensive and hence have limited practical applications. Previous
work using FE model based iterative techniques for inverting conventional eddy
current. data has been reported by Li [60].

In this dissertation, the model based NDT inverse problem that uses the EFG model
for simulating the physical process is proposed. The overall iterative approach using

the forward model is shown in Figure 6.1.

66

  
 
  
 

Initial Defect
Profile d1“) Fomard Model

by EFG Method

 

 

 

 

 

      

 

Update Defect Desired Defect
Profile dl'l Profile (I

Figure 6.1 Schematic representation of the iterative approach of the inverse problem

From the above iteration loop, we can see that this approach consists of two parts:
the development of a forward model of the underlying physical processes, which has
been presented in previous chapters; and a scheme for updating of defect parameters
to derive the desired profile. The different updating schemes are described in this
chapter.

Three approaches are considered in this thesis.
1. Parametric calibration method for 1D problem
2. State space search method that can be used in 2D problem
3. Gradient based method using Adjoint equation for 2D problem
The formulation, implementation and preliminary results of each approach are

presented here.
6.2 1D Problems - Parametric Approach

The case of a 1D tight crack is considered in this section. The NDE inspection uses

a linear induction foil to induce eddy currents and the normal component of

67

associated magnetic flux density is measured using a Giant magneto-resistance (GMR)
sensor. This approach is based on the generation of a calibration curve to characterize
a few crack parameters such as the length and orientation of the crack. The
performance of this procedure depends on (i) definition of defect parameters, (ii)

definition of the cost function and (iii) the defect parameter updating scheme.
6.2.1 Problem Statement and Inversion Scheme

The tight crack is characterized by the parameters r , z and 6 ; where r
corresponds to the radial distance of the crack center from the origin, 2 corresponds
to the length of the crack perpendicular to current direction and 6 corresponds to the
angle made by the crack with the linear current direction as shown in Figure 6.2. We
denote these independent variables by vector X :

X = {r, z, 0}

(6.1)

 

(0.0)

é—induced current

 

 

 

Figure 6.2 Definition of the defect parameters for 1D crack

68

Assuming a given initial value for X = {r,z, 6}, the EFG model solves the forward
problem for simulating the underlying physical phenomenon and generating the GMR
signal. The desired solution is obtained by iteratively updating the parameters to
minimize the error between the given experimental input signal and forward model
prediction.

The EFG model uses the visibility criterion to model the crack. The update
equation and the node discretization of the defect greatly save the computational cost.

The cost function, chosen as the square of the error between the experimental

measurement and model predicted signals is given in equation (6.2).
N 2
F = 2' En "' an l
=1

(6.2)

where B" is the n -th point of the EFG predicted signal representing the magnitude
of the normal component of magnetic flux, an is the n -th point of the measured

(or simulated) signal and N is the number of data points in B and Bm. The

derivative of the cost function is approximated as:

dF _ F(X +AX)-F(X)
dX - AX

 

 

(6.3)

69

6. 2.2 Inversion Results and Discussion

The results of inversion are presented here using simulated data instead of
measured signal. The 2D geometry consisting of an infinite aluminum plate is excited

by an infinite induction foil carrying a 3 kHz sinusoidal current along the x—direction.
A database of . tight cracks, X i = {H ,2i ,6?" ,i =1,2,3,...n} are simulated and the

predicted signals B,- are compared with the measured signal from a true defect Bmi

to get cost function values F,- corresponding to the i -th crack X i. Figures 6.3 (a),

(b) and (0) show the normalized cost function versus the r , z and 6 parameters
independently. The x -axis indicates the difference between predicted and true
parameter values. When the x -axis is zero, the predicted defect parameter is exactly
equal to the true defect parameter value. These curves are used for estimating the

derivatives in equation (6.3).

70

Figure 6.3 Cost function versus (a) r, (h) z and (c) (9

 

0.00 . . . 4 . -

0.05 -

0.04 -
E 0.03 »
0

0.02 ~

0.01

 

 

 

 

-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4
r-rtme (cm)

(a)
0.05 . . . f

 

0.04

0.03 ’

Cost

0.02 .

0.01

 

 

 

 

-0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4
‘r-rtrue (cm)

0))

71

Figure 6.3 continued
0.00 . . .

 

0.05 A
0.04

0.03

Cost

0.02 '
0.01

 

 

 

0 0102030405060700090
eflmkm)

(C)
In order to test the inversion scheme, a straight tight crack of infinite depth, zero

width and length = 1.0 cm is introduced and the corresponding signal is used as the

measured data 3",. The true and predicted cracks using this inversion scheme are

shown in Figure 6.4. The regular grid of dots shows the nodes in the EFG model.
Figure 6.4 (a) shows the initial guess (dotted line), the true defect (solid line), and the
frnal prediction (dashed line) when the straight tight crack is along the direction
perpendicular to the linear current. Figure 6.4 (b) shows the corresponding results

when the crack is at an angle of 55 degrees with the current direction.

72

 

 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

0.5

05

-1

 

 

1 ‘ ‘ ‘ 1 1 1‘ 1 1 ‘ 1 1 1i i‘

(a)

 

 

1 1 1 1 1 1 1 1 1‘ 1 1‘ 1 1 1 1

 

-0.5 '

 

 

1‘ 1 1‘ ‘i‘ 1 1 1 1 1 1 3 ‘ 1 1

0.5

(b)

73

Figure 6.4 Initial guess (dotted line), the true defect (solid line) and the prediction

(dashed line). (a) crack perpendicular to the current (b) crack at an angle of 55

degrees with the current direction

6.2.3 Conclusion

A simple parametric inversion approach for reconstruction of a tight crack has been
validated using EFG model. A calibration curve is applied to characterize the crack
location, length and orientation. For this simplified problem, a simple uniform
discretization in EFG model is sufficient. The inverse model is simple, fast and

efficient.

6.3 Inversion based on State Space Search Techniques for 2D Problem

A model based inversion using EFG model and state space search method is
presented in this section for 2D defect profile reconstruction. The iterative state space
search method using the tree structure is developed for implementing the defect
updating scheme. The inversion scheme is applied to characterize the cross-section of
the defect in terms of both width and depth.

The performance is again evaluated using simulated input signals.

The approach proposed in this section comprises three major parts:

(1) Problem representation — defines the problem geometry and introduces the

defect representation using state space

(2) Search procedure — uses the tree representation to search for the optimal

defect

(3) The dependence of results of cost function

74

6. 3.1 Problem Representation

 

 

 

(a)

 

X-Z Pinned

 

(b)
Figure 6.5 The EFG model Geometry (a) 3D View (b) cross-section view

The geometry for the inverse problem is shown in Figure 6.5. Figure 6.5 (a) shows
the 3D geometry of an infinite conducting plate with a defect. The time-varying
harmonic current in the y-direction is applied using a thin foil as the excitation source.
Since the geometry extends to infinity in the y—direction this can be simplified to a 2D
problem as shown in Figure 6.5 (b), which is the cross-section of the 3D geometry

and defect profile.

75

In order to derive a state space [62]-[63] representation, the cross-section of the
defect is discretized in terms of M discrete cells along the length and N layers of
nodes along the depth into the test specimen. The state X is a MxN matrix with

element 1 or O, which represents presents or absence of a defect respectively in each

cell. The set of all possible states (ZMN) is the state space. Figure 6.6, shows three

example states/defect profiles, where (a) and (b) are the surface defects and (c) is a

sub-surface defect.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(a) (b) (C)
Figure 6.6 Examples of States (Defects)

An initial state describes a possible defect from which the iterative search will start.
An objective state is a defect representation that minimizes the cost function.
The search procedure is performed layer by layer. The best defect profile at each

layer is obtained by minimizing the cost function, which is commonly defined as
N 2
F = XI Bn - 3m |
=1

(6.4)

where B" is the model prediction of imaginary part of the normal component of
magnetic flux, an is the corresponding experimental signal and N is the number

of data points in the signals. The cost function is computed for all the possible values

76

for cells in the current layer and the state, corresponding to the lowest cost function is

selected in each layer.
6. 3.2 Tree Structure

In order to derive the solution to the inverse problem, we search for a simple
best-fit defect state to minimize the difference between the signals corresponding to
true and reconstructed defects. The exhaustive search procedure is summarized using
a tree structure, which is shown in Figure 6.7.

The search is constrained by two assumptions, namely, (a) the defect grows

vertically downwards. (b) the defect grows narrower with depth.

An initial node d( ) is chosen as the imtial state in the tree structure. In each layer,

the nodes are expanded, and the cost function of each node is computed and the node
corresponding to the minimum cost function is selected and retained for expansion to
the next layer. In Figure 6.7, the node inside the square represents the selected
minimum cost nodes in each layer.

The tree expansion continues, until a leaf node is reached. The criten'on to obtain
the leaf node is based on the calculation of the cost function. As the tree expanded to
a next layer, if the minimum cost function begins to increase, then the search stops,
and the selected node in the current layer corresponds to the global minimum cost
function and is labeled as the leaf node. Once a leaf node is obtained, the search
procedure stops and the solution state is reached. In Figure 6.7, the cost function of

the bottom layer gray node is seen to be greater than the minimum cost function of the

77

parent node in the previous layer. The node inside the square is the desired defect

   
 
   
   

 

profile.
6(0) Cost Function CO
0
. . Minimal Cost Function C1<CO
O O
, Minimal Cost Function CZ<C1
. 0 Minimal Cost Function c3<c2
O 0
Minimal Cost Function C4>C3
Figure 6.7 the tree structure

6. 3.3 The Search Procedure

The detailed search procedure is divided into the following steps:

1. Estimation of the initial state (1(0). — involves determining the defect location i.e.
whether the crack is on the top or bottom surface.
(a) Use simple peak detection algorithm to determine the defect location. The

peak-to-peak separation in the input signal is used to estimate the defect length in

the initial layer.

78

 

Figure 6.8 (a) Input signal (b) initial guess of top surface defect (c) initial guess of
bottom surface defect

><1o‘4

 

 

 

 

 

 

 

 

 

 

(b)

79

Figure 6.8 continued

 

 

 

 

 

 

 

 

(C)
(b) In Figure 6.8, (a) represents the input signal, and the peak-to-peak separation of
the signal is from -3 mm to 3 mm, which means the defect length of the initial state is
-3 mm to 3 mm. The location of the defect of length 6 mm could be on the top or
bottom surface as shown in Figures (b) and (c). These two defect geometries are
simulated and the signals are generated. The cost functions are calculated and the

. . . 0
model With lower cost function is chosen as d( ).

2. Profiling the defect depth

Perform exhaustive search in the each layer. Generate the possible allowed defect
shapes in current layer using the defect assumptions and estimate the defect boundary
corresponding to the minimum cost function. As shown in Figures 6.9 (a) through (i),
the possible defects are generated and modeled. Cost fimction associated with each

node is computed.

80

Figure 6.9 Candidate defect shapes used in the search in layer 2

 

(b)

 

 

 

 

(a)

 

 

 

 

 

 

 

 

.::

 

 

 

 

 

 

 

(d)

(C)

 

 

 

(f)

(e)

3. Check for stopping criteria

If the minimum cost function in a layer is larger than the corresponding minimum
cost fimction in the layer above it the search stops. The geometry with the global
minimum cost function provides the desired defect profile or the solution to the

inverse problem.
6.3.4 Initial Results

The inversion strategy is implemented and the results corresponding to two defect
profiles are presented. Linear excitation current of 3 kHz frequency is applied along
the x-axis as shown in Figure 6.5 (a). The upper strip region is the current foil. The
lower strip region represents a section of the aluminum plate. The background dots
are the discretization nodes of the background mesh used in the EFG model. As
before the units are in ms.

The cost function is defined as in equation (6.4).

The measured signal shown in Figure 6.8 (a) is obtained for the true defect profile
in Figure 6.10 (a). The initial guess of a top surface defect is shown in (b) with cost
function 0.1708 and the initial guess of a bottom surface defect is shown in (c) with

d(0)

cost function 0.3198. This leads to the estimation of as a top surface defect. The

reconstructed defect profile is shown in Figure 6.10 ((1).
Similar results for bottom surface defect are shown in Figure 6.11. Figure 6.11 (a)
is the true defect profile. The initial guess of a top surface defect shown in (b) has cost

function 0.3766 and the initial guess of a bottom surface defect shown in (c) has cost

82

 

function 0.1770, so that the initial defect profile is determined as a bottom surface
defect.

When the given input signal is noise free, the reconstructed defect profile is exactly
equal to the true profile.
Figure 6.10 (a) True defect profile and (b) initial guess of a top surface defect (c)

initial guess of a bottom surface defect (d) reconstructed results for top surface
defect

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

83

Figure 6.10 continued

 

 

 

 

 

 

 

 

(C)

 

 

 

 

 

 

 

 

(d)

84

Figure 6.11 (a) True defect profile and (b) initial guess of a top surface defect (c)
initial guess of a bottom surface defect (d) reconstructed result for bottom surface
defect

 

 

 

   

nnnnnnnnnnnnnn

 

 

 

 

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

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

llllllllllllll

--------------

vvvvvvvvvvvvvv

oooooooooooooo

uuuuuuuuuuuu

--------------

 

 

oooooooooooooo
nnnnnnnnnnnnnn
nnnnnnnnnnnnnn
..............

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

 

aaaaaaaaaaaaaa

oooooooooooooo

--------------

nnnnnnnnnnnnnn

--------------

 

 

 

85

Figure 6.11 continued

 

 

 

 

 

 

 

 

 

 

 

 

6.3.5 Choice of Cost Function

The objective of inverse problem is to estimate the “best” profile of the detected
defect, that minimizes the cost function. Hence different cost functions may lead to

different resultant profiles. Hence we need to choose cost function with care.

The cost function acts as the evaluation fimction to describe the error between the
model prediction and input signal, which is obtained from the true defect geometry.

The most commonly used cost function is based on the [Q norm defined as

2
F: an'—an|

1M2

(6.5)

where B" is the model prediction of imaginary part of the normal component of
magnetic flux density, an is the corresponding measured signal and N is the

number of data points in the signals. The search algorithm developed in the above
section is implemented using the cost function defined in Equation (6.4). However
this definition does not yield the best results in all cases. An alternative cost function
that can be employed is the L00 (or Chebychev) norm:
F = max(| Bn - an |),i =1,2,...K
(6.6)
This cost function is defined as the maximum difference between the input signal
and model prediction, among all the data points. The search algorithm is performed
with the new cost function. The results show that the new cost function yields more
accurate defect profiles than the previous cost function when the defect under
consideration is very deep.
Two true defect profiles of larger depth are shown in Figures. 6.12, where Figure
6.12 (a) is top surface defect and 6.12 (b) is bottom surface defect. The reconstructed

top surface defect profile obtained with the Q cost function in equation (6.5) is

87

shown in Figure 6.13 (a), and the corresponding result obtained with the new LOD
cost function in equation (6.6) is shown in Figure 6.13 (b). When the given input
signal is noise free, the reconstructed defect profile is seen to be exactly the same as
the true profile using the new cost function. The old cost function does not converge
to the desired defect profile. Similar results for bottom layer defect are shown in
Figure 6.14 (a) and (b), with the old and new cost functions respectively. The

reconstructed defect profile gives the true profile using the new cost fimction.

 

 

 

 

 

(a)

 

 

 

 

 

 

 

(b)
Figure 6.12 True defect profile (a) top surface (b) bottom surface

88

 

 

 

 

 

(a)

 

 

 

 

 

Figure 6.13 Reconstructed top surface defect profiles (a) with old cost function (b)

with new cost function

03)

89

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

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

--------------

cccccccccccccc

I’

 

 

 

Figure 6.14 Reconstructed bottom surface defect profiles (a) with old cost function
(h) with new cost function

6. 3. 6 Conclusion

An inversion technique using a tree search algorithm is presented and validated for
a 2D eddy current problem with foil excitation. Both the surface and subsurface
defects are reconstructed using the state space search. The reconstructed defect is

shown to be dependent on the choice of cost function minimized.

90

 

Implemented of the technique on experimental signals is currently underway. The
performance and robustness of the inversion technique will then be evaluated on field

signals from aircraft inspection.

6.4 Inversion based on Gradient Search using Adjoint Equation for 2D

Problem

This section presents a defect update approach that makes use of a gradient search
of the multi-dimensional error surface so that the reconstructed defect profile will
result in a signal that matches the given measured signal in the least squares sense. As
shown in section 6.3, currently, the performance is evaluated for the simulated input
signals.

In the gradient based minimization algorithm, the evaluation of the gradient is
essential. To evaluate the gradient of the objective function, each of the partial
derivatives must be evaluated appropriately. The commonly used methods solve the
governing equations each time to calculate each of the partial derivatives. This will
result in a long computational time. To solve this problem, an adjoint equation [64]
based method is proposed to find all the partial derivatives by solving the governing
equations only once. The method reduces the computational time significantly.

When the EFG method is used as forward model, when the defect is updating, the
background discretization nodes remain unchanged and only a few nodes to represent
the defect parameters need to be updated. This Will greatly speed up the forward

computation and reduce the time needed in iterative inversion.

91

This method as well as its robustness is discussed in this section. Preliminary

results validating the approach are presented.
6. 4.1 The Inversion Scheme in 2D

The overall iterative approach is the same as that shown in Figure 6.1. The
geometry for two-dimensional inverse problem is shown in Figure 6.5 along with the
x -, y -, z - directions. The defect is represented in terms of nodal coordinates which
determines the defect profile on the boundary.

The inverse problem solution is the profile that is obtained by minimizing the cost

function is defined as
N 2
F = XI Bn “an I
=1

(6.7)

where B" is the model prediction of imaginary part of the normal component of
magnetic flux, an is the corresponding input signal and N is the number of data

points in the signals. The cost function is computed for each defect profiles and the
best defect profile is selected corresponding to the lowest cost function. The iteration
stops when the cost function is minimized and begins to increase.

A typical iteration scheme consists of four components:

1. a set of state variables (a ;

2. a set of defect parameters X ;

3. an objective function J ((0, X);

92

4. constraint equations F ((0,X )=0.
In the above 2D problem, go represents the magnetic vector potential in the

forward model. The defect parameters X are the nodal coordinates defining the
defect profile. The objective function is the cost function F in equation (6.7).

The constraint equations is the governing partial differential equation in terms of
the parameter set X and unknown to. The matrix equation obtained in the EFG
method is written as

F(¢,X)=K.¢—S=0
(6.8)
where matrix K is the stiffness matrix and S is the source vector.

So the problem can be stated as: iteratively solve the set of equations (6.8) for the

defect parameters X such that the objective function J ((p, X) is minimized,
subject to the constraint equations F ((9, X) = 0.
The commonly used gradient search method is summarized as follows:
1. Assume an initial defect parameter set X (0)
2. Solve the constraint equation (6.8), and calculate the predicted signal B,- , for
i=l,2,. . ...M. Calculate the objective function fi'om equation (6.7).
3. Compute the gradient of the objective fimction VJ over X, and determine
éXm - an updating step.

4. Update the defect parameters X (n+1) = X (n) + 6X(").

93

 

5. Repeat from step 2, until the cost function is minimized and the stopping criterion
is reached.

The details of step 3 are presented below:

Assume X=[x1,x2, ...... ,xm]T and VJ=[fl- 3]— if

6x1,6x2’ ...... ,axm

a .
Foreach xk,k=l,2, ...... ,m 6J —Z.fl.._¢1.=[fl][a_¢]

’ExZ”ja¢j 6xk ago ax
(6.9)
61 w a1 61 T
where [—]=[—,——, ...... ,—] ,and n is the number of the unknown (0.
6¢ a¢l 6¢2 aion

In equation (6.9), [age] can be easily obtained by differentiating equation (6.8),
x

which gives:
6o) 6K
K - — + — = 0
[ ] [51k] [an lb]
(6.10)
We also need to compute [SE] , which is done using the adjoint equation [60]. A
(P

constrained cost function L(¢, X ,A) is defined by (6.11) as
L(¢,X,A) = J (¢,X ) - (F (¢,X ),A)
(6.11)

where A is a set of Lagrangian multipliers. Differentiating equation (6.11) w.r.t. the

unknowns we get:

94

 

-a—L—=O:>F(¢,X)=O
6A

(6.12)
fi=0=[ai].[A]=[fl]
6X 6X 5X
(6.13)
6_L_ r. _ a_J
a¢—0:>[K] [ALB]
(6.14)

By substituting (6.13) and (6.14) into equation (6.9), the derivative can be obtained

as:
aJ 6c) 1‘ 6K
_= . K . _ =_ . __ .
an NT [ 1 [ax] [A] [ax] in]
(6.15)
In the gradient search method,c$X(") = [fl 6.] if , so the updating is

ax] , 5;: , ...... , axm
performed using equation (6.15).

For simplification, we assume the Lagrange multiplier A is a constant number in
each iteration step. Hence equation (6.15) is rewritten as E- = -c - [6_K_] - [p] ,
axk ax]:
where c is a small constant number, chosen to be equal to 0.1.

6. 4.2 Results

The implementation of the above iteration procedure is presented here.

95

 

The true defect profile in the aluminum sample (lower strip region) and the
excitation current foil (upper strip region) are shown in Figure 6.15 (a). The desired
defect is a rectangular defect with the width extending from -0.45 mm to 0.45 mm,
and with height 1.0 mm. The initial guess is a crack-free situation as shown in Figure

6.15 (b). The background nodes are the discretization nodes in the EFG method.

96

 

 

 

 

 

 

 

 

 

 

 

 

 

 

llllllllllllllllllllllllllllllllllllll
oooooooooooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooooooooooooo
llllllllllllllllllllllllllllllllllllll

IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII

nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

 

llllllllllllllllllllllllllllllllllllll
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu
......................................
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
oooooooooooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooooooooooooo
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
--------------------------------------
llllllllllllllllllllllllllllllllllllll
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
oooooooooooooooooooooooooooooooooooo
--------------------------------------
.......................................

nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn

 

 

oooooooooooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooooooooooooo
--------------------------------------
o o I I A o o a u o - I o a n u I n a a v a A u - o I nnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
--------------------------------------
oooooooooooooooooooooooooooooooooooooo

oooooooooooooooooooooooooooooooooooooo

 

oooooooooooooooooooooooooooooooooooooo

 

 

l

2 -1.5 -1 -0.5 0 0.5 1 1.5 2

(b)

Figure 6.15 (a) The desired defect profile. (b) the initial guess

The iteration steps are summarized below:
Determination of defect location and parameter X

(3) Assume the crack is surface-breaking (top surface).

97

I. I—h;

(b) The width of the defect on top surface is determined by the distance between two
peaks in the input signal, which is around -0.45 mm to 0.45 mm.
(c) We assume four nodal points to represents the defect profile. The defect parameter

X is defined as X = {x1,x2,x3,x4,21,22,23,z4}. For simplification of calculation,

the x-coordinates are fixed and equally spaced according to the width of the defect on
the surface. So the x-coordinates in this example are -0.45 mm, -0.15 mm, 0.15 mm,
and 0.45 mm. The z-coordinates are the variables that need to be optimized.

In Figure 6.16 (a), the four black points represents the position of the initial guess
of the defect parameter X .

11. Updating the z -coordinates

(d) From the simplified equation (6.15), 9;]. = _C,I: 321g
k

].[¢] , the coordinates are
62k

updated as Z("+1) =Z(") + 6.] . Where the derivative of the matrix K is

621:")

 

, here we assume

 

calculated numerically as [a—K] = K(zk + Azk ) - K(Zk)

62k (Zk+AZk)-Zk

A2,, = 0.001.

Figures 6.16 (b) to (k) show evolution of four point the defect profile in each
iteration.

111. Check for stopping criteria
(e) In each updating step, the cost function is calculated using equation (6.7). If the

cost function keeps decreasing, the searching procedure continues. When the cost

98

 

function starts to increase, the iteration procedure stops and the optimum defect
profile is obtained corresponding to the minimum cost function.
In Figure 6.16 (l) the optimum defect profile is obtained corresponding to the

minimum cost function.

Figure 6.16 The updating status of the defect profile

 

 

 

 

 

 

 

 

 

 

 

 

 

   

 

OIICOIIDIOOIIDIOIDIOIIUb-IDIOOVUIIO ll ...h..‘.-l \IOIIIOI‘IGI‘IKLIDCUIIQI II
II-OUOIIiii-ll-I-IIIIUIIIIJDIII-IIIDIh ‘CI'ID-I'IIPIIOI'IIG’IOIICCDIUIU‘|IDDI
..OIIIIoA-ullnlI'DIVIIIO‘IIO'VIll... I ...‘I'II-IIIIUI“’-."I'I'lIl'.‘I.'l.'
IOII...‘IJIOIIUIIQDIOIDQOIOII It‘ll... .IIIIQD1I0'IDIQDIIIIODIIIIOIIIIDIII§II
...u..... ....V..V..f.-Y...........yuu 4........ .... .--........-y-
I0!DnlllllIICIDlI-DA'IVOUI Inn-hill... Ont-llIIOtlil‘illlIII-OIADOII-OIIAl ll
"..-."II..II.""I...‘II'..."III..I ..IIIOIIIII\ICII'.'I'I'.I"I.l."..'.'
'QOIIOI‘IUIOI‘IU'11lIIOQIIOICUIIII 0" "-“.I""I"“"'l.l."‘.'.."|".ll
IIIIOIO'<~ODlr'l'lcIvllt-‘htufi-'inIll COI-Clive-IVIOIOID-0IIC.-..IIOICDIO}IU
IOIIIIlll'IID tlelllI-lIl-lrlllllllnu I‘ll-lII)II'lIl-¢OIIAIAAIICIIIO-‘Dlhln
HUG...I.I.IIIu.II-II'IIQIIIDI\.I.II..O ..ill.l~III-..IllI.“.lI...I‘l.l."..'
I-lIOl:IIIOII>ulllO‘IIIJOIIIIIQODCI II Cl)It!l-IlCIIOIIJIIOIIIDIOAIDIIInflblll
.‘IIIIII.I'-|.Il.-‘.'II."-IIDI.I...-. IDDIUIQO‘OQIIV \IllA-lnaiouua
O'IlllllllfillI-ll llllllloll-o-I lllll lb‘IOl'v-utllblI‘ll-IIIIII\IQIIIII"II
IDOIC‘QIIIIOOSVIIODO‘QIDIIDIIOdlfitdbll illbilh‘l‘ll..‘|hl‘l..l'l'."..5...‘..
Ill...I‘dlillIOCQIIIOIII‘UII-alCIIU‘II I‘DIIIIIIIOIlutI‘IIIIIIDIIIIUCIAIOQIII
[DOI’Il‘.".l'.l'lI.‘-IAI‘.I'I."I'.~C cac'vlntnloonuvsunlvo~hllo4r¢-I-‘.1.t.
"‘Utll..ltu—ICIOII‘I"‘IIv..t}l..ll" I‘tldv!4lvl¢IAgloultrotdltfilufilluboil'
..II'III‘.".|'I.-""III.'I"""I..I "’l"lD"...‘l".l"'|.'. IIIIIIIOD"
lonllI-OlbiltlDI:.II.--OIIDI41IIIII|C- vtnllfll‘ll-‘AOI.:...A-ndl-Ollll-vlitl.

 

 

 

 

 

C

99

 

Figure 6.16 continued

 

 

 

 

 

 

(

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

)

(g)

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

')

(i)

l

Figure 6.16 continued

 

oooooooooooooooooooooooooooooooooooooo

.I-t-oo-bI-IA-u-nn nnnnnnnnnnnnnnnnnnnn

.......................................
......................................
--------------------------------------
ailnpitbtlnl ..........................
......................................

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

 

 

 

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

 

oooooooooooooooooooooooooooooooooooo

-------------------------------------

 

 

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

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

 

......................................
......................................
......................................
--------------------------------------
......................................
......................................

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

 

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

---------------------------------------

 

 

 

 

 

6. 4.3 Effect of Noise Level

The above initial results are obtained when the input signal is noise free. While in
realistic problems, the noise level cannot be neglected. The effect of noise on the
performance of the gradient search method is studied to test the robustness of the
inversion scheme. Random white noise varying from 5% to 30% of the peak signal
value was added to the model predicted signal and the corresponding reconstruction
defects were generated.

The true defect profile is shown in Figure 6.17.

101

 

 

 

 

 

 

 

 

 

 

 

-2 -1.5 -1 -0.5 0 0.5 ‘l 1.5 2

Figure 6.17 True defect
Figures 6.18 (a) to (1) show the signals corrupted with 5.0%, 10.0%, 15.0%, 20.0%,

25.0% and 30.0% random noise.

102

Figure 6.18 The signal corrupted with (a) 5.0%. (b) 10.0%. (c) 15.0%. (d) 20.0%. (e)

25.0%. (f) 30.0% random noise

.104

Signal with 0.05 Noise

 

 

l

 

 

-1

Signal with 0.10 Noise

-0.5

0

(a)

 

 

 

 

-0.5

103

240'3

Figure 6.18 continued

«04 Signal with 0.15 Noise

 

 

 

 

 

5
4 r .
3 - -
2 . .
1 . .
o - .
-1 » .
.2 l -
_3 . -
4 . .
'5 -2 41.5 :1 as 0 _ 035 i 1:5 2.104
(C)
5 .104 Signal with 0.20 Noise
4 - .
3 . q
2 . .
1 . .
o - .
.1 - -
.2 ~ -
.3 - -
-4 - .

 

 

 

'5-2 4.5 -1 -o.5 o 0.5 1 1.5 2.104

(d)

Figure 6.18 continued

:104 Signal with 0.25 Noise

 

1

 

 

 

 

(e)
«04 Signal with 0.30 Noise

 

 

 

 

'5-2 -1.5 -1 -0.5 0 0.5 1 1.5 2,104

105

Figures 6.19 (a) to (f) show the reconstructed defect profiles (dashed lines)
corresponding to the signals in Figure 6.18 which are corrupted with 5.0%, 10.0%,
15.0%, 20.0%, 25.0% and 30.0% random noise. The true defect is plotted in the same

figure (solid lines) for comparison.

Figure 6.19 (a) to (i) The reconstructed defect profile (dashed lines) when the signal
is corrupted with (a) 5.0%. (b) 10.0%. (c) 15.0%. (d) 20.0%. (e) 25.0%. (f) 30.0%
random noise

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

106

 

 

Figure 6.19 continued

 

 

LhLLk

 

ALL

 

 

 

 

 

 

(C)

 

 

LL'

L

glL

 

 

 

 

 

 

 

 

(d)

 

 

 

 

 

 

 

'4;

 

 

 

 

(e)

107

Figure 6.19 continued

 

 

 

 

 

 

 

 

 

 

6. 4.4 Conclusion

Model-based inversion technique using the EFG numerical forward model is
presented for a 2D problem. The updating scheme using gradient search in
conjunction with adjoint equation method is developed for the EFG model. In the
iteration procedure only the few coordinates that define the defect profile are updated
while the overall background nodes remain unchanged, which simplifies the solution
procedure and saves computational time. The gradient search method is a fast and

reliable and the results show its robustness in the presence of noise.

108

 

CHAPTER 7. 3D INVERSE PROBLEM USING EFG

METHOD

7.1 Introduction

Chapter 6 discussed the development and implementation of several inversion
schemes for 1D and 2D defect characterization. However in practice, problems are
generally 3D in space, which means that the solution domain and defects are
represented in the x—, y-, and 2- directions. The high-dimensionality is a major
challenge for the inversion scheme [65]-[69]. The defect must be represented by 3D
defect parameters, for example, the length, width, and height, which will make the
defect reconstruction involve large number of unknowns and hence more complex.
The high-dimensional inverse problem is in general ill-posed, since the measured
signals lack the continuous and reliable dependence on the defect profiles and also the
results can be non-unique.

In this chapter, the three dimensional inverse problem is addressed by modifying
the 2D state space search method described in chapter 6. These modifications make it
possible for 3D defect reconstruction. Both surface and subsurface defect

reconstructions are presented for validation and discussion.

109

 

7.2 Problem Statement

Figure 7.1 (a) shows the 3D geometry of an infinite single layer conducting plate
with a defect. The time-varying harmonic current in the y-direction is applied using a
thin foil as the excitation source. Figure 7.1 (b) is the 2D top view of the defect
profile, which is a rectangular shaped surface breaking defect. Figure 7.1 (c) is the
cross-section view of the defect profile. As defined before, a subsurface defect can

also be modeled when the defect initiates from the bottom layer of the plate.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(C)

Figure 7.1 The model geometry (a) 3D view, (b) 2D top view and (c) cross-section
view

The following parameters are used in the forward model:
(1) the thickness of the aluminum plate: 4 mm

(2) the profile of the ideal defect: 3 mm x 3 mm, with depth of 3 mm

110

(3) the excitation current density: 108 A/m2

(4) the excitation current frequency: 3 kHz

(5) the solution domain: [—15,15]mm x[—15,15]mm x {—8, 8] mm
(6) the sensor thickness: 0.2 mm; liftoff: 0.5 mm

The EFG method is used as forward model to obtain the discrete solution domain.
7 .3 3D State Space Representation

In order to reconstruct the 3D defect, the 3D state space representation should be
first developed. This section describes the three dimensional states based on the 2D
definition. The defect is discretized into L layers of nodes along the depth into the test
specimen, which is in z - direction; and in each layer, the surface profile of the defect

is discretized in terms of MxN discrete cells in xy- plane. So a 3D state X is a

MxNxL matrix with element 1 or 0, representing absence or presence of a defect in
the corresponding cell, respectively.

The state space now consists of (ZMNL) possible states.

For simplification we consider a profile of the rectangular defect in Figure 7.1
represented by a 3x3 matrix in xy - plane, and 4 layers of nodes along the depth in

z -direction.
For the 3D defect reconstruction, instead of the line plot, the 2D image of the

normal component of magnetic flux density Bz , which is observed at sensor liftoff

0.5 mm, is used in calculating the cost function. Experimental measurements of the

111

real, and imaginary values of 82 signals can be obtained from the GMR system. In

this thesis however, the input signal from the true defect is simulated using the
forward EFG model.

The cost fiinction is defined to search for the best defect profile, which is:

K1 K2

F: 2 Zlen_Bmmn '2

m=1n=l

(7.1)

where an is the model prediction of the normal component of magnetic flux,
BMW; is the corresponding experimental signal and K 1,K 2 are the number of

signal points in the x- and y-direction in the raster scan. The cost function is

computed for all the possible states of all layers, and the defect corresponding to the

lowest cost function is selected as the objective defect.

7 .4 Search Procedure

The search procedure is based on the same criterion as the tree search explained in

the Figure 6.7. Since the total number of state (2 ) is large, to avoid the tedious

exhaustive search, the search is constrained based on the following assumptions: (a)
the defect grows vertically downwards; (b) the surface defect profile is in a squared
region; (c) the defect grows narrower with depth.

The search procedure is performed layer by layer in z—direction; at each layer, the

3x3 node is expanded to generate the different profiles allowed by the assumption. In

112

order to reduce expansion into all the 29 possible nodes a simplification is introduced

as illustrated in Figure 7.2.

 

 

 

   

 

 

 

 

 

(a) (b) (0)

Figure 7.2 3x3 state space in each layer with three sub-state space
In the first layer, the search starts from the “middle defec ” as shown in Figure 7.2

(a), which can be one of six possible defect profiles as shown in Figure 7.3 from (a) to

(f).The optimum defect profile is selected as the one that minimizes the cost fimction.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 7.3 The candidate state space for sub-state in Figure 7.2 (a)

113

Once the optimum defect from the “middle defect” is estimated, then connected
cells from the “left” and “right” defects corresponding to minimum cost function
value are added to grow the defect in the current layer. When the optimum defect
profile at the current layer is estimated, the search procedure is repeated with the cells
in the next layer, until we reach the global minimum cost function which corresponds
to the leaf node in the tree structure. The desired defect profile is reconstructed and

the search stops.
7 .5 The Detailed Search Procedure

The detailed search procedure is divided into the following steps:
1. Estimation of the initial state (1(0) ;

This step involves determining whether defect location is on the top or bottom
surface.

In the case of the 2D defect a simple peak detection algorithm on the data

(imaginary part of 3,) is used to determine the defect location and the peak-to-peak

separation in the input signal is used to estimate the defect region in the initial layer.

In the 3D inversion scheme, the peak-to-peak separation in the imaginary part of Bz

is seen to vary with frequency whereas the real part is constant for both surface and

subsurface defects. Figure 7.4 (a) and Figure 7.4 (b) show the estimated defect region

114

when the fi'equency is varied from 500 Hz to 5000 Hz for surface and subsurface

defects.

 

 

3.4-

 

S"
N

 

 

the Eat. Defect Leng.&Wdh (mm)
to

 

 

 

 

 

 

 

 

 

 

 

 

 

 

3
2.9
2‘8 —°-ebeo|ute ' 3 c :
2.7I l l l l
0 1000 2000 3000 4000 5000
Operating Frequency
(8)
8 T f I I
-0-rea|
7'5 -+-imag‘nery ‘
5 7_ *ebeolute q
E
g 6.5- -
S: e- .
o
C
3 5.5- .
3 5- T
8
i 4.5" C... T.
m
2 4- -
3.5- .
I
30 1000 2000 3000 4000 5000
Operating Frequency
0))

Figure 7 .4 Plot of estimated defect region using real, imaginary and absolute data,

(a) surface defect, (b) subsurface defect

115

The initial node for the tree search (top or bottom) was therefore determined using

real part of 82. When the operating frequency is 3000 Hz, the defect region of the

initial state is 3.4 mm x 3.4 mm for surface defect and. 3.6 mm x 3.6 mm for the

subsurface defect. Then the estimated defect region is divided into 3x3 state space in
equal size cells. Both the initial defects, (1(0) surface and d(0) subsurface are modeled,

the signals and corresponding cost functions are generated and node with lower cost

C1(0)

function chosen as determines if the defect to be reconstructed is top surface or

bottom surface defect.
2. Profiling the defect depth

Perform the search procedure in each layer. Generate the possible allowed defect
shapes in current layer (expand the tree node) using the defect assumptions and
compute the associated cost function using the forward model prediction and
measured signal. The defect boundary corresponding to the minimum cost function is
selected to be expanded in the next level.
3. Checking for stopping criteria

If the minimum cost function in a layer is larger than the corresponding minimum
cost function in the layer above it the search stops. The geometry with the global
minimum cost function provides the desired defect profile or is the solution to the

inverse problem.

116

7.6 Initial Results

The inversion strategy is implemented on two defect profiles. First, the 3 x 3 x 3 mm
surface defect is assumed to be the true defect as shown in Figure 7.1, and the
operating frequency is 3000 Hz. The measured signal is obtained for the true defect
profile by simulating the EFG forward model and the corresponding real and

imaginary part of 82 are shown in Figure 7.5.

 

(a) (b)
Figure 7.5 The ideal signal of Bz (a) real part (b) imaginary part

The real part of B2 is used to estimate the defect boundary in the initial layer and

the imaginary part of B2 is used to calculate the cost function for minimization.

117

 

...‘CC‘OOOCO

0.0.9.000...

0'. ””00”....

.- eon-oaoeeone

CO eeuvmoeoo

O. 3.” ..”."O

 

u u“..”....
9. “...M....
0. ”Memo...

“oduomcbco
n

(a)

“:4: 1;
Kill?
2t"f1 32

 

Figure 7.6 Initial guess of a top surface defect (a) cross-section view (b) top view

118

 

0. no... 0......

00 eeeeeoeeoeeo

M 00000000.. 00

 

“0.0.0.000.“
00.00....I.COQOOOOOCOOIOOOOCOODOIOIOO.

ne-eeeee-e-oone-oeee-cee-eoeI-eeoeeo~-

 

 

000......OOOWCNOOWOOOOOQOOOO'QCCO...
l

 

 

Figure 7.7 Initial guess of a bottom surface defect (a) cross-section view (b) top view

Top surface defect - Using the peak-to-peak distance to estimate the initial state, a
top surface defect of size 3.4 mm x 3.4 mm and bottom surface defect, of size 3.6 mm
x 3.6 mm as shown in Figure 7.6 (a) and (b) and Fig. 7.7 (a) and (b) respectively were
simulated.

The cost function of top surface defect was 0.2569 whereas the cost function of the

bottom surface defect was 0.5251. This led d(o) to be a top surface defect.

119

The reconstructed 4-layer defect profile is shown in Figure 7.8 and the defect

profiles fi'om top to bottom layers are shown in (a), (b), (c), and ((1), respectively.

Figure 7.8 The reconstructed defect profile from top to bottom layer

 

.u
o
o
.e

 

u e‘o 0.0.. e.

:‘o‘ooo'o-ooo'

 

     

eeloeeeeeeee.

 

    

'eeeoeeeeeOIe-e
eeeeeeeeeeleee

eeoooooeeooecononoee-ee-eI-oweoaeon-ee_

eeooeeeeeeooeeeooeeeeeeeeeeeeeoeeeoee4

..e
O
o
m
e
c
O
o
I
o
0.,
e
a
O
o
e
I
O.
o
5.
O
U

 

 

120

Figure 7.8 continued

3mm,
1mm
~ um

 

The minimum cost function nodes in each layer and the associated cost flmctions

are plotted in Figure 7.9. The x-axis represents the all the intermediate nodes

obtained from considering “middle”, “left” and “right” states shown in Figure 7.2.

121

 

A

min. com fmction
o:

N

 

 

 

01 2 3 4 5 6 7 8 9

the intermediate nodes

Figure 7 .9 The cost function vs. node

Similar results for bottom surface defect are obtained. The true defect is a 3 x 3 x 3
mm subsurface defect. The operating frequency is 3000 Hz. The initial defect profile
is determined as a bottom surface defect using a minimum of two cost functions
corresponding to top surface and bottom surface defects respectively. The defect size
corresponding to the initial node is node is 3.6 mm x 3.6 mm.

The reconstructed 4-layer defect profile is shown in Figure 7.10 and the defect

profiles from bottom to top layers are shown in (a), (b), (c), and ((1), respectively.

122

Figure 7.10 The reconstructed defect profile from bottom to top layer

 

123

Figure 7.10 continued

  

(d)

The minimum cost flmction nodes in each layer and the associated cost functions
are plotted in Figure 7.11. The x-axis represents the all the intermediate nodes

obtained from considering “middle”, “left” and “right” states shown in Figure 7.2.

 

min.coetflnction
.° .° 9 .°
'9 a» a- "I

0
LA

 

3 4 5 6 'r
the intermediate nodes

 

<3

 

a
N

Figure 7.11 The cost function vs. node

7.7 Conclusion

A 3D inversion scheme has been developed based on the state space search
algorithm. The detailed defect representation and search procedure are addressed with
initial validation results. Both the surface and subsurface defects can be reconstructed
with reasonable accuracy. This is the first time the three dimensional defect
reconstruction problem has been attempted that includes both surface and subsurface
defects. However, a lot of work remains to be done to optimize the parameters of

different steps involved in the overall procedure.

125

 

CHAPTER 8. CONCLUSIONS AND FUTURE WORK

8.1 Original Contributions to the Field

An accurate and efficient numerical forward using the element free Galerkin
method has been developed, validated and applied to real 3D airframe geometry. The
forward mode] has been used to understand the underlying physics and visualize the
induced eddy currents in the case of different defect geometries. Development of
model based strategies for signal inversion to reconstruct the defect profile in
electromagnetic NDE is the primary objective of this dissertation.

The research carried out during the course of the dissertation has made the
following contributions to the field:

(1). The development and improvement of the Element-Free Galerkin (EFG) method
for electromagnetic field computations. The improvements include the asymmetric
domain of influence, the orthogonal basis function, and the Penalty boundary
condition. 1D, 2D, and 3D models for Poisson equation and diffusion equation,
describing static or low frequency quasi-static problems, have been presented as
examples to illustrate the validity and effectiveness of EFG method. This method is
shown to be better suited for modeling the complex solution domain compared with
the traditional mesh-based methods.

(2). The use of EFG method in signal inversion has also been investigated and applied

to reconstructing various one and two dimensional defect profiles by means of

126

parametric and non parametric inversion schemes. An iterative parametric algorithm
for one dimensional tight crack reconstruction was demonstrated. Two iterative
nonparametric algorithms are proposed, namely a state space search method and an
adj oint based gradient minimization method. These algorithms been validated and
perform with excellent accuracy in the prediction of the length and depth profile of
defects.

(3). The state space search method is modified and used in the three dimensional
defect reconstruction. These modified state space search scheme are fast and efficient
and avoid the tedious exhaustive search procedure without compromising the
performance. Both the surface and subsurface defects can be reconstructed with

acceptable accuracy.

8.2 Future Work

This dissertation addresses the inverse problem of predicting the defect profiles
using simulated signals using different approaches. However, the application to
inversion of experimental signals needs to be canied out.

Future work includes:

(1). Validation of the inversion schemes using the experimental Giant
Magneto-resistance data. Signal pre-processing may be needed before the data can be
used for inversion. Several experimental parameters need to be considered, such as
the GMR sensor bias, the effect of the scanning liftoff, the GMR sensor

characteristics to map the voltage signal back to the magnetic flux density fields.

127

 

More rigorous study of the cost function must be carried out to control the accuracy
of the solution.

(2). Inverse algorithms in this work are developed for the case of a single defect. The
possibility of multiple defects in practice, which are close to each other, will be more
challenging and needs to be addressed. This is left as a future work to investigate the
physical interaction of fields due to multiple defects.

(3). The inversion scheme needs to be optimized in each stage to improve the

efficiency for online processing after the inspection is completed.

128

 

REFERENCES

[1] D. E. Bray and R. K. Stanley, Nondestructive Evaluation, McGraw Hill, 1989

[2] R. B. Thompson and D. 0. Thompson, “Ultrasonics in Nondestructive
Evaluation”, Proceedings of the IEEE, Vol. 73, No. 12, Dec 1985

[3] D. E. Bray and R. K. Stanley, Nondestructive Evaluation, a Tool for Design,
Manufacturing, and Service, McGraw-Hill, 1989

[4] B. Yang, P. K. Liaw, H. Wang, J. Y. Huang, R. C. Kuo, and J. G. Huang,
Thermography: A New Nondestructive Evaluation Method in Fatigue Damage, JOM,
2003

[5] R. Pradeep, “Neural Network Based Iterative Approaches for Electromagnetic
NDE”, Ph.D Dissertation, Iowa State University, 2001

[6] J. Benson, Automated Analysis of Rotating Probe Eddy Current Data, Interim
Report, Dec. 2003

[7] S. Udpa and L. Udpa, “Eddy Current Nondestructive Evaluation”, Encyclopedia
of Electrical and Electronics Engineering, John G. Webster, Editor, John Wiley, 1999

[8] C. Chen, “Finite Element Modeling of M01 for NDE Applications”, Master
Thesis, Iowa State University, 2000

[9] J. D. Jackson, Classical Electrodynamics, John Wiley & Sons Inc, New York,
1998

[10] L. Udpa and S. S. Udpa, “Application of Signal Processing and Pattern
Recognition Techniques to Inverse Problems in NDE”, International Journal of
Applied Electromagnetics and Mechanics, Vol. 8, pp. 99-117, 1997

[11] M. Yan, M. Afzal, S. S. Udpa, S. Mandayam, Y. Sun, L. Udpa, and P. Sacks,
“Iterative Algorithms for Electromagnetic NDE Signal Inversion”, ENDE’ 97, Reggio
Ca1abria,ltaly, Sep. 1997

[12] S. Hoole, S. Subramaniam, R. Saldanha, J. Coulomb, and J. Sabonnadiere,
“Inverse Problem Methodology and Finite Elements in the Identification of Cracks,

Sources, Materials and their Geometry in Inaccessible Locations”, IEEE Transactions
on Magnetics, Vol. 27, No. 3, 1991

[13] H. A. Sabbagh, L. D. Sabbagh, and T. M. Roberts, “An Eddy-Current Model and

129

 

Algorithm for Three-dimensional Nondestructive Evaluation of Advanced
Composites”, IEEE Transactions on Magnetics, Vol. 24, No. 6, pp. 3201-3212, Nov.
1988

[14] Jacques Hadamard, Princeton University Bulletin, pp. 49-52
[15] Satish Udpa, “Handbook on Electromagnetic NDE Methods”

[16] J. Salon and J. Schneider, “A Comparison of Boundary Integral and Finite
Element Formulations of the Eddy Current Problems”, IEEE Transactions on Power
Systems, Vol. 100, pp. 1473-1479, 1981

[17] Nathan Ida, Numerical Modeling for Electromagnetic Non-Destructive
Evaluation, Chapman & Hall, McGraw Hill, 1984

[18] Oliver Riibenkonig, The Finite Difference Method (FDM)-An Introduction,
Albert Ludwigs University of Freiburg, 2006

[19] J. H. Hwang and W. Lord, “Finite Element Modeling of Magnetic Field/Defect
Interactions”, ASNT Journal of Testing and Evaluation, Vol. 3, No. 1, pp 21-25, 1975

[20] J. H. Hwang and W. Lord, “Magnetic Leakage Field Signatures of Material

Discontinuities”, Proceedings of the 10th Symposium on NDE, San Antonio, Texas,
1975

[21] W. C. Yen, “Finite Element Characterization of Residual Leakage Fields”, Master
Thesis, Colorado State University, 1978

[22] K lshibashi, “Nonlinear Eddy Current Analysis by the Integral Equation
Method”, IEEE Transactions of Magnetics, Vol. 30, No. 5, Sep. 1994

[23] N. Esposito et a1, “Modeling of Three-Dimensional Nonlinear Eddy Current
Problems with Conductors in Motion by an Integral Formulation”, IEEE Transactions
on Magnetics, Vol. 32, No. 3, May 1996

[24] J. R. Bowler, Y. Yoshida, and N. Harfield, “Vector-Potential Boundary-Integral
Evaluation of Eddy-Current Interaction with 3 Crack”, IEEE Transactions on
Magnetics, Vol. 33, No. 5, Sep. 1997

[25] Y. Yoshida and J. R. Bowler, “Vector Potential Integral Formulation for
Eddy-Current Probe Response to Cracks”, IEEE Transactions on Magnetics, Vol. 36,
No. 2, Mar. 2000

[26] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, “Meshless
Methods: An Overview and Recent Developments”, Computer Methods in Applied
Mechanics and Engineering, Vol. 139, Issues 1-4, pp. 3-47, 1996

130

 

[27] R. A. Gingold and J. J. Monaghan, “Smoothed Particle Hydrodynamics: Theory
and Application to Non-spherical Stars”, Royal Astronomical Society, Monthly
Notices, Vol. 181, pp. 375-389, Nov. 1977

[28] B. Nayrolesl , G Touzot, and P. Villon, “Generalizing the Finite Element Method:
Diffuse Approximation and Diffuse Elements”, Computational Mechanics, Vol. 10,
No. 5, pp. 307-318, Sep. 1992

[29] W. K. Liu, Y. Chen, S. Jun, J. S. Chen, T. Belytschko, C. Pan, R. A. Uras, and
C. T. Chang, “Overview and Applications of the Reproducing Kernel Particle
Methods”, Archives of Computational Methods in Engineering, Vol. 3, No. 1, pp.
3-80, Mar. 1996

[30] J. M. Melenk and I. Babuska, “The Partition of Unity Finite Element Method:
Basic Theory and Applications”, Computer Methods in Applied Mechanics and
Engineering, Vol. 139, Issues 1-4, pp. 289-314, Dec. 1996

[31] C. A. M. Duarte and J. T. Oden, “Hp Clouds-an Hp Meshless Method”,
Numerical Methods for Partial Differential Equations, Vol. 12, pp. 673-705, 1996

[32] N. Sukumar, B. Moran, and T. Belytschko, “The Natural Element Method in
Solid Mechanics”, International Journal for Numerical Methods in Engineering, Vol.
43, No. 5, pp. 839-887, Nov. 1998

[33] H. Wendland, “Meshless Galerkin Methods using Radial Basis Functions”,
Mathematics of Computation, Vol. 68, No. 228, pp. 1521-1531, 1999

[34] T. Belytschko, Y. Y. Lu, and L. Gu, “Element-Free Galerkin methods”,
International Journal for Numerical Methods in Engineering, Vol. 37, pp. 229-256,
1994

[35] J. Dolbow and T. Belytschko, “An Introduction to Programming the Meshless
Element Free Galerkin Method”, Archives in Computational Mechanics, Vol. 5, No.
3, pp. 207-241, 1998

[36] V. Cingoski, N. Miyamoto, and H. Yamashita, “Element-free Galerkin Method
for Electromagnetic Field Computations”, IEEE Transactions on Magnetics, Vol. 34,
No. 5, pp. 3236—3239, 1998

[37] S. A. Viana and R. C. Mesquita, “Moving Least Square Reproducing Kernel
Method for Electromagnetic Field Computation”, IEEE Transactions on Magnetics,
Vol. 35, No. 3, pp. 1372—1375, 1999

[38] S. L. Ho, S. Yang, J. M. Machado, and H. C. Wong, “Applications of a Meshless
Method in Electromagnetics”, IEEE Transactions on Magnetics, Vol. 37, No. 5, pp.
3198—3202, 2001

131

 

[39] L. Xuan, Z. Zeng, B. Shanker, and L. Udpa, “Element-fiee Galerkin Method for
Static and Quasi-Static Electromagnetic Field Computation”, IEEE Transactions on
Magnetics, Vol.40, No. 1, pp. 12-20, Jan. 2004

[40] Liang Xuan, “Finite Element and Meshless Methods in NDT Applications”, Ph.
D Dissertation, Iowa State University, 2002

[41] S. Liu, Q. Yang, H, Chen, G. Xu, and F. Liu, “Improvement of the Element -Free
Galerkin Method for Electromagnetic Field Calculation”, IEEE Transactions on
Applied Superconductor, Vol. 14, No. 2, pp. 1866-1869, Jun. 2004

[42] Y. Zhang, K. R. Shao, D. X. Xie, and J. D. Lavers, “Meshless Method Based on
Orthogonal Basis for Computational Electromagnetics”, IEEE Transactions on
Magnetics, Vol. 41, No. 5, pp. 1432-1435, May 2005

[43] Suzhen Liu, Qingxin Yang, Haiyan Chen, Guizhi Xu, and Fugui Liu,
“Improvement of the Element-Free Galerkin Method for Electromagnetic Field
Calculation”, IEEE Transactions on Applied Superconductivity, Vol. 14, No. 2, pp.
1866-1869, Jun. 2004

[44] C. Herault and Y. Marechal, “Boundary and Interface Conditions in Meshless
Methods”, IEEE Transactions on Magnetics, Vol. 35, No. 4, pp. 1450-1453, 1999

[45] T. Belytschko, D. Organ, and Y.Krongauz, “A Coupled Finite Element—Free
Galerkin Method”, Computational Mathematics, Vol. 17, pp. 186-195, 1995

[46] W. M. Xue and S. N. Atluri, “Hybrid and Mixed Finite Element Methods”,
Proceedings of the Winter Annual Meeting, Miami Beach, FL, pp. 91-112, Nov. 1985

[47] T. Belytschko and Y. Krongauz, “Enforcement of Essential Boundary Conditions
in Meshless Approximations using Finite Elements”, Computer Methods in Applied
Mechanics and Engineering, Vol. 131, pp. 133-145, 1996

[48] G. Yagawa and T. Furukawa, “Recent Developments of Free Mesh Method”,
International Journal for Numerical Methods in Engineering, 47, No.2, pp.
1419-1443, 2000

[49] G. L. Fitzpatrick, D. K. Thome, R. L. Skaugset, E. Y. C. Shih, and W. C. L.
Shih, “Magneto-Optic/Eddy Current Imaging of Aging Aircraft: A New NDI
Technique”, Material Evaluation, Vol. 51, No. 12, pp. 1402-1407, Dec. 1993

[50] L. Xuan, B. Shanker, L. Udpa, W. C. L. Shih, and G. L. Fitzpatrick, “Finite
Element Model for MOI Applications Using A-V Formulation”, Review of Progress
in Quantitative Nondestructive Evaluation, AIP, Vol. 557(1), pp. 385-391, Apr. 2001

[51] G. Rubinacci, A. Tamburrino, A. Ventre, F. Villone, L. Udpa, L. Xuan, and Z.

132

 

Zeng, “Numerical Simulation of Magneto-Optic Eddy Current Imaging”, 8th
International Workship on Electromagnetic Nondestructive Evaluations, Saarbrucken,
Germany, Jun. 2001

[52] R. Albanese et al, “A Comparative Study of Finite Element Models for Magneto

Optic Imaging Applications”, 6th International Workshop on Electromagnetic
Nondestructive Evaluations, Budapest, Hungary, Jun. 2000

[53] R. W. Freund, “Transpose-Free Quasi-Minimal Residual Methods for
Non-Hermitian Linear Systems,” Advances in Computer Methods for Partial
Differential Equations, pp. 258-264, IMACS, 1992

[54] L. Udpa, R. Pradeep, J, Benson, and S. Udpa, “Automated Analysis of Eddy

Current Signals in Steam Generator Tube Inspection,” 16th World Conference on
Nondestructive Evaluation

[55] R. Sikora, J. Sikora, E. Cardelli, and T. Chady, “Artificial Neural Network
Application for Material Evaluation by Electromagnetic Methods”, Proceedings of the
International Joint Conference on Neural Networks, Vol. 6, pp. 4027-4032, 1999

[56] W. Qing, S. Xueqin, Y. Qingxin, and Y. Weili, “Using Wavelet Neural Networks
for the Optimal Design of Electromagnetic Devices”, IEEE Transactions on
Magnetics, Vol. 33, No. 2, pp. 1928-1930, Mar. 1997

[57] S. Hoole, “Artificial Neural Networks in the Solution of Inverse Electromagnetic
Field Problems”, IEEE Transactions on Magnetics, Vol. 29, pp. 1931-1934, 1993

[58] R. Popa, and K. Miya, “A Data Processing and Neural Network Approach for
Inverse Problems in ECT”, ENDE, R. Albanese et al. eds., 108 Press, pp. 297-305,
1998

[59] R. Sikora, T. Chady, and J. Sikora, “Neural Network Approach to Crack
Identification”, International Journal of Applied Electromagnetics and Mechanics,
Vol. 9, pp. 391-398, 1997

[60] Yue Li, “Edge Based Finite Element Simulation of Eddy Current Phenomenon
and its Application to Defect Characterization”, Ph. D Dissertation, Iowa State
University, 2001

[61] Guizhong Liu, Yue Li, Y. Sun, P. Sacks, and S. Udpa, “An Iterative Algorithm
for Eddy Current Inversion”, Review of Progress in Quantitative Nondestructive
Evaluation, Vol. 19(A), pp. 497-504

[62] L. Udpa, “An Iterative Numerical Algorithm for the Inversion of Eddy Current
NDT Data”, Electrosoft, Vol. 2, No. 23, 1991

133

 

[63] L. Upda and W. Lord, “A Search-Based Imaging System for Electromagnetic
Nondestructive Testing”, Expert, IEEE, Vol. 4, Issue 4, pp. 18-26, 1989

[64] M. Gunzburger, “Adjoint Equation based Methods for Control Problems in
Incompressible Viscous Flows”, Design Optimal et MDO, Center de Recherche en
Calcul Applique, Montreal, pp. 1-25, 1998

[65] Y. Li, L. Udpa, and S. S. Udpa, “Three-Dimensional Defect Reconstruction from
Eddy Current NDE Signals”, IEEE Transactions on Magnetics, Vol. 40, Issue 2, pp.
410-417, 2004

[66] P. Rarnuhalli, L. Udpa, and S. S. Udpa, “Electromagnetic NDE signal inversion
by function-approximation neural networks”, IEEE Transactions on Magnetics, Vol.
38, Issue 6, pp. 3633-3642, 2002

[67] S. M. Nair and J. H. Rose, “Reconstruction of Three-Dimensional Conductivity
Variations from Eddy Current (Electromagnetic Induction) Data”, Inverse Problems 6,
pp. 1007-1030, 1990

[68] A. Tambunino, R. Fresa, S. S. Udpa, and Y. Tian, “Three—Dimensional Defect
Localization from Time-of-Flight/Eddy Current Testing Data”, IEEE Transactions on
Magnetics, Vol. 40, Issue 2, pp. 1148-1151, 2004

[69] A. Joshi, L. Udpa, and S. S. Udpa, “Use of Radial Basis Function Neural
Network to Predict 3-Dirnensional Depth Profile from Magnetic Flux Leakage

Inspection”, National Seminar on NDT by Indian Society of Nondestructive Testing
(ISNT), Dec. 20

134