Analysis of a Nondestructive Evaluation Technique for Defect Characterization in Magnetic Materials Using Local Magnetic Measurements

Abstract

We study a method for nondestructive testing of laminated strips of a nonlinear magnetic material. Based on local measurements of the magnetic induction at the surface, we are able to reconstruct the proper position of defects inside the material, by solving an inverse problem. This inverse problem is solved by minimizing a suitable cost function using a gradient-based optimization procedure. Calculation of the gradient is done either by the standard method of small perturbations or by solving the sensitivity equation. The latter method yields a significant reduction of the computational time. The validity of the proposed algorithm is confirmed by experimental results.

1. Introduction

Magnetic measurements are indispensable for the identification of magnetic properties of materials. Specifically, electromagnetic non-destructive evaluation (NDE) has gained a lot of interest during the last decade, since it is becoming an important tool in different research topics. The determination of residual stresses in electrical steels, by the aid of the eddy current technique, has been studied [1]. The sizes of surface cracks have been estimated based on a combined numerical “dipole model of a crack”—experimental “Hall element measurements” technique [2]. The magnetic flux leakage technique has been utilized to characterize the defect by solving an inverse problem based on a space mapping methodology [3].

In particular, eddy current testing has been widely used for non-destructive testing such as defects characterization in materials [4]. However, the sensitivity of the characterization process is highly dependent on the probe choice and the operation frequency [5].

In this paper, we aim at developing a new technique for proper defect characterization in magnetic materials. This aim is achieved by solving an inverse problem based on non-destructive local magnetic measurements.

We consider a laminated stack in ferromagnetic material with a defect: an “air gap” in the laminated stack. We construct a finite element model of the measurement setup, which solves the direct problem. For a known position and size of the “air gap” defect, this model returns the value of the magnetic induction on the boundary of the material. In order to solve the inverse problem, we define a cost function that compares these values with the measurements. This cost function depends on the position and size of the air gap and will obtain a minimum if the air gap in the model coincides with the actual defect. We use a gradient-based optimization procedure to compute this minimum, where gradients are evaluated either with the conventional method (using small perturbations in the parameters) or by solving the sensitivity equation. The conventional method requires solving the direct problem for every perturbed parameter, which is very time consuming for non-linear problems. The sensitivity equation has the advantage that it is always linear, but it depends on the solution of the direct problem and can therefore be less accurate, resulting in more iterations. We will compare both techniques for the evaluation of non-linear materials.

In the following section we describe the measurement setup and the mathematical model. The direct problem is defined in Section 3. In Section 4 we derive the sensitivity equation and formulate an algorithm to solve the inverse problem numerically. This algorithm is validated for several test cases in Section 5 and applied to an experimental setting in Section 6. In Section 7 we study the validity of the algorithm for samples with multiple layers. Finally, the conclusions are drawn in Section 8.

2. Measurement Setup and Mathematical Model

The sample under test consists of three electrical steel strips, where the middle one has an air gap, as shown in Figure 1. The magnetic material properties (B-H characteristics) of the strips are known. They can be easily measured using the Epstein frame standard [6].

Figure 1 

Three layers of steel lamination with a defect in the middle one.

The considered NDE experiment is shown in Figure 2. The position 𝛼 and width 𝛽 of the defect are to be determined from the inverse problem. The steel strips are in a single sheet tester (SST) that imposes a horizontal magnetic field in the strips. As the excitation current is known and the reluctance of the flux return path “yoke” is negligible, the magnetic field in the strip can be calculated, according to Ampere's law: 𝐻=𝑁𝑖/𝑙_𝑚, where H, N, I, and 𝑙_𝑚 are the magnetic field in the material, excitation number of turns, excitation current, and magnetic path length in the material, respectively. The magnetic induction is measured locally at many positions along the upper surface of the strips, using the needle probe method (NPM). The NPM is preferable for local magnetic induction measurements as it is a “non-destructive method”.

Figure 2 

Geometry ($xy$-plane) of the non-destructive evaluation setup, consisting of a single-sheet tester and a stack of steel laminations. The $x$ and $y$ axes are symmetry axes of the single-sheet tester. The complete setup is not symmetric any more because of the defect, but the SST itself remains symmetric. The goal is to determine the position $\alpha$ and the width $\beta$ of the defect zone by doing local $B$-measurements on the surface of the stack only.

The needle probe method is based on the fact that when a time-dependent flux Φ is enforced in a steel sheet, eddy currents (𝐽=𝜎𝐸) are introduced, resulting in potential differences measured by the needles [7]. Although the NPM is an efficient technique for measuring the local magnetic induction, it encounters numerous errors, such as errors due to vertical field component and its sensitivity to noise interference [8]. Therefore, we used a modified connection for the needles in order to compensate the air flux; for more details see [9].

The quasi static magnetic measurements are performed at 1 Hz for a sinusoidal current excitation, in order to have a negligible presence of eddy current effects in the magnetic strip.

The results of the measurements are given in Table 1 for several values of the excitation current. These data will be used to determine the position of the air gap.

Table 1 

Measurements of the magnetic induction 𝐵 (in Tesla) at 12 points on the surface of the strip for several values of the imposed current and for a defect with 𝛼=−0.001 m and 𝛽=0.01m.

The problem is modelled in 2D; it is symmetric along the horizontal axis of the middle strip. Our computational domain thus consists of the upper half of the setting and will be denoted by Ω (see Figure 3). The magnetic field is induced by two copper windings carrying a current 𝐼_set and consisting of 100 windings. They are modelled as four regions Ω₁,…,Ω₄. Further, Ω_str denotes the volume occupied by the steel strips, excluding the air gap. The air gap is denoted by Ω_ag and starts at 𝑥=𝛼 with a width of 𝛽.

(a)

(b)

Figure 3 

(a) Computational domain covers the outer rectangle and is denoted by $\mathrm{\Omega }$. Copper windings are denoted by ${\mathrm{\Omega }}_i,\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\hspace{0.17em}\hspace{0.17em}}i=1,\dots ,4$. The current in ${\mathrm{\Omega }}_1$ and ${\mathrm{\Omega }}_4$ flows in the opposite direction as in ${\mathrm{\Omega }}_2$ and ${\mathrm{\Omega }}_3$. The grey area denoted as ${\mathrm{\Omega }}_{\text{str}}$ is one and half steel strip without the air gap. The air gap is denoted by ${\mathrm{\Omega }}_{\text{ag}}$; its left boundary is at $x=\alpha$, and it has a width $\beta$ (see Figure 2). Finally, the yoke is denoted by ${\mathrm{\Omega }}_y$. Short vertical lines mark the positions of local measurements of the magnetic induction. (b) Non-linear material characteristic ${\nu }_{\text{steel}}(B)$ for the steel strips.

The strips in the sample have a thin coating in nonferromagnetic material. Therefore, the numerical model contains thin nonferromagnetic (permeability 𝜇₀) and nonconducting layers between adjacent strips, with a thickness equal to one tenth of the thickness of one strip. If these coating layers are not modelled, the simulated magnetic induction at the surface changes abruptly at 𝑥=𝛼 and at 𝑥=𝛼+𝛽, which is in contrast with the measurements.

We will formulate the model in terms of the magnetic vector potential 𝐀, satisfying ∇×𝐀=𝐁, where 𝐁 denotes the magnetic induction. (The operator ∇× denotes the standard curl operator.) The interaction of the materials with the magnetic field is modelled by the constitutive relation 𝐇=𝜈(|∇×𝐀|,𝐱)𝐁, where the function 𝜈(|∇×𝐀|,𝐱) denotes the magnetic reluctivity of the domain and is a highly non-linear function inside the steel strips. The other materials in the domain are considered linear, such that 𝜈(|∇×𝐀|,𝐱) takes a constant value there. The non-linear behaviour of the steel strips follows from the measured B-H characteristic for a steel sample without defects and is a nonmonotone function 𝜈_steel(|𝐁|) as depicted in Figure 3(b); that is,𝜈(|∇×𝐀|,𝐱)=𝜈_steel(|∇×𝐀(𝐱)|),for𝑥∈Ω_str. (2.1) The non-linear magnetic material characteristic for the steel is rescaled from 𝜈_steel(|∇×𝐀|) to 𝜈_steel(|∇×𝐀|²) for the stability reasons of Newton's method (the Jacobian of |∇×𝐀| contains an unstable term proportional to 1/|∇×𝐀|, which is avoided by taking the square).

The air gap inside the lowest strip is modelled by two parameters 𝛼 and 𝛽. We model the reluctivity of the steel strips as𝜈(|∇×𝐀|²,𝐱)=𝜈₀𝜒_𝛼𝛽(𝐱)+𝜈_steel(|∇×𝐀(𝐱)|²)(1−𝜒_𝛼𝛽(𝐱)), (2.2) where 𝜒_𝛼𝛽 denotes the characteristic function of the air gap, that is, 𝜒_𝛼𝛽(𝐱)=1 for 𝐱∈Ω_ag and is zero elsewhere. Henceforth, we will write the dependence of the magnetic reluctivity on 𝛼 and 𝛽 explicitly as 𝜈_𝛼𝛽(|∇×𝐀|²,𝐱).

The use of the characteristic function 𝜒_𝛼𝛽 has the advantage that we do not need to remesh the domain for every new value of 𝛼 and 𝛽. An alternative method would be to define the position of the air gap explicitly into the geometry of the domain and create a new mesh for every new position of the air gap. Since the steel strips are very thin compared to the other dimensions of the domain, the mesh is very fine there and creating the mesh is time consuming. Moreover, since the inverse problem is solved iteratively, the direct problem has to be solved for several values of 𝛼 and 𝛽. The computation time can therefore be reduced drastically by introducing the characteristic function 𝜒_𝛼𝛽.

The magnetic field is described in terms of the vector potential 𝐀=(0,0,𝐴), which satisfies the magnetostatic equation∇×(𝜈_𝛼𝛽(|∇×𝐀|²,𝐱)∇×𝐀)=𝐉, (2.3) equipped with the Dirichlet condition𝐴=0on𝜕Ω. (2.4) The current density 𝐉=(0,0,𝐽) is defined by−𝐽(𝐱)=100|Ω₁|×𝐼_set,𝐱∈Ω₁∪Ω₄,100|Ω₂|×𝐼_set,𝑥∈Ω₂∪Ω₃,0,elsewhere, (2.5) where 𝐼_set is the excitation current through the windings. These will be chosen from Table 1. The symbol |Ω_𝑖| gives the area of Ω_𝑖.

3. Direct Problem

Since we will use a finite element model, we will reformulate the problem in variational form. The forward problem can then be formulated as follows.

Problem 3.1 (direct problem). For given parameters 𝛼 and 𝛽, find the solution 𝐀_𝛼𝛽 of the non-linear system which satisfies the identity (𝜈_𝛼𝛽(|∇×𝐀_𝛼𝛽|²,𝐱)∇×𝐀_𝛼𝛽,∇×𝝓)=(𝐉,𝝓),𝐴=0,𝑥∈𝜕Ω,𝐱∈Ω, (3.1) for all 𝝓=(0,0,𝜙) with 𝜙∈𝐻¹₀(Ω). We introduce the notation (𝐮,𝐯)_Γ for ∫_Γ𝐮⋅𝐯 and omit the index if Γ is the whole domain. Well-posedness of this problem is guaranteed from the theory of quasi-linear equations [10]. This is of course after the assumption that 𝜈 is greater than a positive constant and has controlled growth. This problem is solved in GetDP [11] using a 2D finite element model. The computational domain (Figure 3(a)) is triangulated into a regular mesh consisting of 24792 triangles with Gmsh [12]. Since the dimensions of the strips are very small compared to the measurement device, the samples are meshed very fine compared to the yoke and the surrounding air. This problem is linearized using Newton's method and solved iteratively, where akima interpolation [13] is applied to interpolate the function 𝜈_steel(𝑠). The characteristic function 𝜒_𝛼𝛽 in (2.2) is approximated by a smooth function.

4. Inverse Problem

To solve the inverse problem, we need to compare the measured values of the magnetic induction as given in Table 1 with the numerical values of |∇×𝐀_𝛼𝛽|. We assume that the magnetic induction 𝐁=∇×𝐀 is horizontal inside the upper strip (this assumption is verified by solving the direct problem). We thus define the cost functional 𝐹 as𝐹(𝛼,𝛽)=_Ωmeas|𝐁_meas−∇×𝐀_𝛼𝛽|², (4.1) where Ω_meas is the domain where the magnetic induction is measured. Assuming that the magnetic induction is uniform over the depth of the upper strip, we can interpret the measured values as a field 𝐁_meas=(𝐵_meas(𝑥),0,0) supported on the upper strip with −34mm≤𝑥≤36mm. The value of 𝐵_meas(𝑥) follows from akima interpolation between the measurements in Table 1. Therefore, the domain Ω_meas is the rectangular part of the upper strip with −34mm≤𝑥≤36mm. The inverse problem is then solved to find the set of parameters {𝛼_min,𝛽_min} which minimizes the cost functional 𝐹(𝛼,𝛽). We will solve this optimization problem iteratively using a conjugate gradients algorithm (CG), which requires the evaluation of the gradient. A straightforward way to compute the gradient is using small perturbations in the parameters𝜕_𝛼𝐹(𝛼,𝛽)≈𝐹(𝛼+𝑡,𝛽)−𝐹(𝛼,𝛽)𝑡,𝜕_𝛽𝐹(𝛼,𝛽)≈𝐹(𝛼,𝛽+𝑡)−𝐹(𝛼,𝛽)𝑡, (4.2) for small 𝑡. The gradient can also be computed by solving the sensitivity equation which describes how 𝐀_𝛼𝛽 changes if the parameters 𝛼 and 𝛽 are modified. It is obtained by differentiation of (3.1) with respect to the parameters. We introduce the following notation for the directional derivative in point 𝜶=(𝛼,𝛽) and direction 𝐡=(ℎ₁,ℎ₂)∶𝛿_𝐡𝐀_𝛼𝛽=lim_𝑡→0𝐀_𝛼𝛽(𝜶+𝑡𝐡)−𝐀_𝛼𝛽(𝜶)𝑡. (4.3) Differentiation of (3.1) with respect to the parameters 𝛼 and 𝛽 leaves (𝜈_𝛼𝛽(|∇×𝐀_𝛼𝛽|²,𝐱)∇×𝛿_𝐡𝐀_𝛼𝛽,∇×𝝓)+(𝛿_𝐡𝜈_𝛼𝛽(|∇×𝐀_𝛼𝛽|²,𝐱)∇×𝐀_𝛼𝛽,∇×𝝓)=0. (4.4)

The derivative 𝛿_𝐡𝜈_𝛼𝛽(|∇×𝐀_𝛼𝛽|²,𝐱) in the second term can be obtained from (2.2) as follows:𝛿_𝐡𝜈_𝛼𝛽(|∇×𝐀_𝛼𝛽|²,𝐱)=(𝜈₀−𝜈_steel(|∇×𝐀_𝛼𝛽|²))𝛿_𝐡𝜒_𝛼𝛽+2𝜈^′_steel(|∇×𝐀_𝛼𝛽|²)(∇×𝐀_𝛼𝛽⋅∇×𝛿_𝐡𝐀_𝛼𝛽)(1−𝜒_𝛼𝛽). (4.5) The derivative 𝛿_𝐡𝜒_𝛼𝛽 of the characteristic function of the air gap can be determined in distributional sense, by evaluating its action on a test function 𝑓(𝐱). Denoting the thickness of one strip by 2𝑑, we obtain (see also Figure 4) _Ω𝛿_𝐡𝜒_𝛼𝛽(𝐱)𝑓(𝐱)d𝐱=lim_𝑡→01𝑡[_Ω𝜒_{𝛼+𝑡ℎ1,𝛽+𝑡ℎ2}(𝐱)𝑓(𝐱)d𝐱−_Ω𝜒_𝛼𝛽(𝐱)𝑓(𝐱)d𝐱]=lim_𝑡→01𝑡[^𝑑₀d𝑦^{𝛼+𝛽+𝑡(ℎ1+ℎ2)}_𝛼+𝛽𝑓(𝑥,𝑦)d𝑥−^𝑑₀d𝑦^{𝛼+𝑡ℎ1}_𝛼𝑓(𝑥,𝑦)d𝑥]=(ℎ₁+ℎ₂)_Γ𝛼+𝛽𝑓(𝑥)d𝑠−ℎ_1Γ𝛼𝑓(𝑥)d𝑠. (4.6) In the last step we used the mean value theorem. The sensitivity equation can thus be formulated as follows.

Figure 4 

Calculation of ${\delta }_{\mathit{h}}{\chi }_{\alpha \beta }$ in the lowest strip with thickness $d=0.25\hspace{0.17em}\text{mm}\text{.}$${\chi }_{\alpha \beta }$ denotes the characteristic function of the air gap between $x=\alpha$ and $x=\alpha +\beta ,$ and ${\chi }_{\alpha +t{h}_1,\beta +t{h}_2}$ is the characteristic function of the perturbed air gap between $x=\alpha +t{h}_1$ and $x=\alpha +\beta +t(h_1+h_2)$. The shaded areas represent the support of both functions.

Problem 4.1 (sensitivity equation). For the given solution 𝐀_𝛼𝛽 to the direct problem, find the directional derivative 𝛿_𝐡𝐀_𝛼𝛽 as the solution of the linear PDE(𝜈_𝛼𝛽(|∇×𝐀_𝛼𝛽|²,𝐱)∇×𝛿_𝐡𝐀_𝛼𝛽,∇×𝝓)+(2𝜈^′_steel(|∇×𝐀_𝛼𝛽|²)(∇×𝐀_𝛼𝛽⋅∇×𝛿_𝐡𝐀_𝛼𝛽)(1−𝜒_𝛼𝛽)∇×𝐀_𝛼𝛽,∇×𝝓)=ℎ1((𝜈₀−𝜈_steel(|∇×𝐀_𝛼𝛽|²))∇×𝐀_𝛼𝛽,∇×𝝓)_Γ𝛼−(ℎ1+ℎ2)((𝜈₀−𝜈_steel(|∇×𝐀_𝛼𝛽|²))∇×𝐀_𝛼𝛽,∇×𝝓)_Γ𝛼+𝛽, (4.7) for all 𝝓=(0,0,𝜙) with 𝜙∈𝐻¹₀(Ω). For a known 𝐀_𝛼𝛽,this equation is linear in 𝛿_𝐡𝐀_𝛼𝛽. To obtain the gradient of 𝐀_𝛼𝛽 with respect to the parameters 𝛼 and 𝛽, we need to solve the sensitivity equation for 𝐡=(ℎ₁,ℎ₂)=(1,0), and (0,1), respectively. The gradient of the cost functional then follows from 𝛿_𝐡𝐹(𝛼,𝛽)=−2_Ωmeas(𝐵_meas−∇×𝐀_𝛼𝛽)⋅∇×𝛿_𝐡𝐀_𝛼𝛽. (4.8)
The solution of the inverse problem is now reduced to the minimization of the cost functional. We use an iterative CG algorithm where gradients are calculated either by applying small perturbations in the parameters (see (4.2)) or by solving the sensitivity equation.

Algorithm 4.2 (non-linear conjugate gradient). (1)Choose 𝜶₁=(𝛼₁,𝛽₁)∈ℝ², 𝜖>0, and set 𝑘=1.(2)Compute 𝐀_{𝛼𝑘𝛽𝑘} from the direct problem. This step contains an inner loop from Newton's method to cope with the nonlinearity.(3)Compute the gradient 𝐠_𝑘 of the cost functional 𝐹 in the point 𝜶_𝑘.(i)Small perturbations: compute 𝐀_{𝛼𝑘+𝑡,𝛽𝑘} and 𝐀_{𝛼𝑘,𝛽𝑘+𝑡} from the direct problem, evaluate 𝐹(𝛼_𝑘+𝑡,𝛽_𝑘) and 𝐹(𝛼_𝑘,𝛽_𝑘+𝑡), then compute the gradient from (4.2). (ii)Sensitivity equation: using the previously calculated 𝐀_{𝛼𝑘𝛽𝑘}, solve the sensitivity equation for ℎ₁=1,ℎ₂=0 and ℎ₁=0,ℎ₂=1, then compute the gradient from (4.8).(4)Calculate 𝐝_𝑘=𝐠_𝑘+𝑐_𝑘𝐝_𝑘−1 with 𝑐_𝑘=(𝐠^𝑇_𝑘𝐠_𝑘)/(𝐠^𝑇_𝑘−1𝐠_𝑘−1) and 𝐝₀=𝐠₀.(5)Calculate the optimal step 𝑡_𝑘, using a line search algorithm.(6)Set 𝜶_𝑘+1=𝜶_𝑘−𝑡_𝑘𝐝_𝑘.(7)If |𝐹(𝜶_𝑘)−𝐹(𝜶_𝑘+1)|<𝜖, then stop;otherwise set 𝑘∶=𝑘+1 and return to step (2).

5. Validation of the Inverse Problem

We compare the conventional inverse algorithm and the inverse algorithm with sensitivity analysis. For validation, we use computed field values at the edge of the lamination for a known air gap as measurement data. The solution algorithms are validated by comparing the numerical result of the inverse problem with the known position of the air gap.

As a first experiment, we suppose that 𝐁_meas is known as a continuous function on the upper strip such that interpolation in the 𝑥-direction is not necessary. We set 𝛼=0m, 𝛽=0.01m and compute the solution of the direct problem. This solution is used as measurement data to solve the inverse problem by the stated algorithm. We repeat this procedure for the 6 different currents given in Table 1. The initial guess is set to 𝜶₁=(0.02,0.01) and stopping criterion 𝜖= 1e-20. The results are stated in Table 2. Both methods (small perturbation and sensitivity) converge sufficiently close to the actual air gap. The sensitivity method converges to the exact position of the air gap, while the perturbation method stops at a value very close to the real air gap.

Table 2 

Results of the inverse problem for simulated measurement data obtained from a sample with air gap 𝛼=(0,0.01)m. The method using sensitivity equation converges to the exact solution of the inverse problem, corresponding to the absolute minimum of the cost functional. For the perturbation method, the algorithm converges to a value very close to the exact solution.

The previous experiment is a good validation of the algorithm but assumes an infinite number of measurements along the surface of the strip and is therefore far from reality. Therefore, we consider now the case where only discrete values of the magnetic induction on the given 𝑥-positions in Table 1 are known. These values are again calculated from the direct problem and thus free of noise. Interpolation between these values is based on akima interpolation. The results of the optimization procedure with discrete measurements are shown in Table 3. The stopping criterion is set to 𝜖=1e-12 and the initial guess is again 𝜶₁=(0.02,0.01). Both methods converge to a local optimal value close to the exact minimizer but are not able to find the global optimum. Even though the measurements were exact, we could not obtain the exact position of the air gap due to the discretization of the measurements. However, taking into account the coarseness of the measurements, our method yields values which are close to the exact air gap and which are consistent for all currents. One cannot expect better results without refining the measurement grid.

Table 3 

Results of the inverse problem for discrete simulated measurement data obtained from a sample with air gap 𝛼=(0,0.01)m. Both methods converge to a local optimum close to the exact solution.

As a last validation of our method, we study the sensitivity of the algorithm with respect to noise in the data. We consider again 12 measurements on the surface of the sample, computed by solving the direct problem. White noise is added to these values by generating a pseudorandom number 𝑟 between −𝑛 and 𝑛 (𝑛=5 or 𝑛=10) and multiplying every value by 1+𝑟/100. These perturbed values are then applied as data for the inverse problem. For a given noise level 𝑛, this procedure is repeated 100 times, that is, the inverse problem is solved 100 times, for every set of noisy data. The difference between the obtained minimizers (𝛼^𝑖_min,𝛽^𝑖_min),𝑖=1,…,100 and the exact position of the air gap (𝛼_ex,𝛽_ex) is plotted in Figure 5. Note that we compare 𝛼_min+𝛽_min to 𝛼_ex+𝛽_ex on the vertical axis, since a negative value of 𝛼_min−𝛼_ex will drive the algorithm to a more positive value of 𝛽_min−𝛽_ex in order to comply with the data. Therefore, we compare the left coordinate 𝑥=𝛼 of the air gap with the right coordinate 𝑥=𝛼+𝛽, instead of the position and width of the air gap.

(a)

(b)

Figure 5 

Result of the inverse problem for 100 different measurements with noise. (a) Noise level 5%, the numerical method was stable for all simulations. (b) Noise level 10%, the method converged only for 69 out of 100 simulations.

For 𝑛=5, the algorithm is stable for every set of perturbed data and converges close to the actual air gap. From Figure 5(a) we conclude that the position of the air gap is determined with an accuracy of 5 mm for almost all simulations. Increasing the noise level to 10% reduces the stability of the inverse problem. Only 69 out of the 100 simulations resulted into a reasonable position of the air gap, and the mean distance of the corresponding points in Figure 5(b) to the origin is larger. The numerical algorithm is thus stable for a noise level up to 5%.

6. Identification of the Air Gap with Real Measurements

We will now apply our algorithm to determine an air gap in a laminated stack by local magnetic measurements on the surface. The evaluation is repeated for several currents on the same sample. In the first part we will focus on the results of the inverse problem and compare the computed minimizers of the cost functional with the actual air gap in the real sample. The comparison between the performance of both methods for gradient computation is done in the second part.

6.1. Results of the Inverse Problem

The measurements in Table 1, obtained from a sample with air gap at 𝛼=−0.001m and 𝛽=0.010m, are now used to build the function 𝐁_meas in the cost functional. We start the optimization algorithm with stopping criterion 𝜖=1e-12 and initial guess 𝜶₁=(−0.02,0.05). The initial guess starts with an air gap which covers the measurement domain almost completely.

In Figure 6, the value of the computed magnetic induction (as result of the inverse problem) at the surface of the upper strip is compared to the measurements for every current. The induction at the initial configuration 𝜶₁ is also plotted. For the larger currents, there is a good correspondence with the measurements. For 𝐼=0.5A and 𝐼=1.0A, the measured values are smaller than the predicted values obtained by the model. At low current values, the output signal of the NPM is very small and very sensitive to the noise, which explains the difference between the model and the measurements.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Comparison between the measured induction (labelled as Meas, together with the akima interpolant) at the surface of the strips with the computed induction in the minimizer of the cost functional ( final in the legend). The dotted curve ( init) represents the computed induction for the initial configuration. (a) $I=3.0\hspace{0.17em}\hspace{0.17em}\text{A}$. (b) $I=2.5\hspace{0.17em}\text{A}$. (c) $I=2.0\hspace{0.17em}\text{A}$. (d) $I=1.5\hspace{0.17em}\text{A}$. (e) $I=1.0\hspace{0.17em}\text{A}$. (f) $I=0.5\hspace{0.17em}\text{A}$.

The numerical values of the theoretically predicted position of the air gap are presented in Table 4. Both methods for gradient evaluation yield a good prediction of the actual air gap for large currents. The position of the predicted air gap is visualized in Figure 7 for all currents.

Table 4 

Results of the inverse problem for the measured data from Table 1, obtained from a sample with air gap at 𝛼=−0.001 m and 𝛽=0.010 m. The results are similar for both methods. We obtain a very good prediction of the position of the air gap for the larger currents.

Figure 7 

Schematic representation of the numerically predicted position of the air gap for all currents. The position of the actual air gap is given by the dashed lines.

For most currents we get a good correspondence between the predicted position of the air gap and its actual position. Only for 𝐼_set=0.5A, the poor agreement between the measurements and the model results into a numerical air gap which is much smaller than the true gap. We conclude that for sufficiently large currents, the position of the air gap can be determined from local boundary measurements on a rather coarse grid of measurements.

6.2. Comparison of Both Gradient Methods

In order to compare the performance of the sensitivity equation for gradient evaluation with the method of small perturbations, we determine the value of the cost function as a function of both the number of function evaluations and the computation time. In Figure 8, the results are presented for 𝐼_set=3.0A. For the other currents, similar results are obtained.

(a)

(b)

Figure 8 

Value of the cost functional in terms of (a) the number of evaluations and (b) the computational time for $I_{\text{set}}=3.0\hspace{0.17em}\text{A}$.

Both methods converge to the same value of the cost function (see also Table 4) and give similar values for 𝛼 and 𝛽. Note that the initial value of the cost functional is already small due to the small integration domain Ω_meas in (4.1). For both methods, we get a significant reduction of the cost functional in the first iterations. From Figure 8(a), we see that the direct method of gradient computation needs less evaluations to reduce the value of the cost functional. The method based on the sensitivity equation is also less stable, as can be seen from the sudden jump in the value of 𝐹. The direct method using small perturbations evaluates the gradient more accurately. Indeed, the coefficients of the sensitivity equation depend on the solution of the direct problem. Small errors on this solution decrease the stability of the sensitivity equation and can lead to large errors on the gradient. This is not the case for the perturbation method where the gradient is calculated by solving the direct problem twice.

In Figure 8(b), the value of the cost functional is given as a function of the time at which it is evaluated. As expected, the method based on the sensitivity equation is much faster than the direct method based on small perturbations. The reason is clear: in every iteration we need to solve three non-linear problems for the direct method, while the sensitivity method needs one non-linear (direct problem) and two linear (sensitivity equation) problems to evaluate the gradient. Even though the sensitivity method needs slightly more evaluations, it is advantageous in yielding a reduction of the computation time by a factor two.

7. Multiple Layers

In this last section, we wish to expand the previous model for samples consisting of multiple layers. We wish to study the possibility to detect air gaps in more complex geometries from a theoretical point of view using artificial measurements. This is only possible if the surface measurements also contain information on the depth of the air gap. We assume that the sample is now built from 5 strips. We denote the highest strip by “strip 1” and assign numbers 2–5 to the strips at lower positions. Again we exploit the symmetry of the problem to model only 2.5 strips. In this stack, we consider two defects: an air gap in strip 2 (and 4, due to symmetry) and an air gap in strip 3. Both gaps are modelled by two parameters, that is, 𝛼 and 𝛽 for the first gap (strip 2) and 𝛾 and 𝛿 for the second gap (strip 3). The inverse problem then consists of the reconstruction of these four parameters.

The solvability of the inverse problem is determined by the uniqueness of the parameters 𝛼,𝛽,𝛾, and 𝛿 for given magnetic induction on the surface. From the previous analysis, we know that the air gap causes a raise in the magnetic induction on the surface of the sample. To be able to distinguish between the two air gaps, they have to produce a detectable difference in the B-values on the surface. Therefore, we study first the direct problem; that is, for the given parameters of the two gaps we compute the value of the magnetic induction on the surface.

The result of two configurations is presented in Figure 9. In Figure 9(a), there is no overlapping between the two air gaps. One can clearly distinguish the effect of both gaps on the local B-values on the surface. The gap in the lower strip 3 causes a smaller rise of the magnetic induction than the gap in strip 2, because the flux in strip 3 will bypass the gap mainly via strip 2, and less via the highest strip 1. Indeed, the magnetic bypass via strip 1 is less favourable, because the flux has to cross the nonmagnetic coating of the strips four times. For the overlapping gaps in Figure 9(b), the position of the gaps is not explicitly presented by the boundary values of the magnetic induction. The position of the upper gap is implicitly present only in the higher values of the magnetic induction. We will now investigate the possibility of our algorithm to reconstruct these configurations.

(a)

(b)

Figure 9 

Numerical evaluation of the direct problem for a sample consisting of 5 strips. The computed magnetic induction on the surface of the stack is plotted as a function of the position. (a) Case of nonoverlapping gaps. There is a significant difference in the values of the magnetic induction on the surface at the position of the two gaps. (b) Partial overlap between the two gaps. The position of the two gaps is now more implicitly contained in the magnetic induction at the surface.

Since the geometry is more complex than in the previous section, more measurement points are needed to reconstruct the geometry of the sample. Therefore, we consider a number of uniformly distributed fictitious measurements in the interval −0.05m≤𝑥≤0.05m. We will study the stability of the method if the number of measurements is reduced. The measurement data are computed by solving the direct problem for a known setting and white noise is added by the same procedure as in Section 5.

As a first experiment, we take 41 measurements and a noise level of 5%. We study two different configurations. For the first configuration, the position of the air gaps is (𝛼_ex,𝛽_ex,𝛾_ex,𝛿_ex)=(−0.03,0.02,0.01,0.02)m, that is, the two gaps do not overlap. After solving the inverse problem, we obtain the following minimizer for the cost function:(𝛼_min,𝛽_min,𝛾_min,𝛿_min)=(−0.030081,0.019847,0.010742,0.019311)m, (7.1) being consistent with the exact configuration. The corresponding induction at the surface is given in Figure 10(a), together with the artificial measurements and the akime interpolant. As a second configuration, we consider overlapping gaps (𝛼_ex,𝛽_ex,𝛾_ex,𝛿_ex)=(−0.01,0.02,0.00,0.02)m. The inverse problem yields(𝛼_min,𝛽_min,𝛾_min,𝛿_min)=(−0.011268,0.021102,0.000812,0.019591)m, (7.2) again in good correspondence to the data as can be seen in Figure 10(b).

(a)

(b)

(c)

(d)

Figure 10 

Results of the inverse problem for a sample with multiple strips. The artificial measurement points are plotted together with the akima interpolant and the surface magnetic induction at the optimal configuration. (a) Nonoverlapping gaps, 41 measurements, 5% noise. (b) Overlapping gaps, 41 measurements, 5% noise. (c) Non-overlapping gaps, 21 measurements, 5% noise. (d) Non-overlapping gaps, 21 measurements, 10% noise.

As a second experiment, we consider again the first configuration and reduce the number of measurement points to 21. The results are shown in Figure 10(c). The numerical values of the air gaps are(𝛼_min,𝛽_min,𝛾_min,𝛿_min)=(−0.029866,0.020295,0.009164,0.021386)m. (7.3) Again we obtain good results, even for this reduced number of measurements. However, if we take less than 21 measurements, the algorithm will not converge to the exact configuration.

Finally, we increased the noise level to 10% on the 21 measurements for the previous configuration. The results are presented in Figure 10(d), and the numerical values of the minimizer are(𝛼_min,𝛽_min,𝛾_min,𝛿_min)=(−0.034322,0.024481,0.009944,0.021653)m. (7.4) The reconstruction of the lower gap is less accurate due to the high level of noise, but we still get an acceptable minimizer from the numerical algorithm.

These experiments indicate that the algorithm can be used for NDE of samples with more strips if the noise on the measurements can be restricted to 5%. However, as the number of strips increases, the effect on the surface measurements of air gaps in the lowest strip decreases and more accurate measurements are needed to detect these gaps.

The previous experiments were repeated for both methods of gradient evaluation. Since the cost function now depends on four parameters, gradient computation by the perturbation method requires solving the direct problem 5 times on every iteration. For the method with sensitivity equation, we need to solve the direct problem once and the sensitivity equation 4 times. In Figure 11, the value of the cost function is given as function of the number of evaluations and as function of the computation time. As for the case with two parameters, the sensitivity method requires slightly more evaluations before the stopping criterion is satisfied. However, the computing time is now significantly reduced by solving the linear sensitivity equation instead of the non-linear direct problem for gradient evaluation.

(a)

(b)

Figure 11 

Comparison of the performance of both methods for gradient computation. (a) Number of function evaluations. For the method based on the sensitivity equation, more evaluations are needed to satisfy the stopping criterion. (b) Value of the cost function as a function of the computing time.

As the number of strips increases, the effect of an air gap inside the sample on the boundary measurements decreases. More accurate measurements are needed to detect the raise in magnetic induction. Moreover, the nonferromagnetic layers between the strips will cause a smoothing effect on the magnetic induction on the boundary, such that the magnetic induction at the surface will change only gradually at the position of the air gap. More strips result also in more unknown parameters. The method based on sensitivity equation can then get very costly, since for 𝑑 parameters, 𝑑 sensitivity equations have to be solved on every iteration. This can be avoided by using the adjoint variable method, where the gradient is computed as the scalar product with a suitable adjoint variable [14].

8. Conclusions

We have described a numerical scheme for the non-destructive evaluation of laminated stacks of a non-linear magnetic material. Based on local measurements of the magnetic induction at the surface of the sample, we were able to reconstruct the position of a crack in the inner layer. The inverse problem was solved by minimizing a suitable cost function using a gradient-based optimization procedure. Calculation of the gradient was done with the standard method of small perturbations or by solving the sensitivity equation. We have shown that a significant reduction of computation time is obtained if the latter method is used.

We have applied this algorithm for several test cases. The algorithm was then used to reconstruct an air gap in a steel sample from real measurements. We obtained a good prediction of the position of the actual air gap.

Finally, we extended the algorithm for more complex samples consisting of more steel layers and indicated the possibility to reconstruct air gaps on different depths if the measurements are sufficiently accurate.

Acknowledgments

Stephane Durand is supported by the BOF-Grant 01D28807 of Ghent University. Ivan Cimrák and Peter Sergeant are supported by the Fund for Scientific Research-Flanders FWO (Belgium). The paper was also supported by the IAP P6/21 project of the Belgian government.