Model of Magnetically Shielded Ferrite-Cored Eddy Current Sensor

Abstract

Computationally fast electromagnetic models of eddy current sensors are required in model-based measurements, machine interpretation approaches or in the sensor design phase. If a sensor geometry allows it, the analytical approach to the modeling has significant advantages in comparison to numerical methods, most notably less demanding implementation and faster computation. In this paper, we studied an eddy current sensor consisting of a transmitter coil with a finitely long I ferrite core, which was screened with a finitely thick magnetic shield. The sensor was placed above a conductive and magnetic half-layer. We used vector magnetic potential formulation of the problem with a truncated region eigenfunction expansion, and obtained expressions for the transmitter coil impedance and magnetic potential in all subdomains. The modeling results are in excellent agreement with the results using the finite element method. The model was also compared with the impedance measurement in the frequency range from 5 kHz to 100 kHz and the agreement is within $3\%$ for the resistance change due to the presence of the half-layer and $1\%$ for the inductance change. The presented model can be used for measurement of properties of metallic objects, sensor lift-off or nonconductive coating thickness.

1. Introduction

Eddy current sensors (probes) respond to changes in the material properties (electrical conductivity, magnetic permeability), geometry or position of the nearby electrically conductive, usually metallic, objects [1,2]. Their applications include nondestructive testing for defects and material changes in objects, such as heat exchangers in nuclear plants, oil well casings, pipelines or aircraft hulls [3,4,5]; measurement of the distance to such objects or their dimensions [6,7]; and material characterization [8,9,10]. Eddy current sensors offer detection of very small surface defects, measurement of dimensions or distances in the micrometer range, insensitivity to nonconductive dirt, such as oil, water or dust, and robustness to the environmental factors such as pressure, temperature or stress. However, absolute measurement using eddy current sensors can be hard to achieve leading to the need for calibration, and human or machine analysis and interpretation of the sensor signals [11,12,13].

Traditionally, an eddy current sensor consists of a transmitter coil that generates a magnetic field and induces the eddy currents in the conductive material [1]. Receiver coils pick up the resultant magnetic field due to the defect and changes in the dimensions and material properties. Receiver coils can be replaced by Hall, magnetoresistive or even super-conducting quantum interference devices (SQUID) sensors [14]. However, the overall characteristics of the eddy current sensors are primarily determined by the excitation, geometry and construction of the transmitter coil.

Fast electromagnetic models of transmitter coils are required in model-based measurements, machine interpretation approach or in parametric and sensitivity analyses in the sensor design phase. When it comes to more complicated sensor geometries, analytical electromagnetic models are not as versatile as numerical ones, e.g., finite-element or boundary-element models. However, if a real sensor and object under test can be approximated by a geometry that is simple enough to be dealt with analytically, then one can expect less demanding numerical implementation and higher computational speed in comparison to fully numerical approaches. More complicated geometries can be treated using hybrid approaches in which the analytical models are used for calculation of the incidence fields and numerical methods for calculation of the fields scattered from the small changes in the tested object [15].

Typically with eddy current sensors, the transmitter coil has a ferrite core to increase the inductance, sensitivity and signal-to-noise ratio or to achieve the desired magnetic field distribution. Since the seminal work by Dodd and Deeds [16], who studied air-cored transmitter coils in planar and cylindrical geometries (layers, rods, tubes), a number of more complicated geometries have been analyzed using the truncated region eigenfunction expansion approach (TREE) [17]. This includes ferrite-cored transmitter coils, which previously required numerical or empirical approaches. In a paper that laid out the modeling path followed by other authors, Theodoulidis analyzed the coil with the cylindrical finite-length ferrite core that was positioned above a half space made of two nonmagnetic and conductive layers [18]. This analysis was extended to a half space made of multiple magnetic and conductive layers by Lu et al. in [19]. Bayani et al., in [20], developed the model of a cup-cored coil above a nonmagnetic conductive double layered half space. In all three studies, the models were verified experimentally and the relative error between the theoretical results and measurements were on average between $0.5\%$ and $2.5\%$. Sakkaki and Bayani, in [21], modeled an E-cored coil above a conductive double layer, and they numerically verified the model by comparing it to a finite element model using Comsol. Tytko and Dziczkowski, in [22,23,24], modeled the I-cored and E-cored coils in the presence of the layered medium with a hole, and they verified their work using a Comsol-based finite element model. Zhu et al., in [25], analyzed an I-cored coil screened with a thin sheet of material with infinite permeability.

In this paper, we studied a coil with a finitely long I core and screened with a finitely thick magnetic shield placed above a conductive and magnetic half-layer using the TREE approach. Figure 1 and Figure 2 depict the geometry in more detail. We made no assumptions regarding the shield thickness and the permeability of the shield and the core. The open upper side of the shield allows for the coil wires and connector towards a measurement instrument. This sensor geometry is typical and practically important, but to the best of our knowledge it was not modeled analytically before. The contributions of this paper are the analytical model of the sensor impedance above a plate made of conductive and magnetic material in Section 2 and accompanied Appendix A, analysis of numerical implementation issues in Section 3, comparison with a finite element model in Section 4, experimental verification in Section 5 and sensitivity analysis with respect to the sensor dimension in Section 6. The immediate application of these contributions could allow model-based measurement of material properties and lift-off (e.g., nonconductive coating thickness) with no or less calibration points (measurement above reference standards with known material and lift-off) compared to the existing measurement procedures [26]. The model can also be used in the sensor design phase where it is often necessary to adjust the sensor properties to the specific application or account for sensor production tolerances. Furthermore, we believe that the modeling approach and discussion on the implementation issues can be useful in modeling other similar induction sensors with ferrites.

2. Analysis

The geometry of the problem is shown in Figure 1 and Figure 2. The coil with a rectangular cross section is wound around the ferrite core with relative magnetic permeability ${\mu }_{\mathrm{C}}$. The coil and its core are screened with the ferrite shield with relative magnetic permeability ${\mu }_{\mathrm{S}}$. The core and the shield are vertically aligned at the plane $z=0$. The core and the shield are assumed to be nonconductive. The sensor is located over the magnetic and conductive half space (plate) with electrical conductivity $\sigma$ and relative magnetic permeability ${\mu }_{\mathrm{r}}$.

Before we derive the vector magnetic potential of the coil, we will first analyze a filamentary sinusoidal current $Iexp\left(j\omega t\right)\delta \left(r-r_0\right)\delta \left(z-z_0\right)\widehat{\mathbf{\phi }}$, where $\delta$ is the Dirac delta function and $\widehat{\mathbf{\phi }}$ is the azimuthal unit vector for the cylindrical coordinate system. The problem is axially symmetric, so there is only $\phi$ component of the vector magnetic potential, i.e., $\mathbf{A}\left(r,z\right)=A\left(r,z\right)\widehat{\mathbf{\phi }}$.

The problem domain is truncated at far enough boundary $r=R$, where we will impose the Dirichlet boundary condition $A\left(r=R,z\right)=0$. The consequence of the truncation is that the final expressions for the potential are in series rather than integral form. The problem domain is divided into horizontal regions indicated generally by index n in the expressions, or specifically by a number or the letter C. If there are multiple materials or the coil in a region, the region is divided into subregions that are indicated generally by index m or specifically by a number following the region designation. For example, $A_1$ is the potential in region 1 (no subregions) and $A_{23}$ in region 2, subregion 3.

Following the separation of variables, the general form of the magnetic potential is

$$\begin{array}{cc}\hfill & A_n\left(r,z\right)=\sum _{i=1}^{\infty }{R}_n\left({\alpha }_{n,i}r\right)Z_n\left({\beta }_{n,i}z\right)=\hfill \\ \hfill & =\sum _{i=1}^{\infty }\left[U_{n,i}{\mathrm{J}}_1\left({\alpha }_{n,i}r\right)+V_{n,i}{\mathrm{Y}}_1\left({\alpha }_{n,i}r\right)\right]\left[C_{n,i}exp\left(-{\beta }_{n,i}z\right)+D_{n,i}exp\left({\beta }_{n,i}z\right)\right],\mathrm{for}n\in \left\{1,5,6\right\},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & A_{nm}\left(r,z\right)=\sum _{i=1}^{\infty }{R}_{nm}\left({\alpha }_{n,i}r\right)Z_n\left({\alpha }_{n,i}z\right)=\hfill \\ \hfill & =\sum _{i=1}^{\infty }\left[U_{nm,i}{\mathrm{J}}_1\left({\alpha }_{n,i}r\right)+V_{nm,i}{\mathrm{Y}}_1\left({\alpha }_{n,i}r\right)\right]\left[C_{n,i}exp\left(-{\alpha }_{n,i}z\right)+D_{n,i}exp\left({\alpha }_{n,i}z\right)\right],\mathrm{for}n\in \left\{2,3,4,\mathrm{C}\right\},\hfill \end{array}$$

(2)

where ${\mathrm{J}}_1$ and ${\mathrm{Y}}_1$ are first-order Bessel functions of the first and second kind, respectively, ${\alpha }_{n,i}$ are eigenvalues arising from the Dirichlet boundary condition at $r=R$ for each region, Section 2.2, and

$${\beta }_{n,i}=\sqrt{{\alpha }_{n,i}^2+j\omega {\mu }_0{\mu }_n{\sigma }_n}.$$

(3)

If a region is homogeneous with relative permeability ${\mu }_n$ and conductivity ${\sigma }_n$, i.e., it has no multiple subregions, the form of the potential is given by (1). For regions with multiple subregions, the general form of the potential in subregion m of region n with relative permeability ${\mu }_{nm}$ is given by (2). It is important for further steps to note that there are no conductive materials in the regions with multiple subregions (only air or ferrites) and that, consequently, both r-related and z-related parts have the same eigenvalues ${\alpha }_{n,i}$. This means that z-related parts of the subregions in a given region are the same in case of (2).

The homogeneous regions $n\in \left\{1,5,6\right\}$ have their corresponding coefficients $V_i=0$ because of the divergence of ${\mathrm{Y}}_1$ at $r=0$, so, consequently, one can freely set $U_i=1$. Coefficients $U_{nm,i}$ and $V_{nm,i}$ in (2) for regions $n\in \left\{2,3,4,\mathrm{C}\right\}$ are determined from the interface conditions along the radial boundaries between the subregions, as described in Section 2.1. The interface conditions in this case are continuity of $B_r$ and $H_z$, i.e., component of $\mathbf{B}$ in $\widehat{\mathbf{r}}$ direction, and component of $\mathbf{H}$ in $\widehat{\mathbf{z}}$ direction.

The unknown coefficients $C_{n,i}$ and $D_{n,i}$ in (1) and (2) are determined from the boundary and interface conditions along the horizontal boundaries between the regions, Section 2.4. The potential must remain finite as $z\to \pm \infty$, so $D_{1,i}=0$ and $C_{6,i}=0$. The interface conditions for the horizontal boundaries are continuity of $B_z$ and $H_r$.

2.1. Continuity of $B_r$ and $H_z$ at Radial Boundaries between Subregions

This condition has to be satisfied for regions 2, 3 and 4. We will analyze a general case of a region n with subregions $m=1,\dots ,M$ and radial interfaces between these subregions at $r_{n1},\dots ,r_{n(M-1)}$. In that case, the radial parts under the sum in (2) for each subregion are

$$\begin{array}{c}{R}_{n1}\left({\alpha }_{n,i}r\right)={\mathrm{J}}_1\left({\alpha }_{n,i}r\right),0\le r\le r_{n1},\end{array}$$

(4)

$$\begin{array}{c}{R}_{n2}\left({\alpha }_{n,i}r\right)=U_{n2,i}{\mathrm{J}}_1\left({\alpha }_{n,i}r\right)+V_{n2,i}{\mathrm{Y}}_1\left({\alpha }_{n,i}r\right),r_{n1}<r\le r_{n2},\\ ⋮\end{array}$$

(5)

$$\begin{array}{c}{R}_{nm}\left({\alpha }_{n,i}r\right)=U_{nm,i}{\mathrm{J}}_1\left({\alpha }_{n,i}r\right)+V_{nm,i}{\mathrm{Y}}_1\left({\alpha }_{n,i}r\right),r_{n(m-1)}<r\le r_{nm},\\ ⋮\end{array}$$

(6)

$$\begin{array}{c}{R}_{nM}\left({\alpha }_{n,i}r\right)=U_{nM,i}{\mathrm{J}}_1\left({\alpha }_{n,i}r\right)+V_{nM,i}{\mathrm{Y}}_1\left({\alpha }_{n,i}r\right),r_{n(M-1)}<r\le R,r_{nM}=R.\end{array}$$

(7)

For the two neighboring subregions m and $m+1$, the radial interface of continuity of $B_r$ and $H_z$ at $r=r_{nm}$ is equivalent to

$$\begin{array}{c}{A}_{nm}\left(r_{nm},z\right)=A_{n(m+1)}\left(r_{nm},z\right),\end{array}$$

(8)

$$\begin{array}{c}\frac{1}{{\mu }_{nm}}\left(\frac{A_{nm}}{r}+\frac{\partial A_{nm}}{\partial r}\right){|}_{r=r_{nm}}=\frac{1}{{\mu }_{n(m+1)}}\left(\frac{A_{n(m+1)}}{r}+\frac{\partial A_{n(m+1)}}{\partial r}\right){|}_{r=r_{nm}}.\end{array}$$

(9)

Since all subregions in region n have the same z-related part, the radial interface conditions can be satisfied by adjusting the coefficients of the r-related part only. The notation will be simpler if we define matrix $\mathbf{R}$ as

$$\mathbf{R}\left({\alpha }_{n,i}{r}_{nm},{\mu }_{nm}\right)=\left[\begin{array}{cc}{\mathrm{J}}_1\left({\alpha }_{n,i}{r}_{nm}\right)& {\mathrm{Y}}_1\left({\alpha }_{n,i}{r}_{nm}\right)\\ {\mathrm{J}}_0\left({\alpha }_{n,i}{r}_{nm}\right)/{\mu }_{nm}& {\mathrm{Y}}_0\left({\alpha }_{n,i}{r}_{nm}\right)/{\mu }_{nm}\end{array}\right].$$

(10)

It can be shown using (2)–(10) that the coefficients in the r-related part of subregion $m+1$ are connected to the coefficients of subregion m as

$$\left[\begin{array}{c}{U}_{n(m+1),i}\\ V_{n(m+1),i}\end{array}\right]={\mathbf{R}}^{-1}\left({\alpha }_{n,i}{r}_{nm},{\mu }_{n(m+1)}\right)\mathbf{R}\left({\alpha }_{n,i}{r}_{nm},{\mu }_{nm}\right)\left[\begin{array}{c}{U}_{nm,i}\\ V_{nm,i}\end{array}\right].$$

(11)

Using recurrence relation (11), and since $V_{n1,i}=0$ and $U_{n1,i}$ can be freely set to 1, we can determine the r-related coefficients for all subregions in a given region:

$$\left[\begin{array}{c}{U}_{n1,i}\\ V_{n1,i}\end{array}\right]=\left[\begin{array}{c}1\\ 0\end{array}\right],\mathrm{for}m=1,0\le r\le r_{n1},$$

(12)

$$\left[\begin{array}{c}{U}_{n2,i}\\ V_{n2,i}\end{array}\right]={\mathbf{R}}^{-1}\left({\alpha }_{n,i}{r}_{n1},{\mu }_{n2}\right)\left[\begin{array}{c}{\mathrm{J}}_1\left({\alpha }_{n,i}{r}_{n1}\right)\\ {\mathrm{J}}_0\left({\alpha }_{n,i}{r}_{n1}\right)/{\mu }_{n1}\end{array}\right],\mathrm{for}m=2,r_{n1}<r\le r_{n2},$$

(13)

$$\begin{gathered} \left[\begin{array}{c}{U}_{nm,i}\\ V_{nm,i}\end{array}\right]=\prod _{k=3}^m\left({\mathbf{R}}^{-1}\left({\alpha }_{n,i}{r}_{n(k-1)},{\mu }_{nk}\right)\mathbf{R}\left({\alpha }_{n,i}{r}_{n(k-1)},{\mu }_{n(k-1)}\right)\right)\left[\begin{array}{c}{U}_{n2,i}\\ V_{n2,i}\end{array}\right], \\ \begin{array}{c}\mathrm{for}m=3,\dots ,M,r_{n(m-1)}<r\le r_{nm},r_{nM}=R.\end{array} \end{gathered}$$

(14)

Because coefficients $U_{nm,i}$ and $V_{nm,i}$ do not depend on spatial coordinates r and z, the radial parts $R_{nm}\left({\alpha }_{n,i}r\right)$ given in (4)–(7) are linear combinations of Bessel functions ${\mathrm{J}}_1$ and ${\mathrm{Y}}_1$ with respect to r. This will be important later on when we make use of the orthogonality property of the Bessel functions. In order to facilitate this, we will add an index to $R_{nm}$ to indicate the order of the involved Bessel functions:

$$R_{\nu ,nm}\left({\alpha }_{n,i}r\right)=U_{nm,i}{\mathrm{J}}_{\nu }\left({\alpha }_{n,i}r\right)+V_{nm,i}{\mathrm{Y}}_{\nu }\left({\alpha }_{n,i}r\right),$$

(15)

where coefficients $U_{nm,i}$ and $V_{nm,i}$ do not depend on $\nu$. It follows from (8) and (9) that

$$\begin{array}{c}{R}_{1,nm}\left({\alpha }_{n,i}{r}_{nm}\right)=R_{1,n(m+1)}\left({\alpha }_{n,i}{r}_{nm}\right),\end{array}$$

(16)

$$\begin{array}{c}\frac{1}{{\mu }_{nm}}{R}_{0,nm}\left({\alpha }_{n,i}{r}_{nm}\right)=\frac{1}{{\mu }_{n(m+1)}}{R}_{0,n(m+1)}\left({\alpha }_{n,i}{r}_{nm}\right).\end{array}$$

(17)

2.2. Eigenvalues

The radial parts of all regions and subregions are defined in the preceding section. The first step in their computation is to find the eigenvalues that satisfy the Dirichlet condition at the truncated boundary of the problem:

$$\begin{array}{c}{A}_n\left(r=R,z\right)=0\Rightarrow {\mathrm{J}}_1\left({\alpha }_{n,i}R\right)=0,\mathrm{for}\mathrm{homogeneous}\mathrm{region},\end{array}$$

(18)

$$\begin{array}{c}{A}_{nM}\left(r=R,z\right)=0\Rightarrow R_{1,nM}\left({\alpha }_{n,i}R\right)=0,\mathrm{for}\mathrm{region}\mathrm{with}M\mathrm{subregions}.\end{array}$$

(19)

The eigenvalues ${\alpha }_{n,i}$ are the real and positive roots of (18) and (19). Homogeneous regions 1, 5 and 6 have identical radial parts and, hence, the same eigenvalues. Similarly, the regions 3, 4 and C, with multiple subregions, have the same radial parts and the same eigenvalues. Region 2 has its own eigenvalues. In order to simplify the notation, which was so far useful for the general case, we will designate the elements of the three sets of eigenvalues as ${\alpha }_i$, $p_i$ and $q_i$:

$$\begin{array}{c}{\alpha }_i={\alpha }_{1,i}={\alpha }_{5,i}={\alpha }_{6,i},\end{array}$$

(20)

$$\begin{array}{c}{p}_i={\alpha }_{2,i},\end{array}$$

(21)

$$\begin{array}{c}{q}_i={\alpha }_{3,i}={\alpha }_{4,i}={\alpha }_{C,i},\end{array}$$

(22)

and for the z-related part in region 6

$${\beta }_i={\beta }_{6,i}=\sqrt{{\alpha }_i^2+j\omega {\mu }_0{\mu }_{\mathrm{r}}\sigma }.$$

(23)

2.3. Final Expressions for Potential in All Regions in Matrix Form

Infinite series in (1) and (2) have to be truncated to N terms in the computation. The potential is

$$A_1\left(r,z\right)=\sum _{i=1}^N{\mathrm{J}}_1\left({\alpha }_{i}r\right)exp\left(-{\alpha }_{i}z\right)C_{1,i}/{\alpha }_i,$$

(24)

$$A_2\left(r,z\right)=\sum _{i=1}^N\left\{\begin{array}{cc}{\mathrm{J}}_1\left(p_{i}r\right),& \mathrm{if}0\le r\le r_{\mathrm{S}1}\hfill \\ R_{1,22}\left(p_{i}r\right),& \mathrm{if}{r}_{\mathrm{S}1}<r\le r_{\mathrm{S}2}\hfill \\ R_{1,23}\left(p_{i}r\right),& \mathrm{if}{r}_{\mathrm{S}2}<r\le R\hfill \end{array}\right\}\left[exp\left(-p_{i}z\right)C_{2,i}+exp\left(p_{i}z\right)D_{2,i}\right]/p_i,$$

(25)

$$A_3\left(r,z\right)=\sum _{i=1}^N\left\{\begin{array}{cc}{\mathrm{J}}_1\left(q_{i}r\right),& \mathrm{if}0\le r\le r_{\mathrm{C}}\hfill \\ R_{1,32}\left(q_{i}r\right),& \mathrm{if}{r}_{\mathrm{C}}<r\le r_{\mathrm{S}1}\hfill \\ R_{1,33}\left(q_{i}r\right),& \mathrm{if}{r}_{\mathrm{S}1}<r\le r_{\mathrm{S}2}\hfill \\ R_{1,34}\left(q_{i}r\right),& \mathrm{if}{r}_{\mathrm{S}2}<r\le R\hfill \end{array}\right\}\left[exp\left(-q_{i}z\right)C_{3,i}+exp\left(q_{i}z\right)D_{3,i}\right]/q_i,$$

(26)

$$A_4\left(r,z\right)=\sum _{i=1}^N\left\{\begin{array}{cc}{\mathrm{J}}_1\left(q_{i}r\right),& \mathrm{if}0\le r\le r_{\mathrm{C}}\hfill \\ R_{1,42}\left(q_{i}r\right),& \mathrm{if}{r}_{\mathrm{C}}<r\le r_{\mathrm{S}1}\hfill \\ R_{1,43}\left(q_{i}r\right),& \mathrm{if}{r}_{\mathrm{S}1}<r\le r_{\mathrm{S}2}\hfill \\ R_{1,44}\left(q_{i}r\right),& \mathrm{if}{r}_{\mathrm{S}2}<r\le R\hfill \end{array}\right\}\left[exp\left(-q_{i}z\right)C_{4,i}+exp\left(q_{i}z\right)D_{4,i}\right]/q_i,$$

(27)

$$A_5\left(r,z\right)=\sum _{i=1}^N{\mathrm{J}}_1\left({\alpha }_{i}r\right)\left[exp\left(-{\alpha }_{i}z\right)C_{5,i}+exp\left({\alpha }_{i}z\right)D_{5,i}\right]/{\alpha }_i,$$

(28)

$$A_5\left(r,z\right)=\sum _{i=1}^N{\mathrm{J}}_1\left({\alpha }_{i}r\right)\left[exp\left(-{\alpha }_{i}z\right)C_{5,i}+exp\left({\alpha }_{i}z\right)D_{5,i}\right]/{\alpha }_i,$$

(29)

$$A_6\left(r,z\right)=\sum _{i=1}^N{\mathrm{J}}_1\left({\alpha }_{i}r\right)exp\left({\beta }_{i}z\right)D_{6,i}/{\beta }_i.$$

(30)

In the above expressions, we divided with the corresponding eigenvalues in order to make later expressions somewhat simpler, which is allowed because it will be canceled out by the final form of $C_{n,i}$ and $D_{n,i}$ [17].

It is convenient to write (24)–(30) in the matrix form. To do that, we will assume that the functions $J_{\nu }\left(\right)$, $R_{\nu ,nm}\left(\right)$, $exp\left(\right)$, as well as the integration and differentiation are applied element wise on a vector or matrix. We will use boldface type to represent matrices (e.g., $\mathbf{\alpha }$ or $exp\left(\mathbf{q}z\right)$), overlined italic type to represent row vectors (e.g., $\overline{\alpha }$ or ${\mathrm{J}}_1\left(\overline{\alpha }r\right)$) and underlined italic type to represent column vectors (e.g., ${\underline{C}}_1$). The dimensions of the matrices are $N\times N$, row vectors $1\times N$ and column vectors $N\times 1$. This notation eases implementation in a programming language such as Matlab or Julia. Using this notation, the matrix form of (24)–(30) is

$$\begin{array}{c}{A}_1\left(r,z\right)={\mathrm{J}}_1\left(\overline{\alpha }r\right){\mathbf{\alpha }}^{-1}exp\left(-\mathbf{\alpha }z\right){\underline{C}}_1,\end{array}$$

(31)

$$\begin{array}{c}{A}_2\left(r,z\right)=\left\{\begin{array}{c}{\mathrm{J}}_1\left(\overline{p}r\right)\\ R_{1,22}\left(\overline{p}r\right)\\ R_{1,23}\left(\overline{p}r\right)\end{array}\right\}{\mathbf{p}}^{-1}\left[exp\left(-\mathbf{p}z\right){\underline{C}}_2+exp\left(\mathbf{p}z\right){\underline{D}}_2\right],\end{array}$$

(32)

$$\begin{array}{c}{A}_3\left(r,z\right)=\left\{\begin{array}{c}{\mathrm{J}}_1\left(\overline{q}r\right)\\ R_{1,32}\left(\overline{q}r\right)\\ R_{1,33}\left(\overline{q}r\right)\\ R_{1,34}\left(\overline{q}r\right)\end{array}\right\}{\mathbf{q}}^{-1}\left[exp\left(-\mathbf{q}z\right){\underline{C}}_3+exp\left(\mathbf{q}z\right){\underline{D}}_3\right],\end{array}$$

(33)

$$\begin{array}{c}{A}_4\left(r,z\right)=\left\{\begin{array}{c}{\mathrm{J}}_1\left(\overline{q}r\right)\\ R_{1,42}\left(\overline{q}r\right)\\ R_{1,43}\left(\overline{q}r\right)\\ R_{1,44}\left(\overline{q}r\right)\end{array}\right\}{\mathbf{q}}^{-1}\left[exp\left(-\mathbf{q}z\right){\underline{C}}_4+exp\left(\mathbf{q}z\right){\underline{D}}_4\right],\end{array}$$

(34)

$$\begin{array}{c}{A}_5\left(r,z\right)={\mathrm{J}}_1\left(\overline{\alpha }r\right){\mathbf{\alpha }}^{-1}\left[exp\left(-\mathbf{\alpha }z\right){\underline{C}}_5+exp\left(\mathbf{\alpha }z\right){\underline{D}}_5\right],\end{array}$$

(35)

$$\begin{array}{c}{A}_6\left(r,z\right)={\mathrm{J}}_1\left(\overline{\alpha }r\right){\mathbf{\beta }}^{-1}exp\left(\mathbf{\beta }z\right){\underline{D}}_6.\end{array}$$

(36)

In the above expressions, matrices $\mathbf{\alpha }$, $\mathbf{p}$ and $\mathbf{q}$ are diagonal with elements of $\overline{\alpha }$, $\overline{p}$ and $\overline{q}$ on their main diagonal, respectively.

2.4. Continuity of $B_z$ and $H_r$ for Horizontal Boundaries between Regions

The continuity of $B_z$ and $H_r$ have to be satisfied along the horizontal interface ($z=\mathrm{const}.$) between two regions. For two neighboring regions n and $n+1$, the horizontal interface condition at $z=z_n$ is equivalent to

$$\begin{array}{c}\left(\frac{A_{n+1}}{r}+\frac{\partial A_{n+1}}{\partial r}\right){|}_{z=z_n}-\left(\frac{A_n}{r}+\frac{\partial A_n}{\partial r}\right){|}_{z=z_n}=0,\end{array}$$

(37)

$$\begin{array}{c}\frac{1}{{\mu }_{n+1}}\frac{\partial A_{n+1}}{\partial z}{|}_{z=z_n}-\frac{1}{{\mu }_n}\frac{\partial A_n}{\partial z}{|}_{z=z_n}=\left\{\begin{array}{cc}{\mu }_{0}I\delta \left(r-r_0\right),\hfill & \mathrm{if}{z}_n=z_0,\hfill \\ \hfill 0,& \mathrm{otherwise}.\hfill \end{array}\right.\end{array}$$

(38)

Relative permeability ${\mu }_n$ in (38) is the piecewise constant function of r for region n if it has multiple subregions, i.e., ${\mu }_n={\mu }_{n1}$ for $0\le r\le r_{n1}$,…, ${\mu }_n={\mu }_{nM}$ for $r_{n(M-1)}<r\le R$.

The interface relations that have to be satisfied for all five horizontal boundaries ($z=L_{\mathrm{S}}$, $z=L_{\mathrm{C}}$, $z=z_0$, $z=0$ and $z=-h$) are given in Appendix A for the sake of clarity. In order for these relations, e.g., (A3), to be satisfied for the entire interval $0\le r\le R$, the procedure followed in Appendix A is to expand the r-dependent parts using a basis formed of continuous piecewise linear combinations of Bessel functions of the first and second kind. The choice of the basis is somewhat arbitrary but it is better, from the standpoint of numerical implementation, to choose these in such a way that as many of the resulting matrices are identical or diagonal. Term-by-term comparison results in a matrix equation for each interface relation involving corresponding vectors $\underline{C}$ and $\underline{D}$. This procedure has to be repeated twice (continuity of $B_z$ and $H_r$) for each of the 5 horizontal interfaces giving rise to 10 algebraic equations with matrix coefficients and 10 unknown vectors $\underline{C}$ and $\underline{D}$.

2.5. Calculation of Coefficients $\underline{C}$ and $\underline{D}$

Coefficient vectors $\underline{C}$ and $\underline{D}$ in (31)–(36) for each of the regions are calculated by solving $10\times 10$ system of algebraic equations derived in Appendix A, i.e., (A5), (A10), (A13), (A16), (A18), (A21), (A25), (A28), (A30) and (A31). First, we define some auxiliary expressions sorted in order of their computational dependencies:

$$\begin{array}{c}{\mathbf{T}}_1=exp\left[-\mathbf{p}\left(L_{\mathrm{S}}-L_{\mathrm{C}}\right)\right]\left[2\mathbf{E}{\left(\mathbf{E}+\mathbf{F}\right)}^{-1}-\mathbf{I}\right]exp\left[-\mathbf{p}\left(L_{\mathrm{S}}-L_{\mathrm{C}}\right)\right],\end{array}$$

(39)

$$\begin{array}{c}\mathsf{\Lambda }=\frac{{\mu }_{\mathrm{r}}\mathbf{\alpha }-\mathbf{\beta }}{{\mu }_{\mathrm{r}}\mathbf{\alpha }+\mathbf{\beta }}exp\left(-2\mathbf{\alpha }h\right),\end{array}$$

(40)

$$\begin{array}{c}{\mathbf{W}}^{\prime }=\left(\mathbf{I}-{\mathbf{T}}_1\right)\mathbf{W},\end{array}$$

(41)

$$\begin{array}{c}{\mathbf{L}}^{\prime }=\left(\mathbf{I}+{\mathbf{T}}_1\right)\mathbf{L},\end{array}$$

(42)

$$\begin{array}{c}{\mathbf{V}}^{\prime }=\mathbf{V}\left(\mathbf{I}+\mathsf{\Lambda }\right),\end{array}$$

(43)

$$\begin{array}{c}{\mathbf{U}}^{\prime }=\mathbf{U}\left(\mathbf{I}-\mathsf{\Lambda }\right),\end{array}$$

(44)

$$\begin{array}{c}{\mathbf{T}}_2=-exp\left(-\mathbf{q}{L}_{\mathrm{C}}\right){\left({\mathbf{W}}^{\prime }+{\mathbf{L}}^{\prime }\right)}^{-1}\left({\mathbf{W}}^{\prime }-{\mathbf{L}}^{\prime }\right)exp\left(-\mathbf{q}{L}_{\mathrm{C}}\right),\end{array}$$

(45)

$$\begin{array}{c}{\mathbf{T}}_3=\frac{1}{2}\left[\left(\mathbf{I}-{\mathbf{T}}_2\right){\mathbf{S}}^{-1}{\mathbf{V}}^{\prime }+\left(\mathbf{I}+{\mathbf{T}}_2\right){\mathbf{S}}^{-1}{\mathbf{U}}^{\prime }\right],\end{array}$$

(46)

$$\begin{array}{c}{\mathbf{T}}_4=\left[\mathbf{W}+\mathbf{L}-\left(\mathbf{W}-\mathbf{L}\right){\left({\mathbf{W}}^{\prime }+{\mathbf{L}}^{\prime }\right)}^{-1}\left({\mathbf{W}}^{\prime }-{\mathbf{L}}^{\prime }\right)\right]exp\left(-\mathbf{q}{L}_{\mathrm{C}}\right).\end{array}$$

(47)

For the filamentary coil at $\left(r_0,z_0\right)$:

$$\begin{array}{c}\underline{X}=\frac{1}{2}{\mu }_{0}I{\mathbf{S}}^{-1}{r}_0{R}_{1,32}\left(\underline{q}{r}_0\right),\end{array}$$

(48)

$$\begin{array}{c}{\mathbf{X}}_{\mathrm{pos}}=exp\left(\mathbf{q}{z}_0\right),\end{array}$$

(49)

$$\begin{array}{c}{\mathbf{X}}_{\mathrm{neg}}=exp\left(-\mathbf{q}{z}_0\right).\end{array}$$

(50)

For the coil with the rectangular cross section:

$$\begin{array}{c}\underline{X}=\frac{1}{2}{\mu }_0{i}_{\mathrm{TX}}{\mathbf{S}}^{-1}{\mathbf{q}}^{-3}int\left(\underline{q}{r}_{\mathrm{T}1},\underline{q}{r}_{\mathrm{T}2}\right),\end{array}$$

(51)

$$\begin{array}{c}{\mathbf{X}}_{\mathrm{pos}}=exp\left(\mathbf{q}{z}_{\mathrm{T}2}\right)-exp\left(\mathbf{q}{z}_{\mathrm{T}1}\right),\end{array}$$

(52)

$$\begin{array}{c}{\mathbf{X}}_{\mathrm{neg}}=-exp\left(-\mathbf{q}{z}_{\mathrm{T}2}\right)+exp\left(\mathbf{q}{z}_{\mathrm{T}1}\right),\end{array}$$

(53)

where

$$int\left(x_1,x_2\right)={\int }_{x_1}^{x_2}x{R}_{1,32}\left(x\right)dx.$$

(54)

Finally, we can write the coefficients $\underline{C}$ and $\underline{D}$ sorted in order of their computational dependencies:

$$\begin{array}{c}{\underline{D}}_5={\mathbf{T}}_3^{-1}\left({\mathbf{T}}_2{\mathbf{X}}_{\mathrm{pos}}+{\mathbf{X}}_{\mathrm{neg}}\right)\underline{X},\end{array}$$

(55)

$$\begin{array}{c}{\underline{C}}_5=\mathsf{\Lambda }{\underline{D}}_5,\end{array}$$

(56)

$$\begin{array}{c}{\underline{D}}_6=2\frac{{\mu }_{\mathrm{r}}\mathbf{\beta }}{{\mu }_{\mathrm{r}}\mathbf{\alpha }+\mathbf{\beta }}exp\left[\left(\mathbf{\beta }-\mathbf{\alpha }\right)h\right]{\underline{D}}_5,\end{array}$$

(57)

$$\begin{array}{c}{\underline{C}}_4=\frac{1}{2}{\mathbf{S}}^{-1}\left({\mathbf{V}}^{\prime }-{\mathbf{U}}^{\prime }\right){\underline{D}}_5,\end{array}$$

(58)

$$\begin{array}{c}{\underline{D}}_4=\frac{1}{2}{\mathbf{S}}^{-1}\left({\mathbf{V}}^{\prime }+{\mathbf{U}}^{\prime }\right){\underline{D}}_5,\end{array}$$

(59)

$$\begin{array}{c}{\underline{C}}_3={\underline{C}}_4+{\mathbf{X}}_{\mathrm{pos}}\underline{X},\end{array}$$

(60)

$$\begin{array}{c}{\underline{D}}_3={\underline{D}}_4-{\mathbf{X}}_{\mathrm{neg}}\underline{X},\end{array}$$

(61)

$$\begin{array}{c}{\underline{C}}_2=\frac{1}{2}exp\left(\mathbf{p}{L}_{\mathrm{C}}\right){\mathbf{K}}^{-1}{\mathbf{T}}_4{\underline{C}}_3,\end{array}$$

(62)

$$\begin{array}{c}{\underline{D}}_2=\frac{1}{2}exp\left(-\mathbf{p}{L}_{\mathrm{S}}\right){\mathbf{K}}^{-1}\left[2\mathbf{E}{\left(\mathbf{E}+\mathbf{F}\right)}^{-1}-\mathbf{I}\right]exp\left[-\mathbf{p}\left(L_{\mathrm{S}}-L_{\mathrm{C}}\right)\right]{\mathbf{T}}_4{\underline{C}}_3,\end{array}$$

(63)

$$\begin{array}{c}{\underline{C}}_1=exp\left(\mathbf{\alpha }{L}_{\mathrm{S}}\right){\left(\mathbf{E}+\mathbf{F}\right)}^{-1}exp\left[-\mathbf{p}\left(L_{\mathrm{S}}-L_{\mathrm{C}}\right)\right]{\mathbf{T}}_4{\underline{C}}_3.\end{array}$$

(64)

2.6. Coil with Rectangular Cross Section

The magnetic potential of a coil with rectangular cross section and $N_{\mathrm{TX}}$ turns, carrying sinusoidal current of amplitude $I_{\mathrm{TX}}$, is calculated by the superposition of a number of filamentary coils. Let a coil with infinitesimal cross section $d{r}_{0}d{z}_0$ carrying current $i_{\mathrm{TX}}d{r}_{0}d{z}_0$ approach a filamentary current I at $\left(r_0,z_0\right)$, where the current density $i_{\mathrm{TX}}$ is

$$i_{\mathrm{TX}}=\frac{N_{\mathrm{TX}}{I}_{\mathrm{TX}}}{\left(r_{\mathrm{T}2}-r_{\mathrm{T}1}\right)\left(z_{\mathrm{T}2}-z_{\mathrm{T}1}\right)}.$$

(65)

In that case, the total magnetic potential for region n is

$$A_n^{\mathrm{tot}}\left(r,z\right)={\int }_{r_{\mathrm{T}1}}^{r_{\mathrm{T}2}}{\int }_{z_{\mathrm{T}1}}^{z_{\mathrm{T}2}}{A}_n\left(r,z,r_0,z_0\right)d{r}_{0}d{z}_0,$$

(66)

where $A_n\left(r,z,r_0,z_0\right)$ is given by (31)–(36) depending on the region. Integration in (66) comes down to integration of ${\mathbf{X}}_{\mathrm{pos}}\underline{X}$ and ${\mathbf{X}}_{\mathrm{neg}}\underline{X}$ in (48)–(50), and it results in (51)–(53). The total potential $A_n^{\mathrm{tot}}$ is given then by the same expressions as $A_n$ in (31)–(36), except one uses (51)–(53) instead of (48)–(50) in (55), and (60) and (61) for ${\underline{D}}_5$, ${\underline{C}}_3$ and ${\underline{D}}_3$, respectively. One should note that superposition for region C ($z_{\mathrm{T}1}\le z\le z_{\mathrm{T}1}$), which is present only in the case of the rectangular cross section coil, requires integration of $A_3\left(r,z,r_0,z_0\right)$ for $z_{\mathrm{T}1}\le z_0\le z$ and $A_4\left(r,z,r_0,z_0\right)$ for $z\le z_0\le z_{\mathrm{T}2}$, i.e.,

$$A_{\mathrm{C}}\left(r,z\right)={\int }_{r_{\mathrm{T}1}}^{r_{\mathrm{T}2}}{\int }_{z_{\mathrm{T}1}}^z{A}_3\left(r,z,r_0,z_0\right)d{r}_{0}d{z}_0+{\int }_{r_{\mathrm{T}1}}^{r_{\mathrm{T}2}}{\int }_z^{z_{\mathrm{T}2}}{A}_4\left(r,z,r_0,z_0\right)d{r}_{0}d{z}_0.$$

(67)

It can be shown that (67) results in

$$\begin{array}{c}{A}_{\mathrm{C}}\left(r,z\right)=A_4^{\mathrm{tot}}\left(r,z\right)+\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left\{\begin{array}{c}{\mathrm{J}}_1\left(\overline{q}r\right)\\ R_{1,42}\left(\overline{q}r\right)\\ R_{1,43}\left(\overline{q}r\right)\\ R_{1,44}\left(\overline{q}r\right)\end{array}\right\}{\mathbf{q}}^{-1}\left[2\mathbf{I}-exp\left(-\mathbf{q}z\right)exp\left(\mathbf{q}{z}_{\mathrm{T}1}\right)-exp\left(\mathbf{q}z\right)exp\left(-\mathbf{q}{z}_{\mathrm{T}1}\right)\right]\underline{X},\hfill \end{array}$$

(68)

where $\underline{X}$ is given by (51). The impedance of the coil, which is the most important output of the model from the measurement point of view, is calculated by the integration of the voltage induced in a single loop over the coil cross section:

$$Z=j\omega \frac{2\pi N_{\mathrm{TX}}}{\left(r_{\mathrm{T}2}-r_{\mathrm{T}1}\right)\left(z_{\mathrm{T}2}-z_{\mathrm{T}1}\right)I_{TX}}{\int }_{r_{\mathrm{T}1}}^{r_{\mathrm{T}2}}{\int }_{z_{\mathrm{T}1}}^{z_{\mathrm{T}2}}r{A}_{\mathrm{C}2}\left(r,z\right)drdz,$$

(69)

where $A_{C2}$ is the potential in region C, subregion 2. Finally, the impedance of the coil is

$$\begin{array}{c}Z=j\omega \frac{2\pi N_{\mathrm{TX}}}{\left(r_{\mathrm{T}2}-r_{\mathrm{T}1}\right)\left(z_{\mathrm{T}2}-z_{\mathrm{T}1}\right)I_{TX}}int\left(\overline{q}{r}_{\mathrm{T}1},\overline{q}{r}_{\mathrm{T}2}\right){\mathbf{q}}^{-4}\hfill \\ \hspace{1em}\hspace{1em} [-\left(exp\left(-\mathbf{q}{z}_{\mathrm{T}2}\right)-exp\left(-\mathbf{q}{z}_{\mathrm{T}1}\right)\right){\underline{C}}_4+\left(exp\left(\mathbf{q}{z}_{\mathrm{T}2}\right)-exp\left(\mathbf{q}{z}_{\mathrm{T}1}\right)\right){\underline{D}}_4+\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left(2\mathbf{q}\left(z_{\mathrm{T}2}-z_{\mathrm{T}1}\right)+exp\left[-\mathbf{q}\left(z_{\mathrm{T}2}-z_{\mathrm{T}1}\right)\right]-exp\left[\mathbf{q}\left(z_{\mathrm{T}2}-z_{\mathrm{T}1}\right)\right]\right)\underline{X}].\hfill \end{array}$$

(70)

The voltage induced in a receiver coil can be calculated using an approach analogous to (69) with the choice of the potential corresponding to the location of the receiver coil.

3. Numerical Implementation

The model can be implemented in a programming language suitable for numerical analysis, such as Matlab or Julia. We have chosen the former. The major functional parts of the implementation are:

• Calculation of coefficients for the radial part $R_{nm}\left({\alpha }_{n,i}r\right)$ given in (4)–(7) using (12)–(14).
• Root finding of (18) and (19) in order to obtain the required eigenvalues.
• Integration of $int\left(x_1,x_2\right)$ in (54).
• Calculation of the matrices in Appendix A ($\mathbf{E}$, $\mathbf{K}$, $\mathbf{F}$, $\mathbf{W}$, $\mathbf{L}$, $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{S}$) to form the set of 10 algebraic equations for 10 unknown coefficients $\underline{C}$ and $\underline{D}$.
• Calculation of the coefficients $\underline{C}$ and $\underline{D}$ using equations in Section 2.5.
• Calculation of the coil impedance Z using (70).

While steps 1, 4 and 6 were implemented straightforwardly following the relevant expressions in the paper, the remaining ones deserve some additional comments.

The functions on the left side of Equations (18) and (19) are real and oscillatory. We found their roots by sampling the functions so that intervals containing exactly one zero crossing can be selected and the roots bracketed. As a root finding algorithm, we used Anderson–Björck modification of the classic regula falsi method [27].

In general, the form of the integral in (54) involving linear combination of Bessel functions of the first and second kind, ${\mathrm{J}}_1$ and ${\mathrm{Y}}_1$, is more easily calculated using Struve functions or expansion in Chebyshev polynomials than by direct numerical integration [18]. We used Struve functions ${\mathbf{H}}_0$ and ${\mathbf{H}}_1$ calculated via their expansion into the series of Bessel functions for which the relevant identities can be found in [28]:

$${\int }_{x_1}^{x_2}t{R}_{1,nm}\left(t\right)dt=\frac{\pi t}{2}\left[{\mathbf{H}}_0\left(t\right)R_{1,nm}\left(t\right)-{\mathbf{H}}_1\left(t\right)R_{0,nm}\left(t\right)\right]{|}_{x_1}^{x_2}.$$

(71)

The expressions for $\underline{C}$ and $\underline{D}$ in Section 2.5 require some matrices to be inverted. Increasing the length N of the series representation increases the eigenvalues. The factors containing exponential functions with positively signed products of an eigenvalue and a spatial dimension can cause a loss of precision and a poor condition number of the matrices that need to be inverted. Thus, care should be taken to arrange the final expressions in such a way that as much as possible of the matrices are diagonal and that the exponential functions have arguments with negative sign.

4. Comparison with FEM

In order to verify the model of the sensor impedance, we compared the model predictions with the results of a finite elements model (FEM) study. In both numerical studies we used the sensor with properties as in

Table 1. Properties of the sensor used in the numerical studies.

Coil
Inner radius	$r_{\mathrm{T}1}$	$2.50$ mm
Outer radius	$r_{\mathrm{T}2}$	$3.15$ mm
Bottom position	$z_{\mathrm{T}1}$	$2.00$ mm
Top position	$z_{\mathrm{T}2}$	$4.00$ mm
Lift-off	h	$0.20$ mm
Number of turns	$N_{\mathrm{TX}}$	128
Core
Radius	$r_{\mathrm{C}}$	$1.75$ mm
Height	$L_{\mathrm{C}}$	$6.00$ mm
Rel. permeability	${\mu }_{\mathrm{C}}$	see the text
Shield
Inner radius	$r_{\mathrm{S}1}$	$3.65$ mm
Outer radius	$r_{\mathrm{S}2}$	$6.05$ mm
Height	$L_{\mathrm{S}}$	$15.50$ mm
Rel. permeability	${\mu }_{\mathrm{S}}$	see the text

. For the FEM study we used free, simple and thoroughly verified 2D solver Finite Element Model Magnetics (FEMM) that was successfully used in a number of studies [29]. The number of the series terms in the presented model of the sensor impedance was $N=140$ and the domain boundary was set to $R=10r_{\mathrm{S}2}$. The identical geometry in FEMM was meshed with approximately $2.4$ million of triangular elements. The computation time for one realization of the problem (i.e., for the selected dimensions, frequency and properties of the core, shield and medium) was around 30 min for FEMM, and less than $0.7$ $\mathrm{s}$ for the calculation of the model, including the eigenvalues computation and without usage of precalculated values.

First, we calculated the impedance of the coil depending on the permeability of the core and shield for the plate with $\sigma =5 \mathrm{M}\mathrm{S}/\mathrm{m}$ and ${\mu }_{\mathrm{r}}=50$ at 60 kHz, Figure 3. For the sake of simplicity, we assumed that ${\mu }_{\mathrm{C}}={\mu }_{\mathrm{S}}$. There is practically no significant further change in the impedance of the sensor for the permeability above 500, which is a finding observed for other cored coils as well [17]. The relative discrepancy between the model and FEM results are bellow $0.07\%$ for the inductance and $0.23\%$ for the resistance. These discrepancies can be further reduced by increasing the density of the FEM mesh.

Second, we calculated the impedance of the coil depending on the conductivity $\sigma$ of the plate for two values of its permeability ${\mu }_{\mathrm{r}}=1$ and ${\mu }_{\mathrm{r}}=50$ at 60 kHz, Figure 4. The core permeability ${\mu }_{\mathrm{C}}=100$ and the shield permeability ${\mu }_{\mathrm{S}}=50$. Typical dependency of the coil impedance on the plate conductivity can be observed, including the maximum point of the resistance [17]. The relative discrepancy between the model and FEM results are bellow $0.07\%$ for the inductance and $0.25\%$ for the resistance.

Finally, we calculated the impedance of the coil depending on the lift-off h of the sensor placed above either a nonmagnetic plate with ${\mu }_{\mathrm{r}}=1$ and $\sigma =25 \mathrm{M}\mathrm{S}/\mathrm{m}$ (properties similar to an aluminum alloy), and a magnetic plate ${\mu }_{\mathrm{r}}=50$ and $\sigma =5\mathrm{M}\mathrm{S}/\mathrm{m}$ (properties similar to a carbon steel), Figure 5. The core permeability was ${\mu }_{\mathrm{C}}=100$ and the shield permeability ${\mu }_{\mathrm{S}}=50$. Again, the relative discrepancy between the model and FEM results are bellow $0.07\%$ for the inductance and $0.25\%$ for the resistance. As expected, the resistance decreases with the lift-off increase for both plates, whereas the inductance in the case of the nonmagnetic plate increases with the lift-off, and decreases in the case of the magnetic plate.

5. Comparison with Experimental Results

We manufactured the sensor according to the geometry in Figure 1 and Figure 2. The sensor properties are given in

Table 2. Properties of the sensor prototype.

Coil
Inner radius	$r_{\mathrm{T}1}$	$5.00$ mm
Outer radius	$r_{\mathrm{T}2}$	$5.09$ mm
Bottom position	$z_{\mathrm{T}1}$	$2.45$ mm
Top position	$z_{\mathrm{T}2}$	$17.45$ mm
Number of turns	$N_{\mathrm{TX}}$	125
Core
Radius	$r_{\mathrm{C}}$	$4.98$ mm
Height	$L_{\mathrm{C}}$	$20.03$ mm
Material		$\mathrm{HF}70$
Rel. permeability	${\mu }_{\mathrm{C}}$	1500
Shield
Inner radius	$r_{\mathrm{S}1}$	$6.45$ mm
Outer radius	$r_{\mathrm{S}2}$	$12.92$ mm
Height	$L_{\mathrm{S}}$	$28.29$ mm
Material		$4\mathrm{W}620$
Rel. permeability	${\mu }_{\mathrm{S}}$	620

. The impedance of the sensor in the air and above an aluminum plate was measured using the precision LCR meter HP 4284A in the frequency range from 5 kHz to 100 kHz and for lift-off values from 0 $\mathsf{\mu }\mathrm{m}$ to 1000 $\mathsf{\mu }\mathrm{m}$ in 100 $\mathsf{\mu }\mathrm{m}$ steps. The plate was made of aluminium alloy EN AW-6082 with a conductivity of $42.2\%\mathrm{IACS}$ or $24.56$ $\mathrm{M}$$\mathrm{S}$/m. The impedance of the coil Z is the sum of the coil impedance in air (without the plate) $Z_0=j\omega L_0$, and the impedance change $\Delta Z=\Delta R+j\omega \Delta L$ due to the proximity of the plate:

$$Z=j\omega L_0+\Delta Z,$$

(72)

where the resistance of the coil in air and the coil parasitic capacitance have been compensated for.

The inductance of the sensor in the air was calculated $1525.67$ $\mathsf{\mu }\mathrm{H}$, while the measured value was from $1524.48$ $\mathsf{\mu }\mathrm{H}$ to $1526.38$ $\mathsf{\mu }\mathrm{H}$ depending on the frequency. The impedance change obtained from the model and the experiment are shown in Figure 6 and Figure 7 for $\Delta R$ and $\Delta L$, respectively. The model and experiment are in practically acceptable agreement with the relative error of $\Delta R$ within $\pm 3\%$ and $\Delta L$ within $\pm 1\%$. These errors are comparable to the ones in other studies on eddy current sensor modeling. In practice, model-based measurements of conductivity and lift-off require calibration procedures with samples of known conductivity and lift-off if one wants to achieve measurement uncertainty better than the reported modeling errors would allow.

6. Sensitivity to Sensor Dimensions

The data in Table 2 was obtained by averaging and rounding to $0.01$ mm the values measured using a Mitutoyo micrometer with resolution of $0.001$ mm and accuracy of $\pm 2 \mathsf{\mu }\mathrm{m}$. We rounded the values because at a higher resolution, the measurement is affected by the applied pressure or measurement position along the measured object, e.g., the ferrite core was not a perfect cylinder at a resolution of $0.001$ mm. We used purposefully overestimated measurement uncertainty of $\pm 0.01$ mm to evaluate its effect on the sensor impedance. With that aim, we calculated the sensitivity of the sensor impedance to the sensor dimensions.

Table 3. Relative uncertainty of sensor impedance in the case of $\pm 0.01$ mm uncertainty in sensor dimensions. The plate has $\sigma = 5 \mathrm{M}\mathrm{S}/\mathrm{m}$ and ${\mu }_{\mathrm{r}}=50$, lift-off is 0 $\mathsf{\mu }\mathrm{m}$ and the frequency is 60 kHz.

Parameter		Value	Uncertainty		Relative Uncertainty of
				${\mathit{L}}_{\mathbf{0}}$	$\Delta \mathit{R}$	$\Delta \mathit{L}$
Coil
Inner radius	$r_{\mathrm{T}1}$	$5.00$ mm	$\pm 0.01$ mm	$\pm 2.6 \mathrm{ppm}$	$\pm 2.6 \mathrm{ppm}$	$\pm 2.6 \mathrm{ppm}$
Outer radius	$r_{\mathrm{T}2}$	$5.09$ mm	$\pm 0.01$ mm	$\pm 2.5 \mathrm{ppm}$	$\pm 2.4 \mathrm{ppm}$	$\pm 2.4 \mathrm{ppm}$
Bottom position	$z_{\mathrm{T}1}$	$2.45$ mm	$\pm 0.01$ mm	$\pm 542 \mathrm{ppm}$	$\mp 872 \mathrm{ppm}$	$\mp 874 \mathrm{ppm}$
Top position	$z_{\mathrm{T}2}$	$17.45$ mm	$\pm 0.01$ mm	$\mp 482 \mathrm{ppm}$	$\mp 856 \mathrm{ppm}$	$\mp 856 \mathrm{ppm}$
Core
Radius	$r_{\mathrm{C}}$	$4.98$ mm	$\pm 0.01$ mm	$\pm 6776 \mathrm{ppm}$	$\pm 5632 \mathrm{ppm}$	$\pm 7646 \mathrm{ppm}$
Height	$L_{\mathrm{C}}$	$20.03$ mm	$\pm 0.01$ mm	$\pm 724 \mathrm{ppm}$	$\pm 879 \mathrm{ppm}$	$\pm 788 \mathrm{ppm}$
Shield
Inner radius	$r_{\mathrm{S}1}$	$6.45$ mm	$\pm 0.01$ mm	$\mp 4858 \mathrm{ppm}$	$\mp 4030 \mathrm{ppm}$	$\mp 6933 \mathrm{ppm}$
Outer radius	$r_{\mathrm{S}2}$	$12.92$ mm	$\pm 0.01$ mm	$\pm 8.9 \mathrm{ppm}$	$\pm 30.4 \mathrm{ppm}$	$\pm 65 \mathrm{ppm}$
Height	$L_{\mathrm{S}}$	$28.29$ mm	$\pm 0.01$ mm	$\pm 0.2 \mathrm{ppm}$	$\pm 0.2 \mathrm{ppm}$	$\pm 0.6 \mathrm{ppm}$

shows the relative changes of the inductance in the air $L_0$ and impedance change $\Delta R$ and $\Delta L$ (for plate with $\sigma =5 \mathrm{M}\mathrm{S}/\mathrm{m}$ and ${\mu }_{\mathrm{r}}=50$ and lift-off 0 $\mathsf{\mu }\mathrm{m}$) if each of the nine sensor dimensions is independently varied within the interval of $\pm 0.01$ mm. The core radius and the shield inner radius uncertainties have the most significant contribution to the impedance uncertainty. This means that the gap between the core and the shield must be precisely controlled during the sensor production.

7. Conclusions

The validity of the impedance model is corroborated by a practically negligible discrepancy in comparison with the FEM approach. We compared the model and FEM approach for a range of practical scenarios: changes of the impedance due to the variations in the core and shield permeabilities, dependence of the sensor impedance on the plate conductivity of magnetic and nonmagnetic plates and impedance dependence on the sensor lift-off above magnetic and nonmagnetic plates. The numerical implementation of the model is very fast in comparison to the FEM approach ($0.7$ $\mathrm{s}$ vs. 30 min), and this can be improved further if the intermediate results are reused for the next iteration. In comparison with the impedance measurement, the relative error better than $3\%$ in the frequency range from 5 kHz to 100 kHz confirms the applicability of the model for eddy current measurement of material properties and lift-off. The model can be extended easily to include multi-layered planar structures.