Comparison of Preisach and Jiles–Atherton models to take into account hysteresis phenomenon for finite element analysis

Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Comparison of Preisach and Jiles–Atherton models to take into account hysteresis phenomenon for finite element analysis A. Benabou, S. Cle! net, F. Piriou* # Laboratoire d’Electrotechnique et d’Electronique de Puissance de Lille, Batiment P2, Universit!e des Sciences et Technologies de Lille, Equipe Mecosyel, 59655 Villeneuve d’Ascq, France Received 21 March 2002; received in revised form 6 November 2002

Abstract In electrical engineering, study and design of electromagnetic systems require more and more accurate models. To improve the accuracy of field calculation code, hysteresis phenomenon has to be taken into account to model ferromagnetic material. This material model has to be accurate and fast. In that context, two macroscopic models are often used: the Preisach and the Jiles–Atherton (J–A) models. In this paper, both models are presented. Field calculation requires a model giving the magnetization M versus either the magnetic field H or the magnetic flux density B: Consequently, from the classical Preisach and J–A, two sub-models MðHÞ and MðBÞ are deduced. Then, we aim at comparing these models in terms of identification procedure facilities, accuracy, numerical implementation and computational effort. This study is carried out for three kinds of materials, which have different magnetic features: ferrites, FeSi sheets and a soft magnetic composite material. Then, the implementation of these models in a finite element code is presented. As example of application, a high-frequency transformer supplied by a rectangular voltage is studied. r 2002 Elsevier Science B.V. All rights reserved. PACS: 75.50.Bb; 75.50.Cc; 75.60.d; 75.60.Ej Keywords: Hysteresis curve; Ferromagnetic material; Preisach model; Jiles–Atherton model; FEM

1. Introduction To design and to study electromagnetic systems, field calculation codes are more and more used. We have at our disposal a virtual prototype in which geometrical dimensions and material characteristics can easily be modified. But, to achieve an efficient design procedure, the used models have *Corresponding author: Tel.: +33-0-320337114; fax: +33-0320436967. E-mail address: [email protected] (F. Piriou).

to be accurate and fast. Satisfying these two criteria simultaneously is not easy, therefore a compromise has to be made. Lots of electrical devices are made up of ferromagnetic materials. To represent the behavior of such materials in field computation codes, nonlinear univoc function are generally used. But, in this case, the hysteresis phenomenon is neglected. In that context, the use of a constitutive relationship which takes into account the hysteresis phenomenon would be more useful to improve the accuracy.

0304-8853/03/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0304-8853(02)01463-4

140

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

The aim of this work is the integration of a simple and few time-consuming hysteresis model in a field computation code based on the finite element method. Different models have been proposed to represent the hysteresis phenomenon and many of them are mathematical models which ignore the underlying physics of the material behavior [1,2]. In our case, we are interested in the use of a model based on physical assumptions. Two macroscopic hysteresis models have been chosen for this study: the Preisach [3,4] and Jiles–Atherton (J–A) [5] models. The Preisach model, based on a phenomenological description of the ferromagnetic materials, is the most used in electrical engineering. The J–A model describes the hysteresis phenomenon as the consequence of a frictional force which opposes to Bloch domain-wall motion. Previous works have shown that both models can be implemented in finite element analysis. Nevertheless, in a 2D field computation code, two formulations are used: the scalar potential formulation (O-formulation) and the vector potential formulation (A-formulation), the most used. To take the material behavior into account, the O-formulation requires a model with the magnetization versus the magnetic field (MðHÞ) whereas the A-formulation requires a model with the magnetization versus the magnetic flux density (MðBÞ). As both original models are presented with H as entry, to implement them in a A-formulation calculation code it is necessary to invert them numerically. For this purpose, dichotomic or Newton numerical schemes are used but it leads to an increase of the computation time. A solution is to use a direct model with B as entry. It is possible to use the J–A model and the Preisach model with B as entry but this is not much studied in the literature. In this paper, we compare the Preisach and J–A models with both H and B as entry in terms of accuracy and time computation. They are applied to model the behavior of three materials. These have been chosen for their use in different areas of electrical engineering: ferrites which are used in power electronics, FeSi sheets for electrical machines at industrial frequencies and iron powder, a soft magnetic composite (SMC), which offers new possibilities for electrical machine design [7]. First

of all, both models are presented with their respective procedures of identification. To evaluate their accuracy, they are tested for different kinds of excitation and using different criterions (hysteresis losses, coercive and remanence fields evolution). They are also compared in terms of computational time. Finally, implementation of the models in a finite element code is presented. As example of application, a high-frequency transformer, made of ferrites, is also studied.

2. Presentation of the models The constitutive relationship of a magnetic material can be described by B ¼ m0 ðH þ MÞ;

ð1Þ

where B is the magnetic flux density, H the magnetic field, M the magnetization and m0 the vacuum permeability. In this section, the Preisach and J–A models are presented. In their original form, they give a relationship between the magnetic field H and the magnetization M: The magnetic flux density B is then obtained from Eq. (1). Nevertheless, models of the magnetization with B as entry, giving the HðBÞ constitutive relationship, can also be deduced from the original form. Both models MðHÞ and MðBÞ are presented in the following. They are restricted to the study of isotropic magnetic materials and in the case of a quasi-static behavior. Dynamic effects (eddy current, dynamic of domain-wall motion, etc.) are not taken into account in the original form of these models. 2.1. The Preisach model and the Everett function 2.1.1. The Preisach model The Preisach model associates to a ferromagnetic material a set of bistable units defined by the function ga;b ¼ 71 (Fig. 1(a)). The switching field couple ða; bÞ characterizing a bistable unit must respect some conditions. If Hsat represents the saturation magnetic field and Msat the corresponding magnetization of the ferromagnetic material, when H > Hsat ; all bistable units are positive and the magnetization is M ¼ Msat : On the opposite

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

141

Fig. 1. (a) Bistable unit, (b) triangle D; and (c) triangle Tðx; yÞ:

side, if Ho  Hsat ; all bistable units are negative and M ¼ Msat : Both previous assumptions lead to the following conditions for the couple ða; bÞ [4]: apHsat bX  Hsat

extrema of the excitation magnetic field Hi ; i.e. some return points of the magnetic field. The memory vector coordinates must verify the following conditions [10]: H0 ¼ 0; h ¼ fH0 ; H1 ; H2 ; H3 ; y; Hn g

As the hysteresis phenomenon is energetically dissipative, we must also have aXb: These three conditions allow us to define a triangle D (Fig. 1), called the Preisach plane. Each couple ða; bÞ characterizing a bistable unit must belong to this plane. A ferromagnetic material is then determined by a statistic distribution pða; bÞ of the switching field couples ða; bÞ belonging to the triangle D: The total magnetization is then given by Z Z M ¼ Msat pða; bÞga;b da db ð2Þ D

The demagnetized state is represented by the equation b ¼ a and the Preisach plane is splitted into two equal surfaces Sþ and S  : Sþ is the surface of couples (a; b) which are such that ga;b ¼ þ1 and S the surface of couples (a; b) such that ga;b ¼ 1: For any other state of the system, the triangle D is splitted into two surfaces Sþ and S separated by a broken line as shown in Fig. 2. Expression (2) can be rearranged as Z Z M ¼ Msat pða; bÞ da db Sþ  Z Z  pða; bÞ da db ð3Þ S

The magnetic state of the system is totally characterized by the broken line. This latter is defined by a memory vector h which includes some

for i ¼ 1; y; n  1 and di ¼ Hi  Hi1 di  diþ1 o0;

ð4Þ

jdiþ1 jojdi j where Hn ; the last component, is the current value of the magnetic field. Using these relations, the memory vector can be easily determined. For example, the magnetic state of Fig. 2 is given by the memory vector h ¼ f0; þHsat ; H1 ; H2 ; H3 g: The knowing of the Preisach density function is sufficient to represent a ferromagnetic material. Several methods for the determination of this function from experimental results are proposed in the literature [8–10]. All these methods require generally numerical derivation and integration, which adds extra numerical errors to experimental ones. So, in contrast of attempts at analytical models, the aim of this work is a fully numerical representation of the Everett function [11] defined by the following equation: Z Z Eðx; yÞ ¼ Msat pða; bÞ da db: ð5Þ Tðx;yÞ

The surface Tðx; yÞ is defined by the right-angled triangle in the Preisach plane (Fig. 1) with ðx; yÞ the vertex coordinates corresponding to the right angle and the hypotenuse is supported by the straight line a ¼ b: The two other legs of the triangle are parallel to a- and b-axis, respectively. Then, if the Everett function is known we can

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

142

Fig. 2. Example of a magnetization process and the corresponding surfaces Sþ and S  :

Fig. 3. Descending part of a centered hysteresis curve and the corresponding Everett function (MðHÞ model).

calculate the magnetization M: In fact, we can show that if H > Hm ; if HoHm ;

MðHÞ ¼ MðHm Þ þ 2EðH; Hm Þ MðHÞ ¼ MðHm Þ  2EðHm ; HÞ

ð6Þ

where Hm is the last return point of the magnetic field (the next to last value of vector h). Then, the Everett function and the magnetization M are linked by a relation which requires no numerical derivation or integration. Experimental determination of the Everett function is presented in the next section. 2.1.2. Identification of the Everett function In the following, we propose a method for the Everett function identification that needs a set of experimental centered minor hysteresis loops. Using these data of MðHÞ; the function EðHm ; HÞ is determined for values of H belonging to ½Hm ; Hm : Fig. 3 gives the descending part of a centered hysteresis loop and the corresponding

function EðHm ; HÞ obtained from Eq. (6) for Hm ¼ 530 A m1 : From n measured centered loops ði ¼ 1; nÞ; a curve set EðHmi ; HÞ supporting the Everett function is obtained (Appendix A). Now, we have to determine this function for an arbitrary point 0 ðHm ; H 0 Þ of the triangle D: This is done by an interpolation method using the previous curve set. This interpolation method must respect the Everett function continuity (this function is a primitive) on the whole studied domain and then for hysteresis curves. The proposed method, already presented in Ref. [12], satisfies this condition. It is based on shape functions used to interpolate the field in 2D finite element method [13]. The expressions of these shape functions are detailed in Appendix A. The Preisach model can also be adapted to obtain a model with B as entry [14]. In this case, the method used for the Everett function identification is detailed in Appendix B. We can note that

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

for each model MðHÞ or MðBÞ an Everett function EðHm ; HÞ or EðBm ; BÞ must be determined. 2.2. The J–A model 2.2.1. The model equation The original J–A model presented in Ref. [5] gives the magnetization M versus the external magnetic field H: This model is based on the magnetic material response without hysteresis losses. This is the anhysteretic behavior which Man ðHÞ curve can be described with a modified Langevin equation:      He a Man ðHÞ ¼ Msat coth ; ð7Þ  He a where He ¼ H þ aM is the effective field experienced by the domains: H is the external applied field and a the mean field parameter representing inter-domain coupling. The constant a is an increasing function of the temperature. The anhysteretic magnetization represents the effects of moment rotation within domains but does not take into account losses induced by domain wall motions. Then, by considering rigid and planar domain walls, the energy dissipated through pinning sites during a domain wall displacement is calculated [5]. The expression of the magnetization energy is obtained under the assumption of a uniform distribution of pinning sites. The magnetization energy is assumed to be the difference between the energy which would be obtained in the anhysteretic case minus the energy due to the losses induced by domain wall motions. Consequently, after some algebraic operations, the differential susceptibility of the irreversible magnetization Mirr can be written as dMirr ðMan  Mirr Þ ¼ kd dHe

ð8Þ

where the constant k is linked to the average pinning site energy. The parameter d takes the value þ1 when dH=dt > 0 and 1 when dH=dto0 with respect to the force which opposes variations of magnetization. However, during the magnetization process, domain walls do not only jump from one pinning site to another: they are flexible and bend when being held on pinning sites. Domain

143

wall bending is associated to reversible changes in the magnetization process. Then, by some physical energy assumptions on the domain wall bending, the obtained reversible magnetization is linearly dependent on Man  Mirr [5]: Mrev ¼ cðMan  Mirr Þ;

ð9Þ

where the reversibility coefficient c belongs to the interval ½0; 1: Assuming that the total magnetization is the sum of the reversible and irreversible components, we have the following expression: M ¼ Mrev þ Mirr

ð10Þ

with Mirr and Mrev defined by Eqs. (8) and (9). Using Eqs. (10) and (9) we can write M ¼ Mirr þ cðMan  Mirr Þ:

ð11Þ

Then, by differentiating this equation with respect to H; the total differential susceptibility of the system is given by the following expression which has already been presented in Ref. [14]: dM ð1  cÞ dMirr =dHe þ c dMan =dHe ¼ : dH 1  ac dMan =dHe  að1  cÞ dMirr =dHe ð12Þ This is the model differential equation that gives the magnetization as a function of the magnetic field H; where Man is given by Eq. (7). The model can also be adapted with B as entry [15]. As for the previous model and using the fact that Be ¼ m0 He ; Eq. (11) is differentiated with respect to B: dM ¼ dB ð1  cÞ dMirr =dBe þ c dMan =dBe : 1 þ m0 ð1  cÞð1  aÞ dMirr dBe þ m0 cð1  aÞ dMan =dBe

ð13Þ In both cases, five parameters a; a; k; c and Ms have to be determined from experimental results. It is important to notice that the J–A parameters are theoretically the same whether B or H is the entry of the model. At the opposite, the Preisach model requires the determination of two independent functions EðHm ; HÞ and EðBm ; BÞ: The physical properties of the five parameters are presented in Table 1.

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

144

Table 1 Physical properties of model parameters a; a; k; c; and Ms Parameter

Physical property

a a k c Msat

Linked to domain interaction Shape parameter for Man Linked to hysteresis losses Reversibility coefficient Saturation magnetization

2.2.2. The parameter identification An identification procedure from experimental data has been described in Ref. [6]. After mathematical developments using Eqs. (7), (8) and (11) for some points of the hysteresis loop, the implicit expressions of the five parameters are obtained in function of experimental data which are: *

* *

*

Hc and Mr the coercive magnetic field and the remanent magnetization, Msat the saturation magnetization, w0ini and w0an the normal and the anhysteretic differential susceptibilities, w0c and w0r the coercive and remanent differential susceptibilities.

Determination of parameters c; a; k and a requires an iterative procedure presented in Ref. [6]. This method is numerically sensitive and does not systematically converge, so we use a slightly different procedure. First of all, we define an objective function Fobj which evaluates the gap between experimental and theoretical results. Z T Fobj ¼ ½Mexp ðtÞ  MJ2A ðtÞ2 dt ð14Þ 0

where Mexp ðtÞ is the experimental magnetization and MJ2A ðtÞ the calculated one, obtained with the same periodic excitation field HðtÞ with a period T: The data ðHðtÞ; Mexp ðtÞÞ corresponds to the major loop, i.e. the centered loop with the highest magnitude. In the first step, we use the following procedure to calculate a first set of the five parameters: (1) systematic choice of ai in the interval [amin ; amax ] by steps Da;

(2) calculation of ai ; ci ; ki from experimental data (Hc ; Mr ; Msat ; w0ini ; etc.), i (3) calculation of the objective function Fobj ; (4) back to (1) until ai ¼ amax ; j (5) determination of (aj ; aj ; cj ; kj Þ using Fobj ¼ i minðFobj Þ: The choice of the interval [amin ; amax ] depends on the materials. But, it can be chosen sufficiently large because this procedure is relatively fast. This first calculation step gives a good estimation of the parameter values but these can be still improved. In this second step, we use an optimization procedure which minimizes the objective function Fobj by modifying independently the five parameters without any constraint. The minimization is carried out considering a set of centered hysteresis loops [16] and not only the major hysteresis loop as for the first step of the identification. Then, the obtained parameters fit well hysteresis loops for small and large magnitudes of excitation. We can note that numerical tests have shown that the second step does not lead to good results without the first step. In fact, this latter enables us to have a set of parameters close to the best solution (in the sense of the chosen objective function) which makes easier the convergence of the optimization procedure. 3. Experimental bench The parameter identification and the comparison of both models can be achieved with experimental data. An experimental bench for magnetic materials characterization is used. This one is described on the synoptic scheme (Fig. 4). The sample is a torus with primary and secondary windings. Excitation is applied using an arbitrary function generator which allows us to impose current or voltage on the primary winding. A computer remotely controls these devices. Using the Ampe" re and Faraday laws, the field H and the magnetic flux density B are calculated from the measurement of the current in the primary winding and the secondary winding voltage. The aim of this study is to compare both models with materials used in different areas of electrical engineering.

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Fig. 4. Synoptic scheme of the experimental bench.

These materials are : * * *

N30 ferrites, used in power electronics; FeSi sheets, used in electrical machines; SMC materials (iron powder) which focuses attention for the new possibilities of electrical machines design that they offer [7].

As the presented models are aimed to be used in the quasi-static case, the measurements were done at the frequency f ¼ 0:5 Hz: At this frequency, we can consider that the three materials are in the quasi-static conditions. In fact, dynamical effects become non-negligible above 100 kHz for the N30 ferrites, 5 Hz for the FeSi sheets (lamination thickness is 0.5 mm) and 1 kHz for the considered SMC material.

cannot be fulfilled and especially for loops with high magnetic field magnitude. Consequently, it is necessary to adjust numerically the measured data in order to center the loops. The more hysteresis loops we use, the more accurate the simulation results. In our case, about 20 centered hysteresis loops have been measured for the three materials. For low excitation fields, it is recommended to have the most important hysteresis loops concentration. In fact, in this area the shape of the hysteresis loop changes a lot (it corresponds to a concavity modification of the first magnetization curve). The Everett functions of the three materials, corresponding to the MðHÞ model, are presented on Figs. 5–7. Their shapes are quite different: the smoother surface is obtained with the SMC material. Indeed, it is characteristic of the hysteresis loop smoothness around the coercive point (Fig. 10). Consequently, it would be difficult to find an analytical expression which enables us to model accurately all these Everett functions. This is one of the advantage of the interpolation method proposed in Section 2.1.2 that does not depend on the studied material. In order to verify the interpolation method robustness, we simulate the same hysteresis loops used for the identification of the Everett function. Simulation results and measured loops are the same. For other magnitudes of excitation, experimental and theoretical results are in good agreements. The gaps between simulation and measurements for the MðHÞ model are shown on Figs. 8–10. These results are similar to those obtained with the MðBÞ model.

4. Comparison of the models 4.1. Identification procedure 4.1.1. Everett function As it was described in Section 2.1.2, the identification of the Everett function needs the measurement of several centered hysteresis loops. The experimental measurement procedure has to be very accurate, especially for high excitation fields where the hysteresis loops have to be rigorously centered otherwise the interpolation method will fail. Experimentally, this condition

145

Fig. 5. Everett function for N30 ferrites.

146

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

4.1.2. Parameters of the J–A model The optimization procedure presented in Section 2.2.2 is applied for the three materials. The parameters are presented in Table 2. The gaps between the calculated major hysteresis loops and

Fig. 6. Everett function for FeSi sheets.

the measurements are shown in Figs. 8–10. One can note that experimental results are well fitted using the optimized parameters. Comparison of the parameters between two different materials is not easy (Table 2). In fact, despite the physical signification of the model parameters, an optimization procedure can lead to good results for the hysteresis loops restitution with parameter values disconnected from any physical signification. Nevertheless, considering the case of N30 ferrites and FeSi sheets, we can note that the parameters lie within the same order of magnitude. Comparison of their parameters can give a good first approach for their physical properties. The parameter a is greater for the FeSi sheets, it is linked to the fact that hysteresis loops for this material have a higher slope at the coercive field. And parameter k shows that hysteresis losses are more important for the FeSi sheets. In the following, we aim at comparing both models for different kinds of excitation wave shapes.

4.2. Sinusoidal excitation

Fig. 7. Everett function for the SMC material.

4.2.1. Hysteresis losses Usually, an hysteresis model is often used to quantify the hysteresis losses, so this criterion was chosen for the comparison. Figs. 11–13 show the evolution of hysteresis losses versus Hmax ; the magnitude of the sinusoidal excitation field. First of all, we can note that the Preisach model is more accurate than the J–A model. The J–A model gives

Fig. 8. N30 ferrites; measured hysteresis loop and gap between both models and measurements for the decreasing part of the loop.

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

147

Fig. 9. FeSi sheets; measured hysteresis loop and gap between both models and measurements for the decreasing part of the loop.

Fig. 10. SMC material; measured hysteresis loop and gap between both models and measurements for the decreasing part of the loop. Table 2 Parameters values for J–A model Parameter a a k c Ms ðA=mÞ

N30 ferrites 5

9:8 10 20 56 0:9 282 100

FeSi sheets 4

1:3 10 59 99 0:55 1 145 500

SMC material 1:8 103 1642 1865 0:8 1 122 600

good results for the three materials even if hysteresis losses are over-estimated for the ferrites and under-estimated for the FeSi sheets. For the SMC material, there is no noticeable difference between both models. In this case, the J–A model gives hysteresis loops in good agreements with experience although the coercive and remanence points evolutions are under-estimated as it will be shown further in Section 4.2.2. This can be

explained by the observation of the hysteresis loop (Fig. 14). Indeed, the over-estimated losses at the hysteresis loop bends (P1 in Fig. 14) are compensated by the under-estimated losses at low field (P2 in Fig. 14). A summary of these results, the maximum relative errors for both models, is presented in Table 3. For the J–A model, the MðHÞ and the MðBÞ models are equivalent, so only one result is given. It can be noted that this latter is well adapted for the SMC material study. At the opposite, both Preisach models (MðHÞ and the MðBÞ) give different results which are nevertheless close. Globally, the Preisach model is the most accurate. 4.2.2. Remanent and coercive point evolutions Another criterion for testing the model robustness is the evolution of two characteristic points of

148

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Fig. 11. N30 ferrites; comparison of hysteresis losses for both models with measurement at 0.5 Hz.

Fig. 12. FeSi sheets; comparison of hysteresis losses for both models with measurement at 0.5 Hz.

Fig. 13. SMC material; comparison of hysteresis losses for both models with measurement at 0.5 Hz.

a hysteresis loop: the remanence and the coercive points. Remanent and coercive points evolution versus Hmax are presented in Figs. 15–20. Simula-

tions for ferrites and FeSi sheets give results close to the measurements. An important aspect of the J–A model is the behavior of the coercive point

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

which becomes constant for low excitation fields whereas this asymptotical evolution appears experimentally at higher excitation fields. As it will be explained in Section 4.4, this behavior comes from the numerical solution of the differential Eq. (12).

Fig. 14. Over- and under-estimated losses for a major hysteresis loop in the SMC material case.

Table 3 Hysteresis losses maximum relative errors (in %) for both models Material

J–A model

Preisach M(H) model

Preisach M(B) model

N30 ferrites FeSi sheets SMC material

+12

4.5

+6.6

+3.7 2

+3.7 +8

9 +2.5

149

4.3. Non-sinusoidal excitation In electrical engineering, systems are not necessarily studied under sinusoidal excitations (pulse width modulation supply or during transient state). To test the capabilities of models to represent the magnetic material behavior under such kinds of excitation, we now consider two types of excitations: Excitation A : HðtÞ ¼ H0 sinðotÞexpðt=tÞ;

ð15Þ

Excitation B : HðtÞ ¼ H0 sinðotÞ þ H00 sinð5ot þ fÞ;

ð16Þ

with t being a decreasing time constant, H0 the magnetic field magnitude, H00 the fifth harmonic magnitude and f its phase lag. The demagnetization-like excitation, denoted excitation A, (Eq. (15)) is often met in the transient state of a system supplied by sinusoidal voltage. However, it allows us to sweep through a great number of magnetic states of the material. In the other hand, as electrical devices are more and more supplied by power electronic converters which induce periodic but non-sinusoidal voltage and current, the second kind of excitation, denoted excitation B, (Eq. (16)) enables us to test the accuracy of both models in such conditions. Results for both excitations are shown in Figs. 21–26. In the case of excitation A, both models are in good agreements with the experimental measurements. These results could be predicted because

Fig. 15. N30 ferrites; comparison with experience of remanent point evolution for both models.

150

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Fig. 16. N30 ferrites; comparison with experience of coercive point evolution for both models.

Fig. 17. FeSi sheets; comparison with experience of remanent point evolution for both models.

Fig. 18. FeSi sheets; comparison with experience of coercive point evolution for both models.

this excitation can be interpreted as a sequence of centered hysteresis loops with a decreasing magnitude. In fact, we have seen in the previous section that both models are accurate in the case of centered loops. Only, a shift between measure-

ments and simulations appears at low excitation fields for FeSi sheets. When applying excitation B, important difference between the two models appears. We can see that the Preisach model is relatively accurate

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

151

Fig. 19. SMC material; comparison with experience of remanent point evolution for both models.

Fig. 20. SMC material; comparison with experience of coercive point evolution for both models.

Fig. 21. N30 ferrites; comparison of theoretical and experimental results using excitation A.

whereas the J–A model gives an incoherent behavior. This behavior is more noticeable in the case of FeSi sheets and SMC material. Figs. 25

Fig. 22. FeSi sheets; comparison of theoretical and experimental results using excitation A.

and 26 include a zoom of the field return-points region where minor loops are not closed. Indeed, this model does not systematically ensure the

152

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Fig. 23. SMC material; comparison of theoretical and experimental results using excitation A.

Fig. 24. N30 ferrites; comparison of theoretical and experimental results using excitation B.

Fig. 25. FeSi sheets; comparison of theoretical and experimental results using excitation B.

Fig. 26. SMC material; comparison of theoretical and experimental results using excitation B.

underlined in numerous papers [17]. Some solutions have been proposed to avoid this problem, but it makes the model more complicated and requires to know a priori the magnitude of the loop. When simulating an electrical device with time-stepping finite elements analysis and using an hysteresis model, the magnitude of each loop on each element of the mesh is not initially known. In this case, the modified J–A model is not adapted. The Preisach model is more adapted in this case thanks to its intrinsic memory properties. In fact, there is a return-point memory aspect of this model which is clearly emphasized by the broken line in the Preisach plane (cf. Section 2.1.1). Each angle of this broken line is a return-point, i.e. an extremum, of the applied field. For this model, the results are quite correct in the case of the N30 ferrites and SMC material. For FeSi sheets, the simulated minor loops are flatter than the experimental ones. Several methods have been proposed to improve the account for minor loops. For example, E. Della Torre has proposed a method called ‘‘moving Preisach model’’ which insures a better restitution for minor loops [18]. An interaction field, that depends on the magnetization, is introduced as a feed-back to the initial model. But, this method increases the complexity of the model and the calculation time. 4.4. Computational effort

enclosure of minor loops. Moreover, the minor loops are very flat in comparison with experimental ones. The phenomenon has already been

To use the J–A model, we have to solve a firstorder differential equation which can be done by

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

the Euler implicit scheme. It has been observed that this model is very sensitive to the chosen magnetic field step dH to solve the equation. The use of a large dH can lead to an important deformation of hysteresis loops. For example, in the case of a sinusoidal excitation, we can see in Fig. 27 a comparison between two dH giving 200 points and 2000 points for a major hysteresis loop. The surface of the loop increases with the magnetic field step dH: Then, it exists a minimal dH beyond which hysteresis loops become stable. Our experience shows that a time step allowing the calculation of about 2000 points for an hysteresis loop is required to have a good compromise between accuracy and computational time for the three materials. As it was presented previously for the minor loops, the numerical sensitivity of this model does not allow calculation of hysteresis loops in the Rayleigh zone. These are not closed and often not centered. It is necessary to calculate several hysteresis curves in order to obtain a satisfying result. On the other hand, the Preisach model is not sensitive to the time step, in field calculation about only 100 points or less are calculated for this model. In this case, considering the same return points for the excitation field, the time step has no effect on the hysteresis loop. In fact, at the opposite of the J–A model, no differential equation is solved (Eq. (6)). Computational times of both models are reported in Table 4. Simulations are made for the FeSi sheets on a major hysteresis loop with time

153

Table 4 Computational times for FeSi sheets (400 MHz Digital WorkStation) Model

Time (ms) for 2000 points

Time (ms) for 100 points

Preisach (Everett) J–A

250 80

100 n.a.

steps giving 2000 points in both cases and 100 points in the Preisach model case. First of all, we can see that computation time for the Preisach model is not proportional to the number of points. In the beginning of the code execution, the curves set supporting the Everett function must be loaded. The time duration of this process is included in the time computation. The Preisach model is more time consuming which is explained by two important points. First, the memory vector management is a heavy numerical method. Moreover, the interpolation method requires a great number of numerical operations. At the opposite, the J–A model requires a reduced number of operations at each calculation step. Despite the number of 2000 points required for this model, it is slightly less time consuming than the Preisach model with 100 points. This slight difference can give advantage to the J–A model when using the finite elements analysis providing no minor-loops are calculated (to avoid the minor-loops non-enclosure). In fact, the magnetization calculation is done on each element and for each time step as it will be presented in the next section. 5. Implementation in field computation In the following, the 2D implementation of an hysteresis model in a field calculation code, with the coupling with external circuit equations, is presented. We consider a domain D bounded by a surface S in the case of magnetostatics. The equations to be solved are: div B ¼ 0;

Fig. 27. Two loops calculated with 200 points and 2000 points per period with the J–A model (N30 ferrites).

curl H ¼ J;

154

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

becomes

with n  B ¼ 0 on Sb ;

curl nFP curl A ¼ J  curl MFP and n 4 H ¼ 0 on Sh ;

ð17Þ

with Sb and Sh being two complementary parts of S; J the current density and n the outward normal vector of S: To take into account the material behavior, the constitutive relationship H ¼ fðBÞ; is added. To model ferromagnetic material, this relationship can be one of the two models previously presented. Nevertheless, these models are scalar models whereas we need a vectorial model in Eqs. (17). Then, to evolve this latter, we assume that the material is isotropic and that B and H are collinear. The magnitude of H is calculated from the one of B and the direction of H is the same as the one of B: All these equations are generally solved using a potential formulation. In the 2D case, the vector potential formulation is generally preferred to the scalar potential one, then we have curl ½fðcurl AÞ ¼ J; with n 4 A ¼ 0 onSb ; and n 4 H ¼ 0 on Sh ;

ð18Þ

where A represents the magnetic vector potential (i.e. B ¼ curl A). One can note that for this formulation a constitutive relationship with B as entry is required, but if a scalar potential formulation (H ¼ grad j) is chosen, a model with H as entry is necessary. The numerical solution of Eq. (18) including hysteresis cannot be done with the same method as the one used with univoc functions (Newton–Raphson scheme for example). We have chosen the fixed-point method already presented in Ref. [19]. The hysteretic constitutive relationship is then rewritten under the form H ¼ fðBÞ ¼ nFP B þ MFP ðBÞ:

ð19Þ

The reluctivity nFP is a constant and must respect some conditions to achieve convergence [20]. The studied hysteretic models assume B and H collinear, consequently the magnetization MFP has the same direction as nFP B: Its magnitude is obtained by calculating MPF ¼ f ðBÞ  nFP B: Finally, the partial differential equation (18)

ð20Þ

The discretization with nodal shape functions for the potential vector of Eq. (20) using the Galerkin method leads to the matrix system ½SFP ½A ¼ ½J  ½MFP ;

ð21Þ

where the vector ½A represents the nodal values of vector potential, ½SFP  a square matrix called stiffness matrix, ½MFP  and ½J the vectors which take into account the magnetization MFP and the current density J. One can note that the matrix ½SFP  is constant because the permeability nFP is constant as well. The non-linearities introduced by ferromagnetic materials are reported in the source term ½MFP  which depends on B (i.e. A). To take into account the coupling with the external circuit of a coil made up of stranded conductors flowed by a current i; a vector ½D is introduced such that ½J ¼ ½Di [21]. Then, we obtain the system " #" # " # " # SFP D A 0 0 d A þ 0 R i Dt 0 dt i " # " # 0 MFP ¼ þ : u 0 This system can be time discretized using an Euler implicit scheme.

6. Example of application To compare the performances of both hysteresis models implemented in finite element analysis, a power electronic transformer, which core is composed of N30 ferrites, has been studied. Conditions of simulation are such that the primary winding is supplied by a square wave voltage V ¼ 400 V at a frequency of 80 kHz: In this case, as the core is made of ferrites, dynamic effects can be neglected. The secondary winding is in series with a resistance R ¼ 10 O: The 2D geometry of the system is given in Fig. 28. Due to the symmetry properties, only a quarter of the structure is studied. The structure mesh is chosen relatively fine with 1801 elements (Fig. 29). We aim at

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

155

comparing local and global quantities obtained from both models. 6.1. Calculation of local and global quantities First, as local quantity, we can compare the hysteresis loops obtained from both models for a given element of the mesh. Both models give close results for different elements of the mesh. As example, Fig. 30 gives the loops obtained on the element located at the point P (see Fig. 28). The global quantities to be compared are the primary and the secondary winding currents. As both models give very close results, only wave shapes of currents and, for more convenience, the difference DI between the currents, given by both models, are reported in Figs. 31 and 32. As we can see, difference between both models is negligible: the maximum gap is about 0.9% of the amplitude for the primary current and about 0.3% for the secondary current. The second global

Fig. 30. Hysteresis loops obtained from both models for an element of the mesh.

quantity to be compared is the evaluation of hysteresis losses. These are very important in electrical engineering for machines design. Hysteresis losses for both models have been determined from the evolution of the density of losses dh versus Bmax : X Se dh ðBemax Þ ð22Þ Ph ¼ l elements

Fig. 28. Geometry of the transformer.

where the sum is over all the elements of the mesh. Se is the surface of the element e; l the depth of the system (in the third dimension) and Bemax the maximum magnitude of the magnetic flux density on the element e. Results of the calculation are reported in Table 5. From these results, current and hysteresis losses evaluation, we can conclude that, in this case, both models give close results for global quantities calculation. Now, it is important to choose the model with the best performances in terms of calculation time. 6.2. Calculation times

Fig. 29. Mesh of the studied system: 1 quarter of the structure.

In Table 6, the ratio of computational times are reported for both models. The J–A model computation time is taken as reference. Calculation have been made for 1 period and 3 periods. In fact, the Preisach model first numerical step is the loading of the experimental Everett functions which are used for the interpolation method. It appears that, after some periods the ratio

156

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Fig. 31. Primary winding current wave shape and difference between the models.

Fig. 32. Secondary winding current wave shape and difference between the models.

Table 5 Hysteresis losses calculated from both models

Hysteresis losses (W)

J–A model

Preisach model

5.69

5.94

Table 6 Ratio of computational times for 1 and 3 periods field calculation (400 MHz Digital WorkStation)

1 period 3 periods

between the models calculation times stabilizes at 2.25 in favor for the J–A model. Moreover, the convergence of non-linear FEM procedure is quite the same. These results show how important is the choice of the user for the ratio (accuracy/calculation time) for the use of one of these models. Of course, this choice has to be made in the case of a sinusoidal excitation. In other cases, the Preisach model is recommended.

J–A model

Preisach model

1 1

2.5 2.25

(J–A model computational time is taken as reference.)

7. Conclusion Main results of the comparison are summarized in Table 7. The Preisach model has been presented with an identification method based on the Everett function determination. This method needs the measurement of rigorously centered hysteresis

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

157

Table 7 Summary of the results Excitation

J–A model

Identification effort

Preisach model (Everett)

Two steps: iterative and optimization procedures for five parameters. Set of experimental centered hysteresis loops. N30 ferrites Sinusoidal ++ Non-sinusoidal + FeSi sheets Sinusoidal + Non-sinusoidal + SMC material Sinusoidal ++ Non-sinusoidal + Computational effort  Require 2000 points per loop but fast  Easy to be numerically implemented

Interpolation method for MðHÞ and MðBÞ models. Experimental loops have to be rigorously centered for the MðBÞ model. ++ ++ ++ + ++ ++  100 points are sufficient but more time consuming  Difficult to implement

++ good results + correct results

loops otherwise gaps between experience and simulation can be observed for the remanent and coercive points, especially in the case of the MðBÞ model. Nevertheless, results for the Preisach model are more accurate whatever the kind of excitation used in this paper. At the opposite, although the J– A model has more physical basis, to obtain good results on a wide range of magnitudes, the parameters are calculated using an optimization procedure based on several centered loops. In terms of identification effort, both models are equivalent. Magnetic materials used for this study were chosen for their different areas of application and for their different hysteretic behavior. Then, results give us a good idea of which model to use in relation to the studied electromagnetic system (Table 7). For example, in the case of a sinusoidal excitation, the use of the J–A model gives results close to those of the Preisach model. Nevertheless, it has been observed that this model is less accurate than the Preisach model for the representation of the N30 ferrites and FeSi sheets. In fact, the Preisach model is well suited for the study of the three materials. But, results show that, in the case of a sinusoidal excitation with a SMC material, the use of the J–A model is recommended. The J–A model main drawback is the non-closure of minor loops when applying excitations with superposi-

tions of higher harmonics. In this case, the Preisach model should be preferred. Concerning the computational aspect, the J–A model is simple to be implemented in a field calculation code and it requires less computational time and memory. In fact, the Preisach model requires the memory vector to be stored in addition to the current magnetic field. When studying a complex electromagnetic system, especially in finite element analysis, the magnetic constitutive relationship calculation must be few time consuming. The example of application with the transformer show how it is important as the Preisach model requires more than 2 times computational time for similar results to those of the J–A model. Nevertheless, this point has to be compensated by the fact that the Preisach model is more accurate than the J–A model in case of nonsinusoidal excitation.

Appendix A. Identification of the Everett function for the M(H) model of Preisach In the following, we present the method to determine the Everett function from n centered hysteresis loops whose magnitudes are between 0 and Hmax : It is then possible to calculate EðHm ; HÞ for ðHm ; HÞA½0; Hmax  ½0; Hmax  with HoHm :

158

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

Only the model with H as entry is detailed. Nevertheless, it is possible to have a similar approach to obtain the Everett function with the magnetic flux density B as entry. In this case, we obtain the magnetization versus the magnetic flux density. If we consider a centered hysteresis loop of magnitude Hm ; from relation (6), the magnetization M along the descending part of the hysteresis loop is: MðHÞ ¼ MðHm Þ  2EðHm ; HÞ:

ðA:1Þ Fig. 33. Curves set supporting the Everett function.

From n measured hysteresis loops ði ¼ 1; nÞ with a magnitude Hmi ; the function EðHmi ; HÞ is calculated for values of H belonging to ½Hmi ; Hmi  and for each loop i: We obtain the set of functions EðHmi ; HÞ supporting the 3D Everett function (Fig. 33) with H between Hmi and Hmi : We can notice that for the particular case where H ¼ Hmi ; the Everett function is equal to zero. In this case, the triangle Tðx; yÞ is a single point. Now, we have to calculate the Everett function for 0 any point ðHm ; H 0 Þ belonging to the triangle D: The proposed method is based on functions used for the magnetic field interpolation in the bidimensional finite element analysis method [13]. Using the symmetry property of the Everett function EðHm ; HÞ ¼ EðH; Hm Þ; it is possible to treat 0 only the case where the value of H is between Hm 0 and þHm : We consider now the values Hmi1 and Hmi which correspond to the magnitudes of the known curves set (Fig. 33) such as 0 Hmi1 oHm oHmi (Fig. 34). In the following, for more convenience, we suppose that i ¼ 2 and the searched value will be interpolated from the two curves with the lowest magnitudes values Hm1 and Hm2 : We have then to distinguish two cases: *

First case: jH 0 jojHm1 j: First, we search for the fields H11 and H12 whose Everett functions values EðHm1 ; H11 Þ and EðHm1 ; H12 Þ are known (i.e. determined from experimental measurements). These fields verify H11 oH 0 oH12 : In the same way, we determine H21 and H22 on the straight line Hm ¼ Hm2 such as H21 oH 0 oH22 : Then, four points a11 ; a12 ; a21 ; a22 are obtained (Fig. 35) which define a quadrilateral with the

Fig. 34. Projection in the D triangle of the Everett functions.

Fig. 35. Interpolation method M1 used to calculate the Everett function in the first case. 0 ; H 0 Þ inside. The value of the Everett point ðHm 0 function EðHm ; H 0 Þ corresponding to this point is obtained by a quadratic interpolation method from the value of the Everett function at points

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

with a ¼ 1 if H > 0; else a ¼ 1;

a11 ; a12 ; a21 and a22 : 0 ðHm  Hm1 Þ 0 ½EðHm2 ; H 0 Þ ; H0Þ ¼ EðHm ðHm2  Hm1 Þ  EðHm1 ; H 0 Þ þ EðHm1 ; H 0 Þ

EðHm2 ; Hc Þ ¼ ðA:2Þ

with EðHmi ; H 0 Þ ¼

ðH 0  Hi1 Þ ½EðHmi ; Hi2 Þ ðHi2  Hi1 Þ  EðHmi ; Hi1 Þ þ EðHmi ; Hi1 Þ ðA:3Þ

*

for i ¼ 1; 2: Second case : jH 0 j > jHm1 j: The previous procedure has to be slightly modified in the case where H 0 is greater than Hm1 : In this case, we consider the straight line crossing the points 0 ðHm1 ; Hm1 Þ and ðHm ; H 0 Þ: The intersection point with the vertical line Hm ¼ Hm2 gives us a point ðHm2 ; Hc Þ (Fig. 36).

Finally, we search for the points H21 and H22 such as H21 oHc oH22 and the Everett function value is then obtained by 0 EðHm ; H 0Þ ¼

159

0 ðHm  Hm1 Þ ½EðHm2 ; Hc Þ ðHm2  Hm1 Þ  EðHm1 ; aHm1 Þ þ EðHm1 ; aHm1 Þ

ðA:4Þ

Fig. 36. Interpolation method M1 used to calculate the Everett function in the second case.

ðHc  H21 Þ ½EðHm2 ; H22 Þ ðH22  H21 Þ  EðHm2 ; H21 Þ þ EðHm2 ; H21 Þ ðA:5Þ

The presented interpolation method, denominated as M1, allow us to insure the continuity of the Everett function on the domain. This characteristic is necessary to obtain continuous MðHÞ loops. Moreover, we observe an exact restitution of the experimental loops used for the interpolation method. Appendix B. Identification of the Everett function for the M(B) model of Preisach The Preisach model can also be adapted with B as entry [14]. In this case, the shape of the Everett function which has to be determined is totally different from the one previously presented in the MðHÞ case. To calculate the Everett function from the MðBÞ model, we can still use the same method as presented for the MðHÞ model. Nevertheless, the interpolation method has to be modified in this case. In fact, the use of the interpolation method M1 gives deformed hysteresis loops at low magnitudes of excitation. Then, we have chosen another interpolation method which is similar to the previous one, and denominated as M2 (Fig. 37). As for the M1 method, the Everett

Fig. 37. Presentation of the interpolation method M2.

160

A. Benabou et al. / Journal of Magnetism and Magnetic Materials 261 (2003) 139–160

0 function is determined for the point ðHm ; H 0 Þ and we have the same definition for the fields Hm1 and Hm2 : Using the straight line d crossing the origin 0 and the point ðHm ; H 0 Þ; we determine the intersections b1 ðH1 ; Hm1 Þ and b2 ðH2 ; Hm2 Þ with the vertical lines H ¼ Hm1 and H ¼ Hm2 : Then, points a11 ; a12 ; a21 and a22 are determined such as presented in Fig. 37. As for the method presented in Appendix A, these four points are used to calculate the 0 Everett function EðHm ; H 0 Þ: The used interpolation function is then given by the Eq. (A.2). It must be noted that, at the opposite to the method M1, method M2 needs the treatment of only one case.

References [1] L.O. Chua, K.A. Stromsmoe, IEEE Trans. Circuit Theory CT-17 (4) (1970) 564. [2] M.L. Hodgdon, IEEE Trans. Magn. 24 (6) (1988) 3120. [3] F. Preisach, Z. Phy. 94 (1935) 277. [4] I.D. Mayergoyz, Mathematical Models of Hysteresis, Springer, New York, 1991. [5] D.C. Jiles, D.L. Atherton, Theory of ferromagnetic hysteresis, J. Magn. Magn. Mater. 61 (1986) 48.

[6] D.C. Jiles, J.B. Thoelke, M.K. Devine, IEEE Trans. Magn. 28 (1992) 27. [7] J. Cros, P. Viarouge, ICEM’02 Bruges, Belgium 2002, CDROM, paper No. 514. [8] H. Debruyne, S. Cl!enet, F. Piriou, Math. Comput. Simulation 46 (1998) 301. [9] G. Bertotti, V. Basso, J. Appl. Phys. 73 (10) (1993) 5827. [10] G. Biorci, D. Pescetti, Il Nuovo Cimento 7 (6) (1958) 829. [11] D. Everett, Trans. Faraday Soc. 51 (1955) 1551. [12] S. Cl!enet, F. Piriou, MGE 2000, France, Lille, December 13–14, 2000, pp. 71–74. [13] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, McGraw-Hill, London, New York, 1989. [14] G.S. Park, S.Y. Hahn, K.S. Lee, H.K. Jung, IEEE Trans. Magn. 29 (2) (1993) 1542. [15] N. Sadowski, N.J. Batistela, J.P.A. Bastos, M. LajoieMazenc, Compumag, France, Evian, Vol. 4, July 2–5, 2001, pp. 246–247. [16] S. Cl!enet, J. Cros, F. Piriou, P. Viarouge, L.P. Lefebvre, Compel 20 (1) (2001) 187. [17] P.I. Koltermann, L.A. Righi, J.P.A. Bastos, R. Carlson, N. Sadowski, N.J. Batistela, Physica B 275 (2000) 233. [18] E. Della Torre, IEEE Trans. Magn. 27 (4) (1991) 3697. [19] O. Bottauscio, D. Chiarabaglio, M. Chiampi, M. Repetto, IEEE Trans. Magn. 31 (6) (1995) 3548. [20] V. Ionita, B. Cranganu-Cretu, B. Iona, IEEE Trans. Magn. 32 (3) (1996) 1128. [21] F. Piriou, A. Razek, IEEE Trans. Magn. 28 (1992) 1295.