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

Authors: Stephane Durand; Ivan Cimrák; Peter Sergeant; Ahmed Abdallh
Journal: Mathematical Problems in Engineering
 (2010)

---

## 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.


---

## Body

## 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 Section3. 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 Figure1. 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 Figure2. 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: H=Ni/lm, where H, N, I, and lm 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 α and the width β 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 (J=σE) 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 Table1 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 inductionB (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.
Needle sensorIset(A)position (mm)0.511.522.53-340.921.241.331.411.451.49-240.891.161.131.31.351.39-140.891.191.221.311.351.35-40.931.281.351.381.41.463.50.981.451.631.681.721.7661.031.381.451.531.561.598.50.991.411.51.561.621.65111.041.351.441.491.551.5813.50.921.271.381.421.481.54160.791.181.321.351.41.51260.81.181.251.31.411.44360.811.171.171.241.421.36The 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 Iset and consisting of 100 windings. They are modelled as four regions Ω1,…,Ω4. Further, Ωstr denotes the volume occupied by the steel strips, excluding the air gap. The air gap is denoted by Ωag and starts at x=α with a width of β.(a) Computational domain covers the outer rectangle and is denoted byΩ. Copper windings are denoted by Ωi,i=1,…,4. The current in Ω1 and Ω4 flows in the opposite direction as in Ω2 and Ω3. The grey area denoted as Ωstr is one and half steel strip without the air gap. The air gap is denoted by Ωag; its left boundary is at x=α, and it has a width β  (see Figure 2). Finally, the yoke is denoted by Ωy. Short vertical lines mark the positions of local measurements of the magnetic induction. (b) Non-linear material characteristic νsteel(B) for the steel strips.
(a)(b)The strips in the sample have a thin coating in nonferromagnetic material. Therefore, the numerical model contains thin nonferromagnetic (permeabilityμ0) 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 x=α and at x=α+β, which is in contrast with the measurements.We will formulate the model in terms of the magnetic vector potentialA, satisfying ∇×A=B, where B 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 H=ν(|∇×A|,x)B, where the function ν(|∇×A|,x) 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 ν(|∇×A|,x) 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(|B|) as depicted in Figure 3(b); that is,(2.1)ν(|∇×A|,x)=νsteel(|∇×A(x)|),forx∈Ωstr.
The non-linear magnetic material characteristic for the steel is rescaled from νsteel(|∇×A|) to νsteel(|∇×A|2) for the stability reasons of Newton's method (the Jacobian of |∇×A| contains an unstable term proportional to 1/|∇×A|, 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(2.2)ν(|∇×A|2,x)=ν0χαβ(x)+νsteel(|∇×A(x)|2)(1-χαβ(x)),
where χαβ denotes the characteristic function of the air gap, that is, χαβ(x)=1 for x∈Ωag and is zero elsewhere. Henceforth, we will write the dependence of the magnetic reluctivity on α and β explicitly as ναβ(|∇×A|2,x).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 potentialA=(0,0,A), which satisfies the magnetostatic equation(2.3)∇×(ναβ(|∇×A|2,x)∇×A)=J,
equipped with the Dirichlet condition(2.4)A=0on∂Ω.
The current density J=(0,0,J) is defined by(2.5)J(x)={-100|Ω1|×Iset,x∈Ω1∪Ω4,100|Ω2|×Iset,x∈Ω2∪Ω3,0,elsewhere,
where Iset is the excitation current through the windings. These will be chosen from Table 1. The symbol |Ωi| gives the area of Ωi.


## 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 Aαβ of the non-linear system which satisfies the identity
(3.1)(ναβ(|∇×Aαβ|2,x)∇×Aαβ,∇×ϕ)=(J,ϕ),A=0,x∈∂Ω,x∈Ω,
for all ϕ=(0,0,ϕ) with ϕ∈H01(Ω). We introduce the notation (u,v)Γ for ∫Γu·v 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(s). 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 Table1 with the numerical values of |∇×Aαβ|. We assume that the magnetic induction B=∇×A is horizontal inside the upper strip (this assumption is verified by solving the direct problem). We thus define the cost functional F as(4.1)F(α,β)=∫Ωmeas|Bmeas-∇×Aαβ|2,
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 Bmeas=(Bmeas(x),0,0) supported on the upper strip with -34mm≤x≤36mm. The value of Bmeas(x) follows from akima interpolation between the measurements in Table 1. Therefore, the domain Ωmeas is the rectangular part of the upper strip with -34mm≤x≤36mm. The inverse problem is then solved to find the set of parameters {αmin,βmin} which minimizes the cost functional F(α,β). 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)∂αF(α,β)≈F(α+t,β)-F(α,β)t,∂βF(α,β)≈F(α,β+t)-F(α,β)t,
for small t. The gradient can also be computed by solving the sensitivity equation which describes how Aαβ 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 h=(h1,h2):(4.3)δhAαβ=limt→0Aαβ(α+th)-Aαβ(α)t.
Differentiation of (3.1) with respect to the parameters α and β leaves (4.4)(ναβ(|∇×Aαβ|2,x)∇×δhAαβ,∇×ϕ)+(δhναβ(|∇×Aαβ|2,x)∇×Aαβ,∇×ϕ)=0.The derivativeδhναβ(|∇×Aαβ|2,x) in the second term can be obtained from (2.2) as follows:(4.5)δhναβ(|∇×Aαβ|2,x)=(ν0-νsteel(|∇×Aαβ|2))δhχαβ+2ν′steel(|∇×Aαβ|2)(∇×Aαβ·∇×δhAαβ)(1-χαβ).
The derivative δhχαβ of the characteristic function of the air gap can be determined in distributional sense, by evaluating its action on a test function f(x). Denoting the thickness of one strip by 2d, we obtain (see also Figure 4) (4.6)∫Ωδhχαβ(x)f(x)dx=limt→01t[∫Ωχα+th1,β+th2(x)f(x)dx-∫Ωχαβ(x)f(x)dx]=limt→01t[∫0ddy∫α+βα+β+t(h1+h2)f(x,y)dx-∫0ddy∫αα+th1f(x,y)dx]=(h1+h2)∫Γα+βf(x)ds-h1∫Γαf(x)ds.
In the last step we used the mean value theorem. The sensitivity equation can thus be formulated as follows.Figure 4
Calculation ofδhχαβ in the lowest strip with thickness d=0.25mm. χαβ denotes the characteristic function of the air gap between x=α and x=α+β, and χα+th1,β+th2 is the characteristic function of the perturbed air gap between x=α+th1 and x=α+β+t(h1+h2). The shaded areas represent the support of both functions.Problem 4.1 (sensitivity equation).
For the given solutionAαβ to the direct problem, find the directional derivative δhAαβ as the solution of the linear PDE(4.7)(ναβ(|∇×Aαβ|2,x)∇×δhAαβ,∇×ϕ)+(2ν′steel(|∇×Aαβ|2)(∇×Aαβ·∇×δhAαβ)(1-χαβ)∇×Aαβ,∇×ϕ)=h1((ν0-νsteel(|∇×Aαβ|2))∇×Aαβ,∇×ϕ)Γα-(h1+h2)((ν0-νsteel(|∇×Aαβ|2))∇×Aαβ,∇×ϕ)Γα+β,
for all ϕ=(0,0,ϕ) with ϕ∈H01(Ω). For a known Aαβ,this equation is linear in δhAαβ. To obtain the gradient of Aαβ with respect to the parameters α and β, we need to solve the sensitivity equation for h=(h1,h2)=(1,0), and (0,1), respectively. The gradient of the cost functional then follows from
(4.8)δhF(α,β)=-2∫Ωmeas(Bmeas-∇×Aαβ)·∇×δhAαβ.
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α1=(α1,β1)∈ℝ2, ϵ>0, and set k=1.(2)
ComputeAαkβk from the direct problem. This step contains an inner loop from Newton's method to cope with the nonlinearity.(3)
Compute the gradientgk of the cost functional F in the point αk.
(i)
Small perturbations: computeAαk+t,βk and Aαk,βk+t from the direct problem, evaluate F(αk+t,βk) and F(αk,βk+t), then compute the gradient from (4.2).(ii)
Sensitivity equation: using the previously calculatedAαkβk, solve the sensitivity equation for h1=1,h2=0 and h1=0,h2=1, then compute the gradient from (4.8).(4)
Calculatedk=gk+ckdk-1 with ck=(gkTgk)/(gk-1Tgk-1) and d0=g0.(5)
Calculate the optimal steptk, using a line search algorithm.(6)
Setαk+1=αk-tkdk.(7)
If|F(αk)-F(αk+1)|<ϵ, then stop;otherwise set k:=k+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 aknown 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 thatBmeas is known as a continuous function on the upper strip such that interpolation in the x-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 α1=(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.
Small perturbationsSensitivityIset (A)αmin (m)βmin (m)Fminαmin (m)βmin (m)Fmin3.0-1.141971e-050.0100181.25e-10-2.171692e-110.0100007.33e-222.5-4.666539e-060.0099992.94e-117.198648e-110.0099993.77e-212.0-5.270462e-060.0100013.31e-11-3.852418e-110.0100001.35e-211.5-4.527130e-060.0099992.91e-112.344620e-110.0099991.43e-211.0-1.880308e-060.0099927.57e-111.920738e-120.0100006.30e-240.5-5.464956e-060.0100013.52e-11-8.964066e-120.0100002.34e-22The 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 givenx-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 α1=(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.
Small perturbationsSensitivityIset (A)αmin (m)βmin (m)Fminαmin (m)βmin (m)Fmin3.00.0015710.0084451.49e-060.0015740.0084571.49e-062.50.0015790.0084331.53e-060.0014950.0085371.53e-062.00.0014360.0085631.57e-060.0013560.0086631.58e-061.50.0012020.0087901.61e-060.0011620.0088381.61e-061.00.0009450.0090311.60e-060.0009340.0090321.60e-060.50.0005600.0094341.51e-060.0005680.0094171.51e-06As 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 numberr between -n and n (n=5 or n=10) and multiplying every value by 1+r/100. These perturbed values are then applied as data for the inverse problem. For a given noise level n, 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 (αmini,βmini),i=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 x=α of the air gap with the right coordinate x=α+β, instead of the position and width of the air gap.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.
(a)(b)Forn=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 Table1, obtained from a sample with air gap at α=-0.001m and β=0.010m, are now used to build the function Bmeas in the cost functional. We start the optimization algorithm with stopping criterion ϵ=1e-12 and initial guess α1=(-0.02,0.05). The initial guess starts with an air gap which covers the measurement domain almost completely.In Figure6, 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 α1 is also plotted. For the larger currents, there is a good correspondence with the measurements. For I=0.5A and I=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.Comparison between the measured induction (labelled asMeas, together with the akima interpolant) at the surface of the strips with the computed induction in the minimizer of the cost functional (finalin the legend). The dotted curve (init) represents the computed induction for the initial configuration. (a) I=3.0A. (b) I=2.5A. (c) I=2.0A. (d) I=1.5A. (e) I=1.0A. (f) I=0.5A.
(a)(b)(c)(d)(e)(f)The numerical values of the theoretically predicted position of the air gap are presented in Table4. 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 Table1, 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.
Small perturbationsSensitivityIset (A)αmin (m)βmin (m)Fminαmin (m)βmin (m)Fmin3.00.0004880.0122931.24e-050.0005220.0125891.24e-052.50.0011400.0088861.26e-050.0012170.0087911.20e-052.00.0006510.0093742.38e-050.0007390.0091792.37e-051.50.0000240.0102034.42e-050.0001530.0097584.42e-051.00.0016890.0063168.37e-050.0016230.0062288.37e-050.50.0067390.0027812.81e-040.0055530.0023192.80e-04Figure 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 forIset=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 Figure8, the results are presented for Iset=3.0A. For the other currents, similar results are obtained.Value of the cost functional in terms of (a) the number of evaluations and (b) the computational time forIset=3.0A.
(a)(b)Both methods converge to the same value of the cost function (see also Table4) 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 F. 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 Figure8(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.


## 6.1. Results of the Inverse Problem

The measurements in Table1, obtained from a sample with air gap at α=-0.001m and β=0.010m, are now used to build the function Bmeas in the cost functional. We start the optimization algorithm with stopping criterion ϵ=1e-12 and initial guess α1=(-0.02,0.05). The initial guess starts with an air gap which covers the measurement domain almost completely.In Figure6, 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 α1 is also plotted. For the larger currents, there is a good correspondence with the measurements. For I=0.5A and I=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.Comparison between the measured induction (labelled asMeas, together with the akima interpolant) at the surface of the strips with the computed induction in the minimizer of the cost functional (finalin the legend). The dotted curve (init) represents the computed induction for the initial configuration. (a) I=3.0A. (b) I=2.5A. (c) I=2.0A. (d) I=1.5A. (e) I=1.0A. (f) I=0.5A.
(a)(b)(c)(d)(e)(f)The numerical values of the theoretically predicted position of the air gap are presented in Table4. 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 Table1, 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.
Small perturbationsSensitivityIset (A)αmin (m)βmin (m)Fminαmin (m)βmin (m)Fmin3.00.0004880.0122931.24e-050.0005220.0125891.24e-052.50.0011400.0088861.26e-050.0012170.0087911.20e-052.00.0006510.0093742.38e-050.0007390.0091792.37e-051.50.0000240.0102034.42e-050.0001530.0097584.42e-051.00.0016890.0063168.37e-050.0016230.0062288.37e-050.50.0067390.0027812.81e-040.0055530.0023192.80e-04Figure 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 forIset=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 Figure8, the results are presented for Iset=3.0A. For the other currents, similar results are obtained.Value of the cost functional in terms of (a) the number of evaluations and (b) the computational time forIset=3.0A.
(a)(b)Both methods converge to the same value of the cost function (see also Table4) 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 F. 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 Figure8(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 Figure9. 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.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.
(a)(b)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≤x≤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:(7.1)(αmin,βmin,γmin,δmin)=(-0.030081,0.019847,0.010742,0.019311)m,
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(7.2)(αmin,βmin,γmin,δmin)=(-0.011268,0.021102,0.000812,0.019591)m,
again in good correspondence to the data as can be seen in Figure 10(b).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.
(a)(b)(c)(d)As a second experiment, we consider again the first configuration and reduce the number of measurement points to 21. The results are shown in Figure10(c). The numerical values of the air gaps are(7.3)(αmin,βmin,γmin,δmin)=(-0.029866,0.020295,0.009164,0.021386)m.
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 Figure10(d), and the numerical values of the minimizer are(7.4)(αmin,βmin,γmin,δmin)=(-0.034322,0.024481,0.009944,0.021653)m.
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 Figure11, 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.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.
(a)(b)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 ford parameters, d 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.



---

*Source: 574153-2010-09-26.xml*