# Matching curvilinear coordinates for the computation of the distribution of eddy currents in a cylindrical tube described by an arbitrary longitudinal internal/external profile

Denis Prémel, G. Granet, F. Caire

► To cite this version:

Denis Prémel, G. Granet, F. Caire. Matching curvilinear coordinates for the computation of the distribution of eddy currents in a cylindrical tube described by an arbitrary longitudinal internal/external profile. European Physical Journal: Applied Physics, 2016, 73, pp.10903. 10.1051/epjap/2015150364. cea-01845385

HAL Id: cea-01845385

https://cea.hal.science/cea-01845385

Submitted on 17 Nov 2023

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.
# Matching curvilinear coordinates for the computation of the distribution of eddy currents in a cylindrical tube described by an arbitrary longitudinal internal/external profile

Denis Prémel<fnref n="1" />, Gérard Granet<fnref n="2" /> and François Caire<fnref n="1" />(*) a

Received: date / Revised version: date

**Abstract.** This paper addresses the problem of the development of a fast numerical model for the computation of the electromagnetic field in the quasi-static regime. A semi-analytical approach is achieved in order to obtain fast simulation tools dedicated to the simulation of complex configurations of Eddy Current (EC) inspection: an EC probe scans a conducting cylindrical tube of complex shape and with varying electromagnetic properties. In this paper, though the cylindrical tube to be inspected presents some angular symmetry, the complexity of the configuration lies in the arbitrary shape of the tube's walls and the dependence of constitutive parameters of the material on the radial or axial coordinates. In this context, no modal approach can be used to expand the components of the fields. This problem is overcome by using a pseudospectral Fourier method. Maxwell's equations are written in a covariant form in order to translate boundary conditions at each interface in an analytical form. Tangential components of the fields with respect to the boundary surface are expanded by using a high order Chebychev polynomials. Spatial derivatives of the components of the fields with respect to the radial coordinate are thus approximated at some collocation points. Numerical validations are discussed in order to show the efficiency of the proposed numerical model.

## Introduction

Most of semi-analytical models implemented into Eddy current modules of CIVA (a software platform dedicated to the simulation of different techniques of Non Destructive Testing (NDT) (Ultrasound, Eddy Current, Radiography, ...)) are based

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

<footnote n="2">² Institut Pascal, Université Blaise Pascal, France</footnote>

<footnote>a F. Caire has transferred his activities at Areva NDE Group Solutions : 4, Rue Thomas Dumorey, Châlon-sur-Saône, 71109 Cedex-France.</footnote>
on Volume Integral Equations and the Green's dyads formalism. This approach is very effective since it leads to accurate and fast numerical models for the Eddy Current virtual inspection of conductive pieces of canonical geometries. Indeed analytical expressions of the dyads are well known for planar or cylindrical geometries and for stratified media. The advantages of such a formalism lies in the fact that boundary conditions are forwardly included in dyadic operators and no mesh is required except for taking into account a volumetric defect in the workpiece. In order to address more complex geometries, it is usual to use purely numerical methods such as the standard Finite Element Method (FEM) despite of a strong degradation of performances. So, in order to avoid any mesh of the complex geometry, the idea consists in introducing a change of coordinates that permits to write easily, at each interface, boundary conditions in an analytical form.

The Curvilinear Coordinate Method or the so called “C Method” is widely used in the community of applied Optics since it has been proposed by Chandezon et al [1,2] in the last two decades. This differential method remains one of the most efficient method for the computation of scattered fields by gratings enlightened by a plane wave [3]. Besides, a great number of ECNDT complex configurations have been solved by using the truncated region eigenfunction expansion (TREE) method. This method is based, in particular, on the second order vector potential formalism which provides, for canonical coordinates systems, advantageous relationships between two scalar potentials and the electromagnetic fields. This is the reason why recent developments have been performed in order to use the co-variant formalism of Maxwell's equations [4] for extending the

second order vector potential formalism in order to obtain any modal expansions of the two scalar potentials in the quasi-static regime [5]. Thus, more complex ECNDT geometries have been efficiently tackled [6–10] thanks to this innovative approach in the domain of eddy current modelling. Unfortunately, even if this formalism can provide quite good solutions for simulating some 2.5D planar geometries with different interfaces of complex shapes (a 3D eddy current probe scanning a conductive piece with a surface varying along a single direction), it fails in the context of this paper for cylindrical surfaces since the profile of the interface may vary arbitrary according to the axial coordinate $z$. In this case, even if the constitutive parameters are constant by considering a homogeneous material, it is not possible to establish any eigenvalue problem to find a modal expansion of the electromagnetic field in the tube’s wall. Moreover, some variations of constitutive properties of the non-homogeneous material are included in this paper. In this case, another approach must be considered in order to obtain another expansion of the fields.

Pseudo-spectral methods are some kind of spectral methods to solve partial differential equations and have been successfully applied for a variety of problems [11–14] and also for analyzing diffraction gratings [15, 16]. This approach has been also successfully applied [10] for simulating EC inspection of planar pieces presenting local complex layers with non-homogeneous properties. However, in the following of the paper, we will see that we will have to introduce some additional fictitious surfaces in order to be able to write modal solutions in air. This differs exactly from the planar case [10]. The Pseudospectral modal method applied in this paper consists in writ-
ing the covariant form of Maxwell's equations [17] with respect with a specific radial coordinate which can fit the variations of the tube's wall. The radial derivatives are calculated using orthogonal functions of high order such as Chebyshev polynomials. Differential matrices are introduced for the differentiation of fields along the radial direction while Fourier basis functions are used for the longitudinal axis. As a consequence, the proposed numerical model is valid for continuous and smooth profiles of the geometry of the tube's wall and for continuous variations of the fields along the radial directions. An abrupt change of the constitutive properties of the material requires to add another layer and to introduce some boundary conditions at the interface. In this paper, some various ECNDE configurations are simulated and numerical results provided by the proposed model are compared to other numerical finite element data in order to confirm numerically the validity of the proposed semi-analytical model.

This paper is organized as follows. First of all, the principle of coordinate matching is presented in order to show that two different coordinate systems must be merged to obtain only one coordinate system in the wall's tube. This last non orthogonal coordinate system leads to introduce the covariant form of Maxwell's equations since a novel metric is used. Then, Maxwell's equations are translated in a differential system applying the natural tangential components of the electromagnetic field with respect to each interface. In the cylindrical coordinate system, since we can define two cylindrical fictitious surfaces, some modal solutions are available in the Fourier space. These modal solutions are combined with pseudospectral modal expansions of the components of the fields in the Fourier space.

The numerical strategy gives some details concerning the discretization of the differential system and how to obtain *in fine* the values of the electromagnetic fields observed on the internal fictitious tube's wall. The response of the EC probe can be deduced from the Auld's formula. In the last section, several numerical results are presented, discussed and compared with those provided by a commercial Finite Element software.

# 1 Formalism

## 1.1 Coordinate Matching and change of coordinate

Fig. 1 displays some kind of example we wish to address in this paper. A cylindrical pipe is inspected by a 2D Eddy Current (EC) probe scanning its longitudinal profile. The complexity in the geometry lies in the arbitrary shape of this profile and moreover, the constitutive electromagnetic parameters of the material may vary according to the radial direction. Let us denote by $\sigma(r)$ and $\mu_r(r)$ respectively the electrical conductivity and the relative magnetic permeability of the material.

<figure><img src="image_1.png" /><figcaption>**Fig. 1.** General configuration considered: a coil inspecting a conductive tube presenting a surface varying along the *z* direction.</figcaption></figure>
In this paper, the cylindrical surface presents a rotational symmetry, the variable $\theta$ is thus omitted in what follows. The coordinate matching consists in connecting two surfaces (interfaces) of coordinates by an analytical transition function denoted by $h(u)$. Fig. 2 presents a slice view of the geometry: the tube's wall is delimited by two surfaces represented by two profiles $a_1(z)$ and $a_2(z)$. In order to be able to write easily

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** A slice view: the baseline of the internal wall is $Rp_1 = 14.5$ mm and those of the external wall is $Rp_2 = 17.5$ mm in this example.</figcaption></figure>

boundary conditions on the two walls on the tube, from the usual cylindrical coordinates system $(r, \theta, z)$, let us introduce a change of coordinates so that the interfaces are aligned along the coordinates lines.

$$ \left\{ \begin{array}{l} r = a(u, z) \\ \theta = v, \\ z = w, \end{array} \right. \qquad (1) $$

where the scalar function $a(u, z)$ depends on a transition function $h(u)$ so that:

$$ a(u, z) = a_1(z) [1 - h(u)] + a_2(z) h(u) \qquad (2) $$

For $h(u) = 0$, the surface fits the internal wall represented by $a_1(z)$ while the surface fits the external wall when $h(u) = 1$.

Between two extremal values $R_{Min}$ and $R_{Max}$, let us consider a set of points $R_{Min} \le u_j \le R_{Max}, j = 0, \dots, N$ and one can choose any transition function:

$$ h(u) = \left( \frac{u - Rp_1}{Rp_2 - Rp_1} \right)^n \qquad (3) $$

with $n$ is an integer. Without any perturbation of the geometry, the internal surface is carried out by the inner radius of the tube denoted by $Rp_1$ while $Rp_2$ stands for the external radius of the tube. By introducing two additional fictitious cylindrical surfaces $r = R_1$ and $r = R_2$, the division of the space between the two fictitious cylindrical surfaces leads to consider three regions (I,II, III) (see Fig. 3 in the following of the text). For each region, a coordinate transformation is chosen:

$$ a(u, z) = \begin{cases} R_1 [1 - h(u)] + a_1(z) h(u), & (I): R_1 \le r \le a_1(z), \\ a_1(z) [1 - h(u)] + a_2(z) h(u), & (II): a_1(z) \le r \le a_2(z), \\ a_2(z) [1 - h(u)] + R_2 h(u), & (III): a_2(z) \le r \le R_2. \end{cases} \qquad (4) $$

Therefore, each change of coordinates leads us to introduce a novel metric tensor:

$$ g_{ij} = \begin{bmatrix} \dot{a}_u^2 & 0 & \dot{a}_u \dot{a}_z \\ 0 & r^2 & 0 \\ \dot{a}_u \dot{a}_z & 0 & 1 + \dot{a}_z^2 \end{bmatrix} \qquad (5) $$

where $\dot{a}_u$ stands for the derivative of $a(u, z)$ with respect to $u$: $\dot{a}_u = \frac{\partial a(u, z)}{\partial u}$ and $\dot{a}_z = \frac{\partial a(u, z)}{\partial z}$. The determinant is denoted by $g = |g_{ij}| = a \dot{a}_u$. By inverting the matrix of covariants components $g_{ij}$, The elements of the metric tensor are
deduced:

$$ \sqrt{g} [g^{ij}] = \begin{bmatrix} \frac{a}{\dot{a}_u} (1 + \dot{a}_z^2) & 0 & -\dot{a}_z a \\ 0 & \frac{\dot{a}_u}{a} & 0 \\ -a \dot{a}_z & 0 & \dot{a}_u a \end{bmatrix} = \begin{bmatrix} g^{11} & 0 & g^{13} \\ 0 & g^{22} & 0 \\ g^{13} & 0 & g^{33} \end{bmatrix} \quad (6) $$

Moreover, taking into account the jacobian matrix J of the coordinate transformation, it is usual to find the relationship between the components of the field in the cylindrical coordinates system with the components of the same field in the transformed coordinates system, so:

$$ \begin{bmatrix} \mathbf{t}_u \\ \mathbf{t}_v \\ \mathbf{t}_w \end{bmatrix} = J \begin{bmatrix} \mathbf{u}_r \\ r \mathbf{u}_\theta \\ \mathbf{u}_z \end{bmatrix} = \begin{bmatrix} \dot{a}_u \mathbf{u}_r \\ a(u, z) \mathbf{u}_\theta \\ \dot{a}_z \mathbf{u}_r + \mathbf{u}_z \end{bmatrix} $$

and

$$ \begin{bmatrix} \mathbf{t}^u \\ \mathbf{t}^v \\ \mathbf{t}^w \end{bmatrix} = g^{ij} \begin{bmatrix} \mathbf{t}_u \\ \mathbf{t}_v \\ \mathbf{t}_w \end{bmatrix} = \begin{bmatrix} (g^{11}\mathbf{t}_u + g^{13}\mathbf{t}_z) \\ g^{22}\mathbf{t}_\theta \\ (g^{13}\mathbf{t}_u + g^{33}\mathbf{t}_z) \end{bmatrix} $$

one can deduce the covariant components and the contravariant components by:

$$ \mathbf{E} = E_u \mathbf{t}^u + E_v \mathbf{t}^v + E_w \mathbf{t}^w \quad (7) $$

$$ = E^u \mathbf{t}_u + E^v \mathbf{t}_v + E^w \mathbf{t}_w \quad (8) $$

Then, Maxwell's equations must be used in its covariant form since we have to use a novel metric depending on the change of coordinates (1).

## 1.2 Covariant Maxwell's equations

In each homogeneous region and without any source term, the so-called covariant Maxwell's equations [4] establish a relationship between the covariant components of the electromagnetic field E and H and its contravariant components. By considering the angular frequency $\omega$ in the harmonic regime $\exp(-i\omega t)$, these equations can be written as:

$$ \nabla \times Z_0 \mathbf{H} = -ik\epsilon_r \sqrt{g} g^{ij} \mathbf{E} = -ik\epsilon_r^{ij} \mathbf{E} \quad (9) $$

$$ \nabla \times \mathbf{E} = ik\mu_r \sqrt{g} g^{ij} Z_0 \mathbf{H} = ik\mu_r^{ij} Z_0 \mathbf{H} $$

where electromagnetic properties of the medium are defined by the relative dielectric permittivity $\epsilon_r$ and the relative magnetic permeability $\mu_r$. Since the medium is conductive, the relative dielectric permittivity includes a loss term due to the high conductivity $\sigma$ of the medium $\epsilon_r = 1 + \frac{i\sigma}{\omega\epsilon_0}$. The dielectric permittivity and the magnetic permeability in the vacuum are also related to the characteristic impedance in the vacuum $Z_0 = \sqrt{\frac{\mu_0}{\epsilon_0}}$. Finally, k is the wavenumber given by $k = \frac{\omega}{\sqrt{\mu_0\epsilon_0}}$. This form of Maxwell equation is usually used in scattering gratings theory rather than for the quasi-static regime. Nevertheless, these notation is very useful in order to be consistent with a great number of papers in the literature associated to gratings theory.

One can note that two matrices $\epsilon_r^{ij}$ and $\mu_r^{ij}$ are defined in order to be used in the following of the paper. It is also important to note that the differential operator $\nabla \times$ is applied in the same manner as that used for the cartesian coordinates system:

$$ \nabla \times Z_0 \mathbf{H} = \begin{bmatrix} 0 & -\partial_w & \partial_v \\ \partial_w & 0 & -\partial_u \\ -\partial_v & \partial_u & 0 \end{bmatrix} \begin{bmatrix} Z_0 H_u \\ Z_0 H_v \\ Z_0 H_w \end{bmatrix} \quad (10) $$

In the particular case of a 2D Eddy current probe, due to the axial symmetry, the components $E_\theta$, $H_u$ and $H_z$ are considered. By some calculations, it is easy to obtain a differential
system which relates the components of interest:

$$ \partial_u \begin{bmatrix} E_v \\ Z_0 H_w \end{bmatrix} = ik \begin{bmatrix} L_{11} & L_{12} \\ L_{21} & L_{22} \end{bmatrix} \begin{bmatrix} E_v \\ Z_0 H_w \end{bmatrix} \quad (11) $$

with:

$$ \begin{align} L_{11} &= \frac{i}{k} \mu_r^{31} [\mu_r^{11}]^{-1} \partial_w \nonumber \\ L_{12} &= \mu_r^{33} - \mu_r^{31} [\mu_r^{11}]^{-1} \mu_r^{13} \tag{12} \\ L_{21} &= \varepsilon_r^{22} + \frac{1}{k^2} \partial_w [\mu_r^{11}]^{-1} \partial_w \nonumber \\ L_{22} &= \frac{i}{k} \partial_w [\mu_r^{11}]^{-1} \mu_r^{13} \nonumber \end{align} $$

The differential system (11) must be satisfied in regions I, II and III. In the rest of the space, especially in air, a classical solution of this differential system can be obtained and a representation of the components of the field can be written by using modal solutions in the cylindrical coordinates system.

## 1.3 Modal solutions in the cylindrical coordinates system

Let us denote by $k_{air}^2$ the wave number in air:

$$ k_{air}^2 = k^2 \mu_0 \epsilon_0 \left( 1 + \frac{i \sigma_1}{\omega \epsilon_0} \right). \quad (13) $$

for numerical reasons, a small conductivity of $\sigma_1 = 10^{-10}$ is assumed in air. For a straight cylindrical tube, $r = a(u, z) = u$ and we obviously obtain: $\mu_r^{12} = 0, \mu_r^{11} = \mu_r^{33} = \mu_r r$ and $\varepsilon_r^{22} = \frac{1}{r}\varepsilon_r$ and finally:

$$ \begin{align} L_{11} &= L_{22} = 0, \quad L_{12} = \mu_r^{33} = \mu_r r \nonumber \\ L_{21} &= \frac{1}{r} \frac{1}{k^2} C_{11}, \quad \text{with } C_{11} = \left[ \frac{\partial}{\partial z} \mu_r^{-1} \frac{\partial}{\partial z} + k^2 \varepsilon_r \right] \tag{14} \end{align} $$

In this particular case, it is easy to find the eigen-solutions of the differential system (11). Let us denote by $\mathcal{L}$ the operator so that $\partial_r V = ik \mathcal{L} V$ with $V = [E_v, Z_0 H_z]^T$. One obtains:

$$ \partial_r \mathcal{L} = \begin{bmatrix} 0 & \mu_r \\ -\frac{1}{r^2} \frac{C_{11}}{k^2} & 0 \end{bmatrix} \quad (15) $$

Moreover, one can shows that: $\partial_r^2 V = ik\partial_r \mathcal{L} - k^2 \mathcal{L}^2 V$ and:

$$ \partial_r^2 V + \mu_r C_{11} V + \frac{1}{r} \partial_r V = \begin{bmatrix} 0 & ik\mu_r \\ 0 & 0 \end{bmatrix} V \quad (16) $$

By considering a Fourier domain along the Z axis of the tube, let us denote par $\alpha$ the spatial frequency corresponding to the z coordinate. In this case, by considering that $\lambda^2 = \alpha^2 - k_{air}^2$, one can obtain $\mu_r C_{11} = \lambda^2$. The second line of equation (16) corresponds to a modified Bessel equation of order 0. Two kind of solutions are thus given by $\hat{H}_z \propto I_0(\lambda r)$ or $\hat{H}_z \propto K_0(\lambda r)$ according to radiation condition at infinity. By taking the second line of the system (11), the component of the electrical field in the Fourier domain is given by $\hat{E}_\theta \propto \frac{1}{\lambda} I'_0(\lambda r)$ or $\hat{E}_\theta \propto \frac{1}{\lambda} K'_0(\lambda r)$.

Finally, the expressions of the components of the electromagnetic field can thus be deduced in air for any profile $a(u, z)$. In the interior of the first fictitious cylindrical surface, i.e. $0 \le r \le R_1$, so:

$$ \hat{E}_{\theta} = \hat{E}_{\theta}^{(Ref)}(r, \alpha) + C(\alpha)I_{1}(\lambda r) \quad (17) $$

$$ Z_0 \hat{H}_z = Z_0 \hat{H}_z^{(Ref)}(r, \alpha) + C(\alpha) \lambda I_0(\lambda r) \quad (18) $$

where $\hat{E}_{\theta}^{(Ref)}(r, \alpha)$ and $\hat{H}_z^{(Ref)}(r, \alpha)$ stand for the components of the field transformed in the Fourier domain and in a reference configuration. What we add in the formulas (17) and (18) corresponds to the perturbation of the components of the electromagnetic field with respect to this reference case. In order to reduce the spectral range of these components, the case of a straight cylindrical tube has been chosen for the reference configuration. Then, out of the exterior of the second fictitious
cylindrical surface, i.e. $r \ge R_2$,

$$ \hat{E}_{\theta} = D(\alpha)K_1(\lambda r) \quad (19) $$

$$ Z_0 \hat{H}_z = -D(\alpha) \lambda K_0(\lambda r) \quad (20) $$

The constants $C(\alpha)$ and $D(\alpha)$ are determined in what's follows by applying boundary conditions at the two fictitious cylindrical surfaces. These conditions are then included in the numerical strategy.

## 1.4 Numerical strategy

In order to solve numerically the differential system (11), we have to introduce some discrete operators to compute a satisfying approximation of the solutions.

### 1.4.1 Discretization of the Fourier space

In the previous section, we have already seen that along the longitudinal axis, the z axis, the solutions are expanded on a Fourier basis. In order to complete the description, we have to discretize firstly the Fourier space. We choose to expand any Fourier transform $\hat{F}(u, \alpha)$ on sinc functions $[\text{sinc}(x) = \sin(\pi x)] / (\pi x)$

since the orthogonality between two sinc functions gives quite advantageous properties for approximating the convolution product of two expanded functions in the Fourier space [18, 19]:

$$ \begin{aligned} a(u, z) F(u, z) &\xrightarrow{\text{FT}} \hat{H}(u, \alpha) = (\hat{a} * \hat{F})(u, \alpha) \\ &\approx \frac{\Delta_{\alpha}}{2\pi} \sum_{m=-\infty}^{m=+\infty} \hat{a}(u, \alpha - \alpha_m) \hat{F}(u, \gamma = \alpha_m) \end{aligned} $$

with $\alpha_m = m\Delta\alpha, m \in \mathbb{Z}$ and $\Delta\alpha$ is the spectral resolution. This is fixed, in practice, by a cut-off integer denoted by $M$ and a maximum value $\alpha_M$ of $\alpha$ depending on the support of

the function $\hat{F}(u, \alpha)$: $\alpha_M = \frac{1}{2} \text{supp}(\hat{F}(u, \alpha))$ in the spectral domain, $M$ stands for the truncation order which gives:

$$ \Delta_{\alpha} = \frac{\alpha_M}{2M + 1}. \quad (21) $$

Finally, by using the point matching method at discrete values $\alpha = \alpha_m$, one obtains:

$$ \hat{H}(u, \alpha) \approx \sum_{p=-M}^{p=+M} h_p \text{sinc} \left( \frac{\alpha - p\Delta_{\alpha}}{\Delta_{\alpha}} \right) $$

with:

$$ h_p = \frac{\Delta_{\alpha}}{2\pi} \sum_{m=-M}^{+M} \hat{a}_p(u) \hat{F}_{p-m}(u) \quad (22) $$

and finally the convolution product $(\hat{a} * \hat{F})(u, \alpha)$ results in the discrete space by a matrix product of the form $[\mathbf{A}] \hat{F}$ where $\hat{F}$ is a column vector such as:

$$ \hat{\mathbf{F}} = [\hat{F}(u, \alpha_1), \hat{F}(u, \alpha_2), \dots, \hat{F}(u, \alpha_{2M+1})]^T \quad (23) $$

and $[\mathbf{A}]$ is a truncated convolution matrix generated by the column vector:

$$ \hat{\mathbf{a}} = [\hat{a}(u, \alpha_1), \hat{a}(u, \alpha_2), \dots, \hat{a}(u, \alpha_{2M+1})]^T. $$

Moreover, the derivative operator $\partial_z$ consists operationally in multiplying the approximated solution by a diagonal matrix:

$$ \partial_z F(u, z) \xrightarrow{\text{FT}} i\alpha \hat{F}(u, \alpha) \approx i[\alpha] \hat{\mathbf{F}} $$

with a diagonal matrix $[\alpha] = \text{diag}(\alpha_1, \alpha_2, \dots, \alpha_{2M+1})$ and the same column vector $\hat{\mathbf{F}}$ of (23).
Considering all these notations, it is possible to write the discrete operators of $\mu_r^{ij}$:

$$ 
\begin{aligned}
\mu_r^{13} &= -[\hat{\mathbf{a}}_z][\hat{\mathbf{a}}] \\
\mu_r^{33} &= [\hat{\mathbf{a}}_u][\hat{\mathbf{a}}] \\
\mu_r^{22} &= \begin{bmatrix} \hat{\mathbf{a}}_u \\ \hat{\mathbf{a}} \end{bmatrix} \\
[\mu_r^{11}]^{-1} &= \mu_r^{22} \begin{bmatrix} \hat{\mathbf{h}} \\ \hat{\mathbf{a}} \end{bmatrix} \\
\text{with } \hat{\mathbf{h}}(u, z) &= \frac{1}{1 + \hat{a}_z \hat{a}_z}
\end{aligned}
\quad (24) 
$$

In these last equations, the notation $\hat{\mathbf{F}}/\hat{\mathbf{G}}$ stands for element-by-element division and $\hat{\mathbf{F}} \cdot \hat{\mathbf{H}}$ element-by-element multiplication. Concerning the computation of last term, the values of the function $\hat{h}(u, z) = \frac{1}{1+\hat{a}_z^2}$ are firstly stored in a column vector $\hat{\mathbf{h}}$ before computing the convolution matrix $[\hat{\mathbf{h}}]$.

All matrix elements $L_{ij}$ can be obviously translated in its discrete form $\mathbf{L}_{ij}$ from (12):

$$ 
\begin{aligned}
\mathbf{L}_{11} &= -\frac{1}{k} \mu_r^{31} [\mu_r^{11}]^{-1} [\alpha] \\
\mathbf{L}_{12} &= \mu_r^{33} - \mu_r^{31} [\mu_r^{11}]^{-1} \mu_r^{13} \\
\mathbf{L}_{21} &= \varepsilon_r^{22} - \frac{1}{k^2} [\alpha] [\mu_r^{11}]^{-1} [\alpha] \\
\mathbf{L}_{22} &= -\frac{1}{k} [\alpha] [\mu_r^{11}]^{-1} \mu_r^{13}
\end{aligned}
\quad (25) 
$$

### 1.4.2 radial direction: Pseudo-spectral approximation

In this subsection, a pseudo-spectral method is described in the case where the components of the electromagnetic field are not discontinuous in any layer. The goal is to obtain an approximation of the quantity $F(u, \alpha)$ for a finite number $N$ of values of $u = u_0, u_2, \dots, u_{N-1}$. Let us consider a set of collocation points $\{u_n\}_{0}^{N-1}$ which are distributed like the Chebyshev points in an interval. The so called Pseudo-Spectral method (Chebyshev collocation method [20,21]) consists in approximating the first order derivative by a differentiation matrix $\mathbf{D}_N$ such that:

$$ \partial_u F(u, z) \approx \mathbf{D}_N \mathbf{F} = \mathbf{D}_N \begin{bmatrix} F(u_0, z) \\ F(u_1, z) \\ \vdots \\ F(u_{N-2}, z) \\ F(u_{N-1}, z) \end{bmatrix} \quad (26) $$

The notation is just a little bit modified in order to define the field $[F]_1^{N-2}$ containing the values of $F(u, z)$ at the discretization points $u_1, \dots, u_{N-2}$ inside the interval $[u_0, u_{N-1}]$. Likewise, $[F]_0$ and $[F]_1^{N-1}$ means the restriction of the values of the field $F(u, z)$ on the limits $u_0$ and $u_{N-1}$ respectively.

$$ \partial_u F(u, z) \approx \mathbf{D}_N \mathbf{F} = \begin{bmatrix} d_1 \mathbf{d}_2 d_3 \\ \mathbf{d}_4 \mathbf{d}_5 \mathbf{d}_6 \\ d_7 \mathbf{d}_8 d_9 \end{bmatrix} \begin{bmatrix} [F]_0 \\ [F]_1^{N-2} \\ [F]_1^{N-1} \end{bmatrix} \quad (27) $$

where the scalars $d_j, j = 1, 3, 7, 9$, the vectors $\mathbf{d}_j, j = 2, 4, 6, 8$ and the matrix $\mathbf{d}_5$ results [22,15] from the partition of the matrix $\mathbf{D}_N$. In particular $\mathbf{d}_2, \mathbf{d}_8$ are row vector of length $N-2$, $\mathbf{d}_4, \mathbf{d}_6$ are column vector of length $N-2$, and $\mathbf{d}_5$ a matrix of dimensions $(N-2) \times (N-2)$. Thus, we can write an approximation of the derivative of $\hat{F}(u, \alpha)$ so that:

$$ 
\begin{aligned}
\partial_u [\hat{F}]_0 &\approx d_1 [\hat{F}]_0 + \mathbf{d}_2 [\hat{F}]_1^{N-2} + d_3 [\hat{F}]^{N-1} \\
\partial_u [\hat{F}]_1^{N-2} &\approx \mathbf{d}_4 [\hat{F}]_0 + \mathbf{d}_5 [\hat{F}]_1^{N-2} + \mathbf{d}_6 [\hat{F}]^{N-1} \\
\partial_u [\hat{F}]_1^{N-1} &\approx d_7 [\hat{F}]_0 + \mathbf{d}_8 [\hat{F}]_1^{N-2} + d_9 [\hat{F}]^{N-1}
\end{aligned}
\quad (28) 
$$

the second variable of $\hat{F}(u, \alpha)$ has been omitted in order to reduce the complexity of the written equations. Since the Fourier space is discretized by $\alpha_1, \dots, \alpha_{2M+1}$, we prefer to introduce $[\hat{\mathbf{F}}_1]_0 = \hat{\mathbf{E}}_v(r = R_1)$ and $[\hat{\mathbf{G}}_1]_0 = Z_0 \hat{\mathbf{H}}_z(r = R_1)$ which stand for the column vectors of the electrical field and the normalized magnetic field in the Fourier space in region (I) at the
first cylindrical fictitious surface $r = R_1$. The field vectors $[\hat{\mathbf{F}}_1]^{N-1}$ and $[\hat{\mathbf{G}}_1]^{N-1}$ mean the values of the same fields observed on the first profile $r = a_1(z)$ in region (I). The values of the electromagnetic field at the interior points of the segment are denoted by $[\hat{\mathbf{F}}_1]_1^{N-2}$ and $[\hat{\mathbf{G}}_1]_1^{N-2}$. In order to achieve the description of all notations which will be necessary in what's follows, let us consider the matrices $\mathbf{D}_q$, $q = 4, 5, 6$ which are computed from the differential operator $\mathbf{d}_q$, $q = 4, 5, 6$ in equation (27) by using a suitable Kronecker product of matrices $\mathbf{D}_q^{(1)} = \mathbf{d}_q^{(1)} \otimes \mathbf{I}_d$ where $\mathbf{I}_d$ denotes the identity matrix of dimension $(2M+1) \times (2M+1)$. An approximation of the first derivative of the field $\hat{F}(u, z)$ becomes:

$$ \partial_u \hat{F}(u, z) \approx \begin{bmatrix} \mathbf{D}_1 & \mathbf{D}_2 & \mathbf{D}_3 \\ \mathbf{D}_4 & \mathbf{D}_5 & \mathbf{D}_6 \\ \mathbf{D}_7 & \mathbf{D}_8 & \mathbf{D}_9 \end{bmatrix} \begin{bmatrix} [\hat{\mathbf{F}}_1]_0 \\ [\hat{\mathbf{F}}_1]_1^{N-2} \\ [\hat{\mathbf{F}}_1]^{N-1} \end{bmatrix} \quad (29) $$

the same notation is also applied for $\hat{G}(u, z)$. The next step consists in approximating the differential equation (11) in each region. It is important to note that the domain between the two fictitious cylindrical surfaces are divided into three regions as it is shown in Fig. 3. Since the discretization of each layer depends on the two extremal values $R_{Min}$ and $R_{Max}$ chosen for each region, it is necessary to write three systems of equations knowing that the differential operators $D_N^{(p)}$, defined for each layer $p = 1, 2, 3$, depend on the thickness $(R_{Max} - R_{Min})$ of each layer (I,II,III).

The first system is obtained by writing that the differential system (11) must be satisfied at the interior discretization

<figure><img src="image_3.png" /><figcaption>Fig. 3. The three regions are discretized with a set of $N+1$ collocation points.</figcaption></figure>

points of the first segment $R_1 \le r \le \max(a_1(z))$:

$$ \begin{align} \mathbf{D}_4^{(1)} [\hat{\mathbf{F}}_1]_0 + \mathbf{D}_5^{(1)} [\hat{\mathbf{F}}_1]_1^{N-2} + \mathbf{D}_6^{(1)} [\hat{\mathbf{F}}_1]_1^{N-1} &= \nonumber \\ & \quad ik [\mathbf{L}_{11}^{(1)} \hat{\mathbf{F}}_1]_1^{N-2} + ik [\mathbf{L}_{12}^{(1)} \hat{\mathbf{G}}_1]_1^{N-2} \tag{30} \\ \mathbf{D}_4^{(1)} [\hat{\mathbf{G}}_1]_0 + \mathbf{D}_5^{(1)} [\hat{\mathbf{G}}_1]_1^{N-2} + \mathbf{D}_6^{(1)} [\hat{\mathbf{G}}_1]_1^{N-2} &= \nonumber \\ & \quad ik [\mathbf{L}_1^{21} \hat{\mathbf{F}}_1]_1^{N-2} + ik \mathbf{L}_{22}^{(1)} [\hat{\mathbf{G}}_1]_1^{N-2} \nonumber \end{align} $$

Let us to introduce now the partitioned matrix of the matrix blocks $\mathbf{L}_{ij}^{(I)}$ coming from the differential operator $\mathcal{L}$ in the first region (I):

$$ \mathbf{L}_{ij}^{(1)} = \begin{bmatrix} L_{1ij}^{(1)} & L_{2ij}^{(1)} & L_{3ij}^{(1)} \\ L_{4ij}^{(1)} & L_{5ij}^{(1)} & L_{6ij}^{(1)} \\ L_{7ij}^{(1)} & L_{8ij}^{(1)} & L_{9ij}^{(1)} \end{bmatrix} \quad (31) $$

Equation (30) becomes:

$$ \begin{bmatrix} \mathbf{M}_5^{(1)} \\ \mathbf{N}_6^{(1)} \end{bmatrix} \begin{bmatrix} [\hat{\mathbf{F}}_1]_1^{N-2} \\ [\hat{\mathbf{G}}_1]_1^{N-2} \\ [\hat{\mathbf{F}}_1]^{N-1} \\ [\hat{\mathbf{G}}_1]^{N-1} \end{bmatrix} = -\mathbf{M}_4^{(1)} \begin{bmatrix} [\hat{\mathbf{F}}_1]_0 \\ [\hat{\mathbf{G}}_1]_0 \end{bmatrix} \quad (32) $$
with

$$ M_k^{(1)} = \begin{bmatrix} (\mathbf{D}_k^{(1)} - ik\mathbf{L}k_{11}^{(1)}) & -ik\mathbf{L}k_{12}^{(1)} \\ -ik\mathbf{L}k_{21}^{(1)} & (\mathbf{D}_k^{(1)} - ik\mathbf{L}k_{22}^{(1)}) \end{bmatrix} \quad (33) \\ \mathbf{N}_k^{(1)} = \begin{bmatrix} \mathbf{D}_k^{(1)} & \mathbf{0} \\ \mathbf{0} & \mathbf{D}_k^{(1)} \end{bmatrix} $$

A similar system of equations can be written for the interior points of the second region (II) in the second segment $\min(a_1(z), a_2(z)) \le r \le \max(a_1(z), a_2(z))$:

$$ \left[ \mathbf{N}_4^{(2)} \quad | \quad \mathbf{M}_5^{(2)} \quad | \quad \mathbf{N}_6^{(2)} \right] \left[ \begin{bmatrix} \hat{\mathbf{F}}_2 \end{bmatrix}_0 \quad \begin{bmatrix} \hat{\mathbf{G}}_2 \end{bmatrix}_0 \quad \begin{bmatrix} \hat{\mathbf{F}}_2 \end{bmatrix}_1^{N-2} \quad \begin{bmatrix} \hat{\mathbf{G}}_2 \end{bmatrix}_1^{N-2} \quad \begin{bmatrix} \hat{\mathbf{F}}_2 \end{bmatrix}_1^{N-1} \quad \begin{bmatrix} \hat{\mathbf{G}}_2 \end{bmatrix}_1^{N-1} \right] = \mathbf{0} \quad (34) $$

where $\mathbf{L}_{ij}^{(2)}$ means the block matrix of $\mathcal{L}$ in the second region (II). The field vectors $[\hat{\mathbf{F}}_2]_0 = \hat{\mathbf{E}}_v(r = a_1(z))$ and $[\hat{\mathbf{G}}_2]_0 = Z_0 \hat{\mathbf{H}}_z(r = a_1(z))$ are observed in region (II) but at the interface $r = a_1(z)$ between the two regions (I) and (II).

Finally, a system of equations similar to (34) can be written for the interior points of the third region (III) for the third segment $\min(a_2(z)) \le r \le R_2$:

$$ \left[ \mathbf{N}_4^{(3)} \quad | \quad \mathbf{M}_5^{(3)} \right] \left[ \begin{bmatrix} \hat{\mathbf{F}}_3 \end{bmatrix}_0 \quad \begin{bmatrix} \hat{\mathbf{G}}_3 \end{bmatrix}_0 \quad \begin{bmatrix} \hat{\mathbf{F}}_3 \end{bmatrix}_1^{N-2} \quad \begin{bmatrix} \hat{\mathbf{G}}_3 \end{bmatrix}_1^{N-2} \right] = -\mathbf{M}_6^{(3)} \left[ \begin{bmatrix} \hat{\mathbf{F}}_3 \end{bmatrix}_1^{N-1} \quad \begin{bmatrix} \hat{\mathbf{G}}_3 \end{bmatrix}_1^{N-1} \right] \quad (35) $$

### 1.4.3 Introducing the source

In these previous equations (32, 35), we assumed that the field vectors at the edges of the fictitious surfaces $[\hat{\mathbf{F}}_1]_0, [\hat{\mathbf{G}}_1]_0$,

$[\hat{\mathbf{F}}_3]^{N-1}$ and $[\hat{\mathbf{G}}_3]^{N-1}$ are known. In reality, these field vectors are obviously unknown and they must be estimated through the numerical process. These vectors are linked to modal solutions previously presented in section 1.3. From the equations (17, 18, 19, 20), one can retrieve the electromagnetic field at the edges of the fictitious surfaces $r = R_1$ and $r = R_2$, so ones obtains:

$$ \begin{align} \hat{E}_{\theta}(R_1, \alpha) &= \hat{E}_{\theta}^{(\text{Ref})}(R_1, \alpha) + C(\alpha)I_1(\lambda R_1) \nonumber \\ Z_0 \hat{H}_z(R_1, \alpha) &= Z_0 \hat{H}_z^{(\text{Ref})}(R_1, \alpha) + C(\alpha)\lambda I_0(\lambda R_1) \tag{36} \\ \hat{E}_{\theta}(R_2, \alpha) &= D(\alpha)K_1(\lambda R_2) \nonumber \\ Z_0 \hat{H}_z(R_2, \alpha) &= -D(\alpha)\lambda K_0(\lambda R_2) \nonumber \end{align} $$

Let us define the unknown vectors **C** and **D** so that:

$$ \begin{align} \mathbf{C} &= [C(\alpha_1), C(\alpha_2), \dots, C(\alpha_{2M+1})]^T \tag{37} \\ \mathbf{D} &= [D(\alpha_1), D(\alpha_2), \dots, D(\alpha_{2M+1})]^T \nonumber \end{align} $$

Likewise, we define the vectors $\hat{E}_{\theta}^{(\text{Ref})}$ and $Z_0 \hat{H}_z^{(\text{Ref})}$ from $\hat{E}_{\theta}^{(\text{Ref})}(R_1, \alpha)$ and $Z_0 \hat{H}_z^{(\text{Ref})}(R_1, \alpha)$ respectively. From (36), we can find finally:

$$ \begin{align} \begin{bmatrix} [\hat{\mathbf{F}}_1]_0 \\ [\hat{\mathbf{G}}_1]_0 \end{bmatrix} &= \begin{bmatrix} \hat{\mathbf{E}}_{\theta}^{(\text{Ref})} \\ Z_0 \hat{\mathbf{H}}_z^{(\text{Ref})} \end{bmatrix} + \begin{bmatrix} \mathbf{\Psi}_{I_1} \\ \mathbf{\Psi}_{\lambda I_0} \end{bmatrix} \mathbf{C} \tag{38} \\ \begin{bmatrix} [\hat{\mathbf{F}}_3]^{N-1} \\ [\hat{\mathbf{G}}_3]^{N-1} \end{bmatrix} &= \begin{bmatrix} \mathbf{\Phi}_{K_1} \\ \mathbf{\Phi}_{\lambda K_0} \end{bmatrix} \mathbf{D} \nonumber \end{align} $$

with some diagonal matrices:

$$ \begin{align} \mathbf{\Psi}_{I_1} &= \text{Diag}(I_1(\lambda R_1)) \nonumber \\ \mathbf{\Psi}_{\lambda I_0} &= \text{Diag}(\lambda I_0(\lambda R_1)) \tag{39} \\ \mathbf{\Phi}_{K_1} &= \text{Diag}(K_1(\lambda R_2)) \nonumber \\ \mathbf{\Phi}_{\lambda K_0} &= \text{Diag}(-\lambda K_0(\lambda R_2)) \nonumber \end{align} $$

The vectors **C** and **D** must be estimated by substituting equations (38) into (32, 35).
### 1.4.4 Supplementary equations

At each interface, the continuity of the fields vectors $\hat{\mathbf{F}}$ and $\hat{\mathbf{G}}$ and also the continuity of their derivatives must be satisfied. For assembling the matrix system to solve, the field vectors $[\hat{\mathbf{F}}_q]_0$ and $[\hat{\mathbf{G}}_q]_0$ are substituted respectively by $[\hat{\mathbf{F}}_{q-1}]^N$ and $[\hat{\mathbf{G}}_{q-1}]^N$. Consequently, the number of unknowns (Matrix Blocks) is about twelve while the number of previous equations is equal to six. So, in order to obtain a square system, we introduce a set of supplementary equations which must be satisfied at different interfaces. At the top of region (I), let us consider that the two components $E_v$ and $Z_0H_w$ of the field vectors must satisfy the differential system (11):

$$ M_1^{(1)} \begin{bmatrix} [\hat{\mathbf{F}}_1]_0 \\ [\hat{\mathbf{G}}_1]_0 \end{bmatrix} + N_2^{(1)} \begin{bmatrix} [\hat{\mathbf{F}}_1]_1^{N-2} \\ [\hat{\mathbf{G}}_1]_1^{N-2} \end{bmatrix} + N_3^{(1)} \begin{bmatrix} [\hat{\mathbf{F}}_1]^{N-1} \\ [\hat{\mathbf{G}}_1]^{N-1} \end{bmatrix} = 0 \quad (40) $$

Similarly, at the bottom of region (III), considering that the two components $E_v$ and $Z_0H_w$ of the field vectors must satisfy the differential system (11), one can obtain another condition at the last interface $r = R_2$:

$$ N_7^{(3)} \begin{bmatrix} [\hat{\mathbf{F}}_3]_0 \\ [\hat{\mathbf{G}}_3]_0 \end{bmatrix} + N_8^{(3)} \begin{bmatrix} [\hat{\mathbf{F}}_3]_1^{N-2} \\ [\hat{\mathbf{G}}_3]_1^{N-2} \end{bmatrix} + M_9^{(3)} \begin{bmatrix} [\hat{\mathbf{F}}_3]^{N-1} \\ [\hat{\mathbf{G}}_3]^{N-1} \end{bmatrix} = 0 \quad (41) $$

In region (II), we can add some conditions at the top and at the bottom of the interval $a_1(z)$ and $a_2(z)$:

$$ \begin{aligned} & (\mathbf{D}_1^{(2)} - ik\mathbf{L}_1^{(2)}_1) [\hat{\mathbf{F}}_2]_0 - ik\mathbf{L}_1^{(2)}_{12} [\hat{\mathbf{G}}_2]_0 \\ & + \mathbf{D}_2^{(2)} [\hat{\mathbf{F}}_2]_1^{N-2} + \mathbf{D}_3^{(2)} [\hat{\mathbf{F}}_2]_1^{N-1} = 0 \end{aligned} \quad (42) $$

$$ \mathbf{D}_7^{(2)} [\hat{\mathbf{F}}_2]_0 + \mathbf{D}_8^{(2)} [\hat{\mathbf{F}}_2]_1^{N-2} + (\mathbf{D}_9^{(2)} - ik\mathbf{L}_9^{(3)}_1) [\hat{\mathbf{F}}_2]_1^{N-1} - ik\mathbf{L}_9^{(3)}_{12} [\hat{\mathbf{G}}_2]_1^{N-1} = 0 $$

We can now construct an unknown augmented vector and propose to assemble the matrix system to solve.

### 1.4.5 Matrix assembly

Considering all previous conditions, we can build up the global matrix in order to evaluate twelve matrix vectors. The unknown matrix blocks are concatenated in a global unknown column vector $\mathbf{X}$ while some source terms are in the second member of the matrix equation $[\mathbf{Z}] \mathbf{X} = \mathbf{Y}$ to solve.

$$ \mathbf{X} = \begin{bmatrix} \mathbf{C} \\ [\hat{\mathbf{F}}_1]_1^{N-2} \\ [\hat{\mathbf{G}}_1]_1^{N-2} \\ [\hat{\mathbf{F}}_1]^{N-1} \\ [\hat{\mathbf{G}}_1]^{N-1} \\ [\hat{\mathbf{F}}_2]_1^{N-2} \\ [\hat{\mathbf{G}}_2]_1^{N-2} \\ [\hat{\mathbf{F}}_2]^{N-1} \\ [\hat{\mathbf{G}}_2]^{N-1} \\ [\hat{\mathbf{F}}_3]_1^{N-2} \\ [\hat{\mathbf{G}}_3]_1^{N-2} \\ \mathbf{D} \end{bmatrix}, \quad \mathbf{Y} = \begin{bmatrix} -\mathbf{M}_4^{(1)} \\ [\hat{\mathbf{E}}_\theta^{(\text{Ref})}] \\ Z_0 \hat{\mathbf{H}}_z^{(\text{Ref})} \\ -\mathbf{M}_1^{(1)} \\ [\hat{\mathbf{E}}_\theta^{(\text{Ref})}] \\ Z_0 \hat{\mathbf{H}}_z^{(\text{Ref})} \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} $$

Even if the structure of the global matrix $[\mathbf{Z}]$ is sparse, we choose to build up a full matrix in order to compute only once
$$
Z_m = \begin{bmatrix}
\mathbf{M}_4^{(1)} \begin{bmatrix} \mathbf{\Psi}_{I_1} \\ \mathbf{\Psi}_{\lambda I_0} \end{bmatrix} & \mathbf{M}_5^{(1)} & \mathbf{N}_6^{(1)} & 0 & 0 & 0 \\
\mathbf{M}_1^{(1)} \begin{bmatrix} \mathbf{\Psi}_{I_1} \\ \mathbf{\Psi}_{\lambda I_0} \end{bmatrix} & \mathbf{N}_2^{(1)} & \mathbf{N}_3^{(1)} & 0 & 0 & 0 \\
0 & 0 & \mathbf{N}_4^{(2)} & \mathbf{M}_5^{(2)} & \mathbf{N}_6^{(2)} & 0 \\
0 & 0 & \begin{bmatrix} (\mathbf{D}_1^{(2)} - ik\mathbf{L}_{11}^{(2)}) - ik\mathbf{L}_{12}^{(2)} \\ \mathbf{D}_7^{(2)} [\hat{\mathbf{F}}_2]_0 \end{bmatrix} & \begin{bmatrix} \mathbf{D}_2^{(2)} & \mathbf{0} \\ \mathbf{D}_8^{(2)} & \mathbf{0} \end{bmatrix} & \begin{bmatrix} \mathbf{D}_3^{(2)} & \mathbf{0} \\ (\mathbf{D}_9^{(2)} - ik\mathbf{L}_{11}^{(3)}) - ik\mathbf{L}_{12}^{(3)} \end{bmatrix} & 0 \\
0 & 0 & 0 & \mathbf{N}_4^{(3)} & \mathbf{M}_5^{(3)} & \mathbf{M}_6^{(3)} \begin{bmatrix} \mathbf{\Phi}_{K_1} \\ \mathbf{\Phi}_{\lambda K_0} \end{bmatrix}
\end{bmatrix}, \quad (43)
$$

its inverse. Thus, even if the values of the second member of
the equation may vary according to the position of the Eddy
Current probe, it is possible to compute all the unknowns for-
wardly. This property leads to reduce the computation time
when it is required to compute the response of the EC probe
due to a longitudinal motion of the probe inside the tube. This
is the main advantage of this semi-analytical approach. More-
over, by using interpolation [23] formulas, it is possible to de-
duce the components of the electromagnetic field anywhere in
the material. Since our applications of interest concern Eddy
Current Non Destructive Testing, it is desirable to compute the
change in coil impedance of the probe. This is the object of the
following section.

## 2 Impedance of the probe

Since, we have considered a configuration of reference for com-
puting the source terms in (38), it is advantageous to compute
the variations of the impedance ΔZ of the EC probe by using
the Lorentz reciprocity principle [24]. According to the for-

mula derived by Auld and Moulder [25,26], the difference of
the coil impedance between a reference state (we have chosen
the coil in a straight conducting tube) and the state of interest
is given by:

$$
I^2 \Delta Z = \iint_{S_F} (\mathbf{E}^{(\text{ref})} \times \mathbf{H} - \mathbf{E} \times \mathbf{H}^{(\text{ref})}) \cdot \mathbf{n} \, ds, \quad (44)
$$

where *I* is the amplitude of the driving current, *S*<sub>*F*</sub> means
an arbitrary closed surface bordering the cylindrical conduc-
tive part but excluding the coil (see Fig. 1), **n** is the unit out-
ward vector, normal to *S*<sub>*F*</sub>. Since, two cylindrical fictitious
surfaces have been considered in the theoretical formalism, it
is natural to choose the internal fictitious surface defined by
*r* = *R*1, −*π* ≤ *θ* ≤ +*π*, −∞ ≤ *z* ≤ +∞ as a part of the
closed surface while the other part corresponds to a cylindrical
surface *S*<sub>∞</sub> with a radius tending to infinity. The superscript
*(ref)* refers to the fields computed along the closed surface *S*<sub>*F*</sub>
in the reference configuration and the fields **E** and **H** stand for
the total fields. Since, the contribution of the surface *S*<sub>∞</sub> tends
towards zero, the formula is simplified as follows:

$$I^2 \Delta Z = \frac{2\pi}{Z_0} \int_{-\infty}^{+\infty} \left( R_1 E_{\theta}^{(ref)} \cdot G_z - F_{\theta} \cdot Z_0 H_z^{(ref)} \right) dz \quad (45)$$

The absolute value of the impedance Z which corresponds to the response of the EC probe is finally obtained by adding the value of the reference impedance $Z^{(Ref)}$ obtained in the reference configuration and the change in coil impedance due to the variation of the geometry $\underline{Z} = Z^{(Ref)} + \Delta Z$. This Eddy Current signal is displayed in figures of results in the next section.

## 3 Numerical validation

In order to provide a numerical validation of the proposed model, a set of numerical experiments have been carried out and this section proposes to compare numerical results obtained from the numerical model, the Pseudo-Spectral Model (PSM) and other numerical data obtained by a commercial Finite Element (FE) solver (COMSOL MULTIPHYSICS [27]). For all numerical experiments, all numerical common parameters are given in Table. 1

### 3.1 First numerical experiment

In this experiment, we consider two different profiles represented in the slice view in Fig. 2. The two analytic forms are expressed as (z in [mm]):

$$a_1(z) = 14.5 + 1.5 \frac{1}{2} [1 + \cos(2\pi z/10)]$$

$$a_2(z) = 17.5 - 1.5 \frac{1}{2} [1 + \cos(2\pi z/10)] [0.141 \exp(10z/10)]$$

<table><caption>**Table 1.** Numerical parameters for the 2D EC probe, the fictitious tube and for the discretization.</caption><thead><tr><th>Parameter</th><th>symbol</th><th>value</th></tr></thead><tbody><tr><td>Inner radius</td><td>r<sub>1</sub></td><td>10 [mm]</td></tr><tr><td>Outer radius</td><td>r<sub>2</sub></td><td>13 [mm]</td></tr><tr><td>With length</td><td>h</td><td>3 [mm]</td></tr><tr><td>Number of turns</td><td>n<sub>t</sub></td><td>200</td></tr><tr><td>Inner radius of fictitious tube</td><td>R<sub>1</sub></td><td>14 [mm]</td></tr><tr><td>Outer radius of fictitious tube</td><td>R<sub>2</sub></td><td>18 [mm]</td></tr><tr><td>Conductivity of the conductive material</td><td>σ</td><td>1 MS/m</td></tr><tr><td>Truncation order (Fourier Space)</td><td>M</td><td>70</td></tr><tr><td>Number of layers (radial direction)</td><td>N</td><td>7</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Real part and imaginary part of the component $E_{\theta}$ for two lines ($r = R_1$ et $r = R_2$).</figcaption></figure>

Figure (4) and 5 display respectively the variations of the angular component $E_{\theta}$ of the electrical field **E** and the axial component $H_z$ of the magnetic field **H**. The components are observed on two lines: $r = R_1$ and $r = R_2$, the eddy current probe being centered at $Z = 0$. The operating frequency is fixed at 100 kHz.
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Real part and imaginary part of the axial component $H_z$ for ($r = R_1$ et $r = R_2$).</figcaption></figure>

The response of the eddy current probe is shown in Fig. 6 when the probe is moving along the axis of the tube with a number of $N_p$ positions. Let us consider some RMS error in % between the two sets of numerical data defined by:

$$ \xi = \sqrt{\frac{\frac{1}{N_p} \sum_{1}^{N_p} |S_{PS} - S_{FE}|^2}{\frac{1}{N_p} \sum_{1}^{N_p} |S_{FE}|^2}}, $$

with $S = R[\Omega]$, or $\Im(\bar{Z}) = X[\Omega]$ or $|\bar{Z}|[\Omega]$. We choose to compute some RMS error on the resistance $\xi_R$, on the reactance $\xi_X$ and on the modulus of the impedance $\xi_{|\bar{Z}|}$. The RMS errors and the computation time ($T_{Cpu}$) for the pseudo-spectral module are reported in Table 2 for this first experiment. For all configurations which were carried out for numerical valuations, the computation time for the reference field remains fixed about 25 s. Concerning the presentation of results, we often prefer to represent the response of the probe in the impedance plane on a Lissajous Curve as it is shown in Fig. 7. This representation will be adopted for all other results in the paper.

<table><caption>**Table 2.** Values of the RMS errors and time computation.</caption><thead><tr><th>ξ<sub>R</sub> [%]</th><th>ξ<sub>X</sub> [%]</th><th>ξ<sub>|Z|</sub> [%]</th><th>T<sub>Cpu</sub> [s]</th></tr></thead><tbody><tr><td>0.15</td><td>0.028</td><td>0.029</td><td>9</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Variations of the resistance and the reactance of the probe according to its position.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Variations of the resistance and the reactance of the probe displayed in a Lissajous Curve in the impedance plane.</figcaption></figure>

## 3.2 Second numerical experiment

In the experiment, we consider that the profile of conductivity of the material may vary along the radial direction. Fig. 8 display the variations of the conductivity between a maximum value of 2 MS/m and a minimum value of 0.5 MS/m. The thick-
<figure><img src="image_8.png" /><figcaption>Fig. 8. The conductivity profile of the material.</figcaption></figure>

ness of the conducting tube is about 2.5 mm while the tube's wall are defined by the same profile (z in [mm]):

$$a_1(z) = 14 + 1.5 \frac{1}{2} [1 + \cos(2\pi z/10)]$$

$$a_2(z) = 16.5 + 1.5 \frac{1}{2} [1 + \cos(2\pi z/10)]$$

<figure><img src="image_9.png" /><figcaption>Fig. 9. A slice view: the baseline of the internal wall is $Rp_1 = 14$ mm and those of the external wall is $Rp_2 = 16.5$ mm.</figcaption></figure>

Fig. 10 gives a representation of the variations in impedance of the probe in the impedance plane. This figure shows a good agreement between the two sets of simulated data. The RMS

<figure><img src="image_10.png" /><figcaption>Fig. 10. Variations of the resistance and the reactance of the probe displayed in a Lissajous Curve in the impedance plane.</figcaption></figure>

<table><caption>Table 3. Values of the RMS errors and time computation.</caption><thead><tr><th>ξ<sub>R</sub> [%]</th><th>ξ<sub>X</sub> [%]</th><th>ξ<sub>|Z|</sub> [%]</th><th>T<sub>Cpu</sub> [s]</th></tr></thead><tbody><tr><td>0.54</td><td>0.024</td><td>0.055</td><td>9</td></tr></tbody></table>

errors and the computation time ($T_{Cpu}$) for the pseudo-spectral module are reported in Table 3.

### 3.3 Third numerical experiment

This experiment considers this time a variation of conductivity along the Z axis. In this kind of example, since the conductivity may vary along the Z axis, it is necessary to substitute the constant parameter $\varepsilon_r^{22}$ by a novel convolution matrix. Since the complex permittivity is related to the conductivity by the relationship $\varepsilon_r = 1 + \frac{i\sigma}{\omega\varepsilon_0}$, we introduce first of all a matrix $\varepsilon_r$ so that:

$$\varepsilon_r = \mathbf{I}_d + i \frac{Z_0}{k} [\hat{\sigma}_z]$$

where $[\hat{\sigma}_z]$ stands for a matrix convolution. Then, we obtain the matrix operator $\Upsilon = \mathbf{I}_r \otimes \varepsilon_r$ where $\mathbf{I}_r$ denotes the identity
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** The conductivity profile of the material according to the Z axis.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** A slice view: the baseline of the internal wall is $R_{p1} = 14$ mm and those of the external wall is $R_{p2} = 18$ mm.</figcaption></figure>

matrix of dimension $M \times M$. So, the operator $L_{21}$ becomes:

$$L_{21} = \epsilon_r - \frac{1}{k^2} [\alpha] [\mu_r^{11}]^{-1} [\alpha] \quad (46)$$

The profile of conductivity is displayed in Fig. 11 while two different profiles are considered (see Fig. 12) for the geometry of the tube. The two analytic forms are expressed as (z in [mm]):

$$a_1(z) = 14 - 0.7 \frac{1}{2} [1 + \cos(2\pi z/10)]$$

$$a_2(z) = 18 - 1.5 \frac{1}{2} [1 + \cos(2\pi z/10)] [0.141 \exp 10z/10]$$

In this numerical experiment, the function $\sigma(z)$ has been chosen in order to be able to simulate edge effects but any function of $\sigma(z)$ may be introduced in the numerical process. Fig. 13 shows the variations in impedance of the probe in the impedance plane. The RMS errors and the computation time ($T_{Cpu}$) for the pseudo-spectral module are reported in Table 4.

<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Variations of the resistance and the reactance of the probe displayed in a Lissajous Curve in the impedance plane.</figcaption></figure>

<table><caption>**Table 4.** Values of the RMS errors and time computation.</caption><thead><tr><td>ξ<sub>R</sub> [%]</td><td>ξ<sub>X</sub> [%]</td><td>ξ<sub>|Z|</sub> [%]</td><td>T<sub>Cpu</sub> [s]</td></tr></thead><tbody><tr><td>0.85</td><td>0.31</td><td>0.32</td><td>9</td></tr></tbody></table>

## 3.4 Fourth numerical experiment

A last numerical experiment has been carried out in order to simulate a weld. The geometry is shown in Fig. 14.
<figure><img src="image_14.png" /><figcaption>Fig. 14. A slice view : the baseline of the internal wall is Rp₁ = 14 mm and those of the external wall is Rp₂ = 18 mm.</figcaption></figure>

The two profiles are defined by their analytical functions (z in [mm]):

$$ 
\begin{aligned}
a_1(z) &= 18 + (z + 5), \quad \text{if } -5 \le z \text{[mm]} \le -1 \\
&= 18 - (z - 5), \quad \text{if } 1 \le z \text{[mm]} \le 5 \\
&= 18, \quad \text{otherwise}
\end{aligned}
 $$

$$ 
\begin{aligned}
a_2(z) &= 18 + 1.5 \sqrt{1 - 4(z/10)^2}, \quad \text{if } -5 \le z \text{[mm]} \le +5 \\
&= 18, \quad \text{otherwise}
\end{aligned}
 $$

Between the first fictitious cylindrical surface (r = R₁ and the first profile a₁(z), the value of the conductivity is fixed at σ₁ = 1.37 MS/m while the conductivity in the weld, i.e between the first profile a₁(z) and the second one a₂(z), the conductivity is equal to σ₂ = 1.64 MS/m. Fig. 15 shows the discretization scheme. The operating frequency is equal to 50 kHz and the number of collocation points has been increased to achieve N + 1 = 14. The variations in impedance of the probe are displayed in Fig. 16, in the impedance plane. The

RMS errors and the computation time (T<sub>Cpu</sub>) are discussed in the next section.

<figure><img src="image_15.png" /><figcaption>Fig. 15. The three regions are discretized each of them with a set of 14 collocation points.</figcaption></figure>

## 3.5 A study of the convergence

The goal of this section is to study the behavior of the numerical model when the number of collocation points (N + 1) may

<figure><img src="image_16.png" /><figcaption>Fig. 16. Variations of the resistance and the reactance of the probe displayed in the impedance plane.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 17.** Errors on the resistance, the reactance and the modulus of the probe impedance and the computation time with the number *N*.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 18.** Structure of the matrix $Z_m$ to invert for $N = 13$.</figcaption></figure>

Fig. 17 displays the variations of errors and the computation time $T_{cpu}$ with an increasing of the number $N$. These results show the consistency of the numerical model. Otherwise, it is clear also that it is preferable to choose an odd integer for $N$.

Fig. 18 displays an example of the structure of the sparse matrix to be inverted for estimating the unknowns in **X**.

We can conclude that all numerical experiments show the validity of the proposed model. The agreement between finite

element data and simulated data provided by the PS module is satisfactory for each configuration and the errors are decreasing with an increasing number of $N$. The time computation remains very low for all cases.

## 4 Conclusion

In this paper, an innovative and efficient semi-analytical method has been proposed in order to be able to carry out fast numerical simulations of an Eddy Current probe moving in a tube with a variation of its geometry or of its constitutive electromagnetic parameters. The profiles of the internal/external tube's walls may vary analytically according to an arbitrary function $a(z)$. The only restriction is that this function must be single-valued and its Fourier Transform must exist. In order to address more general geometries, analytic function $a(z)$ may be also substituted by some expansion on piecewise polynomial functions such as B-spline functions [28,29]. Discontinuities in the geometry could be implemented since factorization rules [30] are satisfied. Since the geometry or the constitutive parameters are $z$ dependant, the diffraction problem cannot be solved by using a modal approach. This problem has been overcome by implementing a pseudospectral modal approach. Chebyshev polynomials have been used to be able to obtain the spatial derivatives of the fields at the Chebyshev collocation points with a high accuracy. Several numerical experiments confirm the validity of the global numerical model and computation times remain very low. This is mainly due to the fact that the source term appears in the second member of an equation written in a matrix form. Moreover, the numerical proposed approach does not require any mesh of the geometry since boundary conditions are
translated analytically thanks to a change of coordinates sys-
tem. This numerical approach can be also extended for solving
a great number of multiple applications in electromagnetic field
computation. Since the advantages of this innovative approach
has bee shown, future work will focus on 2.5D problems con-
sidering a 3D eddy current probe. This crucial step precedes
future research for solving full 3D problems.

References

1. J. Chandezon, D. Maystre, and G. Raoult, “A new theoret-
ical method for diffraction gratings and its numerical ap-
plication,” *Journal of Optics*, vol. 11, p. 235, 1980.

2. L. Li, J. Chandezon, G. Granet, and J. Plumey, “Rigorous and efficient gratings analysis method made easy optical engineers,” *Applied Optics*, vol. 38, no. 2, pp. 304–313, 1999.

3. E. E. Popov, *Gratings: Theory and Numeric Applications*.
Popov, Institut Fresnel, 2012.

4. E. Post, *Formal Structure of Electromagnetics: General Covariance and Electromagnetics*. (Amsterdam: North-Holland): Dover Publications Inc., 1962.

5. D. Prémel, “Generalization of the second order vector po-
tential formulation for arbitrary non-orthogonal curvilinear
coordinates systems from the covariant form of maxwell’s
equations,” *Journal of Electromagnetic Analysis and Ap-
plications*, vol. 4, no. 10, pp. 400–409, 2012.

6. F. Caire, D. Prémel, and G. Granet, “Semi-analytical com-
putation of a quasi-static field induced by an eddy current
probe in a conductor with a rough surface,” *The European
Physical Journal Applied Physics*, vol. 64, 2013.

7. F. Caire, D. Prémel, and G. Granet, “Semi-analytical com-
putation of a quasi-static field induced by a 3d eddy cur-
rent probe scanning a 2d layered conductor with paral-
lel rough interfaces,” *International Journal of Numeri-
cal Modelling: Electronic Networks, Devices and Fields*,
pp. n/a–n/a, 2014.

8. F. Caire, G. Granet, D. Prémel, G. Rubinnacci, A. Tam-
burrino, and S. Ventre, “The simposium project: numeri-
cal modeling of a carbon steel manufacturing benchmark
case,” in *18th workshop on Electromagnetic Nondestructive
Evaluation (ENDE)*, (Bratislava, Slovak republic.), 25-28
of June 2013.

9. F. Caire, D. Prémel, and G. 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*, vol. 51, no. 3, p. 6200604, 2015.

10. F. Caire, *Les équations de Maxwell covariantes pour le calcul rapide des champs diffractés par des conducteurs complexes. Application au Contrôle Non Destructif par courants de Foucault*. PhD thesis, Université Paris 11, 2014.

11. D. Olmos and B. D. Shizgal, “A pseudospectral method of solution of fisher’s equation,” *Journal of Computational and Applied Mathematics*, vol. 193, no. 1, pp. 219–242, 2006.

12. H. H. Yang and B. Shizgal, “Chebyshev pseudospec-
tral multidomain technique for viscous-flow calculation,”
*Computer Methods in Applied Mechanics and Engineer-
ing*, vol. 118, no. 1-2, pp. 47–61, 1994.
13. H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. Teukolsky, "A multidomain spectral method for solving elliptic equations," *Computer Physics Communications*, vol. 152, no. 3, pp. 253–273, 2003.

14. G. H. Rawitscher and I. Koltracht, "An efficient numerical spectral method for solving the schrödinger equation," *Computing in Science and Engineering*, vol. 7, no. 6, pp. 58–66, 2005.

15. D. Song, L. Yuan, and Y. Y. Lu, "Fourier-matching pseudospectral modal method for diffraction gratings," *J. Opt. Soc. Am. A*, vol. 28, pp. 613–620, Apr 2011.

16. Y.-P. Chiou, W.-L. Yeh, and N.-Y. Shih, "Analysis of highly conducting lamellar gratings with multidomain pseudospectral method," *Lightwave Technology, Journal of*, vol. 27, pp. 5151–5159, Nov 2009.

17. K. Varis and A. Baghai-Wadji, "Pseudo-spectral analysis of radially-diagonalized maxwell's equations in cylindrical co-ordinates," *Opt. Express*, vol. 11, pp. 3048–3062, Nov 2003.

18. C. Baudier, R. Dusséaux, K. Edee, and G. Granet, "Scattering of a plane wave by one-dimensionnal dielectric random rough surfaces-study with the curvilinear coordinate method," *Waves Random media*, vol. 14, pp. 61–74, 2004.

19. D. Prémel, "Computation of a quasi-static field induced by two long straight parallel wires in a conductor with a rough surface," *Journal of Physics D: Applied Physics*, vol. 41, no. 245305, 2008.

20. L. Trefethen, *Spectral Methods in MATLAB*. Software, environments, tools, Society for Industrial and Applied Mathematics, 2000.

21. Q. Liu, "A pseudospectral frequency-domain (psfd) method for computational electromagnetics," *Antennas and Wireless Propagation Letters, IEEE*, vol. 1, no. 1, pp. 131-134, 2002.

22. A. Q. C. Canuto, M. Y. Hussaini and T. A. Zang, *Spectral Methods in Fluid Dynamics*. Springer, 1988.

23. B. Costa, "Spectral methods for partial differential equations," *A Mathematical Journal*, vol. 6, no. 6, pp. 1-32, 2004.

24. R. Collin, *Foundations for microwaves engineering*. New York: McGraw-Hill, 1992.

25. J. M. B. Auld, "Review of advances in quantitative eddy current nondestructive evaluation," *Journal of Nondestructive Evaluation*, vol. 18, no. 1, pp. 3-36, 1999.

26. S. Burke and M. Ibrahim, "Mutual impedance of air-cored coils above a conducting plate," *Journal of Physics D: Applied Physics*, vol. 37, no. 13, pp. 1857-1868, 2004.

27. "Comsol multiphysics v4.2 user's guide http://www.comsol.com/," 2011.

28. M. Walz, T. Zebrowski, J. Küchenmeister, and K. Busch, "B-spline modal method: A polynomial approach compared to the fourier modal method," *Opt. Express*, vol. 21, pp. 14683-14697, Jun 2013.

29. G. Granet, "Efficient implementation of b-spline modal method for lamellar gratings," *J. Opt. Soc. Am. A*, vol. 31, pp. 332-337, Feb 2014.

30. L. Li, "Use of Fourier series in the analysis of discontinuities periodic structures," *Journal of Optical Society of America A*, vol. 13, no. 9, pp. 1870-1876, 1996.