# Efficient modeling of ECT signals for realistic cracks in layered half-space

Roberto Miorelli, Christophe Reboud, Theodoros Theodoulidis, Nikolaos

Poulakis, Dominique Lesselier

► To cite this version:

Roberto Miorelli, Christophe Reboud, Theodoros Theodoulidis, Nikolaos Poulakis, Dominique Lesselier. Efficient modeling of ECT signals for realistic cracks in layered half-space. IEEE Transactions on Magnetics, 2013, 49 (6), pp.2886-2892. 10.1109/TMAG.2012.2236102. hal-00767065

HAL Id: hal-00767065

https://centralesupelec.hal.science/hal-00767065

Submitted on 16 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Efficient Modeling of ECT Signals for Realistic Cracks in Layered Half-Space

Roberto Miorelli<fnref n="1" />, Christophe Reboud<fnref n="1" />, Theodoros Theodoulidis<fnref n="2" />, Nikolaos Poulakis<fnref n="3" />, and Dominique Lesselier<fnref n="4" />

Efficient modeling of eddy current testing (ECT) signals is needed in many areas of industry. Design of probes may be improved and interpretation of experimental signals better understood by using dedicated numerical simulation tools, if they are computationally effective and accurate, yet remain simple enough to be applied at an end-user level. A boundary element method (BEM), dedicated to the numerical simulation of ECT signals due to complex narrow cracks within a planar multilayered structure (PMS) and presenting arbitrary orientations, is investigated. The theoretical formulation relies on the calculation of the dyadic Green operator, associated to the PMS, via appropriate vector wave function expansions. Then, the use of the discrete complex image method followed by the application of the generalized pencil of function method is proposed for efficient computation of this operator. Results of the validation of the complete model by comparison with the experimental data acquired in the laboratory-controlled conditions and with data computed by a finite element code are discussed.

*Index Terms*—Boundary element method (BEM), dyadic Green function, eddy current testing (ECT), layered media, narrow cracks.

## I. INTRODUCTION

ECT DDY CURRENT TESTING (ECT) is used in many industrial areas to ensure quality and reliability standards. ECT is widespread, in particular, for inspection of metal structures in power plants and airplanes, for example, where detection and monitoring of mechanical stress and corrosion cracks are major issues. To avoid costly and time-consuming procedures, simulation tools are increasingly used as a complement to hardware experimental design. In addition, real situations often involve complicated layered structures and possibly many topologies of cracks. Therefore, modeling and numerical simulation thereof is more and more employed to better simulate the reality of the ECT signals, which are due to interactions of induced currents with the cracks, and to help their interpretation as well.

Fast and reliable simulations of narrow cracks are hard to address with standard methods such as volume-integral methods (VIMs) and finite element methods (FEMs), in particular since the computational effort required might become very important for complex problems [1], [2]. Indeed, it is easy to understand that a very fine mesh is needed to describe the very small crack volumes due to their tiny openings, while the cracks must still be properly approximated along the other dimensions. The binary element method (BEM) has shown to be an interesting alternative to VIM or FEM for modeling ECT of narrow cracks [3]–[5], [6]. Its main advantages are excellent performance in terms of computational time and a very good accuracy whenever crack openings get close to zero. In addition, this method can be employed in cylindrical structures [7].

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of multilayered structure affected by two arbitrarily oriented straight flaws.</figcaption></figure>

CEA LIST, L2S, and UWM (referring to the acronyms of the authors' laboratories) have collaborated for some time in order to develop a fast and accurate ECT simulation tool. In a recent article [8], the authors proposed a BEM, which is suitable to the case of narrow multiple cracks, be those parallel or orthogonal to one other, within the same layer of a planar multilayered structure (PMS). Validation results showed a good agreement with respect to experimental data, and overall a good computational efficiency was achieved. The aim of the present contribution is to investigate a generalization of this model, strongly needed in view of the targeted applications, by addressing the case of arbitrarily oriented narrow cracks, embedded in any layer of the considered PMS, as sketched in Fig. 1. Here, by arbitrarily oriented cracks (AOCs), narrow cracks arbitrarily positioned and rotated with respect to the vertical axis are those that are meant.

Theory of BEM applied to ECT signal simulations is nowadays rather well established for particular cases and simple topologies. In fundamental works [3], [4], the main assumptions

<footnote n="1">¹CEA, LIST, Laboratoire de Simulation et de Modélisation Électromagnétique, Gif-sur-Yvette 91191, France</footnote>

<footnote n="2">²Department of Mechanical Engineering, University of Western Macedonia, Kozani 50100, Greece</footnote>

<footnote n="3">³Department of Electrical Engineering, Technological Educational Institute of Western Macedonia, Koila 50100, Greece</footnote>

<footnote n="4">⁴Laboratoire des Signaux et Systèmes (L2S) (CNRS-SUPÉLEC-UPS 11) École Supérieure d'Electricité—Plateau de Moulon, Gif-sur-Yvette Cedex 91192, France</footnote>
and validity domain of the BEM are discussed. Recent works have extended the domain of application to a nonideal crack embedded in a plate, with proven efficiency and accuracy [5], [6]. To do so, the accurate derivation of the Green operator associated to the medium undergoing testing, and in particular its singularity, is a key step. The expression of the Green operator is based on its expansion onto a family of vector wave functions [9]–[11] and on generalized reflection and transmission coefficients for transverse electric (TE) and transverse magnetic (TM) modes [12]. It constitutes an important generalization of [8], where the operator was derived using the Hertz vectors.

The Sommerfeld integral (SI) formulation, associated to casting the Green operator via vector wave functions, needs, however, to be properly handled. It is well known that SI kernels are highly oscillatory and that they slowly decay when the distance between source and observation increases. Therefore, their calculation is not easy, and standard adaptive integration techniques might be inaccurate. However, a method based on the discrete complex image method (DCIM), developed in [13], can be used to cast the SI kernels into a form enabling fast and accurate analytical calculations in the spatial domain. The complex coefficients involved in this closed form are obtained by applying an efficient numerical approximation method called the generalized pencil of function (GPOF) [14]–[16].

The present contribution is organized as follows. First, the analytical formulation, describing the general case of cracks embedded in different layers, is introduced as a general extension of the previous work. The derivation of the associated Green operators is then explained in some detail. Subsequently, the treatment of arbitrary crack orientations in PMS is presented. Comparisons of simulation results with experimental data in laboratory-controlled conditions and numerical results obtained with the FEM code Comsol Multiphysics™ v3.5a [17] follow, before conclusions are drawn.

## II. ANALYSIS OF THE PROBLEM

The following work assumes that a given angular frequency $\omega$ and the associated time-harmonic dependence $e^{-i\omega t}$ are omitted inside the formulation. All materials are considered to be linear, isotropic, and nonmagnetic ($\mu = \mu_0$, permeability of air).

### A. Overview of BEM in ECT Signal Simulation

In general, a PMS is made of $N$ layers, where layers 1 and $N$ are both air half-spaces. In our problem, each layer may be affected by the presence of one or more cracks. In the case of Fig. 1, the PMS contains two cracks, located within layers $i$ and $j$. Considering a given couple of cracks, we denote the source crack with subscript $s$ and the observation one with subscript $o$. Their indices span from 1 to $S$, $S$ being the total number of cracks within the PMS.

Let us consider, for simplicity, the case of a single coil inspecting a flawed PMS containing $S$ flaws. The change of the coil impedance $\Delta Z$, due to the crack presence, is computed by using the reciprocity theorem [19] and is given as

$$ \Delta Z = \frac{1}{I^2} \sum_{s=1}^{S} \int_{\Omega_s} \frac{1}{\sigma_s} \left[ \mathbf{n}_s \cdot \mathbf{J}_s^{\text{inc}}(\mathbf{r}) \right] p_s(\mathbf{r}) d\mathbf{r}, \quad \mathbf{r} \in \Omega_s, \quad s \in [1, \dots, S]. \quad (1) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) Typical configuration of interest. (b) Detail of normal dipolar sources inside the crack zone.</figcaption></figure>

In (1), $I$ is the RMS supplied current and $\sigma_s$ is the conductivity associated to the layer affected by the $s$th crack. The incident eddy current density $\mathbf{J}_s^{\text{inc}}(\mathbf{r})$ is emitted by the coil inside the volume $\Omega_s$ when no flaw is present. The crack orientation is identified by its opening direction $\mathbf{n}_s$, $p_s(\mathbf{r})$ represents the horizontal electric dipole (HED) density, oriented toward the direction $\mathbf{n}_s$, corresponding to the opening of the $s$th crack; see Fig. 2(a). This dipole density is the solution of the following system of integral equations:

$$ J_o^{\text{inc}}(\mathbf{r}) = - \sum_{s=1}^{S} f_o(\mathbf{r}) k_o^2 \int_{\Omega_s} G_{os}^{n_o n_s}(\mathbf{r}, \mathbf{r}') p_s(\mathbf{r}') d\mathbf{r}', \quad \mathbf{r}' \in \Omega_s, \quad \mathbf{r} \in \Omega_o, \quad o \in [1, \dots, S]. \quad (2) $$

We define the conductivity $\sigma_o(\mathbf{r})$ of the $o$th flaw and $\sigma_o$ the constant conductivity of the layer hosting the $o$th flaw. Then, in (2), the contrast function $f_o$ is defined as $f_o(\mathbf{r}) = [\sigma_o(\mathbf{r}) - \sigma_o]/\sigma_o$ and is used to obtain the expressions based on fictitious current sources. The wave number $k_o = \sqrt{i\omega\mu_0\sigma_o}$ is associated to the layer hosting the $o$th crack. The electric-electric dyadic Green function $G_{os}^{n_o n_s}(\mathbf{r}, \mathbf{r}')$ generally describes the interaction of crack $s$ on crack $o$. The dyads superscripts $n_o$ and $n_s$ identify the opening orientations of crack $o$ and crack $s$, respectively.

System (2) is solved numerically with the method of moments (MoM) [18] by considering pulse basis and point
test functions, leading to the following linear system of equations:

$$ \begin{bmatrix} J_1^{\text{inc}} \\ \vdots \\ J_S^{\text{inc}} \end{bmatrix} = \begin{bmatrix} G_{11}^{n_1 n_1} & \cdots & G_{1S}^{n_1 n_S} \\ \vdots & \ddots & \vdots \\ G_{S1}^{n_S n_1} & \cdots & G_{SS}^{n_S n_S} \end{bmatrix} \begin{bmatrix} \mathbf{p}_1 \\ \vdots \\ \mathbf{p}_S \end{bmatrix}. \quad (3) $$

Since only one cell is used across the crack gap, the mesh employed is 2-D-like. In (3), $[\mathbf{J}^{\text{inc}}]$ is the known column vector that contains the incident eddy current density, $[\mathbf{G}]$ is the square matrix that contains the Green function elements, and $[\mathbf{p}]$ is the unknown column vector that contains the dipole density values. The accuracy of the method relies on the calculation of the operator $[\mathbf{G}]$, which calls for particular attention in the case of narrow cracks.

## B. Dyadic Green Function in Multilayered Structure: Source and Field Points in Different Layers

In [8], a treatment based on Hertz vector potentials has been applied to find dyadic Green function expressions suitable to handle cases of parallel and orthogonal cracks, buried within the same layer. The approach employed has involved a derivation of the dyadic Green function based on source and reflection parts. In the following, an alternative strategy, based on vector wave functions [9]–[11], is used to describe the effects of a point source onto a field point when they are separated by one or more planar interfaces. This choice is particularly useful in the case of a multilayered structure, when source and observation points are in different layers. It is worth mentioning that the analytical treatment has been kept as general as possible. Therefore, any value of conductivity is admitted, as well as any size of layer thickness. The dyadic Green function is derived starting directly from decompositions into TE and TM modes. Therefore, a direct use of the generalized transmission and reflection coefficients is made. Moreover, the proposed generalization makes the treatments of particular cases easy, like when source and observation points are embedded in the same layer, within a single layer or a half-space. In the following, we only concentrate upon the general case involving dyadic Green function expressions associated to a point source located within a different layer from that of the field point. The Green operator for source and field points located within the same layer has already been addressed in [8]. Let us emphasize that only xx and yx dyadic components are detailed, since only these components are necessary to model flaws interactions with an opening oriented toward the transverse directions (AOCs). This issue of flaw orientation will be dealt with in the next section.

The complete electric–electric dyadic Green function, called hereafter $\underline{\mathbf{G}}_{nm}$, is a solution, in a given layer with wave number $k_m^2 = \omega^2 \varepsilon_n \mu_0$, of the Helmholtz equation

$$ \nabla \times \nabla \times \underline{\mathbf{G}}_{nm}(\mathbf{r}, \mathbf{r}') - k_m^2 \underline{\mathbf{G}}_{nm}(\mathbf{r}, \mathbf{r}') = 0, \quad \mathbf{r}' \in \text{layer } m, \quad \mathbf{r} \in \text{layer } n. \quad (4) $$

Subscripts m and n stand for the layers where source and the field points are located, respectively. Equation (4) needs to be solved in a suitable manner to satisfy the electromagnetic fields behavior on any layer boundary. Therefore, at each layer, $\underline{\mathbf{G}}_{nm}$

must be such that the continuity of the tangential components of electric and magnetic fields at the interfaces is fulfilled, both fields vanishing when $|\mathbf{r}| \to \infty$. Following the works detailed in [9]–[11], we can define the dyadic Green function as

$$ \underline{\mathbf{G}}_{nm}(\mathbf{r}, \mathbf{r}') = k_m^2 (\nabla \times \hat{\mathbf{z}})(\nabla' \times \hat{\mathbf{z}})g^{\text{TE}}(\mathbf{r}, \mathbf{r}') + \frac{k_m^2}{k_{nm}^2} (\nabla \times \nabla' \times \hat{\mathbf{z}})(\nabla' \times \nabla \times \hat{\mathbf{z}})g^{\text{TM}}(\mathbf{r}, \mathbf{r}') \quad (5) $$

where $\hat{\mathbf{z}}$ is the tensor associated to the vertical direction, $k_m^2 = \omega^2 \varepsilon_m \mu_0$, $k_{nm}^2 = \omega^2 \varepsilon_n \mu_0$, and the operators $\nabla(\cdot)$ and $\nabla'(\cdot)$ are applied to the observation and the source coordinates, respectively. Superscripts TE and TM are associated to the transverse electric and magnetic field decomposition. $\underline{\mathbf{G}}(\mathbf{r}, \mathbf{r}')$ in a homogeneous medium reads as

$$ \underline{\mathbf{G}}(\mathbf{r}, \mathbf{r}') = (k_m^2 \mathbf{I} + \nabla \nabla') \frac{e^{ik_m |\mathbf{r} - \mathbf{r}'|}}{4\pi |\mathbf{r} - \mathbf{r}'|}. $$

Expressions associated to the primary and secondary terms for the and TM components are given under the form of SIs as

$$ g^{\text{TE}}(\mathbf{r}, \mathbf{r}') = \frac{i}{4\pi} \int_0^{\infty} \frac{dk_{\rho}}{k_{mz} k_{\rho}} [F_{mn}^{\text{TE}}(k_{\rho}, z, z') J_0(k_{\rho} \rho)] \quad (6) $$

$$ g^{\text{TM}}(\mathbf{r}, \mathbf{r}') = \frac{i}{4\pi} \int_0^{\infty} \frac{dk_{\rho}}{k_{mz} k_{\rho}} [F_{mn}^{\text{TM}}(k_{\rho}, z, z') J_0(k_{\rho} \rho)] . \quad (7) $$

Functions $F_{nm}^{\text{TE}}(k_{\rho}, z, z')$ and $F_{nm}^{\text{TM}}(k_{\rho}, z, z')$ in (6) and (7) can be found in [12]. Then, after substituting these expressions into (5), we obtain

$$ G_{nm}^{xx} = \frac{k_m^2}{4\pi} \times \int_0^{+\infty} \left[ F_{nm}^{\text{TE}}(k_{\rho}, z, z') \sin\left(\frac{k_{mz} \Delta z'}{2}\right) \right] J_0(k_{\rho} \rho) dk_{\rho} \\ + \int_0^{+\infty} \left\{ \left[ \left( \frac{k_{\rho}^2 + k_{mz}^2}{k_{mz}^2} F_{nm}^{\text{TE}}(k_{\rho}, z, z') - \frac{k_m^2 k_{nz}^2}{k_{nm}^2 k_{mz}^2} F_{nm}^{\text{TM}}(k_{\rho}, z, z') \right) \sin\left(\frac{k_{mz} \Delta z'}{2}\right) \right] \times \left( -\frac{\cos^2 \phi}{4\pi} k_{\rho} J_0(k_{\rho} \rho) + \frac{\cos 2\phi}{4\pi \rho} J_1(k_{\rho} \rho) \right) \right\} dk_{\rho}. \quad (8) $$

Superscript xx stands for the preferential directions of the source and the field points, respectively. Subscripts m and n stand for the mth and nth layers, where the source and the observation points are placed, respectively. Furthermore, $J_l$ is the lth-order Bessel function of the first kind, $\rho = \sqrt{x^2 + y^2}$ is the radial distance with respect to the z-axis, and $\phi$ is the angular position of the observation point with respect to the source point. The proposed treatment assumes the relation $k_{iz}^2 = k_i^2 - k_{\rho}^2$ with $i = m, n$

$$ G_{nm}^{yx} = - \frac{\sin \phi \cos \phi}{4\pi} \times \int_0^{\infty} \left[ \left( \frac{k_{\rho}^2 + k_{mz}^2}{k_{mz}^2} F_{nm}^{\text{TE}} - \frac{k_m^2 k_{nz}^2}{k_{nm}^2 k_{mz}^2} F_{nm}^{\text{TM}} \right) \times 2 \sin \left( \frac{k_{mz} \Delta z'}{2} \right) \right] k_{\rho} J_2(k_{\rho} \rho) dk_{\rho}. \quad (9) $$
The evaluation of SIs presented in (8) and (9) has been carried out by employing the DCIM in conjunction with GPOF. Then, the dyadic Green function obtained in closed form has been analytically evaluated [8]. It is worth mentioning that expressions (8) and (9) match exactly those proposed earlier [8], starting from the Hertz vectors, if the layers *m* and *n* are adjacent and have the same properties.

The integration over the discretized source in the vertical interval $[z' - ((\Delta z')/2); z' + ((\Delta z')/2)]$ has been analytically performed. After integration over directions $x'$ and $y'$, in Cartesian coordinates, the dyadic Green functions $G_{nm}^{xx}$ and $G_{nm}^{yx}$ are combined in order to build the matrix block of (3).

### C. Arbitrarily Oriented Crack in PMS

Let us consider now the configuration with two AOCs inside a PMS, sketched in Fig. 2(a). Each crack occupies a volume $\Omega_s$ ($s = 1, 2$), its opening is directed along the direction $\mathbf{n}_s$, and its angle in the transverse plane, perpendicular to the vertical axis, is $\theta_s$. The ECT problem is described by the set of integral state equation (2). Each of these equations is written with respect to the observation domain, their solutions giving the expressions of unknowns $p_s$. In the general case, illustrated in Fig. 2, angles $\theta_1$ and $\theta_2$ can be different from 0° or 90°.

To ensure a good matrix conditioning, and thus a good numerical stability of the solution, each integral equation of system (2) is written in the local coordinate system of its diagonal block. The first equation is written in the frame $(x_1, y_1, z_1)$ and the second one in the frame $(x_2, y_2, z_2)$; see Fig. 2. In each equation, off-diagonal dyadic terms $G_{os}(\mathbf{r}, \mathbf{r}')$ (with $o \neq s$) are combinations of elementary terms (8) and (9), accounting for the relative orientations of the flaws. For instance, in our example, the term $G_{21}$ will be expressed as follows:

$$G_{21}^{n_2 n_1} = G_{21}^{xx}(\mathbf{r} - \mathbf{r}') \cos \theta_{21} + G_{21}^{yx}(\mathbf{r} - \mathbf{r}') \sin \theta_{21} \quad (10)$$

where $\theta_{21} = \theta_2 - \theta_1$. The diagonal terms associated to the crack self-effect ($s = o, \theta_{os} = 0°$) are the same as described in [8].

Then, we can express the numerical system to solve by applying the MoM [18], obtaining the following complex matrix system [1]:

$$\begin{bmatrix} J_1^{\text{inc}} \\ J_2^{\text{inc}} \end{bmatrix} = \begin{bmatrix} G_{11}^{n_1 n_1} & G_{12}^{n_1 n_2} \\ G_{21}^{n_2 n_1} & G_{22}^{n_2 n_2} \end{bmatrix} \begin{bmatrix} p_1 \\ p_2 \end{bmatrix} \quad (11)$$

where diagonal terms on the right-hand side correspond to the contribution of the cracks upon themselves.

In (11), the quantities of the first equation are expressed in the frame $(x_1, y_1, z_1)$ and the quantities of the second one are expressed in the frame $(x_2, y_2, z_2)$. The terms in (10) are the contributions due to one crack on another, calculated with respect to the observation domain. Therefore, each equation in (11) has been written with respect to the frame of the equation. The unknown functions $[p_1]$ and $[p_2]$, associated to cracks 1 and 2, represent the dipole distributions along the normal components $\mathbf{n}_1$ and $\mathbf{n}_2$ of both cracks, respectively. The size of each term of the matrix is $N_{os} \times N_{os}$, with $N_{os}$ representing the total number of fictitious dipoles introduced in the problem. Finally, after solving (11) and applying the reciprocity theorem, the coil response is given by (1) with $S = 2$.

<table><caption>TABLE I
EXPERIMENTAL DATA FOR A MULTILAYERED STRUCTURE BENCHMARK: COIL PARAMETERS, SPECIMEN, AND DEFECT CHARACTERISTICS</caption><thead><tr><th colspan="2">Coil</th><th colspan="2">Crack 1</th></tr></thead><tbody><tr><td>Inner radius r<sub>1</sub></td><td>5.0 mm</td><td>Depth D<sub>1</sub></td><td>1.0 mm</td></tr><tr><td>Outer radius r<sub>2</sub></td><td>9.7 mm</td><td>Length L<sub>1</sub></td><td>29.65 mm</td></tr><tr><td>Length l</td><td>4.0 mm</td><td>Opening O<sub>1</sub></td><td>0.198 mm</td></tr><tr><td>Lift-off lo</td><td>0.1 mm</td><td colspan="2" style="text-align: center;">Crack 2</td></tr><tr><td>Number of turns N</td><td>407</td><td>Depth D<sub>2</sub></td><td>2.0 mm</td></tr><tr><td>Inductance (calc.) L<sub>0c</sub></td><td>2.23 mH</td><td>Length L<sub>2</sub></td><td>29.89 mm</td></tr><tr><td>Inductance (meas.) L<sub>0m</sub></td><td>2.22 mH</td><td>Opening O<sub>2</sub></td><td>0.234 mm</td></tr><tr><td>Frequency f</td><td>1.5 kHz</td><td colspan="2"></td></tr><tr><th>Specimen</th><th>Thickness</th><th colspan="2">Conductivity</th></tr><tr><td>First layer</td><td>1.0 mm</td><td colspan="2">18.72 MS/m</td></tr><tr><td>Second layer</td><td>0.08 mm</td><td colspan="2">0 MS/m</td></tr><tr><td>Third layer</td><td>2.0 mm</td><td colspan="2">17.4 MS/m</td></tr><tr><td>Fourth layer</td><td>half-space</td><td colspan="2">0 MS/m</td></tr></tbody></table>

The choice of the particular coordinate system associated to the crack zone is mandatory. Indeed, in the presence of skewed crack(s), if we express directly the integral equation into an arbitrary coordinate system, different from that of the crack(s), the values associated to the Green dyad will span over an enlarged system of equations that involves the two transverse components instead of the normal component to the crack. Therefore, the linear system in (11) would be artificially larger and the magnitude of the dominant terms would decrease, implying a bad conditioning of the matrix. Then, the solution of the associated complex matrix system would be affected by instabilities and lead to wrong results.

This formulation and its implementation have been tested in several cases that are representative of realistic ECT configurations. Results of validations are presented in the next section.

### III. RESULTS AND VALIDATION

Hereafter, we first present an experimental validation and then a numerical comparison with results obtained with Comsol Multiphysics™ v3.5a [17]. The first case consists of a coil inspecting a stratified planar structure affected by two cracks. The second one is a plate inspection where the test piece is affected by a complex crack.

#### A. Experimental Validation for Multilayered Structure

This validation is based on a set of measurements performed by the MEANDER laboratory, whose experience in precision measurements has been well demonstrated in many experiments [20], [21]. A planar structure made of three layers has been studied. Each one of two metallic plates contains an electrodischarged machined (EDM) through-wall slot. Both these plates are separated by a thin sheet of dielectric insulator. The coil, located above the piece, has been moved by a precision XY-scanner with a spatial resolution of 0.025 mm. Measurements were made using the precision LCR bridge Agilent 4284A and the whole setup has been controlled via PC using Labview. Details of coil, specimens, and cracks parameters are given in Table I.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Experimental setup for AOCs embedded in different layers. Two parallellepiped notches with an angle of 45° are inspected with an absolute coil at 1.5 kHz. Views from both sides are given on the left and above the top view.</figcaption></figure>

Both slabs have been arranged in order to get an angle of 45° between the cracks that share the same center position in the xy-plane (see Fig. 3). The mesh, used during the BEM simulation of each crack, is chosen following the criterion of five elements per skin depth for the crack depth and ten elements per coil outer radius for the crack length.

An inspection of the multilayered structure has been simulated using 60 × 60 positions in steps of 1 mm. To compare the numerical results with experimental measurements, a 1-D signal has been extracted along the axis of crack 1. In Fig. 4(a), the complex plane curves depict the variation of the coil impedance simulated with BEM, compared with experimental measurements. In Fig. 4(b), the normalized real and imaginary parts of the impedance variation $\Delta R/X_0$ and $\Delta X/X_0$, respectively, are also plotted with respect to the coil position. The normalization factor $X_0$ is the reactance of the coil in air where, for the measurement, we have employed factor $X_0 = \omega L_{0,m}$, whereas for the numerical results, $X_0 = \omega L_{0,c}$ is used. Simulation results obtained agree very well with experimental data. Moreover, an excellent performance in terms of computation time has been achieved, since 3600 points of the complete map have been calculated in about one and a half minute on an Intel-Q9550 @2.83 GHz with 8 GB of RAM. Let us emphasize that, due to the very small gap between the two plates, the corresponding FEM calculations do not give satisfactory results at all, due to the extremely large mesh needed to describe the tiny dielectric sheet. The VIM [1] code is also not able to correctly handle a problem like that one, involving two long cracks with very narrow openings. Complete convergence is not achieved due to the lack of memory. Therefore, only BEM results and experimental data are presented in Fig. 4.

## B. Numerical Comparisons for AOCs

In order to validate the theoretical results for configurations involving more complex topologies, we propose here a comparison between Comsol Multiphysics™ v3.5a [17] and the theoretical approach. The geometry studied is shown in Fig. 5, where a succession of five cracks embedded in a plate has been created in order to approximate a hypothetical crack with complex profile. Characteristics of the proposed case are listed in Table II.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Coil impedance change of the probe scanning performed along the axis of crack 1. Comparison between BEM (solid line) and experimental data (circles) in (a) the complex impedance plane and (b) for normalized imaginary and real parts with respect to the coil scan positions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Plate structure affected by a succession of AOCs used to approximate a realistic crack line. Views from the top and both sides are given on the left and above the main plot.</figcaption></figure>

One can notice that in Fig. 5 very tiny gaps (less than 200 µm) are present between every couple of cracks. This choice has been taken in order to simulate as coherently as possible the problem with the two numerical methods. Indeed, differences in mesh methodology are present between FEM and BEM when
<table><caption>TABLE II
ARBITRARILY ORIENTED CRACKS: COIL PARAMETERS,
SPECIMEN, AND DEFECTS CHARACTERISTICS</caption><thead><tr><th colspan="2">Coil</th></tr></thead><tbody><tr><td>Inner radius r<sub>1</sub></td><td>3.0 mm</td></tr><tr><td>Outer radius r<sub>2</sub></td><td>5.0 mm</td></tr><tr><td>Length l</td><td>1.0 mm</td></tr><tr><td>Lift-off lo</td><td>0.2 mm</td></tr><tr><td>Number of turns N</td><td>200</td></tr><tr><td>Inductance L<sub>0</sub></td><td>0.379 mH</td></tr><tr><td>Frequency f</td><td>10 kHz</td></tr></tbody><thead><tr><th colspan="2">Crack 1</th></tr></thead><tbody><tr><td>Depth D<sub>1</sub></td><td>2.0 mm</td></tr><tr><td>Length L<sub>1</sub></td><td>4.75 mm</td></tr><tr><td>Opening O<sub>1</sub></td><td>0.2 mm</td></tr><tr><td>Angle θ<sub>1</sub></td><td>25 deg</td></tr></tbody><thead><tr><th colspan="2">Crack 2</th></tr></thead><tbody><tr><td>Depth D<sub>2</sub></td><td>2.0 mm</td></tr><tr><td>Length L<sub>2</sub></td><td>5.2 mm</td></tr><tr><td>Opening O<sub>2</sub></td><td>0.2 mm</td></tr><tr><td>Angle θ<sub>2</sub></td><td>165 deg</td></tr></tbody><thead><tr><th colspan="2">Crack 3</th></tr></thead><tbody><tr><td>Depth D<sub>3</sub></td><td>2.0 mm</td></tr><tr><td>Length L<sub>3</sub></td><td>2.5 mm</td></tr><tr><td>Opening O<sub>3</sub></td><td>0.2 mm</td></tr><tr><td>Angle θ<sub>3</sub></td><td>150 deg</td></tr></tbody><thead><tr><th colspan="2">Crack 4</th></tr></thead><tbody><tr><td>Depth D<sub>4</sub></td><td>2.0 mm</td></tr><tr><td>Length L<sub>4</sub></td><td>5.0 mm</td></tr><tr><td>Opening O<sub>4</sub></td><td>0.2 mm</td></tr><tr><td>Angle θ<sub>4</sub></td><td>15 deg</td></tr></tbody><thead><tr><th colspan="2">Crack 5</th></tr></thead><tbody><tr><td>Depth D<sub>5</sub></td><td>2.0 mm</td></tr><tr><td>Length L<sub>5</sub></td><td>4 mm</td></tr><tr><td>Opening O<sub>5</sub></td><td>0.2 mm</td></tr><tr><td>Angle θ<sub>5</sub></td><td>170 deg</td></tr></tbody><thead><tr><th colspan="2">Specimen</th></tr></thead><tbody><tr><td>Plate thickness</td><td>2.0 mm</td></tr><tr><td>Conductivity</td><td>17.0 MS/m</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Fig. 6. BEM simulated amplitude map of a complex crack line in a plate structure with 81 × 81 points.</figcaption></figure>

two cracks are partially overlapped. This choice, because of the frequency regime adopted and the specimen conductivity, will give negligible differences compared to a continuous crack line. Furthermore, we have decided to simulate a plate structure with through-wall cracks, because strongest perturbations can give us more trustworthy comparisons of the results.

The amplitude map of the ECT signal, simulated with BEM, is presented in Fig. 6. It corresponds to 81 × 81 positions of the

probe above the PMS, with scanning step of 0.5 mm. The only regions meshed uniformly in 2-D are the flawed regions, using the same criterion as before. To compare such results with those provided by the FEM code, we have carried out three vertical signal extractions at three different positions, as shown also in Fig. 6.

A good agreement between BEM results and FEM simulations can be observed in Figs. 7 and 8, respectively. Moreover, the BEM model has preserved its efficiency since the simulation of the complete map has taken about 12 min on our PC. This time corresponds to the calculation of two coil positions for the FEM code, where the whole electromagnetic scenario has been described with about 71 100 elements. An even better agreement and efficiency has been achieved with BEM for the same configuration at frequencies of 1, 2.5, and 5 kHz.

## IV. CONCLUSION

The generalization of the boundary element model applied in [8] to the ECT signal simulation of a planar multilayered structure has been detailed. This generalization enables to model far more complex crack topologies than before, in terms of orientation and of location inside the PMS. The dyadic Green function (of the layered embedding medium) has been expressed, in a more general way, via vector wave functions instead of

<figure><img src="image_8.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b)</figcaption></figure>

<caption>Fig. 7. Quantitative comparison of BEM (solid lines) and FEM simulations (circles). The change of coil impedance is normalized by its reactance $X_0$ in air. Signal extraction along the crack line shown in Fig. 6: (a) complex plane representation and (b) real and imaginary parts versus coil position.</caption>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Quantitative comparison of BEM (solid lines) and FEM simulations (circles). The change of coil impedance is normalized by its reactance $X_0$ in air. Complex plane representation for signal extraction along the external scan line (a) in -1.5 mm and (b) in 1.5 mm as shown in Fig. 6.</figcaption></figure>

Hertz vectors. Results of the simulations, compared with experimental data and FEM code results, have confirmed the excellent accuracy of the approach, as well as its computational efficiency. This makes the proposed BEM a simulation tool—to be released within the forthcoming version 11 of CIVA software [22]—adapted to applications requiring intensive use of simulation, among which we can cite parametric studies, probe design, database generation, or probability of detection (POD) studies.

Some interesting perspectives of this work are the modeling of cracks with even more complex orientations, and the modeling of configurations involving both narrow and volumetric cracks. Moreover, the precision of MoM can also be improved by the use of higher order projection functions like the Rao-Wilton-Glisson (RWG) basis functions [23].

## ACKNOWLEDGMENT

This work was supported by the CIVAMONT project, aiming at developing scientific collaborations around the nondestructive testing simulation platform CIVA developed at CEA, LIST.

## REFERENCES

[1] C. Reboud, G. Pichenot, S. Paillard, and F. Jenson, "Simulation and POD studies of riveted structures inspected using eddy current techniques," in *Electromagnetic Non-Destructive Evaluation (XIII), Studies in Applied Electromagnetics and Mechanics*, J. Knopp, M. Blodgett, B. Wincheski, and N. Bowler, Eds. Amsterdam, The Netherlands: IOS Press, 2009, pp. 129-136.

[2] G. Rubinacci, A. Tamburrino, S. Ventre, and F. Villone, "A fast 3-D multipole method for eddy-current computation," *IEEE Trans. Magn.*, vol. 40, no. 2, pp. 1290-1293, Mar. 2004.

[3] J. R. Bowler, "Eddy-current interaction with an ideal crack. I. The forward problem," *J. Appl. Phys.*, vol. 75, pp. 8128-8137, Jun. 1994.

[4] J. R. Bowler, "Inversion of open cracks using eddy-current probe impedance," *Rev. Progr. Quant. Nondestruct. Eval.*, vol. 19A, pp. 529-533, 2000.

[5] T. Theodoulidis, "Developments in efficiently modelling eddy current testing of narrow cracks," *NDT&E Int.*, vol. 43, pp. 591-598, Oct. 2010.

[6] T. Theodoulidis, "Efficient computation of eddy current cracks signals," in *Electromagnetic Non-Destructive Evaluation (XIV), Studies in Applied Electromagnetic and Mechanics*, T. Chady, S. Gratkowsky, T. Takagi, and S. S. Udpa, Eds. Amsterdam, The Netherlands: IOS Press, 2010, pp. 96-102.

[7] J. R. Bowler, T. Theodoulidis, H. Xie, and Y. Ji, "Evaluation of eddy current probe signals due to cracks in fastener holes," *IEEE Trans. Magn.*, vol. 48, no. 3, pp. 1159-1170, Mar. 2012.

[8] R. Miorelli, C. Reboud, D. Lesselier, and T. Theodoulidis, "Eddy current modeling of narrow cracks in planar-layered metal structures," *IEEE Trans. Magn.*, vol. 48, no. 10, pp. 2551-2559, Oct. 2012.

[9] W. C. Chew, *Waves and Fields in Inhomogeneous Media*, ser. Electromagnetic Waves. New York, NY, USA: IEEE Press, 1995.

[10] W. C. Chew, J. L. Xiong, and M. A. Saville, "A matrix-friendly formulation of layered medium Green's function," *IEEE Antennas Wireless Propag. Lett.*, vol. 5, pp. 490-494, 2006.

[11] J. L. Xiong and S. C. Chew, "A newly developed formulation suitable for matrix manipulation of layered medium Green's functions," *IEEE Trans. Antennas Propag.*, vol. 58, no. 3, pp. 868-879, Mar. 2010.

[12] W. C. Chew and S. Y. Chen, "Response of a point source embedded in a layered medium," *IEEE Antennas Wireless Propag. Lett.*, vol. 2, pp. 254-258, 2003.

[13] Y. L. Chow, J. J. Yang, D. G. Fang, and G. E. Howard, "A closed-form spatial Green's function for the thick microstrip substrate," *IEEE Trans. Microw. Theory Tech.*, vol. 39, no. 3, pp. 588-592, Mar. 1991.

[14] Y. Hua and T. K. Sarkar, "Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoid in noise," *IEEE Trans. Acoust. Speech Signal Process.*, vol. 38, no. 5, pp. 814-824, May 1990.

[15] M. I. Aksun, "A robust approach for the derivation of closed-form Green's functions," *IEEE Trans. Microw. Theory Tech.*, vol. 44, no. 5, pp. 651-658, May 1996.

[16] Y. Ge and K. P. Esselle, "A new closed-form Green's function for microstrip structures—theory and results," *IEEE Trans. Microw. Theory Tech.*, vol. 50, no. 6, pp. 1556-1560, June 2002.

[17] Comsol [Online]. Available: http://www.comsol.com

[18] R. F. Harrington, *Field Computation by Moment Methods*. New York, NY, USA: MacMillan/Krieger, 1983.

[19] B. A. Auld, F. G. Muennemann, and M. Riaziat, "Quantitative modelling of flaw responses in eddy current testing," in *Nondestructive Testing*. London, U.K.: Academic, 1984, vol. 7, ch. 2, pp. 37-76.

[20] T. Theodoulidis, N. Poulakis, and A. Dragogias, "Rapid computation of eddy current signals from narrow cracks," *NDT & E Int.*, vol. 43, pp. 13-19, Jan. 2010.

[21] J. R. Bowler, T. Theodoulidis, and N. Poulakis, "Eddy current probe signals due to a crack at a right-angled corner," *IEEE Trans. Magn.*, vol. 48, no. 12, pp. 4735-4746, Dec. 2012.

[22] CEA-LIST, "CIVA: Solutions pour les CND," [Online]. Available: www-civa-cea.fr

[23] S. M. Rao, D. R. Wilton, and A. W. Glisson, "Electromagnetic scattering by surface of arbitrary shape," *IEEE Trans. Antennas Propag.*, vol. AP-30, no. 3, pp. 409-418, May 1982.