Thermodynamically Consistent Magnetic Hysteresis Model—Application to Soft and Hard Magnetic Materials Including Minor Loops

A novel thermodynamically consistent macroscopic magnetic hysteresis model is presented. Magnetization is calculated from the reversible part of the magnetic field, while the evolution law of the irreversible part involves physically meaningful material constants: the coercive field and the initial susceptibilities of the hysteretic and anhysteretic curves. It is possible to invert the model through an iterative procedure, allowing either the magnetic field or the flux density as an input to the model. The model is tested on both soft and hard magnetic materials for major and minor loops.


I. INTRODUCTION
T HE main aim of a magnetic hysteresis model is to reproduce as closely as possible experimental magnetic measurements of magnetization m m m as a function of the magnetic field h h h, but computation time and simple parameter identification are also crucial criteria.A hysteresis model must respect two fundamental properties [1].First, the memory property stipulates that the output depends not only on the input value but also on its past evolutions.Second, the rate independence imposes time-scale invariance.One of the most famous models that respect these properties is the Preisach model [2], which allows an infinite set of hysteresis operators.In practice, the weights associated with each elementary operator follow a distribution law described by two mean values and two standard deviations for the coercive and interaction field distributions.Modifications were made by Mayergoyz and Friedman [3] to calculate the energy dissipated at any point during the cycle.Bertotti [4] completed this work by calculating the entropy production at each loading variation.Extensions to the vector case were proposed by Mayergoyz [5] and Vinsintin [6], and Torre et al. [7] take anisotropy into account when considering an elliptical critical surface.However, the use of a large number of simplest hysteresis operators in the Preisach formulation comes with a significant computational burden [8], [9].
Henrotte et al. [10] propose another hysteron formulation based on an energy balance at each time step.The weight of each elementary hysteresis operator is associated with a pinning field density, linking micro-and macroscopic behaviors [11].One can consider an a priori unlimited number of parameters identified in a nontrivial process.To reduce computation time, Tellinen [12] proposed a simple scalar model that considers a density of Barkhausen jumps to modify the permeability.The methodology is based on the ratio between major ascending and descending loops.The model gives good results for symmetrical loops but shows limits for asymmetric loops.Li et al. [13] improved the formulation to make it compatible with dynamic loading and asymmetric cycles.In their work, the constitutive laws are empirical since the formulas depend on the level of induction of the starting minor loop.Similarly, Zirka et al. [14] propose to reconstruct the inner curves of the cycle using splines.The law between the induction and the magnetic field still depends on the level of induction.
Hauser [15] proposed a description at the local scale with dissipated energy that is proportional to the wall motion.Liu and Li [16] have extended their work to make it applicable to different levels of induction.For this purpose, an arbitrary function relates the dissipation coefficient to the induction level.The models are limited to symmetric hysteresis loops and need measurements of several minor loops for identification.
A computationally efficient approach is to consider energy balance at the macroscopic scale rather than at the local scale, as suggested by Jiles and Atherton [17].Their formulation is based on a decomposition of the magnetization into anhysteretic and hysteretic (or irreversible) parts.They also introduce the notion of an effective field, which can be associated with the reversible field.As dissipation was not ensured by the initial model, resulting in a negative slope when unloading that should not occur at low frequencies, Bergqvist [18] proposed a modification of the expression for irreversible magnetization so that magnetization variations follow those of the magnetic field.Nevertheless, this is a posteriori correction that is not intrinsic to the model.Zirka et al. [19] mention that the main limitation of the Jiles-Atherton model is that the so-called irreversible part of the magnetization cannot be directly related to the entropy of the system.This point is also discussed in detail in Crew et al. [20].Harisson [21] has demonstrated that a partition of the magnetic field overcomes this drawback.
A usual approach in mechanics is to start with the strictly positive nature of the dissipation, so that the constitutive laws are restricted to respecting this property by the second principle of thermodynamics, similar to plasticity and damage models [22], [23].It is, therefore, a logical step to use the same starting point to write magnetic evolution laws.Hirsinger et al. [24] propose a model similar to the Jiles-Atherton model with the dual decomposition of the magnetic field into a reversible part h h h r and an irreversible part h h h i , emphasizing the need to ensure dissipation without proving it.
Miehe et al. [25] assume that the magnetization evolves linearly and reversibly until the magnetic field reaches the intrinsic coercive field value H c which corresponds to a threshold above which an irreversible (or remanent) magnetization m m m i appears.The total magnetization is then expressed as the sum of m m m i and a reversible magnetization m m m r .The evolution law of ṁ m m i , denoting the time derivative of m m m i , derives from a dissipation potential.The main limitation of this model is that the magnetization curve up to the intrinsic coercivity is linear in the magnetic field and reversible, which does not follow experimental observations.
We propose here to consider that magnetization depends only on the reversible part of the magnetic field, analogously to mechanical stress which depends only on the reversible part of deformation.The advantages of the Jiles-Atherton formulation, that is, ease of implementation, reduced number of parameters, and low computation time, are retained since the model consists of solving an ordinary differential equation.However, the novel approach proposed in this article ensures dissipation and intrinsically relates the first magnetization to the full hysteresis cycle with a reduced number of parameters.Validations on soft and hard magnetic materials with different loading cases are presented for major and minor loops.An extension to multiphysics coupling is possible to consider thermal and stress effects.

II. METHODS
The notations presented in the introduction have been retained.Where the meaning given to a physical quantity is different, this is specified in the text.

A. New Evolution Law
As mentioned in the introduction, one can consider a partition of the total magnetic field as where h h h r denotes the anhysteretic part and h h h i the irreversible and therefore dissipative part.Following the example of previous authors [21], [24], we consider magnetization to be related to the reversible part of the magnetic field such that it follows the constitutive law where M an describes the evolution of the magnetization as a function of the reversible field.The irreversible field h h h i is an internal state variable.One can postulate the existence of a thermodynamic potential from which the equations of state derive so that a dual variable denoted m m m i is defined.It opposes and causes a delay in the magnetization as the walls of the magnetic domains cross the pinning sites.The Clausius-Duhem inequality leads to the definition of magnetic dissipation [24] To guarantee D ≥ 0, one can consider the evolution law where χ χ χ is a positive definite magnetic susceptibility tensor and λ a positive multiplier.For ferrimagnetic materials, χ χ χ is negative definite, and λ is a negative multiplier to ensure dissipation.The dissipation can then be written as which is always positive.It now remains to specify m m m i and λ .
The variable λ describes the rate of change of the ratio of pinning sites that can be crossed when the domain structure evolves.Assuming that each pinning site requires the same amount of energy to be crossed, the ratio of pinning sites opposing the applied field is 1 − ∥h h h i ∥/H c , where H c corresponds to the coercive field.When there is a change in the direction of the field, the pinning sites previously crossed must be crossed again in the other direction.Consequently with τ 1 being a positive time constant and e e eḣ h h = ḣ h h/∥ ḣ h h∥ the magnetic field evolution direction.This quantity is maximum when no pinning sites have been crossed and tends toward zero when there are no more pinning sites to overcome, that is, the process is reversible.Note that (6) has the nonlocal memory property introduced by Mayergoyz [5]

since future values of h h h i also depend on past values of the input h h h.
The dual variable m m m i results from the application of the irreversible field h h h i .Thus, when the magnetization m m m does not vary, the hysteretic and anhysteretic curves merge so that m m m i = 0.As a first approximation, we can, therefore, assume a linear dependence between m m m i and ṁ m m, so that with τ 2 being another positive time constant.This coefficient is similar to a constant viscosity.In necessary cases, a dependence on other physical quantities (temperature, mechanical stress. ..) or even a time evolution can be included.A change in magnetic properties due to interaction with other physics, particularly in the initial state, can thus be taken into account.Combining ( 4), (6), and (7), we obtain The formula of ( 7) is deliberately chosen so that the variations of h h h i depend on the variations of m m m.It is possible to add a rate-dependent part for the description of the viscose-like phenomena such as excess losses.The full hysteresis model is therefore written as The calculation consists of solving a first-order differential equation.The reciprocal form, which takes the magnetization as input, is easily defined by assuming the inverse anhysteretic function M −1 an is known.If no explicit expression of this function is available, a fixed-point iteration can be used as an approximation.
In many applications, flux density b b b is considered as input.In this case, the implementation of the model is more challenging since the constitutive law relates m m m to h h h.An iterative resolution has to be deployed: for each time instant t, the magnetization and magnetic field are calculated using the previous solutions through the Newton

B. Identification of the Parameters in 1-D
We assume that the field and magnetization measurements are taken in the same direction n n n, parallel to the magnetic field.We thus define Equation ( 8) is then written, in the uniaxial case as with (1/ χ eq ) = (1/ χ )(τ 2 /τ 1 ).In the neighborhood of the demagnetized state, magnetization and reversible magnetic field H r = H − H i are linearly related such that with χ 0 being the initial susceptibility of the anhysteretic part and, according to ( 12) Moreover, close to the demagnetized state, the magnetization is linearly correlated to H such that with χ 1 being the initial susceptibility of the hysteretic part.From ( 13) to (15), one can identify III. VALIDATION

A. Modeling the Anhysteretic Part With a Langevin Function
The choice of M an in (9) is free.Steentjes et al. [26] have shown the benefits of using a double Langevin function to model the anhysteretic magnetization, in particular, for an accurate description of the knee region.We thus define M an as where M si are related to the saturation magnetization M s such that and a i are material parameters that can be related to the initial susceptibility χ 0 (see Fig. 1) by first-order Taylor series expansion, so that In the case where multiphysical coupling effects are to be taken into account, ( 17) is replaced by a multiscale calculation [27].As this method only serves to describe the  anhysteretic part, the influence of couplings on hysteresis can be introduced by adding a dependence on the parameters H c and χ 1 to the different physics.This is a hybrid approach that differs from the fully multi-scale approach proposed by Vanoost et al. [28].In this article, we only consider the use of (17).

B. Modeling a Hypothetical Material
To verify the general performance of the model, various loading cases are investigated in this section for arbitrary values a = 2000 m/A, M s = 1.5 • 10 6 A/m, H c = 2000 A/m, and χ eq = 100.Fig. 2 shows the magnetization obtained for a sawtooth-shaped magnetic field waveform with an amplitude of ±15 kA/m.All loops close without any deviation.The inverse model (10) has very similar behavior, as shown in Fig. 3 where simulations of the magnetic field response to a magnetization sawtooth-shaped signal with an amplitude of 1.3 MA/m.
The hysteresis model also converges to the anhysteretic one when a decreasing field is applied (see Fig. 4).Indeed, for a final target value H = 4000 A/m, the relative error between the anhysteretic and hysteretic magnetizations is about 0.1%.
The simulation of low-amplitude cycles in Fig. 5(a) shows that minor loops require one cycle to converge to a stable position.This phenomenon is known in the literature as accommodation [29] or reputation [30].It is explained by an initial non-symmetry of the remanent magnetization at the start of loading.Denoting M + R (resp.M − R ) the remanence for decreasing (resp.increasing) H values, the cycles are stable once In the case of an initial value , the cycles are intrinsically stable as shown in Fig. 5(b).When the initial state is demagnetized, several cycles are necessary before converging to a stable state, unless the magnetic field amplitude is sufficiently close to saturation, as in Fig. 2. In the latter case, condition (21) is verified in the first loop, since H i reaches the value H c during magnetization.Therefore, major cycles are intrinsically stable.This particularity contrasts with the intrinsic congruency property of the Preisach model that is not verified experimentally when several cycles of accommodation are required before stabilization [5].Kadar and Torre [31], Torre [32], and Zirka et al. [33] have proposed various modifications to Preisach's original model to overcome this inconsistency.Noncongruency is natural here since the predicted magnetic state depends on both the previous magnetization and magnetic field.The number of cycles required before stabilization is related to the value of λ [see (6)].For a given coercive field value, the loops will be more quickly stable if τ 1 is low, that is, if χ eq is low.As a result, simulated hard materials have stable minor cycles faster than soft materials.However, we acknowledge that this behavior is not really a feature of the model that can be controlled but merely a consequence of the formulation itself.

C. Modeling of Different Materials
In this section, the model is tested on different materials with significantly different behaviors.We investigate two soft magnetic materials (galfenol [34] and 10JNEX900 [35]), two permanent magnets (Co 80 Zr 80 B 2 melt-spun ribbons [36] and Vacodym 362 TP) and one hard ferrite (Ferriflex 1 FX3).Measurement data provided in the manufacturer datasheets are used as a reference.The identification of the material parameters is done as follows.
1) M s is the saturation magnetization.When saturation is not reached, the asymptotic value of M is considered.2) H c is defined as H (M = 0).3) χ 0 is approximated using the slope of the major hysteresis loop close to zero magnetization because the anhysteretic curve is not available in the literature for these materials.We compute a using (19).4) When the first magnetization curve is available, χ 1 is identified to compute χ eq using ( 16).Otherwise, it is considered a fitting parameter.For the description of specific materials, a single Langevin function [see (17)] is sufficient, so M s1 = M s2 = M s /2 and a 1 = a 2 = M s /(3 χ 0 ).In the case where two Langevin functions are required, M s2 and a 2 are parameters to be optimized.The values of M s1 and a 1 are then identified using (18) and (19).The different material constants are provided in Table I.
For all the investigated materials, a quasi-dc-field value was imposed, except for galfenol, for which a sinusoidal signal at 0.05 Hz was used.Simulation results of galfenol, 10JNEX900, and Co 80 Zr 80 B 2 melt-spun ribbons are shown in Fig. 6.The relative root mean squared error is given in Table II.The model 1 Registered trademark.seems have no particular problems simulating major loops of materials with widely varying initial susceptibilities and coercive fields since the relative error does not exceed 9%.
As the properties of hard ferrites and permanent magnets are generally given in the form of B(H) curves, both measurements and simulations of induction and magnetization are shown in Fig. 7. Since both are hard magnetic materials, χ 1 has a value close to zero, making it tricky to identify.One can consider χ 1 to be a fitting parameter for both materials.The model seems suitable for describing these cycles since the relative errors in magnetization are 0.7% for Ferriflex and 2.7% for Vacodym 362 TP (see Table II).

D. Minor Loops of Permanent Magnets
We investigated the ability of the model to describe minor loops on permanent magnets.Ruoho et al. [37] carried out several magnetic measurements on a NdFeB magnet (see Table I for parameters).It can be seen in Fig. 8(a) that the model succeeds in reproducing a symmetrical minor loop with a model error of approximately 12%.Regarding the two nonsymmetrical minor loops in Fig. 8(b), a slightly too high demagnetization is computed by the model.This demagnetization is directly linked to χ eq .Indeed, in the case of hard materials, χ eq is approximately equal to χ 1 [see (16)], which is the initial slope of magnetization.The lower this value, the more difficult it is for the material to magnetize and therefore demagnetize.The sharp change in magnetic susceptibility during the second loop is not reproduced at all, which leads to a model deviation exceeding 55% in this region.
Although the model cannot catch exactly the minor loops for positive magnetic field, it models very well the demagnetization behavior for negative magnetic field.This is a very important aspect, especially for the numerical simulation of permanent magnet machines with strong loading capacity, for example, in electric traction systems such as electric vehicles.

E. Experimental Fe-Si Alloy Minor Loops
Quasistatic measurements were carried out on nonoriented Fe-3%Si samples to investigate the ability of the model to reproduce minor loops for soft magnetic materials.A quasi dc-field was imposed by a primary coil positioned around the sample.The field is changing slowly in 24 small steps so that the flow progresses in approximately equal steps, each    integrator.The setup was used by Gürbüz et al. [38] who provided a detailed description.The imposed signal is shown in Fig. 9.
The magnetic response is given in Fig. 10(a).The presence of a remanent magnetization in the initial state is noticeable.By inversion of (20), the model has been initialized with The saturation magnetization of each Langevin function used in (17) and a i parameters are grouped in Table I.Fig. 10(a) shows that the major cycle is correctly reproduced by the model.The focus on minor loops in Fig. 10(b) shows that the evolution and orders of magnitude given by the model are correct.The relative root mean squared error is 5.2%.The slight experimental accommodation is also present in the simulation.We then investigated the loss prediction given by the model for different magnetic field magnitudes.Measurements were performed using the same setup as before.The material parameters of Table I have been retained.Fig. 11 shows a correspondence between experimental and simulated loops for field amplitudes above 0.42 kA/m.However, there is a significant difference at lower amplitudes.Simulations using the original Jiles and Atherton [17] and Bergqvist [18] formulations are presented in Fig. 12.One can see the inherent problem of the initial model already identified in [39] namely the presence of negative slopes in Fig. 12(a).The Bergvist correction in Fig. 12(b) prevents this phenomenon but imposes a non-physical zero slope, which can lead to convergence problems when used in numerical field solvers.
The evolution of experimental field values H (M = 0) with increasing minor loops amplitude H max is given in Fig. 13(a).This appears to be tending toward the H c value identified in Table I.The same applies to saturation magnetization M s in Fig. 13(b).
Measured and predicted losses, as well as the difference between the two, are given in Table III.The last line of the table corresponds to the cycle presented in Fig. 10(a).The loss prediction gap in this case is almost zero since this is the cycle used for parameter identification.The losses predicted by the model saturate after 1.7 T (see Fig. 14).This is consistent with the fact that the material no longer dissipates energy once saturated.
We also considered the sensitivity of the loss prediction to the model parameters on the major loop: we assumed that  the identification error for each parameter could vary between −5% and +5%.Table IV gives the worst-case scenario for a  variation in a pair of parameters (the others are considered to be exactly known).One can see that the highest sensitivity is on the coercive field and saturation magnetization.
Finally, we summarize in Table V the pros and cons of our approach in comparison to the Preisach and Jiles-Atherton models.When a high degree of accuracy is required, the Preisach model remains the most accurate model.On the other hand, to rapidly model a material under complex loading, the proposed method avoids the negative slopes present by the Jiles-Atherton model.

IV. CONCLUSION
A full thermodynamically consistent macroscopic hysteresis magnetic model is provided.It can describe initial magnetization and major and minor loops of a wide range of materials using only easily identifiable physical constants.The anhysteretic part can be modeled freely without any constraints.This work can be extended to couplings by considering the dependence of H c and M an on other physics.The model ( 9) considers a linear evolution of m m m i as a function of ṁ m m.This is an approximation that does not take into account the impact of frequency in excess losses.Implementing it in finite-element solvers and proposing a fully anisotropic evolution law will also be future projects.

Algorithm 1 1 i
-Raphson algorithm to reach the desired value b b b ref as described in Algorithm 1. Computation With b b b ref as Input Input: b b b t ref , m m m t−1 , h h h t−Initialization: k = 0; m m m k = m m m t−1 ; while ∥r r r ∥ ≥ ϵ do use (10) to estimate h h h k i and h h h k ; b b b k

Fig. 10 .
Fig. 10.Minor loops starting from the descending branch.(a) Full magnetic response.(b) Focus on minor loops.

Fig. 11 .
Fig. 11.Experimental and simulated inner loops for various maximum magnetic fields.(a) Measurements.(b) Simulations with the proposed model.

Fig. 13 .
Fig. 13.Evolution of Fe-3%Si properties.(a) Evolution of H (M = 0) as a function of maximum magnetic field magnitude.(b) Evolution of maximum magnetization as a function of maximum magnetic field magnitude.

TABLE III LOSS
PREDICTION Fig. 14.Loss prediction for higher values in magnetic flux.