NOVEL SIMULATION AND DATA PROCESSING ALGORITHMS FOR EDDY CURRENT INSPECTION By Anton Efremov A DISSERTATION Michigan State University in partial fulfillment of the requirements Submitted to for the degree of Electrical Engineering – Doctor of Philosophy 2020 ABSTRACT NOVEL SIMULATION AND DATA PROCESSING ALGORITHMS FOR EDDY CURRENT INSPECTION By Anton Efremov Eddy Current Testing (ECT) is a widely used technique in the area of Nondestructive Evaluation. It offers a cheap, fast, non-contact way for finding surface and subsurface defects in a conductive material. Due to development of new designs of eddy current probe coils and advance of model based solutions to inverse problems in ECT, there is an emerging need for fast and accurate numerical methods for efficient modeling and processing of the data. This work contributes to the two directions of computational ECT: eddy current inspection simulation ("forward problem") and analysis of the measured data for automated defect detection ("inverse problem"). A new approach to simulate low-frequency electromagnetics in 3D is presented, based on a combination of a frequency-domain reduced vector potential formulation with a boundary condition based on Dirichlet-to-Neumann operator. The equations are solved via a Finite Element Method (FEM), and a novel technique for the fast solution of the related linear system is proposed. The performance of the method is analyzed for a few representative ECT problems. The obtained numerical results are validated against analytic solutions, other simulation codes, and experimental data. The inverse problem of interpreting measured ECT data is also a significant challenge in many practical applications. Very often, the defect indication in a measurement is very subtle due to the large contribution from the geometry of the test sample, making defect detection very difficult. This thesis presents a novel approach to address this problem. The developed algorithm is applied to real problems of detecting defects under steel fasteners in aircraft geometry using 2D data obtained from a raster scan of a multilayer structure with a low frequency eddy current excitation and GMR (Giant Magnetoresistive) sensors. The algorithm is also applied to the data obtained from EC inspection of heat exchange tubes in nuclear power plant. ACKNOWLEDGEMENTS I would like to express my sincere appreciation to my primary advisors, Dr. Lalita Udpa and Dr. Antonello Tamburrino, for their guidance and continuous support of my work, attention to the little details, and infinite patience. I’m extremely grateful to the rest of my thesis committee, Dr. Shanker Balasubramaniam and Dr. Mark Iwen, for numerous interesting discussions and conversations, both technical and not, and the constructive criticism. I wish to acknowledge my family: my parents, my sisters, and my uncle, for providing me with their endless support and a keeping me connected to home throughout these years. Special thanks to Andrea Harker, Oleksii Karpenko, Guy Parsey, and Leonid Poroshin. iii TABLE OF CONTENTS . . . . . . . . . . vi . 1.1 Motivation . 1.2 Literature review . . . . . LIST OF TABLES . . . LIST OF FIGURES . . . KEY TO ABBREVIATIONS . CHAPTER 1 . . . . . . INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1 Numerical methods for eddy current simulations . . . . . . . . . . . . . . . 4 1.2.1.1 Methods based on integral formulations . . . . . . . . . . . . . . 5 1.2.1.2 Methods based on differential formulations . . . . . . . . . . . . 7 1.2.1.3 Hybrid approaches . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.1.4 Dirichlet-To-Neumann operator . . . . . . . . . . . . . . . . . . 9 1.2.2 Automatic defect detection algorithms . . . . . . . . . . . . . . . . . . . . 1.2.2.1 NDE of multi-layer aircraft structures . . . . . . . . . . . . . . . 9 1.2.2.2 Data analysis for defect detection . . . . . . . . . . . . . . . . . 10 1.3 Contributions of this work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.4 The structure of the dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 . . CHAPTER 2 THE REDUCED VECTOR POTENTIAL FORMULATION WITH A . . . DIRICHLET-TO-NEUMANN BOUNDARY CONDITION . . . . . . . . . 15 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 Maxwell’s equations . . 17 2.2 Reduced Vector Potential formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Dirichlet-to-Neumann operator on the sphere . 21 2.4 Application of the DtN operator to the boundary term . . . . . . . . . . . . . . . 2.5 Finite Element Method . . 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Application of the boundary condition to the FEM system . . . . . . . . . . . . . . 24 2.7 Stiffness matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.8 Sparsification of the DTN operator . . . . . . . . . . . . . . . . . . . . . . . . . . 25 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.9 Concluding remarks . . . . . . . . . . Implementation details CHAPTER 3 VALIDATION AND APPLICATIONS . . . . . . . . . . . . . . . . . . . . 28 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 The convergence of the numerical DtN operator on the sphere boundary . . . . . . 30 3.3 Conducting sphere in a uniform magnetic field: performances evaluation . . . . . . 34 . 39 3.4 Conducting plate with a coil: numerical validation . . . . . . . . . . . . . . . . . 3.5 Coil inside a tube: experimental validation . . . . . . . . . . . . . . . . . . . . . . 42 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.6 Concluding remarks . . . . CHAPTER 4 GENERALIZED MULTIFREQUENCY FUSION ALGORITHM . . . . . . 47 . 47 4.1 Description of the approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv 4.2 Fusion algorithm . . . . 4.2.1 Secondary image generation . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Optimization algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Defect detection criteria . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.4 Details of numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 . 48 . 50 . 52 . 53 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 . . 4.3 Concluding remarks . 5.1.1 5.1.2 5.1.3 5.1.4 Numerical results . CHAPTER 5 MULTIFREQUENCY FUSION APPLICATIONS . . . . . . . . . . . . . . 54 5.1 Bobbin coil scan of a conducting tube near tube support . . . . . . . . . . . . . . . 54 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Simulation parametrs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Fusion algorithm setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 . 67 5.2.1 EC-GMR probe operating principle . . . . . . . . . . . . . . . . . . . . . 67 5.2.2 Experimental setup for GMR inspection . . . . . . . . . . . . . . . . . . . 69 5.2.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 5.2.4 Configuration of the fusion algorithm . . . . . . . . . . . . . . . . . . . . 73 5.2.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 . 5.2 EC-GMR scan of a riveted structure . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Concluding remarks . Pre-processing stage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 6 CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 82 . 85 APPENDICES . APPENDIX A USEFUL IDENTITIES . . . . . . . . . . . . . . . . . . . . . . . . . 86 APPENDIX B SPHERICAL COORDINATE SYSTEM . . . . . . . . . . . . . . . . 87 APPENDIX C VECTOR SPHERICAL HARMONICS . . . . . . . . . . . . . . . . . 88 APPENDIX D FEM STIFFNESS MATRIX COEFFICIENTS . . . . . . . . . . . . . 89 APPENDIX E FIELD GENERATED BY THE CURRENT RING . . . . . . . . . . . 90 CONDUCTING SPHERE IN A UNIFORM FIELD . . . . . . . . . . 91 APPENDIX F . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 . . BIBLIOGRAPHY . . . . . . . v LIST OF TABLES Table 3.1: Simulation parameters for problem 3.3 . . . . . . . . . . . . . . . . . . . . . . . 34 Table 3.2: Simulation parameters for problem 3.4 . . . . . . . . . . . . . . . . . . . . . . . 40 Table 3.3: Simulation parameters for problem 3.5 . . . . . . . . . . . . . . . . . . . . . . . 42 Table 3.4: Error in the simulated voltage and phase versus experimental data . . . . . . . . 46 Table 5.1: Simulation parameters for problem 5.1. . . . . . . . . . . . . . . . . . . . . . . 57 Table 5.2: Simulated defects for problem 5.1. . . . . . . . . . . . . . . . . . . . . . . . . . 57 Table 5.3: Parameters used in the five investigated test cases . . . . . . . . . . . . . . . . . 76 Table 5.4: Properties of the investigated test cases . . . . . . . . . . . . . . . . . . . . . . 76 Table 5.5: Combination of input parameters corresponding to the highest residue ratio C . . in each test case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 vi LIST OF FIGURES Figure 1.1: FE model showing  at: (a) defect-free aluminum fastener; (b) aluminum fastener site with subsurface notch; (c) defect-free steel fastener; (d) steel fastener site with subsurface notch. Geometries of fastener and defect are marked white. Black contours correspond to 30% of max() for aluminum fasteners, and pink contours correspond to 10% of max() for steel fastener. . 11 Figure 2.1: Schematic of the problem geometry . . . . . . . . . . . . . . . . . . . . . . . . 16 Figure 3.1: The relative error in DtN reconstruction for B field, truncation sphere radius  = 0.3mm vs number of terms in harmonic expansion, for different global quadrature orders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3.2: The relative error in DtN reconstruction for B field, truncation sphere radius  = 0.5mm vs number of terms in harmonic expansion, for different global quadrature orders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 . 32 Figure 3.3: The error in reconstruction of B caused by FEM discretization, sphere  = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Figure 3.4: Sample mesh of the problem 3.3. Darker region is a conducting sphere, brighter - simulated surrounding air. . . . . . . . . . . . . . . . . . . . . . . . 35 Figure 3.5: Imaginary component (dominating part) of the induced current density J in  cross-section at  = 0 on a sample mesh: (a) DtN solution, (b) Error vs analytic formula. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 . . Figure 3.6: (a) Error vs number of degrees of freedom, (b) Error vs number of operations required to compute matrix-vector product. . . . . . . . . . . . . . . . . . . . . 37 Figure 3.7: Effect of the sparsification: (a) Number of operations required to compute matrix-vector product versus number of degrees of freedom, (b) Ratio of operations required to compute the DtN-related part to the total number of operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 . . . . . . Figure 3.8: Sample mesh of the problem 3.4 . . . . . . . . . . . . . . . . . . . . . . . . . 40 Figure 3.9: Impedance of the coil during the linear scan: (a) Real part, (b) Imaginary part, (c) Lissajous plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Figure 3.10: (a) The geometry of the tube with a defect (b) Sample mesh of the defect, . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43  −  plane vii Figure 3.11: Geometry of the problem in  −  plane: 40%TW defect . . . . . . . . . . . . 44 Figure 3.12: Calibrated differential voltage measurement, real versus imaginary part, (a) 140kHz, (b) 280kHz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 Figure 4.1: Flowchart of the multifrequency fusion algorithm. . . . . . . . . . . . . . . . 47 Figure 5.1: (a) The simulated geometry for the 20%TW notch, (b) corresponding mesh. . . 56 Figure 5.2: Lissajous plot of the signal at five simulated frequencies, defect-free geometry . 59 Figure 5.3: Lissajous plot of the signal for five simulated geometries at 140 kHz. Black signal shows the dominant contribution due to support plate. . . . . . . . . . . 60 Figure 5.4: The signal magnitude of the contribution from the smaller defects with 20dB added noise: (a) defect #1, (b) defect #2. . . . . . . . . . . . . . . . . . . . . . 61 Figure 5.5: The fused images (0) and (2) generated with the "reference" algorithm applied to the noisy input with SNR of 20dB, defect #2. . . . . . . . . . . . . . 64 Figure 5.6: The fused images (0) and (2) generated with the "unlocked power" algorithm applied to the noisy input with SNR of 20dB, defect #2. . . . . . . . 64 Figure 5.7: The comparison between truncated fused (cid:48)() images and the base fused . . . . . . . . . . . . . . . Figure 5.8: The comparison between truncated fused images (cid:48)() and the base fused signal (0), SNR=20dB, "reference" algorithm. . 65 signal (0), SNR=20dB, "unlocked power" algorithm. . . . . . . . . . . . . . 65 Figure 5.9: Calculated criteria values  for the varying noise levels: (a) no added noise, (b) SNR 20dB, (c) SNR 17dB. . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Figure 5.10: Designs of excitation coils: (a) top coil (coil 1); (b) bottom coil (coil 2); (c) coil assembly with arrays of GMR sensors. Directions of currents are shown with arrows. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 . . . . . Figure 5.11: Differential GMR array probe with rotating current excitation placed on test sample with a fastener and a defect. . . . . . . . . . . . . . . . . . . . . . . . 68 Figure 5.12: Typical “baseline” field measured using the EC-GMR probe on aluminum plate. The reader is referred to [104] for detailed descriptions of probe configuration, operation and data acquisition. . . . . . . . . . . . . . . . . . . 69 Figure 5.13: Scanning of the sample with the EC-GMR differential array probe. . . . . . . . 70 viii Figure 5.14: EC-GMR scan data after background removal: (a) 100Hz; (b) 700Hz; (c) 1000Hz; (d) schematic of the bottom layer (fastener holes with EDM notches are filled with red color). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 Figure 5.15: Drawing of the top layer of the sample. . . . . . . . . . . . . . . . . . . . . . 72 Figure 5.16: Drawing of the bottom layer of the sample with 8 EDM notches. . . . . . . . . 72 Figure 5.17: Automatic detection of fastener centers (orange dots), and image segmentation using Hough Transform. The blue box corresponds to a segment of image with a defect-free (reference) fastener site. . . . . . . . . . . . . . . . . . . . . 74 Figure 5.18: Typical sizes of ROIs for computing field residues: small radius (1 = 5.4mm, yellow); medium radius (2 = 7.4mm, blue); and large radius (3 = 11.4mm, green). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Figure 5.19: The best fusion result from test case IV. . . . . . . . . . . . . . . . . . . . . . 77 Figure 5.20: Comparisons of normalized field residuals for test cases I-V. The numbers above the bar plots denote fastener site labels. The labels H along the x-axis correspond to defect-free fastener sites. Shown results correspond to ROI with radius 2 = 7.4mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Figure 5.21: Best residue ratios corresponding to fusion test cases I-V. . . . . . . . . . . . . 79 Figure 5.22: Fraction of input parameter combinations that results in successful defect detection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 ix KEY TO ABBREVIATIONS DtN Dirichlet-to-Neumann EC Eddy Current ECT Eddy Current Testing FEM Finite Element Method GMR Giant Magnetoresistive NDE Nondestructive Evaluation NDT Nondestructive Testing POD Probability Of Detection RVP Reduced Vector Potential TSP Tube Support x CHAPTER 1 INTRODUCTION 1.1 Motivation In the modern world Nondestructive Testing (NDT) plays an important role in science and engineering. The main objective of NDT is to assess the structural integrity of critical parts via detection and characterization of flaws, without destroying the usability of the part. NDT is widely used in design, production and maintenance stages in several industries, such as aerospace, nuclear, automotive, etc. Some of the most frequently used methods in NDT include ultrasound, radiographic, elecromagnetic, liquid penetrant testing, etc, and details of these methods can be found in [1]. This thesis is focused on different aspects of electromagnetic NDT methods. The underlying physics of electromagnetic NDT methods is governed by the set of Maxwell’s equations and contains several classes of methods based on different approximations. These methods range from DC (magnetic flux leakage) and magnetoquasistatic limit (eddy currents) to midrange (microwave testing) and high frequency (terahertz imaging, visual inspection). A detailed overview is found in the fundamental work [2]. Eddy Current Testing (ECT) is one of the well-established and widely used electromagnetic NDT techniques. It is used to inspect objects that are made of conductive materials. The physical foundation of the method is the following: an alternating current source generates a time-varying magnetic field, which induces eddy currents in a conducting object. Any imperfections in the test specimen, like cracks or material property variations, will affect the eddy current distribution, and cause the associated induced magnetic flux density to change, This change is detected by external measurements. For example, the induced magnetic flux density variations could be detected directly using GMR sensors [3] or by measuring the impedance of the reception coils [4]. Eddy current testing method is non-contact, which leads to lower inspection time and offers 1 reasonable sensitivity to detect surface and near surface defects. In addition, it does not require complicated equipment compared to other techniques and hence offers a low cost solution. The simplest setup could be build using just an impedance analyzer and a coil excited by an alternating current source. However, eddy currents are limited by the skin effect, whereby the magnitude of the induced current density and magnetic field decreases exponentially with the depth into the conductive media [5]. Further, the decay rate increases as the frequency goes up. It is important to note that in general, the contribution to the eddy current signal due to a defect may be orders of magnitude smaller relative to the background signal in a defect free sample. In addition, it is important to understand the influence of other factors that impact the field distribution, such as geometry and material properties of the inspected specimen, excitation frequency, probe lift-off variations and wobble, and so on. In view of these considerations, it is extremely valuable to have a computational model that can simulate the impact of these factors on the ECT probe measurements. Indeed, such simulation models can help to understand almost all aspects of the testing process, including optimization of sensors and system designs, analysis of dependencies between experimental parameters and resulting signals, generation of defect signatures, visualization, postprocessing and so on. It is impossible to build an analytical model that can take into account all these considerations. Numerical modeling is extremely important in all stages of the nondestructive inspection, from research and development and up to processing and analysis of acquired inspection data. Nowadays computational power is rapidly growing, which allows us to model the inspection of more complex geometries with sophisticated probe designs, conduct the Probability Of Detection (POD) studies, and so on. Consequently, there is always a strong demand for fast and accurate computational models in order to address the challenges imposed by the real-world problems. The present research focuses on two computational aspects where ECT simulations play a key role: the "forward problem" which predicts the signal measured by an ECT probe when it scans over a given specimen, and the "inverse problem" which entails distinguishing between the healthy 2 and damaged regions of the test sample based on information contained in EC NDT measurements. The primary goals of this work could be formulated in the following way: 1. To develop an accurate and computationally efficient numerical method to be used in low- frequency electromagnetics simulations for eddy current testing problems. 2. To develop a signal processing algorithm for automatic defect detection in one- and two- dimensional scan data obtained from the ECT inspection. 1.2 Literature review 1.2.1 Numerical methods for eddy current simulations In classical electrodynamics the behavior of electromagnetic fields is governed by Maxwell’s equations. However, analytic solutions to these equations could be derived in closed form only for a small set of canonical problems. Numerical solutions to Maxwell’s equations gained considerable interest in computational electromagnetics with the growth of computer industry in the middle of 20th century and since then evolved into a broad scientific discipline with numerous practical applications. In certain cases Maxwell’s equations could be simplified while retaining high accuracy. Broadly speaking, the different classes of models that have been developed are based on the nature of the excitation field, electrical size of the investigated object and material properties. Low frequency approximations investigate objects that are small compared to the wavelength. These models are characterized by assuming a quasistatic form of the equations, where time derivative of electric or magnetic fields are neglected, resulting in a set of equations for electric and magnetic fields that are decoupled. In the quasistatic models all aspects of wave propagation are neglected. Eddy current testing is generally applied to detect surface and subsurface defects in a conductive specimen [6]. Due to the skin effect, the frequency range suitable for inspection corresponds to a wavelength which is orders of magnitude larger than the size of the object. For example, an inspection of the conducting alloy tubes with outer diameter of 25mm and a wall thickness of 1mm 3 is typically carried out with excitation frequency in the range of 105Hz, which corresponds to a skin depth of 10−3m and wavelength of 103m. It follows that the inspection process could be effectively modeled by the means of a suitable low frequency model [7]: the magneto-quasistatic approximation [8]. Maxwell’s equations can be formulated in time or frequency domains, which are related by Fourier transform. Time domain formulations allow simulation of an arbitrary excitation pulse, but have to incorporate a time discretization scheme into the numerical model, whereas frequency domain formulations assume a sinusoidal time dependency and does not require a time discretization scheme. Both approaches are widely utilized in eddy current NDT modeling: for example, time domain methods are used to simulate pulsed eddy currents where the source is a step function, whereas frequency domain methods are applied for modeling eddy current inspections where the inspection is done at multiple frequencies. In order to solve Maxwell’s equations numerically, an equivalent "more computational-friendly" formulation of the problem must first be derived. There are three types of formulations applicable to low-frequency case, namely, integral formulations, differential formulations and mixed approaches. 1.2.1.1 Methods based on integral formulations When the investigated partial differential equation (PDE) has a known Green’s function, its exact solution could be written in an integral form [9]. This mathematical concept is the foundation of a family of integral equation methods, (for example, Method of Moments [10]) and it is widely used in computational electromagnetics. Generally, in integral formulation all material regions in the problem are considered as artificial sources. The unknowns are then associated with artificial source current densities (either electric or magnetic) that satisfy proper integral equations which compute the electric and magnetic fields in desired regions. A distinctive feature of integral methods is that the artificial source currents reside in the material regions which are usually bounded domains (except a few cases where infinite canonical shapes like halfspace, slab, cylinder, etc. are applicable). Sometimes it is possible to derive formulations 4 which contain only surface integrals (surface integral equation formulations, SIE), whereas a volume distribution of currents results in volume integral equation formulations (VIE) [11]. Due to their high accuracy, integral equation formulations are frequently utilized to simulate eddy current problems. Various methods have been proposed based on different field and potential formulations. A detailed review of the different methods can be found in literature as follows: general formulations [12, 13], formulations involving source coil [14, 15], and other recent publications [16–20]. There are several challenges related to numerical implementation of integral formulations. For instance, the related integral operators have singular kernels and need to be carefully implemented in order to retain accuracy. Linear systems obtained from VIE might be ill-conditioned [11], and special preconditioning might be required (see [21] for eddy current example). Further, the derived system of linear equations is dense, which imposes a limit on the number of unknowns and leads to increased computation time. However, in some cases an explicit evaluation of the full matrix is avoided. Several fast methods allow to reduce numerical complexity of these problems based on physical or mathematical approximations, such as AIM [22], FMM [23], ACA, and their applications to low-frequency case [24–27]. Compared to differential formulations, it is more difficult to adopt integral methods to model practical eddy current problems which involve nonlinear magnetic materials, though some research on this subject has been conducted in [28–31], and an efficient computational model presented in more recent publications [32–34]. 1.2.1.2 Methods based on differential formulations In this case the problem is formulated as a set of differential equations. The natural way of solving them would be to directly discretize the equations with a numerical scheme, which leads to methods like Finite Difference Method and Finite Volume Method [11]. The more generic and significantly more popular approach is to look for an approximate solution in a certain finite- dimensional subspace of the (infinite-dimensional) solution space – the Finite Element Method [35]. 5 The finite dimensional subspaces provide an approximation of the physical nature of the fields. Typical requirements in low-frequency cases include continuity and discontinuity of field components, curl or div conformance and so on. The set of appropriate subspaces is given and investigated in [36–38] . The first applications of FEM to eddy current problems were reported by Palanisamy and Lord [39–41]. A review of the typical low-frequency formulations in terms of potentials can be found in [42, 43] and its references. One of the widely used formulations is a reduced vector potential formulation [44–46], which allows the decoupling the source field from the "reaction" field, that is the field due to a specimen. This decoupling is advantageous in terms of accuracy, especially when the source field is rather large compared to the reaction field. The solution domain for differential formulations in unbounded domains is usually the entire space. However, to implement such a formulation in terms of a computer program, FEM must operate in a finite domain, where an artificial truncation boundary has to be set and an appropriate boundary condition has to be imposed. Some of the commonly used approaches include Dirichlet or Neumann conditions, Absorbing Boundary Conditions [47] or Infinite Elements [48]. More sophisticated methods are reported in [49] or [50]. Sometimes the system of linear equations arising from FEM is singular, as in the case of first order edge shape functions which curls are not linearly independent. This may or may not be a problem, depending on the particular solver and the gauging techniques applied to enforce uniqueness [51–53]. If multiply connected domains are considered, special treatment is required [54,55]. The related system of linear equations in differential formulations is sparse, which is advantageous from the computational point of view. Further, real-world problems often include complex geometries and materials that require very dense volumetric discretization in order to get the desired accuracy, resulting in large dimensional matrices to be inverted. This issue could be treated by using higher order elements [47] or optimal hp-FEM approcach reported in [56–58]. Another possible option is to utilize domain decomposition methods in order to reduce a large problem into 6 a set of smaller ones [59,60]. It should be noted that the classic FEM approaches (enforcement of tangential/normal continuity) in some practical scenarios could lead to numerical complications: for example, when interfacing two volumes discretized with elements of different shapes, or when using non-conforming meshes typically arising from adaptive mesh refinement techniques. Also, it is somewhat difficult to find a consistent combination of formulations or approximation spaces in different subdomains. One way of dealing with these problems would be to utilize discontinuous finite element subspaces: such approaches are called Discontinuous Galerkin methods [61, 62]. Another solution is offered by so-called meshless or mesh-free methods. A general overview of meshless methods can be found in [63], and some applications to electromagnetics and eddy current testing are presented in [64–68]. 1.2.1.3 Hybrid approaches Hybrid approaches combine different methods in one model in such a way so as to address individual drawbacks of each approach. A good overview of hybrid approaches is found in [11] and [69]. Hybrid techniques in eddy current modeling typically involve the classical FEM inside a domain with a boundary condition based on integral equation [70–72], allowing us to combine the versatility of FEM in modeling complex structures and materials with accuracy of integral methods in unbounded domains. 1.2.1.4 Dirichlet-To-Neumann operator In mathematics, Poincare-Steklov operator establishes a relationship between different boundary data of the solution of a PDE. A rigorous discussion of details and definitions of the approach are presented in [73]. One of its variants, the Dirichlet-to-Neumann operator, maps the Dirichlet data onto the Neumann data, on the boundary of the problem. Most typical applications of this operator in computational electromagnetics involve domain decomposition and utilize the map to build a consistent condition on the interface between domains. 7 There are multiple ways of utilizing the solution to the exterior problem to build a boundary condition. One example is a boundary element method, as in Bossavit and Verite in [74] and further works [75, 76]. In [74], where the H −  formulation is used, the unknowns are magnetic field intensity H inside the conductor (without boundary) and magnetic scalar potential  outside (including boundary). The domain is truncated and the boundary condition is written in terms of normal derivative   , which is approximated from  using boundary element method. The term "Dirichlet-to-Neumann operator" with respect to the boundary conditions was originally introduced by Givoli in [77]. In this work the general form of the DtN operator expressed through a Green’s function is derived for a scalar Laplace equation and then applied to obtain the analytic expression for the operator on a circle and a sphere. Also, the integration with FEM is analyzed. The technique for computing discrete versions of the DtN operator when analytic form of the operator is not known (for example, on domains with an arbitrary shape) was proposed in [78] and extended in [79]. Since then, in computational electromagnetics community the DtN operators have been applied almost exclusively to 2D scattering problems with a circular truncation boundary for which an analytic form of the operator is available [80]. The introduction of the DtN operator for a 3D vector EM problems with spherical boundary is reported in [81]. From a general perspective, DtN entails combining the benefits from two major classes of numerical formulations: differential formulations and integral formulations. The strength of differential formulations is that they give rise to sparse system of equations, whereas strength of integral formulations is that they require discretization of the material regions only. Hence the DtN method results in reduction of the size of the computational domain. Indeed, the DtN boundary condition is exact and, therefore, the boundary of the computational domain can be placed in close proximity of the materials. On the other hand the numerical model is based on a differential formulation, thus yielding an almost sparse system where only the submatrix related to the discretization of the DtN operator is fully populated. This matrix is relatively small because it is related to the unknowns on the boundary of the computational domain only, and it is fully 8 populated because the DtN boundary condition is non-local. The low frequency applications of the DtN operator have not yet been fully investigated, and the present work aims to fill this gap. 1.2.2 Automatic defect detection algorithms A second critical aspect of eddy current NDT is the solution to inverse problem that entails using the information is the probe signal to detect and characterize defects in the test specimen. This problem varies based on a specific application and needs of the industry. Two specific applications addressed in this thesis are described in this section. 1.2.2.1 NDE of multi-layer aircraft structures Fatigue cracks that develop and grow from fastener holes in multi-layer airframe structures pose a serious problem to the structural integrity of aircrafts [82–84]. There is a growing need for robust and rapid NDE techniques capable of detecting defects in bottom layers with the fasteners in place [85]. Traditional methods for locating cracks around fasteners and hidden corrosion include: • ultrasonic testing (UT) [86–88] • X-ray [89] • optical imaging [90] • lock-in thermography [91] • eddy current testing (ECT) [92–94] Low frequency ECT is an indispensable NDE tool for the aviation industry as it provides necessary penetration depth and sensitivity for detection of flaws deep under the surface without using couplants. In conventional ECT, magnetic field around fastener site is sensed by absolute or differential coils, and 2D raster scanning is performed for diagnostic imaging [95]. Recently, 9 giant magnetoresistive (GMR) and tunnel magnetoresistive (TMR) sensors have shown promise to replace field pick-up coils due to their high sensitivity, small size and wide frequency range [96–98]. In order to speed-up the inspection and increase coverage, total magnetic field due to excitation source and sample is measured by an array of GMR or TMR sensors [99, 100]. For instance, eddy current probe for the MAUS NDE system developed by Boeing consists of a planar coil that induces surface current in the test sample [101]. Resulting magnetic field around the fastener site is measured using the GMR array with 16 elements on a line of symmetry of the coil. Such probe is mostly sensitive to subsurface cracks perpendicular to induced current. In this thesis we consider a newly developed probe with rotating eddy current excitation and differential GMR sensor array field pick up (EC-GMR) [102–106]. The sensor array has the capability for scanning along two rows of fasteners simultaneously. Inducing rotating surface current helps to achieve nearly uniform sensitivity to flaws in all directions [104]. 1.2.2.2 Data analysis for defect detection Despite recent advances in eddy current NDT, reliable detection of subsurface defects around fasteners is not trivial. Defect detection from measured magnetic fields is challenging particularly when • scanning over steel fastener sites, or • the probe signal is corrupted by contributions from vertical edges in top or bottom aluminum layers. A simplified finite element model for EC inspection was used to study the challenges in defect detection from data obtained in the inspection of complicated for multilayer riveted geometries. Details of the modeling process, multilayered geometry, excitation coil configuration, experimental implementation and typical signal formation are discussed in the work by Ye from our research group [105]. The data for developing the novel signal processing algorithms was derived from the work [107]. 10 Figure 1.1: FE model showing  at: (a) defect-free aluminum fastener; (b) aluminum fastener site with subsurface notch; (c) defect-free steel fastener; (d) steel fastener site with subsurface notch. Geometries of fastener and defect are marked white. Black contours correspond to 30% of max() for aluminum fasteners, and pink contours correspond to 10% of max() for steel fastener. Figure 1.1a shows a typical image for the aluminum fastener. Field pattern  is circular in defect-free case, and its circularity breaks down with the presence of subsurface flaw (see Figure 1.1b). Hence, contour deformation coefficient was introduced for damage quantification in earlier works of the authors [104]. However, subsequent studies revealed that such damage detection criterion was difficult to apply in case of ferrous fasteners. This is illustrated in Figure 1.1c and Figure 1.1d, when the fastener in the same model presented earlier is made of low carbon steel with relative permeability  = 100. Steel fastener acts as a flux concentrator producing strong  component, but disturbance of the measured field due to defect is small, which means that deep and short cracks may be missed. It should be noted that this problem is not unique to testing the multilayer structures with steel 11 rivets. For instance, similar issues arise in inspection of heat exchanger tubes in nuclear power plants, where coil arrays detect stress-corrosion cracks in a non-magnetic tube that is held by carbon steel support structures [108]. EC-GMR probe used in this work has a relatively large coverage that enables scanning over two rows of fasteners simultaneously. Distortions in the signal are observed even in the absence of defects: when fasteners are too close to each other or the fastener is close to a vertical edge of the structure. Due to the discussed factors, the visual inspection of the acquired data is not sufficient for the defect detection, which requires the development of automated processing algorithms for enhancing defect signals. Multiple signal processing techniques have been designed to solve this problem, including • principal component analysis (PCA) and independent component analysis (ICA) [109,110] • Hilbert-Huang transform (HHT) and support vector machines (SVM) [111] • wavelet analysis and artificial neural networks (ANN) [112] • model-based image processing [113,114] • and frequency fusion [115–117] Dual-frequency fusion reported in [118], provided the most promising and robust results, as it was able to successfully suppress steel fastener field along with amplification of residual field due to defect. One of the objectives of this work is to develop robust signal processing algorithms for detection of subsurface defects, which would extend and generalize the method presented in [118]. 1.3 Contributions of this work There are two main contributions presented in this work. The first contribution is a full analytical derivation of the hybrid DtN-FEM method designed for modeling electromagnetics in the low frequency limit. The physics is governed by Magneto- 12 Quasistatic approximation and is expressed in terms of a Reduced Vector Potential formulation. Equations are derived in a weak form using Galerkin procedure with appropriate test functions for electric scalar and magnetic vector potentials. The Dirichlet-to-Neumann operator is then utilized to build an exact boundary condition on the truncation boundary. The resulting equations are solved by the Finite Element Method. The stiffness matrix in this approach has a dense submatrix related to the DtN operator, which is a nonlocal one. An important contribution of this work is a derivation of an analytical factorization of this submatrix in case of a spherical boundary. This makes the computational overhead due to replacement of a local boundary condition (Dirichlet, Neumann, Robin) with an exact one (DtN) almost negligible, while retaining the gain in term of computational accuracy and, consequently, efficiency. The proposed formulation has been applied to investigate three typical ECT model problems. The first problem is the simulation of a conducting sphere placed in a uniform magnetic field. This problem exhibits an analytical solution and thus it is well-suited to validate and analyze the performance of the method. The second problem simulates the scan of an eddy current probe over a conducting plate with a subsurface volumetric defect. The predicted probe signal is compared to the signal generated by a validated benchmark code. In the third problem, an eddy current probe moves coaxially inside a tube with a volumetric defect in the tube wall, and the results are compared to experimental data. The second contribution is a method for defect detection from multi-frequency EC image data. This work introduces a generalization of the frequency fusion algorithm in two directions: 1. The search for fusion coefficients is formulated as a minimization problem that allows the use of an unlimited number of fusion parameters and significantly improves the overall performance of the algorithm. 2. A nonlinear "power" and "threshold" operations are introduced to add an ability to fuse the signal with the processed versions of itself, effectively extending the method to process data 13 acquired at a single frequency. These improvements lead to higher sensitivity to smaller defects in the noisy inspection data. The proposed algorithm is applied to 1D multi-frequency EC inspection of steam generator tubes with bobbin coils and 2D EC-GMR inspection of riveted plates. 1.4 The structure of the dissertation The dissertation contains two parts and is organized into six chapters. The first Chapter is an introduction, containing the motivation for the work and the literature review. Part I is related to the forward modeling problem and consists of two chapters. Chapter 2 begins with the necessary theoretical background for low frequency electromagnetic modeling and continues with the description and full derivation of the proposed numerical method for eddy current computations. Chapter 3 investigates the applications of the method to a few typical model problems, its performance analysis and validation against experimental data. Part II is devoted to the signal processing methods for defect detection. Chapter 4 introduces the improved algorithm for automatic defect detection applicable to single- and multifrequency ECT data. The applications of this algorithm to 1D and 2D experimental data are presented in Chapter 5. The general overview of the accomplished work, conclusions and the directions for the future work are discussed in Chapter 6. 14 CHAPTER 2 THE REDUCED VECTOR POTENTIAL FORMULATION WITH A DIRICHLET-TO-NEUMANN BOUNDARY CONDITION This chapter contains the overall derivation of the reduced vector potential formulation with the boundary conditions imposed by the means of the Dirichlet-to-Neumann operator defined onto a spherical domain. First, we summarize Maxwell’s equations in the magneto-quasistatic limit and a reduced vector potential formulation in both strong and weak forms. Then, the analytic form of the Dirichlet-to-Neumann operator is derived for a spherical boundary. This latter, when combined with a reduced vector potential formulation in weak form, is used to enforce an exact boundary condition on the spherical boundary, and thus we have a complete boundary-value problem inside the sphere. Eventually, using finite element method to solve this problem in a finite-dimensional subspace, we obtain the system of linear equations whose structure and properties are investigated. Specifically, it was realized that the small but dense submatrix related to the discrete counterpart of the DtN operator affects the efficiency in a non negligible way. Therefore, by taking advantage of the analytical form of the DtN operator, an original sparsification technique was introduced to reduce the matrix-by-vector product for the DtN related submatrix, thus making the overall efficiency extremely competitive. 2.1 Maxwell’s equations In this work, the computational region is  =  ∪  ∪ , where  and  are the conducting and magnetic regions, respectively, and  is the free space, see Figure 2.1. Domain  has a smooth boundary  = Ɖ with an outward-directed normal ˆn. Throughout this work, the symbols D, B, E, H, J , A,  stand for electric and magnetic flux densities, electric and magnetic fields, current density, magnetic vector potential and electric scalar potential. Tensors , ,  represent magnetic permeability, magnetic reluctivity and electrical conductivity, respectively [51]. Subscript "s" stands for "source", 2 = −1. Equations are written 15 Figure 2.1: Schematic of the problem geometry in frequency domain, where the   time dependency with angular frequency  > 0 is assumed. The problem of interest in the overall domain is represented by a diffusion equation in conducting regions  and a magnetostatic one in the complementary domain. Further, the displacement current is neglected (magneto-quasistatic limit), i.e. D ≈ 0. Under this assumption, Maxwell’s equations take the following form: ∇ × E = − B ∇ · B = 0 ∇ × H = J The constitutive relationships are H = B J = E in  0 elsewhere in  in  in  in  (2.1) (2.2) (2.3) (2.4) (2.5) 16 and the continuity equation takes the form ∇ · J = 0 J · ˆn = 0 in  on  (2.6) (2.7) Proper continuity conditions have to be enforced on the boundaries of all material regions, where tensors  and  are discontinuous (like the interface Ɖ between conducting and non-conducting regions): H × ˆn and B · ˆn are continuous on all material boundaries. (2.8) The source current J resides in a source region , which is assumed to be topologically separated from . The magnetic field H and vector potential A are produced by the source current J in free space. They are given by the Biot-Savart law: (cid:90) (cid:90)   H(r) = 1 4 A(r) = 0 4 1 |r(cid:48) − r| dr(cid:48) J(r(cid:48)) × ∇ J(r(cid:48)) |r(cid:48) − r| dr(cid:48) (2.9) (2.10) 2.2 Reduced Vector Potential formulation In general, the magnetic vector potential A and (modified) electric scalar potential  are defined as B = ∇ × A E = − A − ∇ (2.11) (2.12) In the reduced vector potential formulation, the total magnetic vector potential is decomposed as the sum of the known source term A and the "reaction" term A: A = A + A 17 (2.13) Substituting (2.13) in (2.1), we obtain ∇ × ∇ × (A + A) = (− A − A − ∇) + J (2.14) Here we have also decomposed the current density J appearing on the right hand side into a source current density and an induced current density (which is valid since  ∩  = ∅): J = E + J The second equation is obtained from (2.6): ∇ · (− A − A − ∇) = 0 (2.15) (2.16) After simplification, we get the following system (reduced vector potential formulation in strong form): ∇ × ∇ × A + A + ∇ = −∇ × ( − 1)H − A ∇ · ( A + ∇) = −∇ · ( A) (2.17) (2.18) The equations (2.17) and (2.18) are valid in the entire domain  and in the conductive domain , respectively. The reduced magnetic potential A should satisfy the continuity condition on the material boundaries, similar to (2.8): A × ˆn is continuous on all material boundaries. In order to derive a weak form of these equations, let’s define the inner product as (cid:90) (f , g) = f · g  (2.19) (2.20)  Despite the fact that investigated vector functions are complex-valued, the usual complex conjugate operation is not used in this definition. Indeed, the unknown fields are approximated by a linear combinations of real-valued functions with complex coefficients. 18 (cid:90) (cid:90)   (cid:90)  (cid:90) (cid:90) (cid:90)  (cid:90)  +  (cid:90) Ɖ (cid:90) (cid:90) The weak form is derived through a Galerkin procedure by using the test functions N(cid:48) and φ(cid:48) The appropriate subspaces for these functions are defined in Section 2.5. In this case, the equations (2.17, 2.18) can be rewritten as N(cid:48) · [∇ × ∇ × A + A + ∇]  =  (cid:90)  φ(cid:48) · [∇ · ( A + ∇)]  = N(cid:48) · [ − ∇ × ( − 1)H − A]  φ(cid:48) · [ − ∇ · ( A)]  (2.21) (2.22) After applying the divergence theorem and vector identities (A.1)-(A.5) (see Appendix A) and simplifying, equations (2.21) and (2.22) will form the following system (reduced vector potential formulation in weak form [43], [46]): ∇ × N(cid:48) · (∇ × A) d +  N(cid:48) · A d+ N(cid:48) · ∇ d + N(cid:48) · ( ˆn × ∇ × A) d = = − (∇ × N(cid:48)) · ( − 1)H d −  N(cid:48) · A d  ∇φ(cid:48) · A d −   ∇φ(cid:48) · ∇ d =  ∇φ(cid:48) · A d (cid:90) (2.23) (2.24) (cid:90)  −    This formulation is valid in the entire space, with the material properties determined by tensors , , . In numerical FEM implementation these tensors are assumed constant within each element. Without loss of generality, it can be assumed that the magnetic properties of the materials are vanishing in some neighborhood of the outer boundary Ɖ. The resulting quantities will be piecewise-smooth everywhere and discontinuous only on the subdomain boundaries. 2.3 Dirichlet-to-Neumann operator on the sphere In the proposed approach, the Dirichlet-to-Neumann operator ' acts on the reduced vector potential A. Specifically, it maps the tangential component of A on Ɖ onto the tangential 19 component of its curl on Ɖ, but for the exterior problem: ' : (A) → (∇ × A) on Ɖ, (2.25) where the subscript t selects the tangential component of its argument. When the boundary Ɖ has a canonical shape, an explicit form of ' is available. In this work we use a sphere of radius 0 as truncation boundary Ɖ. Using vector spherical harmonics (defined in Appendix C), the following multipole expansion for A in free space can be derived: =1 +∞(cid:88) +∞(cid:88) +∞(cid:88) =1 =1 =− (cid:88) (cid:88) (cid:88) =− =− A = − + + M  [A] M  [A] M  [A] 1 +2 0 1 +2 0 1 +1 0 Y+ U+ V (2.26) Here Y is a radial harmonic, U and V are directed in the tangential space onto the [A] are computed as: [A] and M  surface of the sphere [51]. The expansion coefficients M  [A] = (cid:104)A , U(cid:105) [A] = (cid:104)A , V(cid:105) (cid:104)U, U(cid:105) +2 (cid:104)V, V(cid:105) +1 M  M  0 0 Here (cid:104)f , g(cid:105) = (cid:82) (2.28) f · g∗ d. Using the definition of the curl in spherical coordinates, relevant vector identities and a definition of vector spherical harmonics (see Appendices B, C) , the following analytic expression for the curl of A is obtained as: M  ∇ × A = Y− Ɖ (2.27) =1 +∞(cid:88) +∞(cid:88) +∞(cid:88) =1 =1 (cid:88) (cid:88) (cid:88) =− =− =− − + [A] [A] ( + 1) +2 0  +2 0  +3 0 [A] U+ V M  M  20 From equations (2.26) and (2.29) operator ' can be constructed explicitly. It is important to note that in this case the tangential component of the curl depends only on the tangential component of the original field, in accordance with the uniqueness conditions. (2.29) 2.4 Application of the DtN operator to the boundary term Next, operator ' is used to generate an exact boundary condition. Consider the following term N(cid:48) · ( ˆn × ∇ × A) d (2.30) in (2.23): (cid:90) (cid:90) Ɖ +∞(cid:88) =1 (cid:88) =− +∞(cid:88) (cid:88) =1 (cid:90) +∞(cid:88) +∞(cid:88) =1 (cid:88) (cid:88) =− =−  '(cid:48) = =1 = − 1 0 1 0 =       2+4 0  2+3 0 Substituting equation (2.29) into (2.30) gives us the following: N(cid:48) · ( ˆn × ∇ × A) d = − M  [A] N(cid:48) · ( ˆn × U) d+ Ɖ + M  [A]  +3 0 =− N(cid:48) · ( ˆn × V) d (2.31) The equation (2.31) is valid if  is a scalar in a neighborhood of Ɖ. In our case, Ɖ is topologically separated from other domains, and  = 1 . If Ɖ is a sphere, ˆn = ˆr. By the means 0 of (2.27), (2.28) and formulas (C.6)-(C.11) (see Appendix C), after simplification the equation (2.31) takes the form N(cid:48) · ( ˆn × ∇ × A) d = , (2.32) (cid:90) Ɖ  +2 0 (cid:90) Ɖ (cid:20) '(cid:48) + '(cid:48)(cid:21) where Ɖ '(cid:48) = M  [A] M  [A]  M   M  (cid:104)V, V(cid:105) (cid:104)U, U(cid:105) [N(cid:48)] [N(cid:48)] (2.33) (2.34) (2.35) (2.36) Operators M  [ ] and M  [ ] are defined in equations (2.27) and (2.28) correspondingly. 21 2.5 Finite Element Method The differential formulation defined in (2.17) and (2.18) with an appropriate boundary condition can be written in operator form as Lu = f in  The corresponding weak problem is given by (cid:90) w · (Lu − f)  = 0 for ∀w ∈  (2.37) (2.38)  Here w is a test function that belongs to the subspace . Under certain conditions and appropriate choice of  (Lax-Milgram theorem, [38]), the solution of the weak problem also delivers the solution to the strong one. To translate (2.38) into a computer code, we need to introduce discrete spaces, approximation and discretization. The Galerkin Finite Element Method looks for an approximate solution to the weak problem. It involves the following steps: 1. Select a finite dimensional space:  =  < w1, . . . , w > ⊂ (, ) 2. Approximate the solution as ˜u = (cid:80) 3. Build the system of linear equations(cid:82) =1 w w · (L˜u − f)  = 0,  = 1, . . . ,  4. Solve the system with respect to unknown coefficients   In practical problems it is hard to build functions w that are defined in the entire domain . Typically, the computational domain is discretized into subdomains ("elements"), and  is built using functions that are local, i.e. nonzero only in a few subdomains. 22 From certain physical considerations [38] it follows that magnetic vector potential A and electric scalar potential  belong to the following functional spaces: A ∈ (, )  ∈ (, ) (, ) = {A ∈ 2() | ∇ × A ∈ 2()} (, ) = { ∈ 2() | ∇ ∈ 2()} (2.39) (2.40) (2.41) (2.42) The elementary shape functions have to be elements of the above functional spaces. It can be shown that if we have a polyhedral discretization of , it is possible to build a discrete subspace conformal to the corresponding continuous space by using the Whitney complex (see fundamental book [38]). Specifically, we have to approximate (, ) with nodal functions and (, ) with edge functions. So now, using appropriate basis functions N and φ, our unknown quantities are expressed as (cid:88) (cid:88) =1 =1 A =  =  N  φ , (2.43) (2.44) where  and  represent the number of corresponding basis functions,  and  are unknown If the test functions are chosen to be the same as basis functions (N(cid:48) = N and coefficients. φ(cid:48) = φ) , a square system of linear equations is obtained. The system of equations will be complete when an appropriate boundary condition is imposed on Ɖ, as will be discussed in the next section. 23 2.6 Application of the boundary condition to the FEM system After substituting (2.43) and (2.44) into (2.27), (2.28) we obtain M  [A] = M  [A] = M  [N] M  [N] (2.45) (2.46) Substitution of (2.45), (2.46) into (2.33) and (2.34) leads to the the following final representation of the boundary condition:(cid:90) N · ( ˆn × ∇ × A) d = (cid:20) (cid:88) =1  (cid:21) , '   + '   where Ɖ '   = +∞(cid:88) +∞(cid:88) =1 (cid:88) (cid:88) =− =−  '         = =1 = − 1 0 1 0 = 2+4 0  2+3 0 M  [N] M  [N]  M   M  [N] [N] (cid:104)V, V(cid:105) (cid:104)U, U(cid:105) (cid:88) (cid:88) =1 =1 24 (2.47) (2.48) (2.49) (2.50) (2.51) (2.52) 2.7 Stiffness matrix Next, equations (2.23), (2.24) are combined with (2.47) in order to obtain the complete system of linear equations. In matrix notation it will take the following form:  + ' +       =         Submatrices ,    are defined in Appendix D. Vectors  and  consists of unknown expansion coefficients  and . The elements of the DtN matrix ' are ,   ,   ,    ,   ,   given by this formula: ' = '   + '   (2.53) The linear system of equations has the following properties: 1. The equations are valid for linear anisotropic materials. 2. Matrix ' is symmetric. 3. If  is a symmetric tensor, matrix  is symmetric. 4. If  is a symmetric tensor, matrices ,  are symmetric, and  = (). 5. Matrices , , , ,  are sparse. 6. Matrix ' is square and dense. 7. The system is positive semidefinite. If a gauge is imposed, the system becomes positive definite. 8. Total number of unknowns is  +  When the first order edge and nodal elements are used and no gauge is imposed, the number of unknowns is equal to the number of edges in the mesh plus the number of nodes in a conducting region. The size of matrix ' equals to the number of edges on the spherical domain boundary. 2.8 Sparsification of the DTN operator For practical applications, the infinite series in (2.48) and (2.49) have to be truncated. If all vector spherical harmonics are renumerated with a single index that goes from 1 to  (total number of harmonics used), matrix ' in (2.53) may be rewritten as ' = ' + ', (2.54) 25 where (cid:88) (cid:88) =1 =1 '   = '   = M  [N] M  [N]   M  [N]   M  [N] (2.55) (2.56) If we denote the number of unknowns on the domain boundary as , equation (2.55) takes the form ' =  M M . . . 1 [N M  = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 1 [N1] M 1 [N2] M 2 [N1] 2 [N2] . . . M  . . . M  [N1] [N2] ] M 2 [N ] . . . M  [N (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) ] (2.57) (2.58) (2.59) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171)  =   1 0 . . . 0 0  2  0 . . . 0 . . . 0 0 0 0 . . .    (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) Similar expression for matrix ' can be derived as well. Here matrix  has  rows and  columns, matrix  is a diagonal matrix of order . In the proposed approach, the system of equations (2.52) is sparse, so the preferred solution method is an iterative scheme. Thus, the performance of the approach is determined by the matrix- vector product operation. If a matrix ' is explicitly stored as a dense square matrix, the cost of computing the matrix-vector product  = ' is (2 ). However, if we use the analytical  factorization given by (2.57), we could write ' = (( )), (2.60) so the total number of operations required to compute  = ' in this case would be ( ), and furthermore the assembly and storage of the full matrix ' is avoided. 26 To make a qualitative comparison between the computational costs without and with sparsification, it is worth noting that 1. we are interested in high accuracy inside Ɖ, and 2. the solution is rather smooth on Ɖ (its smoothness increases with 0) The accuracy of the FEM approximation is directly related to the order of shape functions used, as well as the size of the element. When the order of shape functions is fixed, the higher accuracy is achieved through the use of smaller elements, and as a consequence, the increased number of unknowns. The second observation is based on the fact that in the low-frequency limit the discrete approximations of the DtN operator given by (2.55) and (2.56) converge rapidly with , and the desired accuracy could be reached with relatively small values of  (see Section 3.2 for a numerical example). From these two observations we can conclude that in the practical cases of interest the relation  <<  holds, thus the impact of the factorization is relevant and reduces the computational costs from (2 ) to (). In addition, as follows from (2.60), it is sufficent to keep in memory  only the matrix  instead of full matrix ', which leads to similar reduction from (2 ) to ()  in the memory requirements. 2.9 Concluding remarks In this chapter the reduced vector potential formulation with a boundary condition based on a Dirichlet-to-Neumann operator was derived for domains with spherical boudaries. The formulation was coupled with the FEM to build a complete numerical method suitable for modeling low- frequency electromagnetics. A sparsification technique for efficient implementation of the method was introduced. Several applications of the proposed method for the typical eddy current inspection problems will be presented in the next chapter. 27 CHAPTER 3 VALIDATION AND APPLICATIONS This chapter contains an application of the proposed reduced vector potential formulation with a boundary condition based on a Dirichlet-to-Neumann operator to typical problems arising in the field of ECT. The first section contains a description of a numerical implementation of the formulation explained in Chapter 2. The second section shows the convergence behavior of the overall model coupled to the DtN operator for the benchmark problem consisting of a current ring in free space. In particular, the effect of the number of harmonics used to truncate the series expansion representing the DtN operator on a sphere is studied. The first benchmark problem simulates a conducting sphere placed in a uniform magnetic field. This problem can be solved analytically and thus is used to investigate the convergence of the method, as well as to demonstrate the analytic matrix factorization provided in Section 2.8. The quantity of interest is the induced current density J inside the conductor. Along with the analytical solution, the method is compared to a conventional approach based onto the Reduced Vector Potential formulation coupled with Dirichlet boundary condition. The second problem deals with simulation of an eddy current coil probe scanning a conducting plate with a subsurface volumetric defect. In this case, the analytical solution does not exist. The quantity of interest is the voltage induced in the probe coil in every scan position. The numerical results are compared to those obtained with a benchmark code CARIDDI [12,20,24,25]. The third problem simulates an eddy current probe moving coaxially inside a conducting tube with several volumetric defects within the tube wall, such as those encountered during the inspection of a steam generator tube [119]. The quantity of interest is again the voltage induced in the probe coil. The numerical results are compared to the experimental data. 28 3.1 Implementation details This section describes the details of software and codes used in the implementation of the numerical method. The main code was written in Python 3, whereas the computation intensive parts were implemented in C in compliance with C99 standard. Such an approach allowed us to combine the extensive capabilities and code simplicity of Python with the high speed of C to achieve a near optimal computational performance. The code was organized in a modular way, so it could be easily upgraded or adapted to fit different applications. The numerical implementation uses the unstructured tetrahedral meshes generated by Gmsh software [120]. Gmsh has a convenient built-in scripting language that allows to define the geometry and set the desired mesh properties through a semi-automatically generated script and thus integrate a mesh creation step for a given problem into a non-interactive process. The mesh processing and necessary data structures builder was implemented in Python using meshio library. For the numerical tests, first order edge and scalar functions were used. The corresponding finite element module along with the matrix assembler were implemented in C. The computation of spherical harmonics utilizes the GSL (GNU Scientific Library, [121]), and the quadratures for integration were taken from Python quadpy module. The obtained system of linear equations was solved using iterative solvers from PETSc library [122–125]. In order to utilize the analytic factorization technique (see section 2.8), the custom matrix-vector multiplication routine was implemented and integrated into the iterative solver. The unknowns of the linear system are the expansion coefficients of magnetic vector potential A and modified electric scalar potential . The post-processing module that computes the derived physical fields, induced voltages, etc, was written in C. Other post-processing steps were implemented in Python with a use of numpy and scipy libraries. All visualizations were performed with a use of the matplotlib library. 29 3.2 The convergence of the numerical DtN operator on the sphere boundary In this section the following problem is investigated: analytic computation of the DtN operator on the sphere boundary for the field produced by a current carrying loop. The loop is set in  −  plane that is centered at the origin and carries a current of 1A. The analytic formulas for the source potential A  and magnetic flux B  can be found in [126] (see Appendix E for details). The loop has a radius  and the truncation boundary is a sphere with radius . We first investigate the error that is introduced by a numerical implementation of (2.29) – the discrete version of Dirichlet-to-Neumann operator on the sphere boundary. The reconstruction error is calculated with respect to the number of terms in harmonics expansion and the quadrature order used for inner product computation, for different values of truncation boundary radius . The following algorithm was utilized: 1. Calculate the expansion coefficients M  [A  ] on the sphere from equations (2.27) and (2.28) correspondingly. For this calculation, the integral over the sphere is computed through a global Lebedev quadrature of a variable order. [A  ] and M  2. Reconstruct the tangential component of the field B as B = ∇ × A via (2.29). 3. Calculate relative error by the following formula (here (cid:107).(cid:107)2 is an 2 norm on the sphere):  = ∗ 100% (3.1) (cid:107)B  − B(cid:107)2 (cid:107)B  (cid:107)2 The numerical results for a current carrying loop of radius  = 0.2mm are presented in Figures 3.1 and 3.2. It is seen that the DtN reconstruction of the field B is converging to the analytical formula for B  , as expected: for the given source and truncation boundary, we could make the reconstruction error in (3.1) arbitrarily small by selecting an appropriate quadrature and the number of terms in the harmonics expansion. Based on these results, we can make two conclusions. The first conclusion: the reconstruction error rapidly decays with the number of harmonics in the series, even if the truncation boundary is close to the source. This is a key feature because analytical solutions expressed in terms of a series 30 Figure 3.1: The relative error in DtN reconstruction for B field, truncation sphere radius  = 0.3mm vs number of terms in harmonic expansion, for different global quadrature orders. expansion are attractive when rapidly converging. The second conclusion: the quadrature order needs to be matched with the maximum order of harmonics in the expansion. We also investigate a study of the error due to the FEM discretization on the boundary. For calculation of this error, the following algorithm is used: 1. Create a triangulation on the sphere with  edges and define a vector shape function N on every edge  2. Project the vector potential A  to the subspace defined by edge shape functions: (3.2) (3.3) A =  =  N A  · d (cid:88) (cid:90) =1  31 Figure 3.2: The relative error in DtN reconstruction for B field, truncation sphere radius  = 0.5mm vs number of terms in harmonic expansion, for different global quadrature orders. 3. Reconstruct the tangential part of the magnetic flux B using the FEM basis: (cid:88) =1 B = ∇ × N (3.4) 4. Calculate the relative error in reconstruction using the formula (3.1) Using the first order edge shape functions, the numerical results are shown in Figure 3.3 for two configurations of the source loop and the truncation sphere. The numerical error that comes from FEM discretization decays linearly with the number of edges on the boundary, as predicted by theory. It is important to note that the error introduced by DtN operator reconstruction could be made an order of magnitude smaller than the FEM error with a small order of harmonics, even when the truncation boundary is very close to the current loop. The two studies above allow us to find the optimal numerical parameters of the DtN operator for the given problem. They have to be calculated in the following order: 32 Figure 3.3: The error in reconstruction of B caused by FEM discretization, sphere  = 1 1. Set the truncation boundary and create a finite element mesh 2. Determine the FEM discretization error 3. Determine the sufficient number of terms in the harmonic expansion 4. Choose an optimal quadrature order for numerical integration With this procedure we can estimate the optimal parameters of the numerical DtN operator and we are ready to simulate the model problems. 33 3.3 Conducting sphere in a uniform magnetic field: performances evaluation The performances of the proposed methods were first investigated on the canonical problem of a conducting sphere in an uniform magnetic field. This problem allows computation of a closed-form analytical solution [51]. Parameter Sphere  Truncation sphere     (cid:107)(cid:107) Value Unit mm mm 3.774 × 107 S/m kHz mT 3 4 1 1 1 Table 3.1: Simulation parameters for problem 3.3 The simulation parameters are summarized in Table 3.1. The sphere is made of aluminum ( = 3.774×107 S/m,  = 1) and has a radius of 3mm. The sphere is placed in a uniform magnetic field of magnitude 1mT directed along the z axis wih an excitation frequency of  = 1kHz. The skin depth in this case is 2.6mm which is comparable to the radius of the sphere. The mesh for the problem was generated with Gmsh software [120] and is shown in Figure 3.4. It is an unstructured tetrahedral mesh of a uniform density in the entire computational domain. Truncation radius of the computational domain is  = 4mm. The results from the proposed approach are compared with those obtained from analytic formula as well as two other numerical methods: • Conventional model that uses the same formulation and uniform mesh with a Dirichlet boundary condition and computed in a larger domain (sphere with radius of 10mm) • CARIDDI code (see [12, 20, 24, 25]) which is based on an integral formulation and uses second order elements. Non-sparsified version of the code was used. The quantity of interest is the induced current density J inside the conductor, evaluated on a uniform grid inside the sphere, which is illustrated in the Figure 3.5. The error  for all three 34 Figure 3.4: Sample mesh of the problem 3.3. Darker region is a conducting sphere, brighter - simulated surrounding air. numerical methods was computed as the relative 2 error on the grid: (cid:107) −  (cid:107)2 (cid:107) (cid:107)2  = ∗ 100% (3.5) The error plots obtained using several meshes with different resolutions are given in Figure 3.6. It demonstrates that the proposed method significantly outperforms the approach with Dirichlet boundary condition. Moreover, the difference increases when the higher accuracy in terms of smaller  from 3.5 is required. Figure 3.6a represents the error versus number of degrees of freedom in the system (stiffness matrix order). Figure 3.6b shows the error versus number of operations required to perform the matrix-vector product with the stiffness matrix of the system. The number of operations is directly related to the number of nonzero entries in the sparse stiffness matrix, and, in case of the DtN approach, to the number of unknowns on the boundary. This number is directly related to the computational time when an iterative solver is used, as it is the usual case for sparse matrices. In this case (iterative solver), the matrix-by-vector product is the dominating operation. It is worth mentioning that the relative ∞ error in a solution at each evaluation point for the given mesh does not exceed 3%. 35 (a) (b) Figure 3.5: cross-section at  = 0 on a sample mesh: (a) DtN solution, (b) Error vs analytic formula. Imaginary component (dominating part) of the induced current density J in  36 (a) (b) Figure 3.6: (a) Error vs number of degrees of freedom, (b) Error vs number of operations required to compute matrix-vector product. 37 (a) (b) Figure 3.7: Effect of the sparsification: (a) Number of operations required to compute matrix- vector product versus number of degrees of freedom, (b) Ratio of operations required to compute the DtN-related part to the total number of operations 38 The effect of sparsification is demonstrated in Figure 3.7. Two methods (DtN with and without sparsification) are compared on identical meshes, and thus provide exactly the same error in solution. Figure 3.7a is in good match with predicted numerical complexity reduction (from (2 ) to ()). Figure 3.7b shows what part of the operations required to compute matrix-  vector product are taken by computation of the DtN operation. It is seen that with sparsification the time taken by DtN computations approaches zero with increasing mesh density (and, therefore, higher accuracy). Whereas, when sparsification is not applied, the corresponding time quickly reaches 99% of the total computational time. In other words, Figure 3.7 illustrates one of the key features of the investigated approach: with the proposed sparsification, we get an exact non-local boundary condition for the price of a minor increase in computational time. 3.4 Conducting plate with a coil: numerical validation This section describes the application of the proposed method to an application in nondestructive evaluation, namely, an inspection of the conducting plate using an eddy current coil probe [4, 6]. This example is not designed to evaluate the performance of the code, but rather to validate the approach by simulating a realistic structure, so the target quantity would be an impedance of the coil as it scans over the plate with defect. In this case, the solution cannot be written in a closed form. The simulation parameters are summarized in Table 3.2. The plate is made of aluminum ( = 3.774 × 107 S/m,  = 1) and has a thickness of 4mm. Excitation frequency is  = 1kHz. The skin depth in this case equals to 2.6mm, which is comparable to the thickness of the plate. The plate has a subsurface volumetric defect, described by a square cross-sectional area of a size 4mm and depth of 2mm which is located at the bottom of the plate. The coil carries a total current of 1A and has an outer diameter of 20mm and thickness of 5mm placed at a liftoff distance of 0.5mm. The mesh for the geometry is shown in the Figure 3.8. 39 Parameter Plate length Plate width Plate thickness Plate  Plate  Defect length Defect width Defect depth Coil inner radius Coil outer radius Coil height Coil liftoff Total coil current Frequency  Total scan positions Value 70 70 4 Unit mm mm mm 3.774 × 107 S/m mm mm mm mm mm mm mm A Hz 1 4 4 2 5 10 5 0.5 1 1000 31 Table 3.2: Simulation parameters for problem 3.4 Figure 3.8: Sample mesh of the problem 3.4 40 (a) (b) (c) Figure 3.9: Impedance of the coil during the linear scan: (a) Real part, (b) Imaginary part, (c) Lissajous plot 41 The probe scans across the defect, and its complex impedance is computed at each scan position. The simulation results are presented in Figure 3.9. Figures 3.9a and 3.9b present the changes in real and imaginary parts of the coil impedance as the probe scans across the defect. The complex impedance trajectory is plotted in Figure 3.9c. For this problem, the reference results were obtained by CARIDDI code [20]. It is seen that there is a good match in the solution between the two approaches: the maximum relative difference in the real part of the signal (dominating component) is within 5%, and within 10% for the imaginary part. 3.5 Coil inside a tube: experimental validation In this problem we consider a coil moving coaxially inside the conductive tube with a volumetric defect inside the wall. This problem does not admit the closed form solution. The purpose of this example is to validate the proposed approach against experimental data, thus the quantity of interest is the impedance of the coil as it scans along the tube. Parameter Tube outer diameter Tube wall thickness Tube  Tube  Defect diameters Defect depths Coil outer diameter Coil inner diameter Coil thickness Coil liftoff Total coil current Frequency  15.2 12.2 1.3 0.8 1 35, 140, 280, 550 Value 19.0 1.1 8.71 × 105 1 4.8, 4.8, 2.8, 1.3 20, 40, 60, 100 % wall thickness mm Unit mm mm S/m mm mm mm mm A kHz Table 3.3: Simulation parameters for problem 3.5 The simulation parameters are summarized in Table 3.3. The tube is oriented along the z axis. The 3D test geometry of the problem is shown in Figure 3.10a and the x-z cross-section of the tube is displayed in Figure 3.11. The tube has an outer diameter of 19.0mm, the wall thickness of 1.1mm and is made of non-magnetic material with  = 1 characterized by its conductivity of 42 (a) (b) Figure 3.10: (a) The geometry of the tube with a defect (b) Sample mesh of the defect,  −  plane  = 8.71 × 105 S/m (equivalent resistivity value  = 114.8 -cm). Four volumetric defects were simulated: flat bottom holes (FBH) with depths of 20%, 40%, 60% and 100% of the tube wall. The defect diameters were 4.8mm, 4.8mm, 2.8mm, 1.3mm respectively. The defects were positioned on the outer side of the tube. The source field was generated by two geometrically identical coils that were carrying the total current of 1A each and a phase difference of 180◦. This represents the differential measurement [2]: in the section of the tube away from the defect the induced voltages in the coils have equal magnitude (so that the resulting impedance of the two coil system is zero), whereas in proximity of the defect the induced voltages in the two coils are different, which results in a non-zero signal that could be used to detect the presence of the defect. The simulated coil parameters were the following: outer diameter of 15.2mm and cross-section of 1.5mm by 1.3mm. The spacing between the coils is 1.5mm. The coil lift-off in this case equals to 0.8mm. The simulation of the model result in a one-dimensional array of complex voltages evaluated in every scan position. The length of the simulated tube segment was 76mm. The computational domain was a sphere 43 Figure 3.11: Geometry of the problem in  −  plane: 40%TW defect of radius 40mm. It was discretized into unstructured non-uniform tetrahedral mesh generated by Gmsh software [120]. Mesh density was set to be approximately 12 elements per wall in the defect area and 6 elements per wall thickness everywhere else. A sample mesh is shown in Figure 3.10b. In practice, it is impossible to measure the induced voltage in the coil directly: the inspection instrument outputs a response produced by a circuit containing the probe, often with additional hardware signal processing applied. Consequently, in order to correctly compare the simulation results to experimental data, both results need to be calibrated. In this work, the following calibration procedure was applied: the peak magnitude and phase of the signal generated by 100% through-wall defect at each frequency was set to 20V and 140◦ respectively. This was achieved by multiplying the signal by a constant complex-valued scaling factor, which was further applied to other three signals obtained at the same frequency. Thus, the relative magnitude and phase differences between all defects were preserved. The calibrated signals for two frequencies are presented on Figure 3.12. Two metrics were utilized to compare the signals after calibration: the magnitude and phase of the point with the peak magnitude. The results are summarized in the Table 3.4. 44 (a) (b) Figure 3.12: Calibrated differential voltage measurement, real versus imaginary part, (a) 140kHz, (b) 280kHz 45 Frequency, kHz Defect Type Phase difference, deg Peak voltage error, % 35 35 35 35 140 140 140 140 280 280 280 280 550 550 550 550 20% FBH 40% FBH 60% FBH 100% FBH 20% FBH 40% FBH 60% FBH 100% FBH 20% FBH 40% FBH 60% FBH 100% FBH 20% FBH 40% FBH 60% FBH 100% FBH 4 2 -1 0 -4 -3 5 0 -16 -9 -6 0 -5 -11 -9 0 11.00 10.00 3.00 0.00 8.00 7.00 14.00 0.00 14.00 3.00 3.00 0.00 3.00 9.00 20.00 0.00 Table 3.4: Error in the simulated voltage and phase versus experimental data 3.6 Concluding remarks In this chapter, the proposed reduced vector potential formulation with a boundary condition based on a Dirichlet-to-Neumann operator was applied to the typical model problems in the ECT field. The numerical implementation details were discussed. The performance of the method was investigated using two problems that admit analytical solution, and overall performance of the implementation matches the theoretical prediction given in Section 2.8. The method was validated against other simulation software and experimental data. As a conclusion, the proposed formulation could be readily applied for simulation of the low-frequency electromagnetic problems. 46 CHAPTER 4 GENERALIZED MULTIFREQUENCY FUSION ALGORITHM Eddy current inspection consists of three general steps: input data acquisition, data processing and defect characterization. Typical data volume is large, and it is sometimes difficult to distinguish the healthy regions and regions with defects without the use of special methods. In this chapter the generalized multifrequency fusion algorithm for automatic defect detection is introduced. The algorithm takes the acquired data and manually defined "healthy" regions as an input, and finds the potential regions with flaws. The first section of this chapter contains a high-level overview of the method. Then, all necessary mathematical operations for fusing data from different frequencies are defined, and the method is formulated in terms of a minimization problem. Applications of the proposed method are given in Chapter 5. 4.1 Description of the approach The acquisition of test data is discussed in detail in Section 5.2. For the purposes of this chapter it is sufficient to assume the input data in the form of one or more arrays of complex numbers. We also assume the source to have a   time dependency, and thus each array represents the input data acquired with a single excitation frequency. The arrays could be one- or two-dimensional. An example of the 2D input data is shown on Figure 5.14. Figure 4.1: Flowchart of the multifrequency fusion algorithm. 47 The overall schematics of the proposed approach that uses signals acquired with an arbitrary number of frequencies is shown in Figure 4.1. The first row on the diagram represents the input data acquisition and the applied preprocessing steps – backround removal and automatic image segmentation – which are discussed in Section 5.2.3. The second row shows the steps of the fusion algorithm that are discussed in Section 4.2. The algorithm first segments the data into individual fastener images and returns a processed image of the test sample with field residues due to defects. After that, using a simple threshold criterion, the algorithm classifies each fastener site as “defect-free” or “with defect”. The proposed method has the following key features: • generalized linear and nonlinear operations to suppress responses from defect- free regions • robust numerical method for determining optimal fusion coefficients via minimization of an appropriate cost function. 4.2 Fusion algorithm 4.2.1 Secondary image generation In general, the input data is considered as a set of one- or two-dimensional arrays of complex data. Let 1 . . .  be the data measured at frequencies 1 . . . , complex matrices with  rows and  columns. For convenience, these matrices are converted into real matrices with 2 rows and  columns, where the first (second) row is the real (imaginary) part of the corresponding complex matrix stored in row-major order. The data in a matrix form is an input to the fusion algorithm. The generalized transform W is applied to calculate linear combinations of real and imaginary parts of the input data. It could be represented as a product of the input matrix  with a 2 × 2 real matrix : W[] =   = (4.1) 111 + 122 211 + 222  48 Here 1 and 2 are the first and second rows of matrix , i.e. real and imaginary parts of the original image. In the general form, operator W has four degrees of freedom: the real values of 11, 12, 21, 22. Applications of a reduced form of this operator with a smaller number of degrees freedom (one and three) have been studied before in [118]. and form a foundation for the presented method. The variant with one degree of freedom is obtained when  is a rotation matrix defined by a single parameter, rotation angle. The variant with three degrees of freedom adds two real scaling factors 1, 2 and is used as a reference algorithm in this work. It is represented by the following formula: W[] =(cid:169)(cid:173)(cid:173)(cid:171)1 0 0 2 (cid:170)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:171) cos  sin  − sin  cos  1 2  (cid:170)(cid:174)(cid:174)(cid:172) (4.2) The data fusion defined by the linear operator (4.1) and (4.2) is enhanced using a nonlinear operation. In this work we investigate an effect of one possible nonlinear operation: raising the input data to a power. The power operator acting on matrix  is defined by the following formula: (P,[]) = ( − )| − |  (4.3) where the power  and threshold  are two parameters to be estimated. The power value is assumed to be positive:  > 0. Applying an operator P effectively enhances or compresses the amplitudes of real and imaginary components of  in a nonlinear fashion. Such a nonlinear operation is intended to address the effect of field diffusion at low and high inspection frequencies. If  > 1, strong signals are amplified, and weak ones are suppressed. Conversely, 0 <  < 1 reduces the differences between high- and low-amplitude regions of the EC-GMR scan. The threshold parameter  allows to adjust the power operation with respect to the input data magnitude. Similar effect could be achieved by the input data normalization. However, if we use (4.3), parameter  could be used as an extra degree of freedom in the optimization method described below, whereas normalization is a preprocessing technique and is applied before the optimization method. 49 Using the operators W and P defined above, the input data  is transformed into a "secondary image"  in the following way:  = W[P[]] (4.4) Among the input data 1 . . .  we could select one image and call it a reference 0. The image 0 is selected to represent the data acquired at one of the low inspection frequencies, for which the skin depth is large enough to detect all bottom layer defects. We can now define a fused image  as a linear combination of the reference image 0 (without any operations applied) and several transformed versions of input images , potentially including the transformed version of 0 itself: (cid:88) =1 (cid:88) =1  = 0 +   = 0 + W [P [  ]] (4.5) Here  denotes the number of secondary images fused to the reference image, and   are obtained , 1 ≤   ≤ . Note that for different values of  from (4.4) by transforming an input image   operators W  and P  do not have to be the same, and we could apply them to different input images   . The matrix  is a 2 ×  real matrix. (Ò) The fused image  is further segmented into subregions   ()  , which correspond to reference defect-free "healthy" sites and remaining test sites of interest potentially with defects, respectively, and  = 1 . . .   ,  = 1 . . .  −   . Here  is the total number of regions of interest within the input data, and   is the number of "healty" regions, which is typically 1 or 2. and  4.2.2 Optimization algorithm At first, the number of secondary images  is defined, and the input images   are selected in the fusion formula (4.5). The unknown parameters of operators W  and P  for equation (4.5) are calculated through minimizing a cost function Ǝ defined as a sum of the 2 norms of the healthy 50 (Ò) regions   : Ǝ[,  ,   ,  ] = (cid:88) =1 (Ò) (cid:107)   (cid:107)2 Using equations (4.5) and (4.6), the optimization problem is formulated as: min ,  ,  ,  Ǝ[,  ,   ,  ] s.t.   > 0 (4.6) (4.7) An optimal solution to problem (4.7) returns unknown parameters of operators W  and P  that effectively suppress the signals in the healthy regions (Ò). Once these coefficients are estimated (using only the data from defect-free fastener sites (Ò)), they are applied to build the fused subimages () corresponding to subdomains of interest in the test sample, and an entire fused image  from equation (4.5). In a final step, the defect detection criteria are computed in both subdomains (Ò) and () (see Section 4.2.3), and the defects are found. In this work three variants of the optimization method are investigated: 1. Reference: here power and threshold are not applied ( = 1,  = 0), and matrix  is limited to scaling and rotation. In this case optimization space has the dimension of 3( − 1), where  is the number of input frequencies. This version is not applicable if the input data was collected using a single frequency. 2. Locked powers: here the power value  is not an optimization parameter, but rather belongs to a set of predefined powers  = 1 . . . . The threshold is not applied. Optimization space has the dimension of 4( − 1), where  is the number of powers used,  is the number of input frequencies. This version could be applied to a single-frequency input data, as long as more than one power value is used. 3. Unlocked powers: this is the most generic version of the algorithm, where both power  and threshold  are included in the optimization. Optimization space has the dimension of 6( − 1), where  is the number of powers used,  is the number of input frequencies. This version is also applicable to a single-frequency input data. 51 From the equation (4.4) it follows that the residual image  has a continuous dependency on the unknown coefficients: linear operation is continuous, and power operation is essentialy ()|| , which is continuous upon  and  for  > 0. Therefore, if the search domain is a compact set, i.e. closed and bounded, the cost function will have a global minimum in this domain. Thus, the optimization problem in (4.7) should be extended with additional constrains upon input parameters. Also, additional care should be taken in order to prohibit the occurence of the trival : for example, if one of   in (4.5) is the reference image 0, we restrict the corresponding power value   to be different from 1. In this work, the Sequential Least-Square Quadratic Programming (SLSQP, [127]) method was used to solve minimization problem (4.7). SLSQP belongs to a family of well-established numerical methods for multi-dimensional problems with constraints. The SLSQP performs a sequence of quadratic approximations of the cost function (4.6) with imposed constraints. At each step, the quasi-Newton method is applied to find an iterative solution. Numerical experiments show that when the algorithm is applied to the same set of secondary images with an arbitrary starting point, the iterative process always converges to the same solution. 4.2.3 Defect detection criteria In this work, two quantitative criteria were used for classifying regions as “with defect” or “without defect”: (1) = ((cid:107) (cid:48)(cid:107)2) (2) = ((cid:107)((cid:48))(cid:107)2) + ((cid:107)((cid:48))(cid:107)2) (4.8) (4.9) Here (cid:48) denotes a reference subregion (Ò) or a subregion of interest (), and (cid:107).(cid:107)2 stands for a 2 norm of the fused image. The rationale is that both criteria (1) and (2) should be significantly higher in the regions corresponding to fastener sites with defects, and should be close to zero in the healthy regions. 52 4.2.4 Details of numerical implementation As a test example, we considered one of the most computationally intensive case: The numerical implementation of the proposed method was completed in Python 3 and C99. The automatic fastener segmentation subroutine was implemented in Ò using the Hough transform and the OpenCV library [128, 129]. The fusion algorithm was implemented as a library in C for optimized numerical computations, linear operations were done via OpenBLAS library [130], and the multicore C implementation of the SLSQP solver was utilized for the minimization problem. the EC- GMR multifrequency scan data discussed in Section 5.2. The algorithm was applied to the images that were 32 × 250 pixels in size. For the fusion case with the largest optimization space, the average processing time was approximately 0.1s on an 8 core AMD Ryzen 1600 CPU, which is a significantly smaller amount of time than used for the corresponding data acquisition. This example demonstrates that the proposed implementation allows for nearly real-time processing of EC-GMR inspection results. 4.3 Concluding remarks In this chapter a generalized multifrequency fusion algorithm for automatic defect detection was introduced. Input data preprocessing techniques were discussed. The present work significantly extends the previously published algorithm, that could be formulated as as a limiting case of a proposed method. It is worth mentioning that the proposed algorithm could be used with minor changes even if the defect-free regions are unknown or not specified. In this case the fusion is applied to the entire input images instead of just the healthy regions, and the minimization problem in (4.7) is solved to minimize the global norm of the fused image . After the fusion coefficients are found, the fused image  is computed and the regions with flaws are expected to be found near the local peaks in the residual image. Two applications of the proposed fusion algorithm to 1D and 2D problems will be discussed in the next Chapter. 53 CHAPTER 5 MULTIFREQUENCY FUSION APPLICATIONS This chapter presents applications of the proposed multifrequency fusion algorithm to the defect detection in the ECT data obtained from two applications, namely, aerospace and nuclear industries. The first problem investigates the performance of the algorithm on simulated 1D data from a bobbin coil that moves coaxially inside a conductive non-magnetic tube. The tube wall has a circumferential notch with a conductive and magnetic support structure on the outside of the tube. A cross-sectional view of this geometry is shown on Figure 5.1a. The simulation was carried out in COMSOL Multiphysics 5.2 [131] using 2D-axisymmetric geometry and the AC/DC module. Two versions of the fusion algorithm – the "unlocked power" and "reference" versions defined in Chapter 4 – were coputed and compared using the simulated data with artificially added white noise. The second problem investigates the capabilities of the proposed approach for defect detection in a 2D raster scan of a multilayer riveted structure with subsurface defects. The test sample is a lap joint of two conducting panels with several manufactured EDM notches in the bottom layer. The data was collected with an EC-GMR probe at five frequencies. The "reference" and "locked power" versions of the fusion algorithm were applied to the probe measurements in order to investigate the performance of the method for different configurations of algorithm parameters and evaluate the robustness of the approach with respect to the input parameters. 5.1 Bobbin coil scan of a conducting tube near tube support 5.1.1 Problem formulation Steam generator is a heat exchanger in a pressurized-water reactor in the nuclear power plant. It consists of several thousand thin tubes with diameters up to 25mm and lengths of the order of 25m, which are held in place by carbon steel support structures. Due to the harsh operating conditions, they are prone to stress corrosion cracks in the tube wall. Hence, the structural health monitoring 54 of steam generators is crucial for the safe operation of the plant. Eddy current testing with a variety of probe coils is extensively used by the nuclear industry. Differential bobbin probe is one of the widely used ECT probes for the steam generator tube inspection. It consists of two bobbin coils that are positioned coaxially in the tube and carry the opposing currents of equal magnitude. Thus, in the regions without imperfections the induced voltages in the coils cancel each other, but in presence of the defects the differential signal increases, which is used to detect and characterize these defects. The data is collected with several values of excitation frequency. More information about the differential bobbin probe and the sample acquired data is presented in Section 3.5. However, the defects is not the only source that impacts the differential bobbin signal. For example, the tube support structures (TSP) generate a high amplitude signal even when there is no defect in the tube wall. The typical tube material is a corrosion-resistant nonmagnetic alloy, whereas the support structures are made out of a highly conductive and magnetic carbon steel. When the tube scans near the support structure, it generates a signal that is much larger in magnitude than the signal produced by the small defects (see Figures 5.3, 5.4 for an illustration). In this case, it is extremely difficult to separate the defect signal from the support signal without the use of specialized signal processing algorithms. The model problem described in this Section was selected to illustrate the application of the generalized multifrequency fusion algorithm described in Chapter 4 to the representative case of detection of the small defects in the tube near the tube support region, when the input data is corrupted by the tube support signal and noise. 5.1.2 Simulation parametrs The simplified model was designed to simulate an eddy current inspection of the tube with a crack near the tube support. The axisymmetric model geometry consists of four parts: the tube region, the tube support region, the coil and the defect, and is illustrated in the cross-sectional view in Figure 5.1a. The model parameters are summarized in Table 5.1. 55 (a) (b) Figure 5.1: (a) The simulated geometry for the 20%TW notch, (b) corresponding mesh. The tube has an outer diameter of 19mm and the wall thickness of 1.1mm and is made out of non-magnetic alloy with the conductivity value of 9.5 × 105S/m. The simulated tube length is 100mm. The tube support is separated from the tube by an air gap of 2.5mm, has the length of 12.7mm and thickness of 5mm and is made out of carbon steel characterized by the conductivity value of 6×106S/m and relative permeability of 100. Along with a reference signal with no defects, four defect signals were simulated: three 360◦ circumferential notches with the width of 0.1mm and depths of 5, 10, 20% of the tube wall thickness (TW), and one 360◦ groove with the width of 0.8mm and a depth of 5%TW. The z coordinate of the defect position was aligned with the top of the tube support. Thus, every 56 Parameter Tube inner diameter Tube wall thickness Tube length Support inner diameter Support thickness Support length Tube conductivity  Tube relative permeability  Support conductivity  Support relative permeability  Coil inner diameter Coil outer diameter Coil height Coil separation Coil liftoff Total coil current Frequencies  Scan region Number of scan positions Value 19 1.1 100 19.6 5.1 12.7 1 9.5 × 105 6 × 106 100 12.4 15.0 1.5 1.5 0.9 1 [-18, 18] 100 35, 140, 280, 550, 750 Unit mm mm mm mm mm mm S/m S/m mm mm mm mm mm A kHz mm Table 5.1: Simulation parameters for problem 5.1. Defect number #1 #2 #3 #4 Defect site z coordinate Healthy site z coordinate Type 360◦ circumferential notch 360◦ circumferential notch 360◦ circumferential notch 360◦ groove 6.3mm −6.3mm Width, mm Depth, %TW Position Outer Outer Outer Outer 0.1 0.1 0.1 0.8 5 10 15 5 Table 5.2: Simulated defects for problem 5.1. simulated signal will contain the contribution from two potential defect sites: the first site is the tube support edge with a defect in the tube wall that needs to be detected (located at  = 6.3mm), and the second side is the edge without defects in the tube (located at  = −6.3mm) that will be used as a "control group" of healthy sites. Defect properties are listed in the Table 5.2. The differential bobbin probe coil has an outer diameter of 15.0mm, an inner diameter of 12.4mm and height of 1.5mm, which gives a liftoff value of 0.9mm. The coils are separated by 1.5mm and carry the opposing currents of 1A each and are excited at five frequencies: 57 35, 140, 280, 550, 750kHz. The coil scan was emulated by computing the induced voltage on the coils at 100 scan points uniformly distributed in the interval [−18, 18]mm. With the specified parameters, the skin depth inside the tube region varies from 0.6mm at 750kHz to 2.7mm at 35kHz, and the skin depth inside the tube support region varies from 0.02mm at 750kHz to 0.11mm at 35kHz. The simulations were carried out using COMSOL Multiphysics 5.2 software package, with a use of AC/DC Magnetic Fields (mf) module. Since the problem geometry obeys a rotational symmetry about the axis of the tube, the 2D axisymmetric model was selected. The electrical permittivity  is neglected in all domains, so the imposed equations from the mf module correspond to the magneto- quasistatic approximation, which is solved in terms of a generalized magnetic vector potential as an unknown. The computational domain was discretized into an unstructured triangular mesh with the boundary layers near the inner boundary of the tube support domain (see Figure 5.1b). Second order elements were utilized. The scans for all frequencies were simulated with the same mesh, so an appropriate element size distribution was chosen to accomodate the skin depth values for a largest frequency of 750kHz. The obtained systems of linear equations were solved with an iterative solver BiCGStab with a relative tolerance of 10−3. The induced coil voltages for both coils were computed with a built-in interface and then exported into the  files to be processed by the multifrequency fusion code. A dataset of 25 signals were obtained at 5 frequencies for the 5 selected sample and defect geometries. The simulated signals at five frequencies for the defect-free domain with no additional noise are displayed on the Figure 5.2. The signals from five simulated geometries (defect-free and four defects) simulated at 140kHz are shown in Figure 5.3. These figures demonstrate that the defect contribution to the total induced voltage is significantly smaller than the signal due to the tube support. A defect contribution could be approximated as a difference between the signal obtained from a geometry with a defect and a signal obtained from a corresponding defect-free geometry. This contribution for defect #2 is illustrated on the Figure 5.4. Unfortunately, in the real world the 58 Figure 5.2: Lissajous plot of the signal at five simulated frequencies, defect-free geometry pure signal from the corresponding defect-free geometry is not available, hence the objective of the fusion algorithm is to provide similar level of "base signal removal" as shown in Figure 5.4. Lastly, every data acquisition system introduces a certain level of the input noise due to the various sources such as the inspection instrument electronics, the off-center coil position inside the tube, probe tilting or probe wobble during motion, etc. In order to analyze the robustness of the proposed method with respect to the input noise levels, the artificial uniform (white) noise was introduced. In addition to the simulation without noise, two values of the noise magnitude were chosen to approximately match the 10% and 20% of the expected defect signal magnitude for a differential measurement, which translates into a Signal to Noise Ratio (SNR) values of 20 and 17dB, correspondignly. 5.1.3 Fusion algorithm setup For this study, the portions of the tube with different defects were simulated independently, so the fusion algorithm was applied to the entire signal. We consider the simulated signal from the 59 Figure 5.3: Lissajous plot of the signal for five simulated geometries at 140 kHz. Black signal shows the dominant contribution due to support plate. defect-free geometry as "healthy" reference sites for computation of the cost function (4.6). Two versions of the fusion algorithm were used: the "reference" version (4.2) and the "unlocked power" version (4.3). The computations were implemented in the following order: 1. First, we select one of the 3 input noise levels and apply it to the defect-free data 2. Then, as explained in Section 4.2, we calculate the fusion coefficients from the noisy defect- free data by solving a minimization problem in (4.7). 3. The calculated coefficients are substituted into (4.5) to obtain the "base" image (0)– the fused image for a case without the defect. 4. After that, the same coefficients are applied to calculate the remaining 4 fused images (), for each defect case. For this study, the reference image 0 in equation (4.5) was chosen to represent the data acquired with the highest frequency of 750kHz, which could be explained as follows. The investigated 60 (a) (b) Figure 5.4: The signal magnitude of the contribution from the smaller defects with 20dB added noise: (a) defect #1, (b) defect #2. 61 geometry has a defect outside the tube wall, in a close proximity of the tube support. Moreover, the contribution from the tube support is a dominant part of the signal regardless of the frequency, but due to the skin effect the ratio between the defect and TSP contributions will be larger for the higher-frequency data. Consequently, we are interested in filtering out the TSP signal from the data representing the highest frequency available. In (4.5) we set  = 4, and the images 1, 2, 3, 4 represent the data acquired at 35, 140, 280, 550 kHz, respectively. The resulting dataset for two algorithm variations, three noise levels and five investigated sites contains a total of 30 fused images with dimensions of 2 × 100. In case of the "reference" algorithm, the optimization space has a dimension of 12, and the In parameters are searched inside the box  case of the "unlocked powers" algorithm, the optimization space has a dimension of 24, and ( ) 22 , ∈ [−10, 10], ( ) ∈ [0.1, 10], the parameters are searched inside the box  ( ) ∈ [− max | |, max | |],  = 1, 2, 3, 4. ∈ [−10, 10], ( ) ∈ [0, 2],  = 1, 2, 3, 4. ( ) 12 ,  ( ) 21 ,  ( ) 1 ,  ( ) 2 ( ) 11 ,  The elements of the "base" fused image (0) are expected to have the close to zero magnitudes in the entire domain, with respect to the simulation tolerance and an introduced artificial noise (0)  |. Thus, we detect the defects level. We denote the maximum magnitude of (0) as  = max  in the fused images () by finding regions with magnitudes larger than . It is convenient to use the truncated set of images (cid:48)(),  = 1, 2, 3, 4, that will contain only the data from the potential regions of interest: |  (cid:48)()  = | ()  | if | ()  |>=   otherwise (5.1) (5.2) Now the defect detection criteria is computed for the set (cid:48)() in the following way: 100(cid:88) =1  = |(cid:48)()  − | The defect is successfully detected if ||>= . The regions of interest for the fused image () are determined as the coordinates  of the points where the magnitude of the fused signal is 62 larger than the maximum of the "base" image : |(cid:48)()  |>  5.1.4 Numerical results (5.3) The fused images (0) and (2) for the case with noise level of 20dB are shown in Figure 5.5 for the "reference" algorithm and in Figure 5.6 for the "unlocked powers" version. The filled area marks the region that is used for the residual computation. As expected, the "base" images (0) are uniformely low in the entire domain for both algorithms, and the fused images (2) correctly indicate the site with a defect (which is at  = 0.25in) and suppress the signal at the healthy site ( = −0.25in). The comparison between the four truncated fused images (cid:48)(),  = 1, 2, 3, 4, and the "base" image (0) for the applied noise level of 20dB is demonstrated in Figures 5.7 and 5.8 for the "reference" and "unlocked powers" methods, correspondingly. In this case, the "reference" method detected 3 defects and missed the smallest one (defect #1), whereas the "unlocked power" method successfully detected all four defects. The criteria  computed by (5.2) and the corresponding detection threshold level  from (5.3) for the three noise levels are plotted in Figure 5.9. Without additive noise, both methods detect all four defects. With a noise level of 20dB the "unlocked powers" version detects all defects and the "reference" version misses the smallest one. With the noise level of 17db the "unlocked powers" algorithm detects three defects, and the "reference" version detects only the two largest defects. In all cases, the criteria values  for the defect-free sites is exactly zero, as expected. The fusion algorithm is expected to show the best results in the noise-free case, and the detection threshold level goes up with the noise level, which is also depicted in Figure 5.9. The "unlocked powers" algorithm shows a better performance compared to the "reference" version: it suppresses the defect-free sites better than "reference" when the noise is present, and the detection threshold grows slower. 63 Figure 5.5: The fused images (0) and (2) generated with the "reference" algorithm applied to the noisy input with SNR of 20dB, defect #2. Figure 5.6: The fused images (0) and (2) generated with the "unlocked power" algorithm applied to the noisy input with SNR of 20dB, defect #2. 64 Figure 5.7: The comparison between truncated fused (cid:48)() images and the base fused signal (0), SNR=20dB, "reference" algorithm. Figure 5.8: The comparison between truncated fused images (cid:48)() and the base fused signal (0), SNR=20dB, "unlocked power" algorithm. 65 (a) (b) (c) Figure 5.9: Calculated criteria values  for the varying noise levels: (a) no added noise, (b) SNR 20dB, (c) SNR 17dB. 66 5.2 EC-GMR scan of a riveted structure 5.2.1 EC-GMR probe operating principle A brief description of the EC-GMR probe used in this work is provided here to make the reader familiar with its design and operating principle [104]. The EC-GMR probe has two unidirectional coils oriented in orthogonal directions and carrying currents with 90 degree phase shift (see Figure 5.10). The rotating current produced by this coil configuration is represented as: J = ˆJ0 +  ˆJ0, (5.4) √−1. The two orthogonal coils are positioned where J0 is current amplitude in both coils, and  = 5mm above the surface of the test sample, they are 130mm long and share the same geometric center. The induced eddy currents rotate in x-y plane. Figure 5.10: Designs of excitation coils: (a) top coil (coil 1); (b) bottom coil (coil 2); (c) coil assembly with arrays of GMR sensors. Directions of currents are shown with arrows. EC-GMR probe has two 32-element GMR arrays that are positioned symmetrically 4mm above and below the excitation coils (see Figure 5.11). EC-GMR probe uses GF708 sensors from Sensitec GmbH. The sensor pitch is 1.6mm, and the total length of the array is 51.2mm. Line scans are performed perpendicular to the sensor array along the x direction. 67 Figure 5.11: Differential GMR array probe with rotating current excitation placed on test sample with a fastener and a defect. The output from the single channel of the differential probe can be expressed using the following equation:  = 2 − 1 = [B(2) − B(1)], (5.5) where 2 and 1 are phasor outputs of sensor from GMR array 2 and GMR array 1, respectively; 2 and 1 are the corresponding observation point vectors; B is the normal component of magnetic flux density, and  is the sensor gain constant. Using the superposition property, the magnetic field can be decomposed into three contributions (1) generated by the  directed Coil 1 in free space, the field () related to the induced eddy currents in (2) generated by  directed Coil 2, and field B at every point: the field vectors B vectors B the sample: B = B (1) + B (2) + B () (5.6) GMR sensor arrays are placed on the line of symmetry of Coil 1 and are perpendicular to the line of symmetry of Coil 2. Hence, in the ideal case when the probe is strictly parallel to the specimen, the specimen is homogeneous, and there is no edge signal, we can state that there is no background due to Coil 1 at all sensor locations – equation (5.7), and there is a strong background field B contribution is offset due to Coil 2 at sensors away from the line of symmetry. The B field B (1)  (2)  (2)  68 Figure 5.12: Typical “baseline” field measured using the EC-GMR probe on aluminum plate. The reader is referred to [104] for detailed descriptions of probe configuration, operation and data acquisition. by the differential measurement – equation (5.8). (1)  (2)  ()  (2) = B (2) − B (2) − B (1)  (2)  ()  B B B (1) = 0 (1) = 0 (1) (cid:54)= 0 (5.7) (5.8) (5.9) Equations (5.7), (5.8), (5.9) show that the output of the probe is proportional to the field B ()  due to the eddy currents in the aluminum plate. Distribution of B is non-uniform along the y-axis and is referred to as a “baseline” signal (see Figure 5.12). The baseline needs to be subtracted from measured data. ()  5.2.2 Experimental setup for GMR inspection The experimental setup for acquiring EC-GMR measurements is shown in Figure 5.13. The test sample was a two-layer aluminum panel made of high strength alloy Al-7075 with conductivity of 33.4% IACS. The top and bottom layers are fastened with fourteen fasteners HI-LOK HL-525-100 manufactured as per Mil-S-5000 from low alloy steel 4340 with conductivity of 6.95% IACS. Prior to the inspection, the fasteners were demagnetized by placing them inside of a demagnetizing coil. Eight of the fourteen fastener sites had second layer notches of different lengths machined through the thickness of the bottom layer. The notches were realized by using the electric discharge 69 Figure 5.13: Scanning of the sample with the EC-GMR differential array probe. machine (EDM) and have a width of 0.15mm. Locations and sizes of all defects are shown in Figure 5.16. The sample with 14 fasteners was placed in a high-precision XYZ gantry, and inspection data were acquired at excitation frequencies of 100Hz, 200Hz, 400Hz, 700Hz and 1000Hz. These excitation frequencies 1 = 100Hz to 5 = 1000Hz correspond to a range of skin depths of 1 = 11.51mm to 5 = 3.64mm, respectively. Hence, the high frequency field largely characterized the fastener head, while the low frequency data contained information regarding the bottom layer with defects. In all experiments, step size along x-direction was 1mm. EC-GMR data acquired at frequencies 100Hz, 700Hz and 1000Hz are shown in Figures 5.14a, 5.14b, 5.14c, respectively. Figure 5.14d illustrates fastener numbering and schematics of the bottom layer of the test sample. Note that fastener sites 10, 13 and 14 with flaws can be identified by visual inspection (see Figure 5.14a). For instance, fastener 10 is highlighted in the 100 Hz EC-GMR scan very well, because the eddy current flow between the hole and nearby bottom edge is blocked completely. Another important observation is that defects at fastener sites 11 and 12 are barely visible in Figure 5.14a, even though they are larger than the notch at fastener site 10. These artifacts can be contributed to the fact that the defect signal is dominated by the strong signal from the steel rivet, making defect detection very challenging. 70 Figure 5.14: EC-GMR scan data after background removal: (a) 100Hz; (b) 700Hz; (c) 1000Hz; (d) schematic of the bottom layer (fastener holes with EDM notches are filled with red color). 71 Figure 5.15: Drawing of the top layer of the sample. Figure 5.16: Drawing of the bottom layer of the sample with 8 EDM notches. 72 5.2.3 Pre-processing stage In this application, the pre-processing consists of two steps: background removal and automatic image segmentation to indicate the fastener sites. In this section, we define background as the low spatial frequency trends in the EC-GMR data due to (i) baseline signal from the EC-GMR probe as explained in Section 5.2.1, and (ii) signals from vertical edges in the top or bottom layers of the test sample. If vertical edges are parallel to the scan direction, their contributions can be removed by subtracting the edge signal from a calibration sample without fasteners. This simple approach was applied in the present work. More sophisticated methods that exploit differences in spatial frequencies between fastener and edge fields (e.g. polynomial fitting or robust sparse coding [114]), can be applied for background removal on more complicated samples. Automated image segmentation is critical for accurate computation of localized field residuals when classifying fastener sites as “defect-free” or “with defect”. At higher inspection frequencies, the magnetic fields are more focused around the fastener thereby preserving the circular shape of the fastener image. Further, since the skin depth is also small at high frequency, the image data is not affected by defects in the bottom layer. Hence, the Hough transform is applied to the high frequency signal for circle detection. This step is followed by Canny edge detection to get the outline of the circle. Once the circles are detected and their centers are identified, an appropriate square region of interest (ROI) is assigned to each fastener. 5.2.4 Configuration of the fusion algorithm Automatic image segmentation algorithm was applied to the 1000Hz EC-GMR data, and typical results are presented in Figure 5.17. All fastener sites were successfully identified as highlighted by the green circles. The centers of fasteners were detected with 1.5mm accuracy owing to slight field distortions. Figure 5.18 shows circular ROIs for computing field residues in the optimization problem. Three sizes of ROIs, namely, small (1 = 5.4mm), medium (2 = 7.4mm) and large (3 = 11.4mm) 73 Figure 5.17: Automatic detection of fastener centers (orange dots), and image segmentation using Hough Transform. The blue box corresponds to a segment of image with a defect-free (reference) fastener site. Figure 5.18: Typical sizes of ROIs for computing field residues: small radius (1 = 5.4mm, yellow); medium radius (2 = 7.4mm, blue); and large radius (3 = 11.4mm, green). were investigated to understand which regions of measured magnetic field had more information about bottom layer defects. As illustrated in Figure 5.18, the smallest ROIs fields measured strictly on top of fastener heads, and larger ROIs had more measurement points outside the fastener heads. The maximum radius was selected so as to avoid effects of field coupling between the neighboring fasteners. 74 For this study, the "locked powers" version of the fusion algorithm was used (see Section 4.2.2 for details). Powers defined in (4.3) were assigned fixed values in order to reduce dimensionality of the optimization space. An exhaustive search was applied to determine the best set of powers that would maximize detectability of all defects. The domain P with powers used in exhaustive search is specified in Table 5.3. It was observed that raising the image data to powers very close in value produced nearly identical secondary images. In order to avoid this effect, the powers used in the fusion algorithm were separated by at least 0.5. The defect-free fastener sites 1-6 were used as references for estimating the optimal weights in (4.5), which were then used to compute the residues at each segmented fastener image. As shown in Table 5.3, the total of 14 combinations of defect-free fasteners were used as reference regions in the fusion algorithm: six single fasteners - numbers [1, 2, 3, 4, 5, 6], and eight pairs of fasteners - [(1,5), (1,6), (2,5), (2,6), (3,5), (3,6), (4,5), (4,6)]. The results of the exhaustive search were grouped in five test cases that are summarized in Table 5.3. In each test case, exhaustive evaluation of results was performed by generating the fused images corresponding to all possible inputs – see Table 5.4. Total number of combinations that were obtained during exhaustive search for all five test cases was 236040. Test case I shows the fusion of single frequency data with its secondary images. The purpose of this case is to investigate the applicability of the algorithm to practical cases where the multifrequency data is not available or is hard to acquire. The fusion algorithm was applied to two secondary images obtained by raising 100Hz data in two power values (see Table 5.3). Test cases II and IV demonstrate the performance of the method when the power operation is disabled: no secondary images are generated, and the two or three frequency inputs are fused. The frequency combinations for these test cases were selected in such way that 100Hz or 200Hz data was always used as a reference image. It is worth noting that the test case II is a more generic version of the mixing algorithm presented in the previous work [118], which was used as a baseline for evaluating all other test cases. Test cases III and V represent the most generic algorithm with a power operation. For example, 75 List of frequency combinations for test case I List of frequency combinations for test case II List of frequency combinations for test case III List of frequency combinations for test case IV List of frequency combinations for test case V List of powers (P) for test case I List of powers (P) for test cases III and V List of reference fastener sets List of ROI radii, mm (100) (100, 400), (100, 700), (100, 1000), (200, 400), (200, 700), (200,1000) (100, 200), (100, 400), (100, 700), (100, 1000), (200, 400), (200,700), (200, 1000) (100, 400, 700), (100, 400, 1000), (100, 700, 1000), (200, 400, 700), (200, 400, 1000), (200, 700, 1000 ) (100, 200, 400), (100, 200, 700), (100, 200, 1000), (100, 400,700), (100, 400, 1000), (100, 700, 1000), (200, 400, 700), (200,400, 1000), (200, 700, 1000 ) 0.05, 0.07, 0.1, 0.12, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.8, 2.0, 3.0 0.02, 0.03, 0.04, 0.05, 0.07, 0.1, 0.12, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.8, 2.0, 3.0 1, 2, 3, 4, 5, 6, (1,5), (1,6), (2,5), (2,6), (3,5), (3,6), (4,5), (4,6) 5.4, 7.4, 11.4 Table 5.3: Parameters used in the five investigated test cases Test case I II III IV V input Number of frequency combinations 1 6 7 6 9 Number of power combinations Number ROI radii of 264 – 334 – 334 3 3 3 3 3 of Number reference fastener combinations 14 14 14 14 14 Total: Total number of input parameter combinations 11088 252 98196 252 126252 236040 Table 5.4: Properties of the investigated test cases 76 in test case III we first select two input images, then select two frequencies from the list P, then generate four secondary images by raising two input images in two selected powers, and then fuse them together, using one of the combinations for reference fasteners and ROI size. In test case V, three inputs and two powers are selected, and six secondary images are fused. 5.2.5 Numerical results Fusion of multi-frequency data was performed for all test cases as per Table 5.3, and the results were quantified using the criterion established in Section 4.2.3. A typical output of the fusion algorithm is presented in Figure 5.19. It corresponds to the test case IV, in which EC-GMR data acquired at 100Hz, 700Hz and 1000Hz were fused together without applying power operations. Figure 5.19: The best fusion result from test case IV. Figure 5.20 shows comparisons of field residuals for all five mixing cases. The bar plots for each case were normalized by the maximum residual among six defect-free fastener sites. The first six groups of bar plots labeled H on the x-axis corresponded to defect-free (healthy) fastener sites 1-6 shown in Figure 5.14d. The remaining eight groups of bar plots were sorted by the defect sizes from smallest to largest. Bar plots in Figure 5.20 also represented the best results for each mixing case. The criterion for these selections was the following ratio, that denotes the relative separation between the residuals computed within healthy domains and domains with defects (which could be computed since we 77 Figure 5.20: Comparisons of normalized field residuals for test cases I-V. The numbers above the bar plots denote fastener site labels. The labels H along the x-axis correspond to defect-free fastener sites. Shown results correspond to ROI with radius 2 = 7.4mm. know the exact locations of the defects):  =(cid:104) min  ()  () max   − 1(cid:105) × 100%, (5.10) where () is the set of residuals corresponding to all fastener sites 7-14 with defects, and () – the residuals corresponding to all defect-free fastener sites 1-6. For instance, in test case II, 252 sets of residuals are obtained based on the combinations of input parameters from Table 5.3. However, the best result corresponds to a single combination of input parameters (see Table 5.5) with the highest ratio  in (5.10), that provides the best visual separation between defective (7-14) and defect-free fastener (1-6) sites. Numerical computations revealed that in all five test cases the highest criterion values were obtained for the medium radius of circular ROIs (2 = 7.4 mm, see Figure 5.18 and Table 5.5). Hence, this radius was used to compute bar plots in Figure 5.20. Results shown in Figure 5.20 demonstrate that field residuals were higher for fastener sites with defects, which provides a numerical validation for the multi-frequency fusion approach. For 78 mixing cases with two and three frequencies, the residuals grew monotonically with increased defect sizes. However, the trend was non-linear, since notches at fastener sites 10, 13 and 14 were easily detectable. This could be explained by significant distortions of eddy currents in the bottom layer of the test sample. For instance, a long defect connecting fastener holes 13 and 14 would completely blocks currents orthogonal to it. Another useful observation is that all defects were also detected by mixing a low frequency data with its secondary versions. This is particularly important when multi-frequency data is not available or hard to acquire. However, in this case the difference between residuals at defective and defect-free fastener sites is smaller compared to that obtained with multiple frequencies. Figure 5.21 shows the ratios  from (5.10) computed for all fusion test cases, and Table 5.5 presents the corresponding input parameters. Quantitative results demonstrated that adding more frequencies and powers helped better suppress signals from defect-free fastener sites as well as to enhance defect indications. Figure 5.21: Best residue ratios corresponding to fusion test cases I-V. It is important to investigate the robustness of defect detection with respect to the choice of input parameters (frequencies, powers, sizes of residue evaluation region and indices of reference fasteners). For this purpose, the results of exhaustive search were utilized in the following way: for every input parameter combination (total of 236040, see Table 5.4) two sets of residuals, namely } (fastener sites 7-14 with Ò = { ()  } (defect-free or “healthy” fastener sites 1-6) and  = { ()  79 Test case residue Best ratio C, % I II III IV V 6.7 27.6 35.8 38.6 62.3 List powers of List frequencies, Hz 100 100, 400 100, 700 100, 400, 700 1 100, 400, 700 1, 3 0.05, 1.1, 2 1 0.8, 1.5 of ROI, mm 7.4 7.4 7.4 7.4 7.4 Reference fastener number (4,5) 1 4 3 4 Table 5.5: Combination of input parameters corresponding to the highest residue ratio C in each test case. ∈  that was defects), were computed using the formula in Section 4.2.3. Each residual  smaller than the largest residual in Ò was indicating a missed defect at the region with index , otherwise it was counted as successfully detected defect: ()  > max Ò: successfully detected defect ()  ()  ≤ max Ò: missed defect   (5.11) (5.12) Figure 5.22: Fraction of input parameter combinations that results in successful defect detection. The number of missed and successfully detected notches was computed for each input combination of test cases I-V using (5.12). After this, the input combinations were grouped by the number of 80 successfully detected defects, and obtained results are summarized in Figure 5.22. Since the numbers of input parameter combinations were different for test cases I-V, the vertical axis of Figure 5.22 was normalized by the corresponding values from the right column of Table 5.4. Note that for each test case, almost any selection of input parameters leads to detection of at least 3 defects – these are the same defects that could be detected visually from unprocessed images. More than 90% of the input combinations lead to successfull detection of 5 defects for all five test cases. From this we can conclude that the algorithm provides the output as accurate as visual inspection with almost any input parameters, and the that the algorithm delivers the acceptable results even if the non-optimal input parameters are selected. 5.3 Concluding remarks In this Chapter an application of the generalized multifrequency fusion algorithm for automatic defect detecton in 1D and 2D ECT data was demonstrated. The first model problem represents a detection of cracks within a 1D simulated data of the tube scan with a small defect located by tube support region. The fusion method effectively suppressed the tube support signal and was able to reveal the detect sites even for the smaller defects in presence of the higher levels of input noise. The second problem shows the application of the fusion method to the 2D experimental data, which is a 2D raster image of EC-GMR probe scan over the riveted structure with the bottom-layer defects. The performance of the several different configurations of the method were investigated, and multiple combinations of the input parameters were studied. An ability to apply the method to a single frequency data was shown, however with a lower sensitivity compared to the multifrequency fusion. In all investigated model cases the proposed algorithm was able to successfully detect all defects within the input data and demonstrated a significantly better performance compared to the classic method. 81 CHAPTER 6 CONCLUSIONS AND FUTURE WORK The contribution of this dissertation to the field of eddy current testing is two-fold. The first contribution is a novel low-frequency numerical formulation that solves the forward problems with high accuracy and efficiency. The second contribution is an improved multifrequency fusion algorithm that solves the inverse problem of defect detection and characterization from experimental data. The new numerical formulation combines a reduced vector potential formulation with an exact boundary condition imposed by the means of a Dirichlet-to-Neumann operator. The equations are derived in weak form and solved using a finite element method. The corresponding numerical model is sparse, apart from a "small" square block corresponding to the discrete couterpart of the Dirichlet-to-Neumann operator. This block is of the size of the number of unknowns on the boundary of the computational domain where the boundary condition is imposed exactly by means of the Dirichlet-to-Neumann operator, as aforementioned. This block is fully populated because of the non-local nature of the Dirichlet-to-Neumann operator. Despite its "small" size, this fully populated block affect the efficiency of the method, so a sparsification was developed to reduce the cost of the matrix-by-vector product. Moreover, the sparsification was carried out analitically for the Dirichlet-to-Neumann operator onto a sphere. This sparsification allows to embed the exact boundary condition based onto the Dirichlet-to-Neumann operator with a negligible overhead, if compared to simpler (local) boundary conditions. The numerical performance of the method has been investigated for a model problem where an analytical solution was available: a conducting sphere in a uniform magnetic field. Other than this analytical benchmark case, the new numerical model was compared to other simulation codes as well. Specifically, this latter comparison was carried out for a practical eddy current testing problem consisting in simulating the C-scan of a coil over a conducting plate with a subsurface volumetric defect. Another validation was carried out against experimental data acquired with a 82 bobbin coil scanning a conducting tube. The proposed formulation could be readily adopted to simulate more complicated ECT problems, such as domains with nonlinear magnetic materials or eddy current probes with complex designs. However, there are multiple directions of the future improvements of this method. The spherical coordinates used in this work is a good starting point for two reasons: it is possible to derive a simple expression for the analytic DTN map on the boundary and it is easy to investigate the performances of the method. On the other hand, most geometries of interest have either one (tubular) or two (planar) dominating dimensions, where an ellipsoidal boundary might be more computationally efficient. Helmholtz equation allows separation of variables in both oblate and prolate spheroidal coordinates, so in principle the proposed approach is applicable in these cases. It is an interesting subject to research, though the mathematical derivation could be a more complex. In present work the boundary condition is based on the exterior DtN map, where the analytic solution in (unbounded) region exterior to computational domain is mapped to the truncation boundary, and the solution in interior region is computed numerically. Under certain conditions it is possible to use a similar approach and build an interior DtN for bounded domains. It also could be combined with domain decomposition FEM tecnhiques and further improve the performances of the method. With the current implementation, the obtained system of linear equations is solved by a classical iterative solver with a custom matrix-vector multiplication routine, and is capable of utilizing multiple cores on a single compute node. This could be improved by implementing a more generic solver that is scalable to several compute nodes via MPI and utilizes GPUs through CUDA or OpenCL. One of the possible applications would be an integration into the Steam Generator Tube Inspection Simulation (SGTSIM) software [119], developed in the Nondestructive Evaluation Laboratory at Michigan State University. Another contribution of this work is an improved fusion algorithm for defect detection and characterization. The method introduces a generic set of fusion coefficients for a single- and multifrequency data, as well as a parametrized nonlinear operation. These coefficients are found 83 by solving a constrained optimization problem and then applied to compute a fused image, which allows to detect the sites with potential defects by computing an appropriate cost function in the corresponding subdomains. The method has been applied to 1D and 2D data and shows a better performance in comparison to conventional multifrequency approaches. The practical importance of the method is the ability to detect smaller defects in presence of noise, where a conventional method fails. Another benefit is an ability to use it for the data acquired on a single frequency, and an efficient implementation for an almost real-time operation in the data acquisition systems. In certain practically important cases (highly magnetic materials, higher levels of noise) the smallest defects could still be sometimes missed. Consequently, the further improvements of the method must aim towards a better performance in these circumstances, potentially through defining a different nonlinear operation or a more elaborated cost function. 84 APPENDICES 85 APPENDIX A USEFUL IDENTITIES (cid:90) ∇ · A  = (cid:90) A · ˆ    A · (B × C) = B · (C × A) ∇ · (A × B) = (∇ × A) · B − A · (∇ × B) ∇ · (∇) = ∇2 + ∇ · ∇ ∇ · (A) = ∇ · A + ∇ · A (A.1) (A.2) (A.3) (A.4) (A.5) 86 APPENDIX B SPHERICAL COORDINATE SYSTEM The standard spherical coordiante system (, , ) is used:  – the distance from origin  2  ∈ [−  ] – zenith angle 2 ,  ∈ [0, 2) – polar angle The gradient and curl in spherical coordinates are: ∇  =   ˆr + ∇ × A = 1   1     ˆθ + 1   ˆφ  (() −  )ˆr + (B.1) (B.2) (B.3) (B.4) 1  ( 1    − ( )) ˆθ + 1  (( ) −  ) ˆφ (B.5) 87 APPENDIX C VECTOR SPHERICAL HARMONICS Let’s use the following definition of scalar spherical harmonics: (cos )  ,  = 0, 1, . . . ,  = −, . . . ,  (cid:115)( − ||)! ( + ||)! || 4 2 + 1 (cid:48)(cid:48)  (cid:90) (, ) = ∗ (cid:48)(cid:48)  = The vector spherical harmonics are defined as Y =  ˆr U = ∇ V = ∇ × ˆr Using the curl representation, we obtain V 1 ∇ × Y =  ∇ × U = −1 ∇ × V = ( + 1) V   Y + 1  U If  represents a scalar function, these relations are valid: 1   V ∇ × (  Y) = ∇ × (  U) = −(  +  ∇ × (  V) = ( + 1)  1  )V   Y + (   + 1   )U 88 (C.1) (C.2) (C.3) (C.4) (C.5) (C.6) (C.7) (C.8) (C.9) (C.10) (C.11) (D.1) (D.2) (D.3) (D.4) (D.5) (D.6) (D.7) APPENDIX D FEM STIFFNESS MATRIX COEFFICIENTS Submatrix elements: (cid:90)  = (∇ × N) · (∇ × N)  Right hand side:  = −   (cid:90) (cid:90)    =    ∇ × N · ( − 1)H d −  (cid:90)  N · A d   =    =   N · N  N · ∇φ   = −    = −   ∇φ · N  ∇φ · ∇φ  (cid:90) (cid:90)   (cid:90) (cid:90)   ∇φ · A d 89 APPENDIX E FIELD GENERATED BY THE CURRENT RING This section contains the solution to the model problem: computation of a magnetic vector potential  and magnetic flux density  of a field generated by a current loop (see [126]). Source is a circular current loop of radius , carrying current , located in  −  plane at origin. Analytic expressions of the field are given by (cid:112)2 + 2 + 2 sin   (2) (2 − 2)(2) − 2(2) 2 [(2 + 2 cos 2)(2) − 2(2)]  = 0   = 2 cos  2   = 22 sin  2 = 2 + 2 − 2 sin  2 = 2 + 2 + 2 sin  2 = 4 sin  2 + 2 + 2 sin  = 1 − 2 2  = 0  (E.1) (E.2) (E.3) (E.4) (E.5) (E.6) (E.7) Here () represent the elliptic integral of the first kind, () represent the elliptic integral of the second kind. 90 APPENDIX F CONDUCTING SPHERE IN A UNIFORM FIELD This section contains an analytic solution to the model problem: conducting sphere in a uniform magnetic field (see [132]). The sphere occupies domain  = {||≤ } and has constant  and . The source is a uniform  along  axis, excitation frequency  is assumed. In this case,  = ˆ1 2 sin . Coloumb gauge ∇ ·  = 0 is enforced. Constants:  =  1 2   = (  )  = () 3 2  =  = ( − 0)−1 2 (2 + 0)−1 2 ( − 0)−1 2 3 + (0(1 + 2) − )1 2 − (0(1 + 2) + 2)1 2 + (0(1 + 2) − )1 2 Total magnetic vector potential inside and outside the sphere:  = ˆ  = ˆ Total magnetic flux density: 1 2 1 2 ) sin  1 2−1 ((  ) 2 3 2 ( + −2) sin  , = −(1 −  23 ) sin  , = (1 +  3 ) cos  2[−1 , = −1 2 (  ) , = −3 2 3 2 ((  ) ((  ) 1 2 1 2 1 2 ) cos  91 1 2 ) − −3 2 3 2 ((  ) 1 2 )] sin  (F.1) (F.2) (F.3) (F.4) (F.5) (F.6) (F.7) (F.8) (F.9) (F.10) (F.11) BIBLIOGRAPHY 92 BIBLIOGRAPHY [1] [2] [3] [4] [5] [6] [7] [13] [14] IEEE Press Series on P. Mix, Introduction to Nondestructive Testing: A Training Guide. Wiley, 2005. P. Moore and S. Udpa, Nondestructive Testing Handbook: Electromagnetic Testing. Nondestructive testing handbook, American Society for Nondestructive Testing, 2004. C. Ye, Y. Huang, L. Udpa, and S. S. Udpa, “Novel rotating current probe with GMR array sensors for steam generate tube inspection,” IEEE Sensors Journal, vol. 16, no. 12, pp. 4995–5002, 2016. J. Xin, N. Lei, L. Udpa, and S. S. Udpa, “Rotating field eddy current probe with bobbin pickup coil for steam generator tubes inspection,” NDT&E International, vol. 54, pp. 45 – 55, 2013. R. Harrington, Time-Harmonic Electromagnetic Fields. Electromagnetic Wave Theory, Wiley, 2001. D. J. Hagemaier and D. J. Hagemaier, Fundamentals of eddy current testing. American Society for Nondestructive Testing, 1990. A. Buffa, H. Ammari, and J. C. Nedelec, “A justification of eddy currents model for the Maxwell equations,” SIAM Journal on Applied Mathematics, vol. 60, no. 5, pp. 1805–1823, 2000. H. Haus and J. Melcher, Electromagnetic Fields and Energy. Prentice Hall, 1989. J. Stratton, Electromagnetic Theory. An IEEE Press classic reissue, Wiley, 2007. [8] [9] [10] R. F. Harrington, Field computation by moment methods. Wiley-IEEE Press, 1993. [11] [12] R. Albanese and G. Rubinacci, “Solution of three dimensional eddy current problems by integral and differential methods,” IEEE Transactions on Magnetics, vol. 24, no. 1, pp. 98– 101, 1988. J. R. Bowler, “Eddy-current interaction with an ideal crack. I. The forward problem,” Journal of Applied Physics, vol. 75, no. 12, pp. 8128–8137, 1994. J. Bowler, S. Jenkins, L. Sabbagh, and H. Sabbagh, “Eddy-current probe impedance due to a volumetric flaw,” Journal of Applied Physics, vol. 70, no. 3, pp. 1107–1114, 1991. J. Jin, Theory and Computation of Electromagnetic Fields. Wiley, 2011. [15] H. Sabbagh, “A model of eddy-current probes with ferrite cores,” IEEE Transactions on Magnetics, vol. 23, no. 3, pp. 1888–1904, 1987. [16] R. Miorelli, C. Reboud, T. Theodoulidis, N. Poulakis, and D. Lesselier, “Efficient modeling of ECT signals for realistic cracks in layered half-space,” IEEE Transactions on Magnetics, vol. 49, no. 6, pp. 2886–2892, 2013. 93 [17] J. R. Bowler and T. Theodoulidis, “Boundary element calculation of eddy currents in cylindrical structures containing cracks,” IEEE Transactions on Magnetics, vol. 45, no. 3, pp. 1012–1015, 2009. [18] R. K. Murphy, H. A. Sabbagh, E. H. Sabbagh, L. Zhou, W. Bernacchi, J. C. Aldrin, D. Forsyth, and E. Lindgren, “Nondestructive damage characterization of complex aircraft structures by inverse methods: advances in multiscale models,” in AIP Conference Proceedings, vol. 1706, p. 090007, AIP Publishing, 2016. [19] H. A. Sabbagh, R. K. Murphy, E. H. Sabbagh, J. C. Aldrin, and J. S. Knopp, “Computational electromagnetics and model-based inversion,” New York, NY, USA: Springer-Verlag doi, vol. 10, pp. 978–1, 2013. [20] M. Morozov, G. Rubinacci, A. Tamburrino, and S. Ventre, “Numerical models with experimental validation of volumetric insulating cracks in eddy current testing,” IEEE Transactions on Magnetics, vol. 42, no. 5, pp. 1568–1576, 2006. [21] Schmidlin, G. , Fischer, U. , Andjelic, Z. and Schwab, C., “Preconditioning of the second- kind boundary integral equations for 3D eddy current problems,” International Journal for Numerical Methods in Engineering, vol. 51, no. 9, pp. 1009–1031, 2001. [22] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, “AIM: Adaptive Integral Method for solving large-scale electromagnetic scattering and radiation problems,” Radio Science, vol. 31, no. 5, pp. 1225–1251, 1996. [23] L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” Journal of Computational Physics, vol. 73, no. 2, pp. 325 – 348, 1987. [24] G. Rubinacci, A. Tamburrino, S. Ventre, and F. Villone, “Fast computational methods for large-scale eddy-current computation,” IEEE Transactions on Magnetics, vol. 38, pp. 529– 532, March 2002. [25] G. Rubinacci, A. Tamburrino, S. Ventre, and F. Villone, “A fast 3-D multipole method for eddy-current computation,” IEEE Transactions on Magnetics, vol. 40, pp. 1290–1293, March 2004. [26] J. Smajic, Z. Andjelic, and M. Bebendorf, “Fast BEM for eddy-current problems using h-matrices and adaptive cross approximation,” IEEE Transactions on Magnetics, vol. 43, pp. 1269–1272, April 2007. [27] Y. Bao, Z. Liu, and J. Song, “Adaptive cross approximation algorithm for accelerating BEM in eddy current nondestructive evaluation,” Journal of Nondestructive Evaluation, vol. 37, p. 68, Aug 2018. [28] K. Ishibashi, “Nonlinear eddy current analysis by the integral equation method,” IEEE Transactions on Magnetics, vol. 30, pp. 3020–3023, Sep. 1994. [29] M. H. Lean and D. S. Bloomberg, “Nonlinear boundary element method for two-dimensional magnetostatics,” Journal of Applied Physics, vol. 55, no. 6, pp. 2195–2197, 1984. 94 [30] K. Ishibashi, “Nonlinear eddy current analysis of induction heating by boundary element method,” Electrical Engineering in Japan, vol. 105, no. 1, pp. 63–71, 1985. [31] K. Ishibashi, “Nonlinear eddy current analysis by volume integral equation method,” IEEE Transactions on Magnetics, vol. 23, pp. 3038–3040, Sep. 1987. [32] M. d’Aquino, G. Rubinacci, A. Tamburrino, and S. Ventre, “Efficient numerical solution of magnetic field problems in presence of hysteretic media for nondestructive evaluation,” IEEE transactions on magnetics, vol. 49, no. 7, pp. 3167–3170, 2013. [33] M. d’Aquino, G. Rubinacci, A. Tamburrino, and S. Ventre, “Three-dimensional computation of magnetic fields in hysteretic media with time-periodic sources,” IEEE transactions on magnetics, vol. 50, no. 2, pp. 53–56, 2014. [34] M. d’Aquino, S. Minucci, C. Petrarca, G. Rubinacci, A. Tamburrino, and S. Ventre, “3D efficient simulation of a magnetic probe for characterization of ferromagnetic specimens,” Studies in Applied Electromagnetics and Mechanics, vol. 40, p. 11, 2015. [35] O. Zienkiewicz and R. Taylor, Finite Element Method. Butterworth-Heinemann, 2000. [36] J. C. Nedelec, “Mixed finite elements in R3,” Numerische Mathematik, vol. 35, no. 3, pp. 315–341, 1980. J. C. Nedelec, “A new family of mixed finite elements in R3,” Numerische Mathematik, vol. 50, no. 1, pp. 57–81, 1986. [37] [38] A. Bossavit and I. Mayergoyz, Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements. Electromagnetism, Elsevier Science, 1998. [39] R. Palanisamy and W. Lord, “Finite element modeling of electromagnetic NDT phenomena,” IEEE Transactions on Magnetics, vol. 15, no. 6, pp. 1479–1481, 1979. [40] R. Palanisamy and W. Lord, “Prediction of eddy current probe signal trajectories,” IEEE Transactions on Magnetics, vol. 16, no. 5, pp. 1083–1085, 1980. [41] W. Lord and R. Palanisamy, “Development of theoretical models for nondestructive testing eddy-current phenomena,” in Eddy-Current Characterization of Materials and Structures, ASTM International, 1981. [42] O. Biro, K. Preis, W. Renhart, K. R. Richter, and G. Vrisk, “Performance of different vector potential formulations in solving multiply connected 3-D eddy current problems,” IEEE Transactions on Magnetics, vol. 26, no. 2, pp. 438–441, 1990. [43] O. Biro, “Edge element formulations of eddy current problems,” Computer Methods in Applied Mechanics and Engineering, vol. 169, no. 3, pp. 391–405, 1999. [44] E. X. Xu and J. Simkin, “Total and reduced magnetic vector potentials and electrical scalar potential for eddy current calculation,” IEEE Transactions on Magnetics, vol. 40, no. 2, pp. 938–940, 2004. 95 [45] Z. Zeng, X. Liu, Y. Deng, L. Udpa, and S. Udpa, “Reduced magnetic vector potential and electric scalar potential formulation for eddy current modeling,” Przeglad Elektrotechniczny, vol. 83, no. 6, pp. 35–37, 2007. [46] Z. Zeng, L. Udpa, S. S. Udpa, and M. S. C. Chan, “Reduced magnetic vector potential formulation in the finite element analysis of eddy current nondestructive testing,” IEEE Transactions on Magnetics, vol. 45, no. 3, pp. 964–967, 2009. J. Jin, The Finite Element Method in Electromagnetics. John Wiley & Sons Canada, Limited, 2003. [47] [48] Y. Kagawa, T. Yamabuchi, and S. Kitagami, “The infinite boundary element method and its application to a combined finite boundary element technique for unbounded field problems,” COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, vol. 2, no. 4, pp. 179–193, 1983. [49] L. Nannen, T. Hohage, A. Schadle, and J. Schoberl, “Exact sequences of high order Hardy space infinite elements for exterior Maxwell problems,” SIAM Journal on Scientific Computing, vol. 35, no. 2, pp. A.1024–A1048, 2013. [50] D. Givoli, “High-order local non-reflecting boundary conditions: a review,” Wave Motion, vol. 39, no. 4, pp. 319–326, 2004. J. Jackson, Classical Electrodynamics. John Wiley & Sons, 1999. [51] [52] R. Albanese and G. Rubinacci, “Magnetostatic field computations in terms of two-component vector potentials,” International Journal for Numerical Methods in Engineering, vol. 29, no. 3, pp. 515–532, 1990. J. B. Manges and Z. J. Cendes, “A generalized tree-cotree gauge for magnetic field computation,” IEEE Transactions on Magnetics, vol. 31, no. 3, pp. 1342–1347, 1995. [53] [54] Z. Ren and A. Razek, “Boundary edge elements and spanning tree technique in three- for Numerical dimensional electromagnetic field computation,” International Journal Methods in Engineering, vol. 36, no. 17, pp. 2877–2893, 1993. [55] A. A. Rodriguez, P. Fernandes, and A. Valli, “Weak and strong formulations for the time- harmonic eddy-current problem in general multi-connected domains,” European Journal of Applied Mathematics, vol. 14, no. 4, pp. 387–406, 2003. I. Babuska and B. Guo, “The h, p and h-p version of the finite element method; basis theory and applications,” Advances in Engineering Software, vol. 15, no. 3, pp. 159–174, 1992. [56] [57] P. D. Ledger and S. Zaglmayr, “Hp-Finite element simulation of three-dimensional eddy current problems on multiply connected domains,” Computer Methods in Applied Mechanics and Engineering, vol. 199, no. 49-52, pp. 3386–3401, 2010. [58] H. Xu, C. D. Cantwell, C. Monteserin, C. Eskilsson, A. P. Engsig-Karup, and S. J. Sherwin, “Spectral/hp element methods: Recent developments, applications, and perspectives,” Journal of Hydrodynamics, vol. 30, no. 1, pp. 1–22, 2018. 96 [59] W. Yao, J.-M. Jin, and P. T. Krein, “A highly efficient domain decomposition method applied to 3-D finite-element analysis of electromechanical and electric machine problems,” IEEE Transactions on Energy Conversion, vol. 27, no. 4, pp. 1078–1086, 2012. [60] W. Yao, J.-M. Jin, and P. T. Krein, “A 3D finite element analysis of large-scale nonlinear dynamic electromagnetic problems by harmonic balancing and domain decomposition: FETI-DP on nonlinear dynamic EM problems,” International Journal of Numerical Modelling: Electronic Networks, Devices and Fields, vol. 29, no. 2, pp. 166–180, 2016. [62] [61] S. Auberhofer, O. Biro, and K. Preis, “Discontinuous Galerkin formulation for eddy-current problems,” COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, vol. 28, no. 4, pp. 1081–1090, 2009. I. Perugia and D. Schotzau, “The hp-local discontinuous Galerkin method for low- frequency time-harmonic Maxwell equations,” Mathematics of Computation of the American Mathematical Society, vol. 72, no. 243, pp. 1179–1214, 2003;2002;. I. Babuska, U. Banerjee, and J. E. Osborn, “Survey of meshless and generalized finite element methods: A unified approach,” Acta Numerica, vol. 12, pp. 1–125, 2003. [63] [64] T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994. [65] 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. [66] L. Xuan, Z. Zeng, B. Shanker, and L. Udpa, “Element-free Galerkin method for static and quasi-static electromagnetic field computation,” IEEE Transactions on Magnetics, vol. 40, no. 1, pp. 12–20, 2004. [67] Y. Zhang, K. R. Shao, Y. Guo, and J. D. Lavers, “A boundary meshless method for transient eddy current problems,” IEEE Transactions on Magnetics, vol. 41, no. 10, pp. 4090–4092, 2005. [68] S. Li and X. Cui, “An edge-based smoothed finite element method for nonlinear magnetostatic and eddy current analysis,” Applied Mathematical Modelling, vol. 62, pp. 287–302, 2018. J. L. Volakis, K. Sertel, and B. C. Usner, Frequency domain hybrid finite element methods for electromagnetics, vol. 10;10.;. San Rafael, Calif. (1537 Fourth Street, San Rafael, CA 94901 USA): Morgan & Claypool Publishers, first ed., 2006. [69] [70] F. Matsuoka and A. Kameari, “Calculation of three dimensional eddy current by FEM-BEM coupling method,” IEEE Transactions on Magnetics, vol. 24, no. 1, pp. 182–185, 1988. [71] Z. Ren, F. Bouillault, A. Razek, A. Bossavit, and J. C. Verite, “A new hybrid model using electric field formulation for 3-D eddy current problems,” IEEE Transactions on Magnetics, vol. 26, no. 2, pp. 470–473, 1990. 97 [72] Q. Chen, A. Konrad, and P. P. Biringer, “A finite element-Green’s function method for the solution of unbounded three-dimensional eddy current problems,” IEEE Transactions on Magnetics, vol. 30, no. 5, pp. 3048–3051, 1994. [73] M. S. Agranovich, Sobolev Spaces, Their Generalizations and Elliptic Problems in Smooth and Lipschitz Domains. Cham: Springer International Publishing, 2015 ed., 2015. [74] A. Bossavit and J. Verite, “The "TRIFOU" code: Solving the 3-D eddy-currents problem by using H as state variable,” IEEE Transactions on Magnetics, vol. 19, no. 6, pp. 2465–2470, 1983. [75] J. C. VERITE, “Application of a 3-D eddy current code (TRIFOU) to non-destructive testing,” COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, vol. 3, no. 3, pp. 167–178, 1984. [76] A. Bossavit, “Coupling finite elements and boundary elements in eddy-current computations,” in Advanced Boundary Element Methods, pp. 49–54, Springer, 1988. [77] D. Givoli and J. B. Keller, “A finite element method for large domains,” Computer Methods in Applied Mechanics and Engineering, vol. 76, no. 1, pp. 41–66, 1989. [78] A. S. Deakin and J. R. Dryden, “Numerically derived boundary conditions on artificial boundaries,” Journal of Computational and Applied Mathematics, vol. 58, no. 1, pp. 1–16, 1995. [79] D. Givoli, I. Patlashenko, and J. B. Keller, “Discrete Dirichlet-to-Neumann maps for unbounded domains,” Computer Methods in Applied Mechanics and Engineering, vol. 164, no. 1, pp. 173–185, 1998. [80] M. J. Grote and C. Kirsch, “Dirichlet-to-Neumann boundary conditions for multiple scattering problems,” Journal of Computational Physics, vol. 201, no. 2, pp. 630–650, 2004. [81] M. J. Grote, “Local nonreflecting boundary condition for Maxwell’s equations,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 29-32, p. 3691, 2006. [82] N. R. Council et al., Aging of US Air Force Aircraft, vol. 488. National Academies Press, 1997. J. Schijve, “Fatigue damage in aircraft structures, not wanted, but tolerated?,” International Journal of Fatigue, vol. 31, no. 6, pp. 998–1011, 2009. [83] [84] X. Yu, “On the fatigue crack growth analysis of spliced aircraft wing panels under sequential axial and shear loads,” Engineering Fracture Mechanics, vol. 123, pp. 116–125, 2014. [85] E. Lindgren, J. Knopp, J. Aldrin, G. Steffes, and C. Buynak, “Aging aircraft NDE: capabilities, challenges, and opportunities,” in AIP Conference Proceedings, vol. 894, pp. 1731–1738, American Institute of Physics, 2007. 98 [86] M. Moles, O. Dupuis, F. Cancre, P. Herzog, J. T. Miller, and J. Hatmaker, “Phased array ultrasonic NDT system for fastener inspections,” feb 2008. US Patent 7,328,619. [87] H. Chan, B. Masserey, and P. Fromme, “High frequency guided ultrasonic waves for hidden fatigue crack growth monitoring in multi-layer model aerospace structures,” Smart Materials and Structures, vol. 24, no. 2, p. 025037, 2015. [88] D. Hopkins, M. Datuin, J. Aldrin, M. Warchol, L. Warchol, and D. Forsyth, “Localization and characterization of fatigue cracks around fastener holes using spherically focused ultrasonic probes,” in AIP Conference Proceedings, vol. 1806, p. 090007, AIP Publishing LLC, 2017. J. Xu, T. Liu, X. Yin, B. Wong, and S. B. Hassan, “Automatic x-ray crack inspection for aircraft wing fastener holes,” in 2nd International Symposium on NDT in Aerospace, vol. 1, 2010. [89] [90] T. Hutt and P. Cawley, “Feasibility of digital image correlation for detection of cracks at fastener holes,” NDT & e International, vol. 42, no. 2, pp. 141–149, 2009. [91] A. Gleiter, G. Riegert, T. Zweschper, G. Busse, et al., “Ultrasound lock-in thermography for advanced depth resolved defect selective imaging,” Insight-Non-Destructive Testing and Condition Monitoring, vol. 49, no. 5, p. 272, 2007. [92] P. Underhill and T. Krause, “Eddy current analysis of mid-bore and corner cracks in bolt holes,” NDT & E International, vol. 44, no. 6, pp. 513–518, 2011. [93] G. L. Burkhardt, J. L. Fisher, J. S. Stolte, S. R. Kramer, and K. L. Cobble, “NDE of aging aircraft structure using orthogonal-axis eddy current probes,” in Review of Progress in Quantitative Nondestructive Evaluation, pp. 1021–1027, Springer, 1997. [94] Y. He, F. Luo, M. Pan, F. Weng, X. Hu, J. Gao, and B. Liu, “Pulsed eddy current technique for defect detection in aircraft riveted structures,” NDT & e International, vol. 43, no. 2, pp. 176–181, 2010. [95] P. McIntire, Electromagnetic Testing: Eddy Current, Flux Leakage and Microwave Nondestructive Testing, vol. 4 of Nondestructive Testing Handbook. American Society for Nondestructive Testing, 1986. [96] C. Dolabdjian, G. Wach, and L. Perez, “Improvement in the detection of subsurface fatigue cracks under airframe fasteners using improved rotating giant magneto-resistance magnetometer head,” Insight-Non-Destructive Testing and Condition Monitoring, vol. 49, no. 3, pp. 133–136, 2007. [97] D. Rifai, A. N. Abdalla, K. Ali, and R. Razali, “Giant magnetoresistance sensors: A review on structures and non-destructive eddy current testing applications,” Sensors, vol. 16, no. 3, p. 298, 2016. [98] T. Dogaru, C. Smith, R. Schneider, and S. Smith, “Deep crack detection around fastener holes in airplane multi-layered structures using GMR-based eddy current probes,” in AIP Conference Proceedings, vol. 700, pp. 398–405, American Institute of Physics, 2004. 99 [99] A. Yashan, W. Bisle, and T. Meier, “Inspection of hidden defects in metal-metal joints of aircraft structures using eddy current technique with GMR sensor array,” Proc. 9th ECNDT, Berlin, 2006. [100] C. Ye, Y. Wang, and Y. Tao, “High-density large-scale TMR sensor array for magnetic field imaging,” IEEE Transactions on Instrumentation and Measurement, vol. 68, no. 7, pp. 2594–2601, 2018. [101] D. Motes, J. C. Aldrin, M. Keiser, G. Steffes, and D. S. Forsyth, “GMR sensing array technique validation study for the inspection of multi-layer metallic structures,” in AIP Conference Proceedings, vol. 1511, pp. 1562–1569, American Institute of Physics, 2013. [102] G. Dib, G. Yang, C. Ye, A. Tamburrino, L. Udpa, and S. S. Udpa, “EC-GMR array with rotating current excitation for multilayered riveted structures inspection,” in AIP Conference Proceedings, vol. 1650, pp. 368–376, American Institute of Physics, 2015. [103] G. Yang, G. Dib, L. Udpa, A. Tamburrino, and S. S. Udpa, “Rotating field EC-GMR sensor for crack detection at fastener site in layered structures,” IEEE Sensors Journal, vol. 15, no. 1, pp. 463–470, 2014. [104] C. Ye, Y. Huang, L. Udpa, and S. S. Udpa, “Differential sensor measurement with rotating current excitation for evaluating multilayer structures,” IEEE Sensors Journal, vol. 16, no. 3, pp. 782–789, 2015. [105] C. Ye, L. Udpa, and S. Udpa, “Optimization and validation of rotating current excitation with GMR array sensors for riveted structures inspection,” Sensors, vol. 16, no. 9, p. 1512, 2016. [106] Z. Su, C. Ye, A. Tamburrino, L. Udpa, and S. Udpa, “Optimization of coil design for eddy current testing of multi-layer structures,” International Journal of Applied Electromagnetics and Mechanics, vol. 52, no. 1-2, pp. 315–322, 2016. [107] O. Karpenko, A. Efremov, C. Ye, and L. Udpa, “Multi-frequency fusion algorithm for detection of defects under fasteners with EC-GMR probe data,” NDT & E International, vol. 110, p. 102227, 2020. [108] H. Huang and T. Takagi, “Crack shape reconstruction from noisy signals in ECT of steam generator tube,” in Industrial Electronics Society (IECON) 2000 26th Annual Conference, vol. 4, pp. 2507–2512, IEEE, 2000. [109] P. Horan, P. Underhill, and T. Krause, “Pulsed eddy current detection of cracks in F/A-18 inner wing spar without wing skin removal using modified principal component analysis,” NDT & E International, vol. 55, pp. 21–27, 2013. [110] J. Kim, G. Yang, L. Udpa, and S. Udpa, “Classification of pulsed eddy current GMR data on aircraft structures,” Ndt & E International, vol. 43, no. 2, pp. 141–144, 2010. [111] G. Yang, G. Dib, J. Kim, L. Zhang, J. Xin, and L. Udpa, “Nonlinear, non-stationary image processing technique for eddy current NDE,” in AIP conference proceedings, vol. 1430, pp. 689–696, American Institute of Physics, 2012. 100 [112] F. Lingvall and T. Stepinski, “Automatic detecting and classifying defects during eddy current inspection of riveted lap-joints,” NDT & E International, vol. 33, no. 1, pp. 47–55, 2000. [113] J. S. Knopp, J. C. Aldrin, and K. V. Jata, “Computational methods in eddy current crack detection at fastener sites in multi-layer structures,” Nondestructive Testing and Evaluation, vol. 24, no. 1-2, pp. 103–120, 2009. [114] S. M. Safdarnejad, Z. Su, C. Ye, L. Udpa, and S. Udpa, “Analysis of EC-GMR data for detection of cracks under fasteners (CUFs),” International Journal of Applied Electromagnetics and Mechanics, vol. 52, no. 1-2, pp. 415–423, 2016. [115] B. van den Bos, S. Sahlen, and J. Andersson, “Automatic scanning with multi-frequency eddy current on multi-layered structures,” Aircraft Engineering and Aerospace Technology, 2003. [116] M. S. Safdernejad, O. Karpenko, C. Ye, L. Udpa, and S. Udpa, “A robust multi-frequency mixing algorithm for suppression of rivet signal in GMR inspection of riveted structures,” in AIP Conference Proceedings, vol. 1706, p. 180005, AIP Publishing LLC, 2016. [117] B. A. Lepine, “Eddy current imaging using multi-frequency mixing methods for aircraft structural integrity verification,” in Review of Progress in Quantitative Nondestructive Evaluation, pp. 355–362, Springer, 1998. [118] O. Karpenko, C. Ye, and L. Udpa, “Dual frequency fusion for defect signal enhancement in EC-GMR inspecton of riveted multilayer structures,” NDT & E International, vol. 92, pp. 97–103, 2017. [119] “Steam Generator Tube Inspection Simulation Software (SGTSIM) v4.” https://www. epri.com/#/pages/product/000000003002011673/?lang=en-US, 2017. Online; accessed May 13, 2020. [120] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-D finite element mesh generator with built- in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, pp. 1309 – 1331, 2009. [121] B. Gough, GNU Scientific Library Reference Manual - Third Edition. Network Theory Ltd., 3rd ed., 2009. [122] S. Abhyankar, J. Brown, E. M. Constantinescu, D. Ghosh, B. F. Smith, and H. Zhang, “PETSc/TS: A modern scalable ode/dae solver library,” arXiv preprint arXiv:1806.01437, 2018. [123] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, and H. Zhang, “PETSc Web page.” https://www.mcs.anl.gov/petsc. Online; accessed May 13, 2020. 101 [124] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev, D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, and H. Zhang, “PETSc users manual,” Tech. Rep. ANL-95/11 - Revision 3.13, Argonne National Laboratory, 2020. Online; accessed May 13, 2020. [125] S. Balay, W. D. Gropp, L. C. McInnes, and B. F. Smith, “Efficient management of parallelism in object oriented numerical software libraries,” in Modern Software Tools in Scientific Computing (E. Arge, A. M. Bruaset, and H. P. Langtangen, eds.), pp. 163–202, Birkhäuser Press, 1997. [126] J. C. Simpson, J. E. Lane, C. D. Immer, and R. C. Youngquist, “Simple analytic expressions for the magnetic field of a circular current loop,” 2001. [127] D. Kraft, “A software package for sequential quadratic programming,” Forschungsbericht- Deutsche Forschungs- und Versuchsanstalt fur Luft- und Raumfahrt, 1988. [128] G. Bradski et al., “The opencv library (2000),” Dr. Dobb’s Journal of Software Tools, 2000. [129] “OpenCV.” https://opencv.org/about.html. Online; accessed May 13, 2020. [130] “OpenBLAS.” http://www.openblas.net. Online; accessed May 13, 2020. [131] S. COMSOL AB, Stockholm, “Comsol multiphysics v. 5.2.” www.comsol.com. Online; accessed May 13, 2020. [132] W. Smythe, Static and Dynamic Electricity. International series in pure and applied physics, McGraw-Hill, 1967. [133] C. Balanis, Advanced Engineering Electromagnetics, 2nd Edition. Wiley, 2012. [134] Z. Zeng, L. Udpa, and S. S. Udpa, “An efficient finite element method for modeling ferrite- core eddy current probe,” International Journal of Applied Electromagnetics and Mechanics, vol. 33, no. 1-2, pp. 481–486, 2010. [135] S. Meddahi and V. Selgas, “A mixed–FEM and BEM coupling for a three-dimensional eddy current problem,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 37, no. 2, pp. 291–318, 2003. [136] L. Lehti, J. Keranen, S. Suuriniemi, and L. Kettunen, “Coil winding losses: Decomposition strategy,” IEEE Transactions on Magnetics, vol. 52, no. 1, pp. 1–6, 2016. [137] J. Rappaz and R. Touzani, Mathematical Models for Eddy Currents and Magnetostatics : With Selected Applications. Dordrecht: Springer, 2014 ed., 2013;2014;. 102