Fast computation of the response of any 3D Eddy Current probe scanning a 3D stratified conductor characterized by a set of arbitrary interfaces of complex shape

Abstract

This paper concerns the development of a semi-analytical model dedicated to the fast computation of the response of a 3D Eddy Current (EC) probe scanning a planar conductor of complex shape. The workpiece is characterized by a finite number of complex interfaces k, 1 ≤ k ≤ N, each one being defined by any 2D arbitrary surface a _k(x, y). A previous semi-analytical model based on the Curvilinear Coordinate Method (CCM) has been presented for the computation of quasi-static fields induced by any 3D EC probe scanning a half-space of complex shape. This paper gives a natural extension of this preliminary work. CCM consists in introducing a change of coordinates in order to be able to write analytically and easily boundary conditions separating two media. No mesh is needed, nevertheless, the covariant form of Maxwell’s equations is required due to a novel generalized metric space. Boundary conditions at multiple interfaces are efficiently implemented thanks to the S-matrix algorithm. Finally, some numerical experiments show the validity of the numerical model by some comparison between simulated data obtained by a FE commercial code and those provided by the proposed numerical model.

Keywords

Covariant form of Maxwell’s equations eddy current simulation Eddy current nondestructive testing

1. Introduction

In order to answer to industrial needs, the fast computation of a 3D eddy current (EC) probe scanning a conductor of complex shape is often required. Though some efficient semi-analytical models, based on the Green’s dyad formalism, have been developed and largely implemented notably into EC modules themselves integrated into the CIVA platform, these numerical models can address only canonical geometries such as planar stratified media or tubes of finite thickness. In order to avoid the use of purely numerical methods for some complex geometries, the Curvilinear Coordinate Method (CCM) [1] which is widely used in the optical society for solving rigorously some scattering problems for crossed gratings and periodic structures has been evaluated and transferred from the high frequency range to the low frequency range. This efficient and original method based on the covariant form of Maxwell’s equations has been applied recently for Eddy Current calculations in the planar case [2–4] for 2.5D configurations characterized by a 3D eddy current probe scanning a 2D layered stratified conducting media. In this case, the geometry of the workpiece is described by a set of 1D analytical profiles a _k(x),1 ≤ k ≤ N. The extension to a 3D problem has been dealt with by applying the formalism to a 3D half-space conducting media [5] characterized by a single arbitrary surface a(x, y) depending on the two directions X and Y. The main advantage of this formalism comes from the fact that no mesh is required since all boundary conditions which must be satisfied at each interface can be analytically written by some obvious equalies. Thanks to a smart change of variables, by using the covariant form of Maxwell’s equations, it is possible to compute efficiently a modal expansion of the tangential components of the electromagnetic fields taking into account the complex geometry of the interface. In this paper, the properties of the stratified media are included into the so called S-matrix algorithm [6,7] in order to compute efficiently the response of a 3D eddy current probe scanning a 3D multilayered complex structures. Nevertheless, the main difficulty to overcome comes from the numerical computation of eigenmodes characterizing the structure since the discretization of the problem requires some truncation of the fields into the Fourier domain, so a finite number of modes can limit some kinds of analytical shapes of profiles. This paper is organized as follows: firstly, a new curvilinear coordinate system is associated to each layer, a new tensor metric is therefore introduced. A modal expansion of the tangential components of the field is obtained from the modal expansion of the longitudinal components of the electromagnetic field. The S-Matrix algorithm is then described for parallel interfaces and non-parallel interfaces implying in this case another transformations. To validate the new numerical model, some numerical experiments are performed and some comparison to other simulated data confirm the validity of the model. In the last section, some discussion is carried on concerning some limitations of the approach and some comments are given for new developments in the future.

2. A summary of the formalism

2.1. Principle of the change of coordinate system

In this section, the 3D formalism of the Curvilinear Coordinate Method (CCM) is recalled and translated in our context of ECNDT, in the low frequency range. Lets us consider a stratified conducting media characterized by a set of interfaces (see Fig. 1). Each interface separating two consecutive media is described by an analytic function a _p(x, y) depending on the two directions x and y. The air domain, denoted by 0, contains an EC current probe. The subscript p stands for other interfaces so that 1 ≤ p ≤ N. The last media is assumed to be infinite. Each interface is associated to a non-orthogonal coordinates system: x ¹ = x, x ² = y and $x_{p}^{3}=z-a_{p}(x,y)$ . Therefore, boundary conditions at the interface becomes easier to write since the covariant coordinate $x_{p}^{3}$ is equal to a constant. Due to a change of coordinates, a new metric tensor must be introduced corresponding to a map between the initial cartesian coordinate system $\{x,y,z\}$ and the new coordinates system $\{x^{1},x^{2},x_{p}^{3}\}$ . Since the partial derivative forms denoted by ${\dot{a}}_{1}=\partial _{x}a_{p}$ and ${\dot{a}}_{2}=\partial _{y}a_{p}$ with respect to the analytical function a _p(x, y) = a _p(x ¹, x ²) exist depending on the two variables x and y, the tensor metric g _ij and the conjugate tensor metric g ^ij = [g _ij]⁻¹ are given by: $\begin{eqnarray}\displaystyle g_{ij}=\left[\begin{array}{@{}ccc@{}}1+{\dot{a}}_{1}^{2} & {\dot{a}}_{1}{\dot{a}}_{2} & {\dot{a}}_{1}\\[3.0pt] {\dot{a}}_{1}{\dot{a}}_{2} & 1+{\dot{a}}_{2}^{2} & {\dot{a}}_{2}\\[3.0pt] {\dot{a}}_{1} & {\dot{a}}_{2} & 1\end{array}\right],\quad [g^{ij}]=\left[\begin{array}{@{}ccc@{}}1 & 0 & -{\dot{a}}_{1}\\[3.0pt] 0 & 1 & -{\dot{a}}_{2}\\[3.0pt] -{\dot{a}}_{1} & -{\dot{a}}_{2} & 1+{\dot{a}}_{1}^{2}+{\dot{a}}_{2}^{2}\end{array}\right]=\left[\begin{array}{@{}ccc@{}}1 & 0 & g^{13}\\[3.0pt] 0 & 1 & g^{23}\\[3.0pt] g^{13} & g^{23} & g^{33}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (1) The determinant g = |g _ij| = 1. The covariant components of the electrical (magnetic) field E ₁, E ₂ (H ₁, H ₂), in this natural coordinate system, correspond to the tangential components of the electrical (magnetic) field at the interface p.

Fig. 1.

A 3D air-core probe scanning a 3D stratified conductor media. Only a few smooth interfaces are represented.

2.2. Modal expansion of the tangential components

In what follows, let us assume some notations. let denoting by $k={\omega}\sqrt{{\varepsilon}_{0}{\mu}_{0}}$ the wavenumber in vacuum while in each layer l, the wavenumber is given by $k_{l}^{2}=k^{2}{\mu}_{l}{\varepsilon}_{l}$ where ϵ₀ and 𝜇₀ stand for the vacuum magnetic permeability and the vacuum dielectric permittivity and ϵ_l and 𝜇_l stand for the relative magnetic permeability and the relative dielectric permittivity of the medium l. Denoting by σ_l the conductivity of the medium, the relative permittivity is a complex number given by ${\varepsilon}_{l}=∼1-\frac{i{\sigma}_{l}}{{\omega}{\varepsilon}_{0}}$ (assuming the time dependance e ^iωt). Moreover the vacuum wave impedance is denoted by $Z_{0}=\sqrt{\frac{{\mu}_{0}}{{\varepsilon}_{0}}}$ . Now, keeping in mind that the contravariant metric tensor g ^ij does not depend on the variable $x_{p}^{3}$ , (i.e., $\partial_3$ g ^ij ≡ 0), it is possible [5,8], from the covariant form of Maxwell equations [9], to deduce, at any interface, the tangential components of the electromagnetic field E ₁, E ₂, H ₁, H ₂ from two longitudinal components E ₃ and Z ₀ H ₃. In summary, the application of the principle of separation of variables and the assumption that an exponential dependance according to the covariant coordinate x ³ can be applied so that ${\Phi}(x^{1},x^{2},x^{3})={\phi}(x^{1},x^{2})∼e^{-i{\gamma}_{p}x_{p}^{3}}$ . The two longitudinal components E ₃ and H ₃ must satisfy the same Helmholtz’s equation translated in the new coordinates system [8]: $\begin{eqnarray}\displaystyle \left[g^{33}\partial _{3}^{2}+(\partial _{1}g^{13}+g^{13}\partial _{1})\partial _{3}+(\partial _{2}g^{23}+g^{23}\partial _{2})\partial _{3}+\partial _{1}^{2}+\partial _{2}^{2}+k_{l}^{2}\right]{\phi}_{p}=0 & & \displaystyle\end{eqnarray}$ (2) where 𝜙_p ≡ E ₃ or H ₃ and the partial derivative operator $\partial _{3}^{2}$ is substituted by $-{\gamma}_{p}^{2}$ . By using an auxiliary variable $i{\phi}_{p}^{\prime }={\gamma}_{p}{\phi}_{p}$ , a differential system of first order is given by: $\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}{\mathcal{L}}_{p}\end{array}\right]\left[\begin{array}{@{}c@{}}i{\phi}_{p}^{\prime }\\[3.0pt] {\phi}_{p}\end{array}\right]={\gamma}_{p}\left[\begin{array}{@{}c@{}}i{\phi}_{p}^{\prime }\\[3.0pt] {\phi}_{p}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (3) Finally, in order to obtain an algebraic version of this differential equation, a discrete version of a 2D-Fourier transform along the covariant coordinates x ¹ and x ² is considered. The spatial frequencies are denoted by 𝛼 and 𝛽 respectively. Therefore, the potential ${\phi}_{p}(x^{1},x^{2},x^{3})={\mathcal{T}}{\mathcal{F}}^{-1}[\hat{{\phi}}_{p}({\alpha},{\beta},x^{3})]$ may be reconstructed in the spatial domain from the potential $\hat{{\phi}}_{p}$ defined in the spectral domain. Nevertheless, a finite number of modes leads to a truncation of the spectral domain ([−M, …, +M] for 𝛼 and [−N, …, +N] for 𝛽. Thus, the translation of the differential system (3) in the truncated Fourier domain can give a modal expansion of the two longitudinal components. The longitudinal components of the fields are expressed as a linear combination of a TE field (E ₃ = 0) and a TM field (H ₃ = 0): $\begin{eqnarray}\displaystyle \hskip -1pc\left[\begin{array}{@{}c@{}}\widehat{E}_{3}^{(p)}\\[3.0pt] Z_{0}\widehat{H}_{3}^{(p)}\end{array}\right]=\left[\begin{array}{@{}cc@{}}\boldsymbol{{\psi}}_{p,q}^{TE(+)} & 0\\[3.0pt] 0 & \boldsymbol{{\psi}}_{p,q}^{TM(+)}\end{array}\right]\left[\begin{array}{@{}cc@{}}\boldsymbol{d}_{p}^{+} & 0\\[4.0pt] 0 & \boldsymbol{d}_{p}^{+}\end{array}\right]\left[\begin{array}{@{}c@{}}b_{p,q}^{TE(+)}\\[4.0pt] b_{p,q}^{TM(+)}\end{array}\right]+\left[\begin{array}{@{}cc@{}}\boldsymbol{{\psi}}_{p,q}^{TE(-)} & 0\\[4.0pt] 0 & \boldsymbol{{\psi}}_{p,q}^{TM(-)}\end{array}\right]\left[\begin{array}{@{}cc@{}}\boldsymbol{d}_{p}^{-} & 0\\[4.0pt] 0 & \boldsymbol{d}_{p}^{-}\end{array}\right]\left[\begin{array}{@{}c@{}}b_{p,q}^{TE(-)}\\[4.0pt] b_{p,q}^{TM(-)}\end{array}\right] & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$ (4) where $\boldsymbol{{\psi}}_{p,q}^{\pm }$ denotes the eigenvectors associated to the eigenvalues ${\gamma}_{p,q}^{\pm }$ and $\boldsymbol{d}_{p}^{\pm }$ is a diagonal matrix diag[e ^{i𝜆_{p, −M}}, …, e ^{i𝜆_{p, +M}}]. The coefficients $b_{p,q}^{TE\pm }$ and $b_{p,q}^{TM\pm }$ must be determined by applying boundary conditions. The forward/backward contributions of the field differ from the sign ^± of the exponents and non-physical solutions (growing towards infinity) are suppressed. Let us consider $\boldsymbol{b}_{p}^{\pm }=[b_{p,q}^{TE(\pm )},b_{p,q}^{TM(\pm )}]^{T}$ . Finally a modal expansion of the tangential components of the fields is summarized as a global form: $\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}\widehat{\mathbf{E}}_{2}^{(p)}\\[4.0pt] Z_{0}\widehat{\mathbf{H}}_{1}^{(p)}\\[4.0pt] Z_{0}\widehat{\mathbf{H}}_{2}^{(p)}\\[4.0pt] \widehat{\mathbf{E}}_{1}^{(p)}\end{array}\right]=\left[\begin{array}{@{}cc@{}}\boldsymbol{{\psi}}_{p,q}^{(+)} & \boldsymbol{{\psi}}_{p,q}^{(-)}\end{array}\right]\left[\begin{array}{@{}cc@{}}\boldsymbol{D}_{p}^{+} & 0\\[4.0pt] 0 & \boldsymbol{D}_{p}^{-}\end{array}\right]\left[\begin{array}{@{}c@{}}\boldsymbol{b}_{p}^{(+)}\\[4.0pt] \boldsymbol{ b}_{p}^{(-)}\end{array}\right]. & & \displaystyle\end{eqnarray}$ (5)

2.3. Multiple interfaces: the S-Matrix algorithm

The S-Matrix algorithm described in [3,7] can lead to a relationship between the amplitudes $\mathbf{b}_{1}^{+}$ and $\mathbf{b}_{1}^{-}$ in the first conducting layer relatively to the amplitudes $\mathbf{b}_{N}^{+}$ and $\mathbf{b}_{N}^{-}$ in the last layer by using: $\left[\begin{array}{@{}c@{}}\mathbf{b}_{1}^{+}\\[3.0pt] \mathbf{b}_{N}^{-}\end{array}\right]=S(N,1)\left[\begin{array}{@{}c@{}}\mathbf{b}_{1}^{-}\\[3.0pt] \mathbf{b}_{N}^{+}\end{array}\right]$ where S (N,1) results from a recursive scheme implying the S-interface matrix linking two successive layers: $\left[\begin{array}{@{}c@{}}\mathbf{b}_{j-1}^{+}\\[3.0pt] \mathbf{b}_{j}^{-}\end{array}\right]=S(j,j-1)\left[\begin{array}{@{}c@{}}\mathbf{b}_{j-1}^{-}\\[3.0pt] \mathbf{b}_{j}^{+}\end{array}\right]$ . The amplitudes $\mathbf{b}_{1}^{+}$ and $\mathbf{b}_{1}^{-}$ must be combined to another condition for taking into account the first interface air-conductor: $\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}\boldsymbol{b}_{0}^{+}\\[3.0pt] \boldsymbol{ b}_{1}^{-}\end{array}\right]\stackrel{x_{0}^{3}=0}{=}\left[\begin{array}{@{}c@{}}\boldsymbol{{\Psi}}_{0}^{+},-\boldsymbol{{\Psi}}_{1}^{-}\end{array}\right]^{-1}\boldsymbol{{\Psi}}_{1}^{+}\boldsymbol{b}_{1}^{+}-\left[\begin{array}{@{}c@{}}\boldsymbol{{\Psi}}_{0}^{+},-\boldsymbol{{\Psi}}_{1}^{-}\end{array}\right]^{-1}\boldsymbol{{\Psi}}_{Reference} & & \displaystyle\end{eqnarray}$ (6) where the reference field 𝛹_Reference is due to the presence of any 3D air-core eddy current probe over a planar conducting slab. The reader can refer to [5]. This reference field is computed in the air region, satisfying the condition $x_{0}^{3}=0$ on the first surface z = a ₀(x, y). The planar interface of the conducting slab concides with the lowest point of the first interface a ₀(x, y). At the end of the computations, the unknown amplitude $\mathbf{b}_{0}^{+}$ is determined. Now, in order to take into account a more general geometry where the interfaces are non necessary parallel, another transformation, the S-matrix layer must be introduced since two non-orthogonal coordinates systems may coexist in some layer. The principle of this transformation is described in [4] and it is extended here to the 3D case. In summary, considering that two coordinates systems: the upper $x_{p}^{3}=z-a_{p}(x,y)$ and the lower $x_{q}^{3}=z-a_{q}(x,y)$ may coexist in one layer, the amplitudes of the forward/backward waves are given by $\mathbf{b}_{p}^{+}=\widetilde{S}_{q\rightarrow p}∼\mathbf{b}_{q}^{+}$ and $\mathbf{b}_{q}^{-}=\widetilde{S}_{p\rightarrow q}∼\mathbf{b}_{p}^{-}$ . These transformations are incorporated into the previous global algorithm in order to be able to compute the global S-matrix S (N,1) for a general stratified media characterized by N layers. Then, after translating the Auld’s formula in the new coordinates system, on the first interface, the change in the impedance of the probe is deduced.

Fig. 2.

A 3D air-core probe scanning a 3D mockup of finite thickness. The two interfaces are rather arbitrary smooth.

Fig. 3.

Variations of the real part R and the imaginary part X of the probe impedance, compared to FE data.

3. Numerical results

In order to show the numerical validity of the model, a test configuration has been chosen with two smooth interfaces. The smoothness can drastically reduce the number of harmonics which are necessary in the spectral domain. The two surfaces are characterized by some analytical function a _p(x, y) = h _p a _p(x − s _p)a _p(y), $p\in \{0,1\}$ with $a(w)=\frac{1}{2}[1+\cos (2{\pi}\frac{w}{L_{w}})]$ , ∀w ∈ [−L _w∕2, L _w∕2], w = x or y. h _p corresponds to the height of the surface since the function a (⋅) is normalized to 1. s _p represents some spatial shift along the X direction. The height of the first surface is h ₀ = −0.75 mm while the length of the profile along X is L _x = 10 mm. The profile is shifted on the left side with s ₀ = −2 mm. For the second interface, the lower one, h ₁ = +0.25 mm and the profile is shifted on the right side with s ₁ = +2 mm. The functions a ₀(y) and a ₁(y) remains symmetrical with respect to the Y axis in order to be able to carry on FE simulations with only a half-geometry. Let us consider L _y = 10 mm too. In the other hand, the CCM numerical model does not require any symmetry property to be able to perform some calculations. The thickness between the two surfaces is fixed at 1 mm. Other parameters concerning the probe are described in [5] and the parameters presented in Table 1 are not modified. The geometry is displayed in Fig. 2 while Fig. 3 shows the variations of the real part and the imaginary part (given in 𝛺) of the impedance when the probe is scanning the sample. FEM data have been obtained thanks to a commercial version of a Finite Element software. The agreement between simulated data is quantified by the two separated errors [5] ϵ_R < 1.04% and ϵ_X < 0.05%. Moreover, the longitudinal component H ₃ which is a critical component in this configuration is also represented in Fig. 4 while Fig. 5 displays two sliceviews. The agreement between simulated data is quite satisfying since these results have been obtained with some parameters M = 23 and N = 21. The convergence is this method is reached thanks to the smoothness of the profiles. An increasing number of modes requires an improvement of the algorithms in order to be able to address other severe geometric profiles.

Fig. 4.

H ₃(z = a ₀(x, y)), the probe is fixed at x = 0.

Fig. 5.

A sliceview of the H ₃ component.

4. Conclusions and future works

This contribution extends to the full 3D case previous significative works concerning the development of a semi-analytical model for simulating some ECNDT configurations when a 3D EC probe is scanning a 2D conducting stratified media. The numerical model can address now a set of multiples interfaces of complex shape. Due to numerical and memory space limitations, the profile of the surfaces is rather smooth in order to reduce the number of harmonics along the two covariant axis x ¹ and x ² for each interface. First improvements concern the implementation of dedicated algorithms to enlarge the sizes of the matrices to deal with and other works lie on the parallelization of the S-matrix algorithm. The time computation is relatively reasonable (about a few minutes for all the probe positions) in comparison to other simulating methods due to the complexity of the geometry but some attention must be payed to the problem of solving the eigenvalue problem due to the size of the matrix. Moreover, the convergence of the numerical method can be strongly difficult to reach for full 3D problems. On the other hand, there are no really practical limitations for 2D workpieces on usual computers.

References

Chandezon

, A new theoretical method for diffraction gratings and its numerical application, Journal of Optics 11 (1980), 235.

Caire

Prémel

and Granet

, Semi-analytical computation of a quasi-static field induced by an eddy current probe in a conductor with a rough surface, The European Physical Journal Applied Physics 64(2) (2013), 24511, 9 p.

Caire

Prémel

and Granet

, Semi-analytical computation of a quasi-static field induced by a 3D eddy current probe scanning a 2D layered conductor with parallel rough interfaces, International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 27(3) (2014), 600–613.

Caire

Prémel

and Granet

, Fast computation of the fields diffracted by a multi-layered conductor with non-parallel rough interfaces. Application to eddy-current non-destructive testing simulation, IEEE Transaction on Magnetics 51(3) (2015), 6200604.

Prémel

and Granet

, Development of the curvilinear coordinate method for the computation of quasi-static fields induced by an eddy current probe scanning a 3d conductor of complex shape characterized by an arbitrary 2D surface, International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 31 (2018), pp. e2219–n/a.

Cotter

N.P.K.

, Scattering-matrix approach to multilayer diffraction, J. Opt. Soc. Am. 12 (1995), 1097–1103.

, Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings, J. Opt. Soc. Am. A 13(5) (1996).

Granet

, Analysis of diffraction by surface-relief crossed gratings with use of the Chandezon method: application to multilayer crossed gratings, J. Opt. Soc. Am. A 15(5) (1998), 1121–1131.

Post

, Formal Structure of Electromagnetics: General Covariance and Electromagnetics, North-Holland, Amsterdam, 1962.