Impedance Variation of a Coaxial Coil Encircling a Metal Tube Adapter

1. Introduction

A tube adapter is a component connecting two tubes of different diameters. The standard analytical method of Dodd and Deeds [1] is unable to investigate the interaction of an induction coil with a metal tube adapter, due to the end effects involved in this problem. The truncated region eigenfunction expansion (TREE) method, pioneered by Hannakam and Tepe [2], and developed by Theodoulidis, Kriezis and Bowler [3-9] for modeling of the eddy current nondestructive testing (EC NDT), is capable of analyzing the end effects and establishing analytical models. However, successful implementation of TREE for the model of end effects depends on the solution of relevant eigenvalue equations, which are transcendental and complex roots should be determined. Conventionally, the Newton-Raphson algorithm [10-13] or contour integral based on the Cauchy’s theorem [14-17] are applied to solve the eigenvalue equations. A novel method based on the Sturm-Liouville theory and Galerkin approach has been proposed recently [18, 19], which greatly simplifies the process of locating the complex eigenvalues.

Unfortunately, the TREE cannot be used immediately for the problem of this work, even with the novel approach. The bottleneck lies in the numerical issues for the regions composed of multiple air-metal subdomains, more specifically, the numerical overflow occurs constantly when the symbolic eigenfunctions are employed for the n-subdomain ( ) regions. Admittedly, for certain three-subdomain problems, the overflow can be suppressed by splitting the model into even and odd parts [6, 8, 20-22], but the technique is merely available for the problems with additional symmetry.

More recently, the difficulties of overflow have been overcome in [23] by an innovative approach based on the FEM solution of the Sturm-Liouville equation. By this method, no numerical overflow arises in the solving process, therefore, many new models can be dealt with the TREE. According to [23], the Clenshaw-Curtis quadrature rule is appropriate for the evaluations of the integrals involving the numerical eigenfunctions, to obtain an efficient algorithm. The investigation of this work will be based on the approach.

In Section 2, the TREE solution is given for the metal tube adapter surrounding by a coaxial coil. The permeability of the metal is not restricted to μ0 in comparison with [23]. In Section 3, the numerical eigenfunctions are obtained by the FEM solution of the Sturm-Liouville equations and normalized, and the Clenshaw-Curtis quadrature rule is applied to the computation of the integrals involving the numerical eigenfunctions. In Section 4, the TREE results are compared with those from the FEM simulation.

2. Formulation

A metal tube adapter of conductivity σ and permeability $\mu ={\mu }_r{\mu }_0$ (μ_r is supposed to be constant), is encircled by a coaxial induction coil excited by a time harmonic current of frequency ω and amplitude I (See Figure 1). The geometry of the coil and tube adapter is shown in Figure 2. Perfect electric boundary is imposed at z=0 and z=b to discretize the eigenvalues of this boundary value problem (BVP).

The solution domain is divided into 5 regions along the r-axis (See Figure 2). The vector potentials A₁ to A₅ satisfy the Laplace or Helmholtz equations in the corresponding regions:

$${\nabla }^2{\mathit{A}}_{1,5}=0$$

(1)

$${\nabla }^2{\mathit{A}}_{2,3,4}=k^2{\mathit{A}}_{2,3,4}$$

(2)

where $k=\sqrt{i\omega \sigma {\mu }_0{\mu }_r}$ is the wavenumber of the metal.

Only the φ-component of the vector potential exists, due to the axisymmetry of the BVP, i.e. $\mathit{A}=A{\mathit{e}}_{\phi }$, and the vector Laplacian of equations (1) and (2) is reduced to

$${\nabla }_{\phi }^2=\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\frac{\partial }{\partial r}-\frac{1}{r^2}+\frac{{\partial }^2}{\partial z^2}$$

(3)

2.1. Vector Potential of the Source Coil

The formulation of source vector potential can be obtained by the source expansion of the Poisson equation [24, 25]. The vector potential of the coil can be written in the form (See Figure 3)

$$A_I\left(r,z\right)={\mathbf{S}}^{\mathrm{T}}\left(z\right){\mathbf{I}}_1\left(\mathsf{\alpha }r\right){\mathbf{C}}_1^{\left(e\right)}$$

(4a)

$$A_{II}\left(r,z\right)={\mathbf{S}}^{\mathrm{T}}\left(z\right)\left[{\mathbf{I}}_1\left(\mathsf{\alpha }r\right)C_2^{\left(e\right)}+{\mathbf{K}}_1\left(\mathsf{\alpha }r\right){\mathbf{D}}_2^{\left(e\right)}+\mathbf{V}\left(r\right)\right]$$

(4b)

$$A_{III}\left(r,z\right)={\mathbf{S}}^{\mathrm{T}}\left(z\right){\mathbf{K}}_1\left(\mathsf{\alpha }r\right){\mathbf{D}}_3^{\left(e\right)}$$

(4c)

where the source vector V(r) is

$$\mathbf{V}\left(r\right)=\left[\begin{array}{c}{v}_1\left(r\right)\\ v_2\left(r\right)\\ ⋮\end{array}\right]$$

with the elements

$$v_i\left(r\right)={\kappa }_i{\mathbf{L}}_1\left({\alpha }_{i}r\right)$$

(5)

where ${\alpha }_i=i\pi /b$, and Ln(x) is the modified Struve function of order n, and

$${\kappa }_i=\frac{2{\mu }_{0}J}{i{\alpha }_i^2}\sin \left[\frac{{\alpha }_i}{2}\left(z_1-z_2\right)\right]\sin \left[\frac{{\alpha }_i}{2}\left(z_1+z_2\right)\right]$$

(6)

Other matrices and vectors in (4a)-(4c) are

$$\mathsf{\alpha }=\left[\begin{array}{ccc}{\alpha }_1& 0& \cdots \\ 0& {\alpha }_2& \cdots \\ ⋮& ⋮& \ddots \end{array}\right],{\mathbf{I}}_1\left(\mathsf{\alpha }r\right)=\left[\begin{array}{ccc}{I}_1\left({\alpha }_{1}r\right)& 0& \cdots \\ 0& I_1\left({\alpha }_{2}r\right)& \cdots \\ ⋮& ⋮& \ddots \end{array}\right],{\mathbf{K}}_1\left(\mathsf{\alpha }r\right)=\left[\begin{array}{ccc}{K}_1\left({\alpha }_{1}r\right)& 0& \cdots \\ 0& K_1\left({\alpha }_{2}r\right)& \cdots \\ ⋮& ⋮& \ddots \end{array}\right],\mathbf{S}\left(z\right)=\left[\begin{array}{c}\sin \left({\alpha }_{1}z\right)\\ \sin \left({\alpha }_{2}z\right)\\ ⋮\end{array}\right],$$

where I_n(x) and K_n(x) are the modified Bessel functions of the first and second kinds of order n, respectively, and ${\mathbf{C}}_1^{\left(e\right)}$，$C_2^{\left(e\right)}$，$D_2^{\left(e\right)}$，${\mathbf{D}}_3^{\left(e\right)}$ are coefficients to be determined. With the interface conditions of B_r and H_z at r=r₁ and r=r₂, the coefficients can be found:

$$C_{1,i}^{\left(e\right)}={\kappa }_i\left[\chi \left({\alpha }_i{r}_1\right)-\chi \left({\alpha }_i{r}_2\right)\right]$$

(7a)

$$C_{2,i}^{\left(e\right)}=-{\kappa }_i\chi \left({\alpha }_i{r}_2\right)$$

(7b)

$$D_{2,i}^{\left(e\right)}={\kappa }_i\eta \left({\alpha }_i{r}_1\right)$$

(7c)

$$D_{3,i}^{\left(e\right)}={\kappa }_i\left[\eta \left({\alpha }_i{r}_1\right)-\eta \left({\alpha }_i{r}_2\right)\right]$$

(7d)

where

$$\chi \left(x\right)=x\left[K_1\left(x\right){\mathbf{L}}_0\left(x\right)+K_0\left(x\right){\mathbf{L}}_1\left(x\right)\right]$$

(8a)

$$\eta \left(x\right)=x\left[I_1\left(x\right){\mathbf{L}}_0\left(x\right)-I_0\left(x\right){\mathbf{L}}_1\left(x\right)\right]$$

(8b)

For the function χ(x) used for the subsequent analysis, it is advisable to adopt an alternative form for the practical evaluations, namely

$$\chi \left(x\right)=\{\begin{array}{l}\frac{x^2}{2}{\displaystyle \sum _{m=0}^{m_0}\frac{{\left(x/2\right)}^{2m}}{\Gamma \left(m+3/2\right)}\left[\frac{K_1\left(x\right)}{\Gamma \left(m+3/2\right)}+\frac{x{K}_0\left(x\right)}{2\Gamma \left(m+5/2\right)}\right],x<15}\\ 1+\frac{x}{{\pi }^2}{\displaystyle \sum _{m=0}^{m_1}\frac{{\Gamma }^2\left(m+1/2\right)}{{\left(x/2\right)}^{2m}}\left[\frac{K_0\left(x\right)}{m-1/2}-\frac{K_1\left(x\right)}{x/2}\right],x\ge 15}\end{array}$$

(9)

Expression (9) is obtained by the Maclaurin and asymptotic expansions of L_n(x) [26], and high accuracy can be achieved by setting m₀=23 and m₁=10, respectively.

2.2. Impedance Change of the Coil Encircling the Metal Tube Adapter

The vector potentials in the 5 regions of Figure 2 are expansible by the separation of variables

$$A_1\left(r,z\right)=S^{\mathrm{T}}\left(z\right)I_1\left(\mathsf{\alpha }r\right)C_1$$

(10a)

$$A_2\left(r,z\right)=F^{\mathrm{T}}\left(z\right)\left[I_1\left(P_{1}r\right)C_2+K_1\left(P_{1}r\right)D_2\right]$$

(10b)

$$A_3\left(r,z\right)=G^{\mathrm{T}}\left(z\right)\left[I_1\left(P_{2}r\right)C_3+K_1\left(P_{2}r\right)D_3\right]$$

(10c)

$$A_4\left(r,z\right)=H^{\mathrm{T}}\left(z\right)\left[I_1\left(P_{3}r\right)C_4+K_1\left(P_{3}r\right)D_4\right]$$

(10d)

$$A_5\left(r,z\right)=S^{\mathrm{T}}\left(z\right)\left[I_1\left(\mathsf{\alpha }r\right)C_1^{\left(e\right)}+K_1\left(\mathsf{\alpha }r\right)D^{\left(s\right)}\right]$$

(10e)

where P₁, P₂ and P₃ are the eigenvalue matrices of the Regions 2, 3 and 4, respectively,

$$P_1=\left[\begin{array}{ccc}{p}_{1,1}& 0& \cdots \\ 0& p_{1,2}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right],{\mathbf{P}}_2=\left[\begin{array}{ccc}{p}_{2,1}& 0& \cdots \\ 0& p_{2,2}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right],{\mathbf{P}}_3=\left[\begin{array}{ccc}{p}_{3,1}& 0& \cdots \\ 0& p_{3,2}& \cdots \\ ⋮& ⋮& \ddots \end{array}\right],$$

and

$$F\left(z\right)=\left[\begin{array}{c}{f}_1\left(p_{1,1},z\right)\\ f_2\left(p_{1,2},z\right)\\ ⋮\end{array}\right],G\left(z\right)=\left[\begin{array}{c}{g}_1\left(p_{2,1},z\right)\\ g_2\left(p_{2,2},z\right)\\ ⋮\end{array}\right],\mathbf{H}\left(z\right)=\left[\begin{array}{c}{h}_1\left(p_{3,1},z\right)\\ h_2\left(p_{3,2},z\right)\\ ⋮\end{array}\right]$$

are the axial eigenfunctions satisfying the relevant Sturm-Liouville equations:

$$\frac{d^2{f}_i\left(z\right)}{d{z}^2}-k_1^2\left(z\right)f_i\left(z\right)=-p_{1,i}^2{f}_i\left(z\right),f_i\left(0\right)=f_i\left(b\right)=0$$

(11a)

$$\frac{d^2{g}_i\left(z\right)}{d{z}^2}-k_2^2\left(z\right)g_i\left(z\right)=-p_{2,i}^2{g}_i\left(z\right),g_i\left(0\right)=g_i\left(b\right)=0$$

(11b)

and

$$\frac{d^2{h}_i\left(z\right)}{d{z}^2}-k_3^2\left(z\right)h_i\left(z\right)=-p_{3,i}^2{h}_i\left(z\right),h_i\left(0\right)=h_i\left(b\right)=0$$

(11c)

with

$$k_1\left(z\right)=\{\begin{array}{l}k,b_1\le z\le b_3\\ 0,\mathrm{others}\end{array}$$

(12a)

$$k_2\left(z\right)=\{\begin{array}{l}k,b_2\le z\le b_3\\ 0,\mathrm{others}\end{array}$$

(12b)

and

$$k_3\left(z\right)=\{\begin{array}{l}k,b_2\le z\le b_4\\ 0,\mathrm{others}\end{array}$$

(12c)

Taking account of the interface conditions of Br and Hz at r=a₁, r=a₂, r=a₃ and r=a₄, the following equations for the coefficients $C_2$, $C_3$, $C_4$, $D_1$, $D_2$ and $D_3$ can be derived

$$\frac{b}{2}{I}_1\left(\alpha a_1\right)C_1=T_1\left[I_1\left(P_1{a}_1\right)C_2+K_1\left(P_1{a}_1\right)D_1\right]$$

(13a)

$$T_1^{\mathrm{T}}\alpha I_0\left(\alpha a_1\right)C_1=P_1\left[I_0\left(P_1{a}_1\right)C_2-K_0\left(P_1{a}_1\right)D_2\right]$$

(13b)

$$I_1\left(P_1{a}_2\right)C_2+K_1\left(P_1{a}_2\right)D_2=T_2\left[I_1\left(P_2{a}_2\right)C_3+K_1\left(P_2{a}_2\right)D_3\right]$$

(13c)

$$T_2^{\mathrm{T}}{P}_1\left[I_0\left(P_1{a}_2\right)C_2-K_0\left(P_1{a}_2\right)D_2\right]=P_2\left[I_0\left(P_2{a}_2\right)C_3-K_0\left(P_2{a}_2\right)D_3\right]$$

(13d)

$$I_1\left(P_2{a}_3\right)C_3+K_1\left(P_2{a}_3\right)D_3=T_3\left[I_1\left(P_3{a}_3\right)C_4+K_1\left(P_3{a}_3\right)D_4\right]$$

(13e)

$$T_3^{\mathrm{T}}{P}_2\left[I_0\left(P_2{a}_3\right)C_3-K_0\left(P_2{a}_3\right)D_3\right]=P_3\left[I_0\left(P_3{a}_3\right)C_4-K_0\left(P_3{a}_3\right)D_4\right]$$

(13f)

$$T_4\left[I_1\left(P_3{a}_4\right)C_4+K_1\left(P_3{a}_4\right)D_4\right]=\frac{b}{2}\left[I_1\left(\mathsf{\alpha }{a}_4\right)C_1^{\left(e\right)}+K_1\left(\mathsf{\alpha }{a}_4\right)D^{\left(s\right)}\right]$$

(13g)

$$P_3\left[I_0\left(P_3{a}_4\right)C_4-K_0\left(P_3{a}_4\right)D_4\right]=T_4^{\mathrm{T}}\mathsf{\alpha }\left[I_0\left(\mathsf{\alpha }{a}_4\right)C_1^{\left(e\right)}-K_0\left(\mathsf{\alpha }{a}_4\right)D^{\left(s\right)}\right]$$

(13h)

where

$$T_1={\displaystyle {\int }_0^{b}S\left(z\right)F^{\mathrm{T}}\left(z\right)}dz$$

(14)

$$T_2={\displaystyle {\int }_0^b\frac{1}{{\mu }_r^{\left(1\right)}\left(z\right)}F\left(z\right)G^{\mathrm{T}}\left(z\right)}dz$$

(15)

$$T_3={\displaystyle {\int }_0^b\frac{1}{{\mu }_r^{\left(2\right)}\left(z\right)}G\left(z\right)H^{\mathrm{T}}\left(z\right)}dz$$

(16)

$$T_4={\displaystyle {\int }_0^{b}S\left(z\right)H^{\mathrm{T}}\left(z\right)}dz$$

(17)

In (13a)-(13h), the orthogonalities of the eigenfunctions

$${\displaystyle {\int }_0^{b}S\left(z\right)S^{\mathrm{T}}\left(z\right)}dz=\frac{b}{2}I$$

(18a)

$${\displaystyle {\int }_0^b\frac{1}{{\mu }_r^{\left(1\right)}\left(z\right)}F\left(z\right)F^{\mathrm{T}}\left(z\right)}dz=I$$

(18b)

$${\displaystyle {\int }_0^b\frac{1}{{\mu }_r^{\left(2\right)}\left(z\right)}G\left(z\right)G^{\mathrm{T}}\left(z\right)}dz=I$$

(18c)

$${\displaystyle {\int }_0^b\frac{1}{{\mu }_r^{\left(3\right)}\left(z\right)}H\left(z\right)H^{\mathrm{T}}\left(z\right)}dz=I$$

(18d)

have been adopted, where I is the identity matrix, and

$${\mu }_r^{\left(1\right)}\left(z\right)=\{\begin{array}{l}{\mu }_r,b_1\le z\le b_3\\ 1,\mathrm{others}\end{array}$$

(19a)

$${\mu }_r^{\left(2\right)}\left(z\right)=\{\begin{array}{l}{\mu }_r,b_2\le z\le b_3\\ 1,\mathrm{others}\end{array}$$

(19b)

$${\mu }_r^{\left(3\right)}\left(z\right)=\{\begin{array}{l}{\mu }_r,b_2\le z\le b_4\\ 1,\mathrm{others}\end{array}$$

(19c)

The orthonormalization relations of (18b)-(18d) will be expounded in Section III.

Matrix algebra of (13a)-(13h) yields the equation system

$$\left[\begin{array}{cccc}{A}_{11}& A_{12}& \mathbf{0}& \mathbf{0}\\ A_{21}& A_{22}& A_{23}& A_{42}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ \mathbf{0}& \mathbf{0}& A_{43}& A_{44}\end{array}\right]\left[\begin{array}{c}{C}_2\\ D_2\\ C_4\\ D_4\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ E\end{array}\right]$$

(20)

where

$$A_{11}=U_1{I}_1\left(P_1{a}_1\right)-P_1{I}_0\left(P_1{a}_1\right)$$

(21a)

$$A_{12}=U_1{K}_1\left(P_1{a}_1\right)+P_1{K}_0\left(P_1{a}_1\right)$$

(21b)

$$A_{21}=I_1\left(P_2{a}_3\right)M_3+K_1\left(P_2{a}_3\right)M_1$$

(21c)

$$A_{22}=I_1\left(P_2{a}_3\right)M_4+K_1\left(P_2{a}_3\right)M_2$$

(21d)

$$A_{23}=-T_3{I}_1\left(P_3{a}_3\right)$$

(21e)

$$A_{24}=-T_3{K}_1\left(P_3{a}_3\right)$$

(21f)

$$A_{31}=T_3^{\mathrm{T}}{P}_2\left[I_0\left(P_2{a}_3\right)M_3-K_0\left(P_2{a}_3\right)M_1\right]$$

(21g)

$$A_{32}=T_3^{\mathrm{T}}{P}_2\left[I_0\left(P_2{a}_3\right)M_4-K_0\left(P_2{a}_3\right)M_2\right]$$

(21h)

$$A_{33}=-P_3{I}_0\left(P_3{a}_3\right)$$

(21i)

$$A_{34}=P_3{K}_0\left(P_3{a}_3\right)$$

(21j)

$$A_{43}=U_2{I}_1\left(P_3{a}_4\right)+P_3{I}_0\left(P_3{a}_4\right)$$

(21k)

$$A_{44}=U_2{K}_1\left(P_3{a}_4\right)-P_3{K}_0\left(P_3{a}_4\right)$$

(21l)

$$E=T_4^{\mathrm{T}}\mathsf{\alpha }\left[I_0\left(\mathsf{\alpha }{a}_4\right)+K_0\left(\mathsf{\alpha }{a}_4\right)K_1^{-1}\left(\mathsf{\alpha }{a}_4\right)I_1\left(\mathsf{\alpha }{a}_4\right)\right]C_1^{\left(e\right)}$$

(21m)

with

$$U_1=\frac{2}{b}{T}_1^{\mathrm{T}}\mathsf{\alpha }{I}_0\left(\mathsf{\alpha }{a}_1\right)I_1^{-1}\left(\mathsf{\alpha }{a}_1\right)T_1$$

(22a)

$$U_2=\frac{2}{b}{T}_4^{\mathrm{T}}\mathsf{\alpha }{K}_0\left(\mathsf{\alpha }{a}_4\right)K_1^{-1}\left(\mathsf{\alpha }{a}_4\right)T_4$$

(22b)

$$M_1=a_2\left[X_1{I}_1\left(P_1{a}_2\right)-X_2{I}_0\left(P_1{a}_2\right)\right]$$

(22c)

$$M_2=a_2\left[X_1{K}_1\left(P_1{a}_2\right)+X_2{K}_0\left(P_1{a}_2\right)\right]$$

(22d)

$$M_3=a_2\left[X_3{I}_1\left(P_1{a}_2\right)+X_4{I}_0\left(P_1{a}_2\right)\right]$$

(22e)

$$M_4=a_2\left[X_3{K}_1\left(P_1{a}_2\right)-X_4{K}_0\left(P_1{a}_2\right)\right]$$

(22f)

$$X_1=P_2{I}_0\left(P_2{a}_2\right)T_2^{-1}$$

(22g)

$$X_2=I_1\left(P_2{a}_2\right)T_2^{\mathrm{T}}{P}_1$$

(22h)

$$X_3=P_2{K}_0\left(P_2{a}_2\right)T_2^{-1}$$

(22i)

$$X_4=K_1\left(P_2{a}_2\right)T_2^{\mathrm{T}}{P}_1$$

(22j)

Solving equation (20) will give the coefficients C₂, D₂, C₄ and D₄, and other coefficients can be found by

$$C_1=\frac{2}{b}{I}_1^{-1}\left(\alpha a_1\right)T_1\left[I_1\left(P_1{a}_1\right)C_2+K_1\left(P_1{a}_1\right)D_2\right]$$

(23)

$$C_3=M_3{C}_2+M_4{D}_2$$

(24)

$$D_3=M_1{C}_2+M_2{D}_2$$

(25)

The coefficient required for calculation of ΔZ is

$$D^{\left(s\right)}=K_1^{-1}\left(\mathsf{\alpha }{a}_4\right)\left\{-I_1\left(\mathsf{\alpha }{a}_4\right)C_1^{\left(e\right)}+\frac{2}{b}{T}_4\left[I_1\left(P_3{a}_4\right)C_4+K_1\left(P_3{a}_4\right)D_4\right]\right\}$$

(26)

Accordingly, the coil impedance variation is given by

$$\begin{array}{l}\Delta Z=\frac{i\omega }{I^2}{\displaystyle {\int }_V{\mathit{A}}^{\left(s\right)}\cdot \mathit{J}dV}\\ =\frac{{\pi }^{2}i\omega N^2}{{\left(r_2-r_1\right)}^2{\left(z_2-z_1\right)}^2}{\displaystyle \sum _{n=1}^{\infty }\frac{\cos \left({\alpha }_n{z}_1\right)-\cos \left({\alpha }_n{z}_2\right)}{{\alpha }_n^3}}\left[\chi \left({\alpha }_n{r}_2\right)-\chi \left({\alpha }_n{r}_1\right)\right]d_n^{\left(s\right)}\end{array}$$

(27)

where the current density J has been omitted (letting J=1) to simplify the expression.

3. Eigenfunctions and the Associated Integrals of the Multi-Subdomain Regions

The conventional TREE method is limited to the two-subdomain problems, apart from certain three-subdomain problems which can be split in the even and odd parts. In general, it is impossible to evade the numerical overflow of the symbolic eigenfunctions of Regions 2, 3 and 4. Accordingly, the approach of [23] will be applied to cope with F(z), G(z) and H(z). In accord with [18, 19], the eigenvalues of (11a) can be obtained by the solution of a generalized eigenvalue equation

$$\mathbf{K}{\mathbf{U}}_i=p_{1,i}^2\mathbf{W}{\mathbf{U}}_i$$

(28)

where K is the stiffness matrix with the elements

$${\mathbf{K}}_{mn}={\displaystyle \underset{0}{\overset{b}{\int }}\frac{1}{{\mu }_r^{\left(1\right)}\left(z\right)}[\frac{d{\phi }_m\left(z\right)}{dz}\frac{d{\phi }_n\left(z\right)}{dz}}+k_1^2\left(z\right){\phi }_m\left(z\right){\phi }_n\left(z\right)]dz$$

(29)

and W is the damping matrix of the elements

$${\mathbf{W}}_{mn}={\displaystyle \underset{0}{\overset{b}{\int }}\frac{{\phi }_m\left(z\right){\phi }_n\left(z\right)}{{\mu }_r^{\left(1\right)}\left(z\right)}dz}$$

(30)

where φ_m and φ_n are FEM functions consisting of the Lagrange polynomials defined on the reference interval $-1\le \xi \le 1$ (the shape functions).

Sparse matrix K will be generated from the FEM basis. Hence, equation (28) can be solved by the efficient algorithm such as Arnoldi iteration [27]. The solution provides both the eigenvalues p_1,i and the eigenvectors U_i, which are the discrete eigenfunctions $f_i\left(z\right)$. Moreover, denoting

$$\mathbf{U}={\left[{\mathbf{U}}_1,{\mathbf{U}}_2,\mathrm{...}\right]}^{\mathrm{T}}$$

(31)

and by virtue of the vector normalization

$${\mathbf{U}}^{\prime }=\frac{\mathbf{U}}{\sqrt{\mathrm{diag}\left(\mathbf{U}\mathbf{W}{\mathbf{U}}^{\mathrm{T}}\right)}}$$

(32)

the eigenfunction normalization

$${\displaystyle \underset{0}{\overset{b}{\int }}\frac{1}{{\mu }_r^{\left(1\right)}}{f}_i^2\left(z\right)dz}=1$$

(33)

can be established automatically. Equations (32) and (33) can be validated by inspecting the diagonal entries of $\mathbf{U}\mathbf{W}{\mathbf{U}}^{\mathrm{T}}$, and taking equation (30) into account. Consequently, the orthonormality of (18b)-(18d) can be established.

The requirement of the accurate and efficient algorithm leads to the choice of high order Lagrange polynomials for the FEM basis. Here, we choose the cubic Lagrange polynomials

$$\{\begin{array}{l}{N}_0\left(\xi \right)=-\frac{1}{16}\left(\xi -1\right)\left(3\xi -1\right)\left(3\xi +1\right)\\ N_1\left(\xi \right)=\frac{9}{16}\left(\xi -1\right)\left(\xi +1\right)\left(3\xi -1\right)\\ N_2\left(\xi \right)=-\frac{9}{16}\left(\xi -1\right)\left(\xi +1\right)\left(3\xi +1\right)\\ N_3\left(\xi \right)=\frac{1}{16}\left(\xi +1\right)\left(3\xi -1\right)\left(3\xi +1\right)\end{array}$$

(34)

The cubic interpolation of the eigenfunction is

$$f_i\left(z\right)={\displaystyle \sum _{e=0}^3{u^{\prime }}_{l+e}{N}_e\left(z\right)}$$

(35)

where ${u^{\prime }}_{l+e}$ are successive 4 entries of ${{\mathbf{U}}^{\prime }}_i$, N_e(z) is obtained by (34) with the change of variable

$$\xi =\frac{2z-z_a-z_b}{z_b-z_a}$$

(36)

where z_a and z_b are the mesh nodes corresponding to the reference interval. The numerical overflow of $f_i\left(z\right)$ is eliminated by this procedure. They are consequently well adapted for the subsequent integral computation. Furthermore, it appears to be very effective to evaluate directly the integrals (14)-(17) with the Clenshaw-Curtis quadrature rule, which is quoted here for completeness [28-30]

$${\displaystyle {\int }_{-1}^{1}f\left(x\right)}\mathrm{d}x={\displaystyle \sum _{k=0}^n{w}_{k}f\left(x_k\right)}$$

(37)

where the weights w_k are given by

$$w_k=\frac{g_k}{n}\left(1-{\displaystyle \sum _{j=1}^{\lfloor n/2\rfloor }\frac{b_j}{4j^2-1}\cos \left(2jk\pi /n\right)}\right)$$

(38)

and the quadrature nodes are

$$x_k=\cos \left(\frac{k\pi }{n}\right),k=0,1,\mathrm{...},n$$

(39)

with

$$g_k=\{\begin{array}{l}1,k=0,n\\ 2,\mathrm{otherwise}\end{array}$$

(40)

$$b_j=\{\begin{array}{l}1,j=\frac{1}{2}n\\ 2,\mathrm{otherwise}\end{array}$$

(41)

It follows from equations (37)-(41) that the matrix elements of T₁ can be computed by

$$\begin{array}{l}{\mathbf{T}}_{1,ij}={\displaystyle \underset{0}{\overset{b}{\int }}\sin \left({\alpha }_{i}z\right)f_j\left(z\right)dz}\\ =\frac{b}{2}{\displaystyle {\int }_{-1}^1\sin \left[\frac{b}{2}{\alpha }_i\left(1+x\right)\right]f_j\left[\frac{b}{2}\left(1+x\right)\right]}dx\\ =\frac{b}{2}{\displaystyle \sum _{k=0}^n{w}_k\sin \left({\alpha }_i{z}_k\right)f_j\left(z_k\right)}\end{array}$$

(42)

where

$$z_k=\frac{b}{2}\left(1+x_k\right)$$

(43)

The matrix elements of T₂ are likewise given by

$$\begin{array}{l}{\mathbf{T}}_{2,ij}={\displaystyle \underset{0}{\overset{b}{\int }}\frac{1}{{\mu }_r^{\left(1\right)}}{f}_i\left(z\right)g_j\left(z\right)dz}\\ ={\displaystyle \underset{0}{\overset{b_1}{\int }}{f}_i\left(z\right)g_j\left(z\right)dz}+\frac{1}{{\mu }_r}{\displaystyle \underset{b_1}{\overset{b_3}{\int }}{f}_i\left(z\right)g_j\left(z\right)dz}+{\displaystyle \underset{b_3}{\overset{b}{\int }}{f}_i\left(z\right)g_j\left(z\right)dz}\\ =\frac{b_1}{2}{\displaystyle \sum _{k=0}^n{w}_k{f}_i\left(z_k^{\left(1\right)}\right)g_j\left(z_k^{\left(1\right)}\right)}+\frac{b_3-b_1}{2}{\displaystyle \sum _{k=0}^n{w}_k{f}_i\left(z_k^{\left(2\right)}\right)g_j\left(z_k^{\left(2\right)}\right)}+\frac{b-b_3}{2}{\displaystyle \sum _{k=0}^n{w}_k{f}_i\left(z_k^{\left(3\right)}\right)g_j\left(z_k^{\left(3\right)}\right)}\end{array}$$

(44)

where

$$z_k^{\left(1\right)}=\frac{b_1}{2}\left(1+x_k\right),z_k^{\left(2\right)}=\frac{b_3+b_1}{2}+\frac{b_3-b_1}{2}{x}_k,z_k^{\left(3\right)}=\frac{b+b_3}{2}+\frac{b-b_3}{2}{x}_k$$

(45)

The same analysis is also applicable to the matrix elements of T₃ and T₄.

4. Numerical Validation

The proposed method will be verified with the parameters of the metal tube adapter and the induction coil given in Table 1, Table 2 and Table 3. The nonmagnetic alloy UNS (Unified Numbering System) C96400 (70-30 Copper-Nickel), and the magnetic stainless steels S31600 (austenitic) and S32760 (super duplex) [31] are used for the numerical validation. The coil impedance variations are calculated and plotted for these metal materials, with different coil positions. The TREE results are compared with those from the FEM simulation of Comsol Multiphysics^®, shown in Figure 4, where the theoretical and FEM data are denoted by solid lines and circles, respectively. The reactance of the isolated induction coil is X₀=ωL₀, with L₀=4.104132mH, which can be found by the method such as in [32, 33].

Further calculations are carried out for the coil impedance variation with respect to the frequencies. For the alloys of lower μ_r (C96400 and S31600), the calculation frequency ranges from 1kHz to 100kHz, for higher μ_r (S32760) the frequency interval [100Hz, 10kHz] is chosen. The results are shown in Figure 5 and Figure 6, where the TREE data are plotted by solid lines, in connection with the circles representing the data of FEM simulation. Other parameters are referred to the Table 2 and 3.

Very good agreement is obtained between the TREE and FEM results in the numerical comparisons. The calculations were implemented on a personal computer of a 4.2 GHz processor (Intel® Core i7-7700K) and 16 GB RAM. Additional algorithm details are shown in Table 4, where the frequencies, summation terms (matrix size), mesh elements and quadrature nodes used in the computation are listed. The execution time of the eigenvalue and eigenfunction computation, and the total execution time of TREE evaluation are also provided. No more than 1.5 seconds (including the time consumed by the calculation of eigenvalues and eigenfunctions) are consumed for a TREE evaluation. The satisfactory algorithm efficiency is in evidence.