Incorporation of a 3-D Energy-Based Vector Hysteresis Model Into the Finite Element Method Using a Reduced Scalar Potential Formulation

In this work, the reduced scalar potential is utilized to derive a finite element formulation capable of handling hysteretic nonlinear material laws. The Fréchet derivative is employed to deduce a quasi-Newton scheme in the weak form for solving the nonlinear partial differential equation that describes the magnetostatic field. Methods for evaluating the Jacobian are presented, and their performance is compared on a 3-D domain under uniaxial and rotational excitation. The numerical results demonstrate the necessity of a flexible approximation to overcome the non-uniqueness of the Jacobian at reversal points, which naturally occurs in hysteresis loops. Consequently, an exact or excessively localized evaluation would give rise to difficulties in states of material magnetization characterized by these reversal points.


Incorporation of a 3-D Energy-Based Vector Hysteresis Model
Into the Finite Element Method Using a Reduced Scalar Potential Formulation Lukas Daniel Domenig , Klaus Roppert , and Manfred Kaltenbacher Institute of Fundamentals and Theory in Electrical Engineering, TU Graz, 8010 Graz, Austria In this work, the reduced scalar potential is utilized to derive a finite element formulation capable of handling hysteretic nonlinear material laws.The Fréchet derivative is employed to deduce a quasi-Newton scheme in the weak form for solving the nonlinear partial differential equation that describes the magnetostatic field.Methods for evaluating the Jacobian are presented, and their performance is compared on a 3-D domain under uniaxial and rotational excitation.The numerical results demonstrate the necessity of a flexible approximation to overcome the non-uniqueness of the Jacobian at reversal points, which naturally occurs in hysteresis loops.Consequently, an exact or excessively localized evaluation would give rise to difficulties in states of material magnetization characterized by these reversal points.

I. INTRODUCTION
T HE consideration of local magnetic properties is a crucial aspect when simulating highly optimized electric motors and transformers.One challenge in this area is the efficient and robust incorporation of vector-based magnetic material laws into the finite element method (FEM), which must also account for hysteretic effects.The Jiles-Atherton (JA) model [1] has been implemented by various authors using a differential permeability approach described in [2] for the reduced scalar potential formulation employing the Newton method [3], [4].A similar approach to incorporate the JA model has also been done for the vector potential formulation in [5], [6], and [7].However, the JA model originates from a differential equation that describes the differential permeability/reluctivity and provides, therefore, an elegant way to use it in the finite element formulation.For energy-based hysteresis models, which are used in this work to describe hysteretic effects, this is not the case anymore.Nevertheless, in the literature, there exist approaches that use this kind of hysteresis model and use it in the FEM for vector and reduced scalar potential formulations [8], [9].However, they lack a description of the evaluation of the Jacobian matrix for the Newton scheme and its convergence behavior.In addition, other models have been incorporated into the FEM, such as the inverse G model in [10] and the parametric algebraic model (PAM) in [11].Nonetheless, again, the Jacobian matrix can be either represented analytically or the description of the numerical evaluation scheme is not discussed in detail.Furthermore, since the differential permeability is a tensor, there has to be more than one method to compute it.In [12], the convergence of the Newton method is analyzed, and the authors employ automatic differentiation to evaluate the differential reluctivity tensor.This may be a way that works.However, it gives little insight into the actual behavior, and it is not compared to other methods.Moreover, automatic differentiation canot be used in every programming language or every FEM software package and is therefore less accessible.For this reason, this article aims to show different methods to compute the Jacobian and compare them.In addition, this article provides a derivation of the nonlinear reduced scalar potential formulation for the Newton method by using the Fréchet derivative.The article is structured as follows.Section II provides a description of the employed hysteresis model, presented in operator notation as a history-dependent constitutive law.Also, the techniques used to identify the parameters in this work are described.Section III investigates the nonlinear reduced scalar potential formulation and presents the derivation of the Newton formulation, utilizing the Fréchet derivative.The subsequent section provides a more detailed description of the numerical implementation of the formulation.This includes the evaluation of the Jacobian using quasi-Newton methods and the determination of the line search parameter required to modify the nonlinear scheme.In Section IV, two problems are solved using the derived scheme, and all methods for evaluating the Jacobian are compared in terms of their ability to accurately approximate the Jacobian and reduce the residual within an efficient number of nonlinear iterations.These two problems involve a 3-D domain subjected to the presence of a uniaxial and a rotating magnetic field, as well as the 3-D example of the TEAM problem 32.

II. ENERGY-BASED VECTOR HYSTERESIS MODEL
The conservation of energy within the context of magnetic fields has been demonstrated in previous studies [13], [14] and is where u(M) is the internal energy, d denotes a dissipation functional, H corresponds to the magnetic field intensity, and B represents the magnetic flux density.To accurately capture ferromagnetic behavior, it is essential to appropriately model the functionals u(M) and d.A detailed modeling of these functionals can be found in [15].There it is stated that the reversible part of the magnetic field is governed by the anhysteretic curve, while the dissipation functional d is designed to model the overcoming of pinning forces that arise due to non-magnetic inclusions within the material, hindering the movement of magnetic domains and, consequently, impeding magnetization changes [16] using a dry-friction model.The parameters of the dry-friction model are a set of weighted pinning forces (ω i , κ i , where i = 1, 2, 3, . . ., N ), which attributes for the force that needs to be overcome in order to change the magnetization.The solution of the resulting governing equation is obtained by reformulating it as an optimization problem.This can then be efficiently solved by using a Newton-Raphson scheme, as described in [17].
In the operator notation, the operator M (H) is defined, which outputs a magnetization depending on the magnetic field.
Using the general constitutive law of magnetics results in the operator B(H) that takes the vector H as an input and yields the corresponding magnetic flux density B according to the underlying hysteresis model This operator is used in the finite element formulation as a law that describes the material involving hysteresis and saturation effects.

A. Parameter Identification
To fully identify the parameter of the hysteresis model, the parameter of the used anhysteresis curve and the set of pinning forces κ and weights ω are required.In this study, the double Langevin function [18] is employed as the anhysteretic curve, which is characterized by four parameters where is the so-called Langevin function.Uniaxial measurements of the major loop are used to calculate the measured anhysteretic curve, as shown in Fig. 1.
A comparison is then made between the calculated anhysteretic curve and the simulated one using the least squares error, and the parameters are identified by minimizing this error.To determine the values for the set, (ω i , κ i ), a method developed in [19] is employed.This method enables the derivation of a continuous distribution function ω(κ), which is subsequently discretized optimally to obtain the pinning forces and weights.

III. REDUCED SCALAR POTENTIAL FORMULATION
To solve the magnetostatic field by means of a scalar potential approach, the magnetic field strength is decomposed into a rotational and an irrotational component [20] where H s represents the rotational part that satisfies Ampere's law (∇ × H s = J s ).By substituting (5) into Gauss law for magnetic fields, the partial differential equation to be solved is where H s is a source magnetic field produced by the applied current density and φ denotes the reduced scalar potential.
To derive the weak form, ( 6) is multiplied by appropriate test functions φ ′ , integrated over the entire domain, and integration by parts is applied.
For a given H s , find u ∈

A. Generation of a Source Magnetic Field Strength
In general, the source magnetic field strength can be computed by any method that fulfills Ampere's law [21].However, when subtracting the source field from the gradient field of the reduced scalar potential, cancellation errors may occur.To mitigate these errors, it is recommended to employ edge elements for representing H s .Hence, in this study, the proposed approach involves solving the magnetostatic Aformulation, where the Nédélec element [22] is the natural choice for discretizing the computational domain and solving the problem.The strong form of the formulation is given by where A s is the source magnetic vector potential, ν 0 is the vacuum reluctivity, and J s is the source current density.The weak form is given as follows. Find Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
for all A ′ ∈ H 0 (curl, ).The parameter σ is an artificial conductivity to avoid a gauging condition, which is necessary to guarantee the uniqueness of the solution [23].The source magnetic field is then

B. Nonlinear Scheme
To resolve the nonlinear and hysteretic behavior of the material, a Newton-Raphson method is applied, that is, where η is a line search parameter and s is the Newton direction.The residual F in the case of the φ-formulation using the material law (2) and the hysteresis operator with the magnetization as output M (H) is The Jacobian ∂F is derived by the use of the Fréchet derivative [24], [25], a generalization of the derivative of a real-valued function with only a single real variable to functions that are vector-valued and have multiple variables.This enables the approximation of the differentiation of the residual Substituting the residual ( 13) into ( 14) results in where Term 1 is the gradient of the Newton direction and Term 2 is approximated as follows: In doing so, the Fréchet derivative is where ∂B(H)/∂ H is the Jacobian that uses the magnetic flux density as the output of the hysteresis operator.Therefore, the Newton scheme reads as follows.
For a given H s , find s ∈ V ( ) such that for all φ ′ ∈ W ( ) until the convergence criterion ∥F ∥ < ϵ is fulfilled, where and the Newton step is damped by a line search parameter η IV. IMPLEMENTATION A. Evaluation of the Jacobian Matrix ∂B/∂ H One of the main challenges with the proposed nonlinear scheme is the evaluation of the Jacobian in (19), which involves the material tensor ∂B(H)/∂ H.This tensor is defined as the partial derivative of the nonlinear operator B(H) with respect to its input.The literature offers several methods to approximate this material tensor using secantbased approaches, which results in quasi-Newton methods.The simplest approach is to use the last and the current point of B(H) and H to obtain a secant, which represents a simplified version of the method of finite differences.This numerical approximation (simpleFD) of the Jacobian is here defined as where To avoid a division by zero, there is an offset of 1 × 10 −8 A m −1 added to H.This offset is small enough to not affect the results of the simulation.On the other hand, the real definition of finite differences can be used in the form of central differences (FD), that is where ϵ is a small increment.Here, we chose ϵ = 10 −13 A m −1 , and e j is the unit vector with i, j = x, y, z.
Another method is the Broyden method [26] that is of the family of update formulas that use the Jacobian from the last nonlinear iteration as information for the current one based on the secant equation The Broyden method minimizes a Frobenius norm (Euclidian norm for matrices) of the form which results in the following updated formula of rank 1: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The Davidon-Fletcher-Powell (DFP) [27] formula is rank 2 and preserves symmetry and positive-definiteness and, therefore, finds a solution to The update formula then is Since both update ( 29) and ( 26) require the Jacobian from the previous step, in the initial step, the Jacobian can be approximated using the finite difference approaches FD and simpleFD, which serve as good starting points.

B. Evaluation of the Line Search Parameter η
Generally, the quasi-Newton method does not exhibit global convergence [28] and requires modification.Therefore, a line search parameter η is introduced, which must be selected to achieve global convergence.To achieve this, there are several methods available.The fundamental concept is to modify the length of the Newton step in the Newton direction to meet specific criteria.The simplest and most straightforward approach is decreasing the residual monotonically in every Newton iteration known as simple decrease To achieve this behavior, a specific set of values for η is used to calculate the residual.For example, η l = 1/τ l , where τ controls the rate of decrease of η.Usually, τ is chosen to be (1.6, 2 or 10).This method works for most cases; however, it may have limitations.There are scenarios, where oscillations of the residual can be problematic as there may not be a significant decrease compared to its step length.To address this issue, one approach is to enforce a sufficient decrease in the residual, known as Armijo's rule where α is a tuning parameter that should be chosen rather small (e.g., α = 10 −4 ) [29].Using this technique, a good value for the line search parameter η can be determined without having high computational costs.However, the parameters τ and α must be chosen appropriately to achieve a satisfactory line search.This may not be a robust solution for an arbitrarily chosen geometry, excitation, and material parameter configuration.A very robust, yet less efficient method is to perform an exact line search where, in every Newton iteration, the energy functional E , which minimizer is the solution of the problem, is minimized along the Newton direction s.The energy functional E is derived with respect to η, which gives a first-order optimality condition [30] The root of the resulting equation, which just consists of the residual multiplied with the current Newton direction, is sought by using Brent's method [31].This is a robust and fast 1-D root-finding algorithm.The root is always the optimal step length η opt , which ensures the convergence of the nonlinear scheme.Finding the root is sufficient with Brent's method, even if the termination criterion is set to be inaccurate.In the experiments, accuracy to the second decimal was always sufficient, which makes it competitive with the other mentioned methods.

V. NUMERICAL RESULTS
This section presents two examples aimed at testing the formulation in 3-D.The objective of these examples is to compare the methods described for evaluating the Jacobian and the line search technique discussed in Section IV-A.The focus of this analysis is on determining the number of iterations required to solve the nonlinear system of ( 13) using the quasi-Newton scheme ( 19)-( 21) depending on the chosen method to evaluate the Jacobian and the line search method.

A. Step-Lap Joint Transformer
This 3-D setup presents a realistic scenario in magnetics involving a magnetic field with a non-negligible z-component.The core of the transformer is formed by step-lap arranged steel sheets, causing the magnetic field to avoid the air gap.In Fig. 2, the field can be observed in a 45 • cut through the L-joint, demonstrating the switch to a different lamination of the core.In addition, a 3-D view is provided, which highlights the arrangement of the steel sheets, which causes the air gap.The edges of the sheets that are seen from this perspective Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II PARAMETERS OF THE USED EB MODEL FOR THE TP32 REGARDING THE SET OF PINNING FORCES AND WEIGHTS
are marked with a black line.The surrounding air is colored gray.For the sake of simplicity, the air gap length and steel sheet width were chosen to not reflect reality.This is because resolving the realistic width of a steel sheet is computationally costly and is not within the scope of this article.Fig. 3 shows the mesh used, which includes the air domain, core domain, and three coil domains for each coil of the three phases (u, v, and w).This enables an observation of alternating and rotational fields.To enforce an alternating field, coils u and w are connected in series and supplied by a sinusoidal current of the form i u,alt = i w,alt = (IN/A coil ) min(t/(7π ), 1) sin(t) (33) where IN/A = 5 × 10 4 A m −2 and coil v is not supplied.
In addition, the amplitude is ramped to see the effect of minor loops.For the rotational field, the coils (u, v, and w) are supplied by sinusoidal input currents that are 120 • phaseshifted where IN/A coil = 10 × 10 4 A m −2 , IN represents the ampere-turns, and A coil represents the cross section of the coil.For the alternating excitation, the measurement point P meas,alternating is placed in the vicinity of the step-lap joint to obtain a good understanding of the 3-D characteristics of this geometry setup.In the case of a rotating excitation, a point in the neighborhood of the T -joint is investigated such that the rotating field can be observed.Also, for simplicity, the anhysteretic curve for the magnetization is chosen to be a simple arctan() function instead of the abovementioned double Langevin function where M s is the saturation magnetization and A is a shape parameter.They are set to M s = 1.23 × 10 6 A m −1 and A = 38 A m −1 .For the parameter that attributes for the hysteresis, ten linear increasing κ i 's are chosen, with a minimum of 0 A m −1 and a maximum of 140 A m −1 .For the weights ω i , an evenly distributed set is chosen ω i = 1/10 ∀i. 1) Alternating Excitation: In the L-joint, the shape of a typical hysteresis loop can be observed in all three directions, as displayed in Fig. 4. Furthermore, it can be seen that the magnetization saturates in all directions, despite the saturation magnetization is not reached in a single direction.In the hysteretic case, there are reversal points due to the tapering of the B-H curve, which makes the evaluation of the Jacobian matrix difficult.Also, it may change suddenly, as seen in the reversal point in Fig. 4. For this reason, an overly precise evaluation of the Jacobian may cause difficulties when approaching such a situation.This is illustrated in Fig. 5, where the black vertical lines indicate a reversal point.When evaluating the Jacobian matrix using the finite difference technique, it may require more Newton iterations in the process.This issue can be avoided by using the Broyden formula.However, the DFP formula fails for higher amplitudes of the input current, as the nonlinear scheme does not converge properly.This problem is described in [32], where it is noted that the ratio between Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the input and output increments should not be too small or too large.In this scenario, the increment of the magnetic field intensity H can be significantly larger compared to the increment of the flux density M due to the saturation effect.The mean number of iterations using the Broyden method in this example is 4.50, while the finite difference method needs 9.99.Regarding the selection of the line search method to choose η, it appears that the simple decrease [see (30)] is not robust, as the residual begins to oscillate.Armijo's rule works well, but the user has to adjust the α parameter until it works, which makes it less robust.The energy minimization method [see (32)] is able to solve the problem without requiring a user adjustment and, in addition, converges fast.
2) Rotating Excitation: For the rotating excitation, the B loci are displayed in Fig. 6.Again, the contribution of the zcomponent of the field is not negligible anymore compared to the xand y-components near the step lap of the core.It can be seen that the field rotates not only in the x y plane but also in the zx and zy planes.Due to the rotating field, there are higher harmonics in the resulting B field.This might yield more reversal points and, thus, difficulties in evaluating the Jacobian using finite differences, leading to a large number of needed Newton iterations.Furthermore, again, the DFP approach fails for higher amplitudes of the input current (see Fig. 7).The mean number of Newton iterations using finite differences is 7.44, while, for Broyden, it is 4.69.Moreover, the energy minimization line search method remains the best choice in terms of robustness and efficiency.

B. TEAM Problem 32
The TEAM problem 32 (TP32) is a benchmark used to validate and compare the accuracy and efficiency of various field analysis techniques [33].The considered geometry is a three-limb transformer with two windings that can be seen in Fig. 8.The parameter identification is done as described in Section II-A and yields the parameter shown in Table I for the parameter of the anhysteretic curve and in Table II for the set of pinning forces and weights.The key aspects of the evaluation of the Jacobian and the use of a line search method are already presented in Section V-A.For this reason, the TP32 serves to validate the accuracy of the implementation, if applied to a practical problem, where comparative measurements are available.Important to note is that this work only focuses on magnetostatics.Therefore, the measured current is used as the input of the model.For this reason, case 4 of the TP32 is omitted.The type of excitation for the transformer depends on three specific cases that are explained in Sections V-B1-V-B3.
1) Case 1 (Uniaxial Flux Density): In case 1, the two windings are connected in series, and a sinusoidal voltage source is applied.This configuration results in uniaxial flux density in each of the three limbs, which is measured by Bcoils C5 and C6.The results regarding the flux density are shown in Fig. 9 for C5 and C6, respectively, using the Broyden method to evaluate the Jacobian.The results demonstrate a good agreement in the B-coil C5 but do not match with measurements for the B-coil C6.This inconsistency is discussed Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.    in detail in [33].The iteration profile is depicted in Fig. 10, where, due to the reversal points, the FD method needs the most iterations.
2) Case 2 (Uniaxial Flux Density With Higher Harmonics): In the second case, the windings are also connected in series but are excited by a voltage signal with higher harmonics.This results in minor loops appearing in the B-H curve.The comparison with measurements can be seen in Fig. 11.Minor discrepancies are present, which may be attributed to an imperfect hysteresis model or imperfect parameter identification.The expected behavior, regarding the necessary iterations, can be observed in Fig. 12, which is due to the higher harmonics even more clear, since more reversal points appear.
3) Case 3 (Rotational Flux Density): In case three, the behavior in the presence of a rotating field is investigated.To do this, the windings are connected in parallel and are  fed by voltage sources that are 90 • phase-shifted.Fig. 13 provides a comparison between the measurements and simulation results.Minor deviations are present, which can be attributed to the assumption that the simulated model behaves isotropically, while the used steel sheets are not.In this configuration of the input current, there are fewer reversal points, resulting in a more effective finite difference method.However, the Broyden method still converges faster than all the other methods tested, as can be seen in Fig. 14.On average, the Broyden method requires 2.76 iterations, while FD requires 3.34 and simpleFD 3.52.

VI. CONCLUSION
This article presents a finite element formulation utilizing the reduced scalar potential to incorporate a vector hysteresis model as a history-dependent constitutive law in the magnetostatic case.The energy-based serves as the material law in this study.To solve the resulting nonlinear system of equations, the modified quasi-Newton method with a line search approach is employed.The Jacobian of the residual, necessary for determining the Newton direction, is obtained using the Fréchet derivative, which generalizes the ordinary derivative of functions.Subsequently, methods for evaluating the Jacobian in the current iteration steps are presented.These methods, except for the finite differences method, only require information from the last Newton step, ensuring computational efficiency.The formulation is then tested on two 3-D problems with a particular emphasis on evaluating the performance of Jacobian determination methods, as it significantly influences the number of nonlinear iterations required.The numerical results demonstrate that all methods effectively decrease the residual, except for the DFP method, as explained in Section V-A.Among the other methods, Broyden consistently performs the best.This is attributed to the presence of reversal points where the hysteresis loop tapers off, causing methods such as FD to struggle and necessitate a higher number of iterations.The authors also suggest that an exact derivative would encounter similar challenges in this region of the hysteresis loop.Thus, methods approximating the Jacobian based on the previous time step are more suitable.

Fig. 3 .
Fig. 3. Mesh and corresponding measurement points of the step-lap joint transformer.

Fig. 7 .
Fig.7.Needed Newton iterations using the different methods to evaluate the Jacobian matrix for rotating fields.The cross denotes a failure of the Newton scheme using the DFP formula.

Fig. 9 .
Fig. 9. Comparison between the measurement and simulation result for search coil C5 (top) and C6 (bottom) for case 1.

Fig. 10 .
Fig.10.Needed Newton iterations for case 1 for all used methods to evaluate the Jacobian.

Fig. 11 .
Fig. 11.Comparison between the measurement and simulation result for search coil C6 for case 2.

Fig. 12 .
Fig.12.Needed Newton iterations for case 2 for all used methods to evaluate the Jacobian.

Fig. 14 .
Fig.14.Needed Newton iterations for case 3 for all used methods to evaluate the Jacobian.

TABLE I PARAMETERS
OF THE USED EB MODEL FOR THE TP32 REGARDING THE ANHYSTERETIC FUNCTION Fig. 2. 3-D view and magnetic field plot of the L-joint from a 45 • cut through the L-joint.