A Numerical Comparison between Preisach, J-A and D-D-D Hysteresis Models in Computational Electromagnetics

Abstract

The incorporation of hysteresis models in computational electromagnetic software is of paramount importance for the accurate prediction of the ferromagnetic devices’ performance. The Preisach and Jiles-Atherton (J-A) models are frequently used for this purpose. The former is more accurate and can represent a broad range of magnetic materials, but it is computationally expensive. The latter is more efficient but can accurately model only soft ferromagnetic materials. In this paper, a recently proposed hysteresis model, referred to as the D’Aloia-Di Francesco-De Santis (D-D-D) model, is shown to have the best trade-off between accuracy and computational burden. For the first time, a numerical comparison between the Preisach, J-A and D-D-D models is provided for a large class of hysteresis loops including soft, semi-hard and hard ferromagnetic materials.

1. Introduction

Ferromagnetic materials, both soft and hard, are widely used as magnetic cores or permanent magnets in electrical devices, such as inductors, transformers, generators and motors, due to their high permeabilities or magnetization properties [1]. When designing such devices, the knowledge and characterization of their hysteretic behavior are therefore of paramount importance, and having a good mathematical model capable of reproducing this behavior is important as well [2,3]. However, to achieve an efficient design procedure, the used models have to be accurate and fast. Satisfying these two criteria simultaneously is not easy; therefore, a compromise has to be made.

During the past century, a large number of hysteresis models have been proposed and applied in computational electromagnetics [4,5,6,7,8,9,10]. Based on their application level, they can be classified into scalar and vector models, depending on whether or not the field remains collinear at all times, as well as between static and dynamic models, depending on whether or not the constitutive relationships are indifferent to rates. Another classification is based on their accuracy and distinguishes between mathematical or physical models. The former are generally faster but less accurate [4,5], as they ignore the underlying physics of the material behavior, whereas the latter can describe the complex magnetization process more faithfully and are thus more applicable to engineering problems. The physical models can be further classified into physics-based and phenomenological models, depending on whether they are based or not on the actual physical processes occurring in the matter subjected to magnetic fields.

Among the former are the Stoner-Wohlfarth (S-W) [6] and Jiles-Atherton (J-A) [7] models, while the Preisach [8], Play [9] and neural network-based [10] models are usually considered as phenomenological ones. The implementation of the S-W model requires the knowledge of the distribution of the particles in the ferromagnetic material. However, the statistical distribution function is difficult to determine. In order to achieve sufficient accuracy, a complex identification process is needed based on a vast number of experimental measurements, which limits its practical application [11]. The J-A model is therefore the only physics-based hysteretic model that can be easily implemented into the finite element method (FEM) and hence finds a wide use in the design of electromagnetic devices. The original J-A model was a static scalar model, but an extension to dynamic and vector hysteresis models has been proposed in [12,13], respectively.

Phenomenological models are widely used and can be incorporated into FEM simulations as well. The classical Preisach model (CPM) was born as a static scalar model; however, many extensions and modification attempts have been made over the years [14,15,16,17]. For instance, Della Torre included the dynamic effect and proposed a moving Preisach model (MPM) [14], while Cardelli provided an extension to the vector Preisach model (VPM) [15]. Attempts to reduce the computational efforts for either CPM and VPM have been made in [16,17], respectively. Even the Play model and neural network-based models have been improved to account for the vectorization or numerical speed-up of the hysteresis models [18,19,20,21,22,23,24,25].

In this study, only the CPM and original J-A models have been considered, as they are the most popular ones. A numerical comparison between these two models can be found in [26,27]. The static J-A model is considered to be more efficient in terms of computational time and memory requirements. On the other hand, the CPM model is not only more accurate but can represent a wider class of materials because of the arbitrariness of the material function [27]. Besides these, a novel phenomenological model has been proposed by the authors in [28], referred to as the D’Aloia-Di Francesco-De Santis (D-D-D) model. In that paper, only the mathematical background of a scalar model for hard magnetic materials was described, while in this paper, a numerical comparison with the CPM and J-A models will be provided for a large class of ferromagnetic materials. From this comparison, it will be shown that the D-D-D model is the best solution, being as accurate and generic as the CPM and even faster than the J-A model.

2. Materials and Methods

2.1. Materials

A broad range of ferromagnetic materials have been considered for the numerical and experimental comparison in order to challenge the fitting capabilities of the several hysteresis models. Specifically, the measurements of MnZn ferrite provided in [29] and of Fe-Si taken from the TEAM benchmark 32 [30] have been selected as representative of soft materials, while the NdFeB at two different temperatures (80 °C and 27 °C) [28] has been deemed to be representative of semi-hard and hard materials, respectively. Note that only major loops were treated in the fitting procedure, as the interest has been in finding challenging shapes rather than the physical meaning of the hysteresis models.

2.2. Hysteresis Models

In this section, the CPM, J–A and D-D-D models are presented. In their original form, they give a relationship between the magnetic field vector H and the magnetization M of the ferromagnetic material. The magnetic flux density B is then obtained from the constitutive equation given by:

$$\mathit{B}\left(\mathit{H}\right)={\mu }_0\left(\mathit{H}+\mathit{M}\right)$$

(1)

where μ₀ is the vacuum permeability.

2.2.1. CPM Model

The CPM is a phenomenological model that presumes a ferromagnetic material to be composed of a large number of rectangular switches, similar to a Barkhausen jump [8]. The scalar output M(t) of the CPM is calculated from the scalar input H(t) as follows [3]:

$$M\left(t\right)={{\displaystyle \iint }}_T^{}\mu \left(h_1,h_2\right)m\left(h_1,h_2,H\left(t\right)\right)d{h}_{1}d{h}_2$$

(2)

where T denotes the Preisach triangle, μ(h₁, h₂) is the Preisach distribution function and m(h₁, h₂, H(t)) represents the elementary relay operators which can take the values +1 or −1. Depending on the history of the applied magnetic field, two domains can be distinguished on the Preisach triangle: one where the operators are switched down (h₁) and one where the operators are switched up (h₂). The boundary between domains is called staircase line.

Given its vertices, the magnetization can be determined at any instant with the aid of the Everett function:

$$E\left(x,y\right)={{\displaystyle \int }}_x^y{{\displaystyle \int }}_x^{h_2}\mu \left(h_1,h_2\right)d{h}_{1}d{h}_2$$

(3)

where x, y is a vertex of the staircase line. The Everett function can be measured; however, it is common to approximate the Preisach function analytically and perform the double integration (3) numerically. Compared to the original Preisach model, this way can greatly accelerate the computation by avoiding the double integration in (2). A closed form of the Everett function (3) was provided in [31]:

$$E\left(x,y\right)=\frac{a^2}{b^2}\frac{\left(c^2-1\right)\left(-{\mathrm{e}}^{bx}+{\mathrm{e}}^{by}\right)-\left(c+{\mathrm{e}}^{by}\right)\left(1+c{\mathrm{e}}^{bx}\right)\log \frac{\left(1+c{\mathrm{e}}^{by}\right)\left(c+{\mathrm{e}}^{bx}\right)}{\left(1+c{\mathrm{e}}^{bx}\right)\left(c+{\mathrm{e}}^{by}\right)}}{{\left(c^2-1\right)}^2\left(c+{\mathrm{e}}^{by}\right)\left(1+c{\mathrm{e}}^{bx}\right)}$$

(4)

where a, b and c are the parameters of the inverse cosine hyperbolic type Preisach distribution function, as explained in [32].

The reversible component can be considered as M_r = dH, where d is a constant and added to the magnetization produced by the Preisach model. It should be noted that an improved closed form with additional polynomials was proposed in [32], but for the sake of comparison, the original version was hereby adopted to reduce the parameter numbers.

The Matlab implementation of the above-described CPM is provided by Szabó in [33], where a nonlinear least-squares (NLS) algorithm based on the Levenberg–Marquardt method was used for the identification of the parameters.

2.2.2. J-A Model

The scalar J-A model [7] is one of the most famous physics-based hysteresis models that takes domain wall motion into account. The two modes of domain wall transitions (both its bending and translational motions) result in a reversible and an irreversible component of magnetization, respectively. The mutual coupling of the domains, i.e., the pinning of domain walls at impurity sites during the motion and flexibility of the domain walls, is also taken into account. The energy dissipated in overcoming these pinning sites contributes to hysteresis loss. The total magnetization inside a material is computed using an ordinary differential equation (ODE) and is given by [34]:

$$\frac{dM\left(t\right)}{dH}=\frac{c_J}{1+c_J}\frac{d{M}_{an}\left(M_s,a_J\right)}{dH}+\frac{1}{1+c_J}\frac{M_{an}\left(M_s,a_J\right)-M\left(t\right)}{\mathrm{sign}\left(\frac{dH}{dt}\right)k-{\alpha }_J\left(M_{an}\left(M_s,a_J\right)-M\left(t\right)\right)}$$

(5)

where M_an is the anhysteretic magnetization that is computed using Langevin’s polynomial [7], while M_s, α_J, a_J, k and c_J are the five J-A model parameters, which represent the saturation magnetization, the inter-domain coupling coefficient, the effective domain density, the energy-dissipative features in the microstructure and the reversibility coefficient, respectively [35].

The MATLAB code of the J-A model is provided by Szewczyk in [29], where the ode23 function was used to solve (5) and an NLS method based on the Nelder and Mead Simplex algorithm was employed to identify the five parameters.

2.2.3. D-D-D Model

The D-D-D method is a phenomenological model based on the analogy between the change in magnetization caused by a single hysteron and the change in the velocity of disk-shaped solids elastically colliding with each other [28]. By assuming a scalar and time-harmonic applied field H with an angular frequency ω, a simple ODE is obtained [28]:

$$\frac{dM\left(t\right)}{dt}=\gamma sin\left(\omega t+\beta \right) e^{-{\alpha }_{D}co{s}^2\left(\omega t+\beta \right)}$$

(6)

where γ, α_D and β are the three D-D-D parameters. In particular, the former affects the amplitude of the magnetization, whereas α and β influence the slope and width of the obtained hysteresis loop, respectively. It should be noted that the proposed D-D-D model is conceptually similar to the Preisach model but mathematically closer to the J-A model, with an even simpler ODE to be solved and two parameters less to be identified. Since the latter are straightforward to determine, a tentative approach has been undertaken.

3. Results

The curve fitting of the several hysteresis loops for the three considered models is reported in Figure 1, while the coefficients obtained by the identification parameters are summarized in

Table 1. Parameters of the considered hysteresis models for the several materials.

Material	Model	Parameter Number 1
		1	2	3	4	5
MnZn ferrite	CPM	0.049	−8.73	17.495	0.95∙10−3	-
	J-A	3.178∙105	1.099∙10−7	12.649	12.448	0.844
	D-D-D	2.74∙105	5.5	0.116	-	-
Fe-Si loop 11	CPM	0.044	69.12	23.212	4∙10−4	-
	J-A	1.23∙106	1∙10−4	47	66	0.99
	D-D-D	1.194∙106	300	0.07	-	-
NdFeB at 80 °C	CPM	6.6∙10−3	4.56∙105	1.49∙105	0.5∙10−3	-
	J-A	1.18∙106	0.46	1.25∙105	5.24∙105	0.05
	D-D-D	1.13∙106	12.5	0.3	-	-
NdFeB at 27 °C	CPM	0.016	1.24∙106	5.44∙104	0.2∙10−4	-
	J-A	0.954∙106	1.1	1.45∙105	1∙106	1∙10−6
	D-D-D	0.909∙106	20	1.02	-	-

¹ Parameter number reads as follows: CPM: a, b, c and d; J-A: M_s, α_J, a_J, k and c_J; D-D-D: γ, α_D and β.

. As can be observed, the CPM and D-D-D models can fit hard materials and MnZn ferrite well but suffer to reproduce the near saturation region of the major loop 11 in the rolling direction of the Fe-Si. It should be noted that a better fitting for this material has been obtained with the CPM in [32], but a higher order (up to 3) of the polynomial used to discretize (3) has been taken, bringing the number of parameters to be identified up to nine, plus the d coefficient for the reversibility. On the other hand, the J-A model is well suited for soft materials [29] but suffers to fit semi-hard and hard materials. The main reason for this is the natural sigmoid-shaped hysteresis curve coming from Langevin’s polynomial used for the anhysteretic magnetization [35], which is not always suitable for hard materials. Therefore, modifications to the original J-A model have been provided to overcome this issue but complicate the model and the identification procedure, with eight parameters to be identified [36].

3.1. Accuracy

To quantify the accuracy of the proposed models, two metrics are usually employed. The first one is related to the core losses obtained by the area integral A of the hysteresis loops, where the percentage error can be evaluated as:

$$e=\frac{\left|A_m-A_c\right|}{A_m}·100$$

(7)

A_m and A_c are the measured and calculated areas, respectively. The second one is presented to elaborate the goodness of fit of the proposed models (i.e., the Pearson correlation coefficient), which is expressed as [24]:

$$r^2=1-\frac{{{\displaystyle \sum }}_{k=1}^{N_m}{\left(B_{m,k}-B_{c,k}\right)}^2}{{{\displaystyle \sum }}_{k=1}^{N_m}{\left(B_{m,k}-B_{m,avg}\right)}^2}$$

(8)

where N_m is the total number of sampled measurement points, B_m,k is the k-th measured point, B_c,k is the k-th calculated point and B_m,avg is the average value of all measured points. The value range of r² is [0, 1], and the closer it is to 1, the more accurate the fitting is.

Table 2. Numerical comparison of the considered models. Bolded values are the most performant metrics.

Material	Model	e (%)	r2	CT (s)
MnZn ferrite	CPM	4.05629	0.99785	0.09182
	J-A	1.51375	0.99948	0.15134
	D-D-D	4.53686	0.99718	0.00685
Fe-Si loop 11	CPM	15.24575	0.91456	0.08678
	J-A	5.0458 3	0.99212	0.22523
	D-D-D	4.91653	0.97532	0.00976
NdFeB at 80 °C	CPM	0.19651	0.99927	0.07080
	J-A	9.62853	0.93845	0.26254
	D-D-D	2.03716	0.99921	0.00766
NdFeB at 27 °C	CPM	2.18167	0.99584	0.07430
	J-A	31.31527	0.79565	0.35150
	D-D-D	2.16800	0.99872	0.00660

summarizes the accuracy values obtained by the several models for the considered materials. As can be noted, the J-A model results are the most accurate for soft materials but suffer for semi-hard and hard materials, and vice versa for the CPM. Instead, the proposed D-D-D model shows good accuracy for all materials with an error in the core losses always below 5% and a Pearson coefficient higher than 0.97.

3.2. Computational Time

Table 2 also shows the computational time (CT) obtained by the several hysteresis models for the considered materials. They have been implemented in MATLAB R2021a (©The MathWorks Inc.) and executed on a PC with an Intel Core i7-1165G7 CPU running at a clock speed of 2.80 GHz. As can be observed, the D-D-D model is the faster method, followed by the CPM and the J-A model.

4. Discussion

Since we used open-source codes, CTs for all the hysteresis models are evaluated for the single iteration (the final one), without accounting for the identification procedure, which could depend on the adopted strategy affecting the comparison. Moreover, for the CPM, the number of sampling points for the ascending and descending branches affects not only the accuracy but also the CT [27,31,32]. Therefore, a number of N_asc = 100 samples for each branch has been deemed as a good trade-off between accuracy and CT.

As for the CPM, the number of sampling points N_asc is shown to affect the accuracy, even for the D-D-D model, the number of points N_c = 2 N_asc used to discretize the phase shift vector ωt would affect the accuracy. Indeed, it has been experienced that a number N_c < 100 yielded unclosed and asymmetric loops. This could happen with coarse measurement points (N_m < 100) in evaluating the Pearson correlation coefficient (i.e., N_c = N_m). In such a case, an oversampling of both measured and simulated points will avoid this issue.

To further improve the accuracy of the D-D-D model, it is also desirable to refine the sampling step where the gradient of the magnetization is steepest and relax it where it is smoothest, as expected from (6). Once again, this is somewhat not controllable during measurement campaigns but could be overcome in the post-processing phase. Another improvement that could be made is to add a reversible term, as for the CPM, which will give a slope to the magnetization after saturation when needed, such as for Fe-Si loop 11.

Finally, we noted that a tentative approach has been undertaken for the parameter identification of the D-D-D model, but in the future, an automated optimization procedure could be easily afforded.

5. Conclusions

The incorporation of hysteresis models in finite element analysis (FEA) is of utmost importance for the accurate predictions of the ferromagnetic devices’ performance. The Preisach and Jiles-Atherton (J-A) models are typically used for this purpose. The former is more accurate and can represent a broad range of magnetic materials, but it is computationally expensive unless closed forms of the Everett function are used. The latter is more efficient but can accurately model only soft ferromagnetic materials [27].

In this paper, a hysteresis model previously proposed by the authors, referred to as the D’Aloia-Di Francesco-De Santis (D-D-D) model, is compared with the Preisach and J-A models for a large class of hysteresis loops including soft, semi-hard and hard ferromagnetic materials. The results of this comparison have shown that the proposed model represents the best solution, as it is as accurate as the Preisach model and even faster than the J-A model. This makes the D-D-D model the perfect candidate in FEA.