Inverse Scheme to Locally Determine Nonlinear Magnetic Material Properties: Numerical Case Study

: We are interested in the determination of the local nonlinear magnetic material behaviour in electrical steel sheets due to cutting and punching effects. For this purpose, the inverse problem has to be solved, where the objective function, which penalises the difference between the measured and the simulated magnetic flux density, has to be minimised under a constraint defined according to the corresponding partial differential equation model. We use the adjoint method to efficiently obtain the gradients of the objective function with respect to the material parameters. The optimisation algorithm is low-memory Broyden–Fletcher–Goldfarb–Shanno (BFGS), the forward and adjoint formulations are solved using the finite element (FE) method and the ill-posedness is handled via Tikhonov regularisation, in combination with the discrepancy principle. Realistic numerical case studies show promising results.


Introduction
For manufacturers of magnetic materials, measurement techniques are used to determine macroscopic properties in terms of characteristics (losses, magnetisation, coercive and remanent fields, etc.) and curves (virgin curve, hysteresis loop, commutation curve, etc.).In addition, these measurement techniques are used for the systematic study and development of materials.Many of the common test methods are specified in IEC 60404 ff.Standard techniques include the hysteresis graph, the coercive meter, the Epstein frame, the single sheet tester and the double-C yoke (see, for example, [1]).
In order to better understand and model the various loss mechanisms that occur in electrical machines during operation, the actual degradation of electrical steel sheets due to cutting, punching, welding etc. is of paramount importance.A finite-element (FE) model has been presented in [2] in order to consider the local degradation effect.The approach introduces a certain number of boundary layers along the punched edges, where the relative permeability decreases according to a certain degradation profile.Thereby, the degradation profile has been assumed to be with an initial value γ at the cut edge and the degradation skin depth δ. Figure 1 displays such a degradation profile for γ = 0.6 and δ = 1.5 mm evaluated in [2].In [3], the effects of three different cutting techniques, namely punching, waterjet and laser cutting, on the magnetic properties of grain-oriented electrical steel sheets have been investigated, and a loss model based on measurements (using an Epstein frame) has been developed.The result showed that the relative increase in core losses was in the range of 180-120%, 50-40% and 6% for laser, punching and waterjet cutting, respectively.In addition to such approaches, measurements based on sensors that evaluate the magnetic field locally, e.g., in [4,5], have been developed.In [4], the needle probe method was used in combination with a scalar Preisach hysteresis operator in 2D FE simulations (based on the diffusion equation for the magnetic field intensity), and the Preisach distribution function was identified via least squares minimisation using the sequential programming method.In [5], the hysteresis properties of electrical steel sheets in the centre of the test material could be determined based on a single U-shaped yoke electromagnet, two Hall sensors (measuring the tangential surface fields) and a sensing coil around the yoke.Impressive research has recently been presented in [6], where pressure needle probes and micrometric giant magneto-resistance (GMR) sensors were mounted on each of the steel sheets in a laminated core.This allowed measurements to be made of the magnetic hysteresis curves for all the individual steel sheets (the averaged magnetic property for each individual steel sheet; no consideration of the cutting-edge effect).In conclusion, there is currently no method available to locally determine the degradation of electric steel sheets due to punching and/or cutting processes so that material models can be fitted and numerical simulations can be performed with sufficient accuracy.It should be noted that the problem at hand is quite complex, as one can only measure the magnetic field quantities outside the steel sheets, i.e., the stray field.Therefore, the local determination of the magnetic material behaviour requires the solution of an inverse problem.
In a previous publication [7], we demonstrated the successful identification of the nonlinear BH-curve (relation between magnetic induction, B, and magnetic field strength, H) based on Epstein frame measurements using (1) a combined Newton-type method with a full multigrid method as a linear solver (i.e., iterative regularisation via coarse discretisation) and (2) a nonlinear full multigrid method.The aim of our current research is to develop a combined method based on measurements, numerical simulations and inverse schemes to determine a sequence of nonlinear BH-curves to locally characterise the magnetic properties of electrical steel sheets and, thus, be able to quantify the cutting-edge effects.Recently, we have proposed an inverse scheme capable of locally determining the change in the linear magnetic property (the linear relationship between B and H via magnetic reluctance) due to the cutting-edge effects.The approach is based on a quasi-Newton method using Broyden's update formula to compute the gradients of the objective function with respect to the local linear magnetic reluctivity [8].In order to advance the practical relevance, this research paper extends the linear magnetic material model to a nonlinear one and uses the adjoint method (see, e.g., [9]) to compute the gradients of the searched-for parameters, which are the coefficients in all BH-curves.We would like to emphasise that the developed inverse scheme is capable of determining local material parameters in general as soon as the measurement setup can be accurately simulated using a physical model based on partial differential equations.
The remainder of the article is structured as follows.In Section 2, we define the problem and discuss the physical model of a sensor-actuator system with which the measurements will be made (in our tests, these measurements will be provided through numerical simulations).Next, Section 3 describes in detail the inverse scheme, including the formulation of the forward, as well as adjoint, problems, the computation of the gradient and the regularization.Section 4 presents the numerical results, and their detailed discussion is summarised in Section 5. Finally, we conclude and provide an outlook for the next research steps.

Problem Definition
We assumed an experimental setup displayed in Figures 2 and 3.The goal is to identify the nonlinear BH-curve in each subdomain Ω l , especially to find out how the magnetic properties change in the direction of the cutting edges.In doing so, we assumed N s sensors and N p positions of the sensor-actuator system, and for each position, we excited the magnetic field from a minimum to a maximum over N e excitation levels (the amplitude of the current density g).Therefore, the abstract problem may be defined as follows In ( 2), the model operator A is defined according to the partial differential equation (PDE) of magnetostatics, which has to be solved in order to get the magnetic vector potentials, u, and in a postprocessing step, the magnetic flux densities, B (for details, see Section 3.1).In doing so, the objective functional is defined according to where Ω sens s is the domain of sensor s, B p,e,s = ∇× u p,e,s is the computed magnetic flux density and B meas p,e,s the measured one, u p,e,s the computed magnetic vector potential, θ = (θ 1 , . . ., θ N sd ) the vector of parameters of the BH-curves, θ ref the vector of reference parameters of the BH-curves and D a diagonal weighting matrix.Since the entries in the parameter vector θ are of different sizes, the diagonal matrix D scales them to a value of about 1.0.Furthermore, A : V × R n θ → W is the model operator (defined according to the weak form of the PDE, including boundary conditions) with appropriate function spaces, V and W.

Inverse Scheme
Let us introduce the parameter-to-state map, which is implicitly defined via the model equation The Lagrange function is introduced as where < •, • > W denotes the dual pairing between some function space, W, and its dual W * .In the next step, the reduced objective function can be defined as follows Assuming that the parameter to state map S is well defined, and A, J, S are differentiable, we can express the partial derivatives of the reduced functional with reference to each parameter θ i as for any w ∈ W. Using the identity containing the adjoint operator defined according to we obtain Choosing w such that it solves the adjoint equation allows us to rewrite (8) as For the nonlinear BH-curve, we make the following nonlinear function ansatz: with constants c 0 , c 1 and c 2 yet to be determined.The reluctivity is computed via and its derivative via Since we assume a possibly different ν(B) curve in each of the N sd subdomains (see Figure 3), we have n θ = 3 N sd parameters to identify.In addition, the model operator describes the magnetostatic field The weak form reads as with the function spaces W = V ⊆ H 0 (curl) [10].The boundary condition included in the definition of the function spaces W, V is u × n = 0, resulting in closed field lines corresponding to a solenoidal vector field.

Formulation of the Forward Problem
The sensors in the experimental setup shown in Figure 2 are Hall sensors capable of measuring static magnetic fields being generated via DC currents in the coils at different excitation levels.To solve the magnetic field problem according to (16), Newton's scheme in the weak form for each sensor actuator position, p, and excitation level, e, reads as [10] where k is the iteration counter, and η is some line search parameter determined using Armijo's rule [11].Edge finite elements are applied using Nédélec basis functions and hierarchical polynomials from [12] to discretise the continuous H 0 (curl) function space.
To ensure coercivity, a small regularisation parameter is used to scale the artificial mass matrix, whose entries are six orders of magnitude smaller than those of the stiffness matrix (for details, see [10]).

Formulation of the Adjoint Problem
Since A is self-adjoint, the left-hand side of (10) for each position, p, of the sensoractuator system and excitation level, e, of the current density, g, in the coils is (see (18)) For the right-hand side, we explore the Gâteaux derivative and obtain for each position, p, and excitation level, e,  Note that the adjoint equation is linear and, therefore, extremely fast to solve compared to the forward simulation, which has to solve a nonlinear system of equations for all unknown vector potentials in the computational mesh (see ( 18)).

Computation of the Gradient
For the evaluation of the gradients, the first term in (11), according to (3), computes as The second term in (11) with the solution u of the forward Problem (18) (note that B = ∇× u) and the solution w of the adjoint Equation ( 21) is calculated according to The derivative of the magnetic reluctivity ν with reference to the coefficients in ( 23) is computed for each subdomain Ω l according to the chosen ansatz (13) via (24)

Overall Inverse Scheme
According to the derivations in Sections 3.1-3.3,all components are available for the inverse scheme.Since the parameter vector θ is improved iteratively, a stopping criterion and a strategy for the regularisation parameter α in the objective function (3) are needed.To this end, the following relative L 2 error norm is defined: In addition, the choice of the regularisation parameter is crucial in obtaining an optimal solution during the iterative process.If the regularisation parameter is set too high, the minimisation will prioritise the regularisation term, resulting in a strong bias, while setting it too low can lead to instability and the divergence of the iterative process.Depending on the accuracy and resolution of the sensors, there is an a priori upper bound β for the error norm defined in (25) with B p,e,s = B exact p,e,s .Using the discrepancy principle [13], an initial regularisation parameter α init is chosen, and it is reduced at each iteration step via until (25) with B p,e,s = B comp p,e,s is smaller than β = ε stop .In all of the following calculations, the values a = 0.9 and α init = 10 were chosen.

Numerical Results
Cutting electric steel sheets causes a severe degradation of the magnetic properties in the immediate vicinity of the cut edges.For this reason, the 0.5 mm-thick electrical steel sheets are divided into five subregions (see Figure 3), and the density of the regions is significantly higher due to the strong material gradient in the immediate vicinity of the cut edges.The affected region is assumed to have a total length, x CE , of 3 mm, and the lengths, ∆x l , of the respective subregions, Ω l , are summarised in Table 1.This distribution of subdomains means that Ω 5 represents the bulk material (unaffected by cutting-edge effects), and it is assumed that the material data for this subdomain are already known from corresponding measurements with an Epstein frame or a single-sheet tester (SST).The sensor-actuator system is equipped with N s = 7 Hall sensors evenly distributed along the straight line [(−1.5, 0.4, 0), (1.5, 0.4, 0)] in mm.The Hall sensors are capable of directly measuring all three spatial components of the stray field (B x , B y , B z ).The electrical steel sheet is measured at N p = 5 different positions so that the Hall sensor s 4 (see Figure 2) is once directly over the cutting edge and once over the centre of each of the subdomains.Since the material behaviour of the two sheets is symmetrical, only the subdomains in the positive x-domain are measured.The resulting positions are x p = [0, 0.125, 0.5, 1.125, 2.25] in mm, where x p is the x-position of sensor s 4 .At each position N p , the sheets are magnetised at N e = 15 excitation levels by varying the current density, g, of the excitation coils.The different current density levels, g(e), are defined according to the following relations where |g init | = 10, 000A/m 2 is the initial current density and e = 0.1, 0.2, ..., 1.5 the excitation levels.
Each subdomain is associated with a nonlinear magnetic material behaviour according to (12), where the parameters c 0 , c 1 and c 2 are the unknowns to be determined in the inverse scheme.Since Ω 5 is the already-known bulk material, 12 unknown parameters must be determined θ = [θ l,i ], where l = 1, 2, 3, 4 is the corresponding subdomain, and i = 0, 1, 2 is the parameter index.For example, θ 1,0 is the material parameter, c 0 , of the subdomain Ω 1 .For the numerical study, the parameters of the nonlinear material function, ν l , are classified into three cases, ν exact   2. The excitation levels are marked with black, solid, vertical lines.
The measured data, B meas , are generated by artificially applying the forward simulation with the exact material data, ν exact l , as listed in Table 2, and the defined excitation levels (see ( 27)).The measurements are taken at the previously defined positions.Thereby, measurement noise (Gaussian white noise, N (0, σ 2 ) with σ = 10%) is superimposed on the resulting data.To avoid an inverse crime, a finer computational grid of the sensor-actuator system was used to generate the measurement data compared to the model for the inverse scheme (see Figure 5).The stopping criterion, ε stop , for the iterative procedure can be determined a priori.To do this, the FE model for the inverse method (Figure 5b) is used to calculate the magnetic flux density, B, based on the exact material data, ν exact l .Here, we carried out a mesh independence study to obtain reliable magnetic flux density data at the sensor position, resulting in a mesh of approximately 100,000 unknowns (82,000 hexahedral and wedge elements).Based on these results and the previously simulated measurement data B meas on the fine grid (see Figure 5a) superimposed by Gaussian noise with σ = 10%, the stopping value according to (25) gives ε stop = 0.1023.The optimisation algorithm is a low-storage BFGS from the NLopt software library [14], to which the computed gradient (11) and the relative L 2 error norm (25) are passed at each iteration step.In addition to the stopping criterion, a relative error |∆θ k |/|θ k | < θ tol_rel = 0.001 and a relative residual ∆ε k /ε k < ε tol_rel = 0.001 are specified to allow the early termination of the optimisation if no further improvements are achieved.
Based on the given data, the sought parameter vector, θ, is calculated using the proposed method, and the results are presented below.The optimisation was stopped after 48 iterations with a remaining relative error of ε = 0.1048.Since the specified stopping value, ε stop = 0.1023, could not be reached, the additionally defined tolerances, θ tol_rel and ε tol_rel , were applied.The relative error of the unknown parameter vector, θ, in each iteration step is shown in Figure 6 and the general convergence behaviour of the relative L 2 error ε in Figure 7.  |, are listed in Table 3.The corresponding BH-curves are compared in Figure 8 to better illustrate the difference between the optimised parameters, θ opt , and the exact parameters, θ exact .8. BH-curves for θ exact (solid line) and θ opt (dashed line) for the different subdomains Ω l , based on the colour code as defined in Table 2.
Finally, we would like to emphasise that the computational grid used for the inverse scheme resulted in about 39,000 unknowns (30,000 hexahedral and wedge elements), and the total time for obtaining the results (all 48 iterations, see Figure 7) was 2 h on a PC with 64 GB of memory (10 GB were used for the computation) and five cores.The computation for one FE solution of the nonlinear magnetostatic PDE needed about 1.5 min.This clearly demonstrates the efficiency of the developed inverse scheme.In addition, the scheme does not rely on the low-storage BFGS , and other optimisers are being tested.The computation of the forward and adjoint problem was performed with the open-source FE solver software openCFS [15].

Discussion
The convergence behaviour of ε in Figure 7 shows two characteristic peaks at iteration steps 20 and 38, preceded by a region of barely noticeable change from the residual.These anomalies are also reflected in the behaviour of the relative parameter error, θ rel l,i , in Figure 6.For the iteration range between 12 and 19, there were no more significant changes in the parameters, although the parameters θ 1,1 and θ 4,1 still showed pronounced errors.The same can be observed for the iteration range of 30 to 37, with the exception that only the parameter θ 4,1 shows increased deviations.Due to the underlying nonlinear problem, these anomalies indicate local minima where the optimisation procedure gets stuck and, therefore, there is no improvement in the relative L 2 error or the relative error of the parameters θ rel l,i .To overcome these local minima, the optimiser makes a significant change to the parameters, which is reflected in the subsequent deviation of the parameters.Thus, after iteration step 20, the relative error of the parameters θ 1,1 and θ 4,1 can be significantly improved.After iteration step 38, however, this did not lead to success, so the optimisation was stopped according to the additionally defined stopping criteria, although the predefined ε stop was not reached.
Comparing the BH-curves of the exact and optimised parameters, as shown in Figure 8, it can be seen that the curves based on the optimised parameter are lower than those based on the exact parameter.This behaviour is mainly due to the following points.Firstly, the measurement data used are superimposed with Gaussian white noise, which allows a solution to be found only in the vicinity of the exact solution.Secondly, due to the illposedness of the inverse problem and the consequent need to use a regularisation (Tikhonov regularisation), the regularisation term and, thus, the a priori information are included in the calculated parameters.This results in the remaining deviations listed in Table 3.However, the increased relative errors, θ rel l,i , for the subdomain Ω 4 are noteworthy.This could be due to the chosen excitation levels or the positioning of the sensor-actuator system in the sense that the sensitivity to the parameters in this subdomain is less pronounced compared to the others.These parameters were chosen arbitrarily.In order to use the optimal excitation levels and positioning of the sensor-actuator system with respect to the determination of the unknown material parameter vector θ, a design of experiment (DOE) can significantly improve the current convergence behaviour.

Conclusions
In this research, we have presented an inverse scheme capable of determining the nonlinear magnetic material properties locally in the vicinity of the cut edges of electrical steel sheets.The approach was based on FE simulations of the magnetostatic field using the magnetic vector potential, u, and it locally identifies the magnetic reluctance as a function of the magnetic flux density, B. The objective function, which penalises the difference between the measured and simulated magnetic flux densities, is minimised under a constraint defined according to the corresponding partial differential equation for the magnetic vector potential.To efficiently obtain the gradients of the objective function with respect to the material parameters, the adjoint method is applied.The optimisation algorithm is a lowmemory BFGS from the NLopt software library, the forward and adjoint formulations are solved using the open source FE solver openCFS [15], and the ill-posedness is handled via Tikhonov regularisation, in combination with the discrepancy principle.The numerical tests show promising results, and in future work, we will build the sensor-actuator system to perform real measurements and apply the developed inverse scheme to these data.

Figure 1 .
Figure 1.Degradation of magnetic property (in this case, magnetic permeability) towards the cutting edges.

Figure 2 .
Figure 2. Sensor-actuator system with two electrical steel sheets, denoted as sample 1 and sample 2.

Figure 3 .
Figure 3. Electrical steel sheet discretisation into N sd subdomains Ω l with l = 1, 2, ..., N sd (colourcoded); each subdomain assigned with reluctivity ν l , being a nonlinear function of the magnetic flux density |B|.x CE , is the affected area due to cutting, and ∆x l is the length of the subdomains.

l
, ν init l and ν ref l .Thereby, ν exact l defines the theoretically real material behaviour and serves as a reference in the sense of calculating the relative errors between the optimised ν opt l and the exact parameters.Furthermore, ν init l is the initial configuration for the iterative inverse scheme, and it has a 50% deviation from the real configuration, ν exact l .In addition, ν ref l is the reference configuration used for the Tikhonov regularisation, and it is here chosen as the arithmetic mean between ν init l and ν exact l .The parameters for the different cases are listed in Table2.The representative BH-curve for ν exact l and ν init l is shown in Figure4.

Figure 5 .
Zoomed FE model of the sensor-actuator system for (a) computing the measurement data and (b) used for the inverse scheme.

Figure 6 .
Figure 6.Relative error of the searched-for parameters θ.

Table 1 .
Length ∆x l for each subdomain Ω l .