Three-Dimensional Magnetic Hysteresis Modeling Based on Vector Hysteresis Operator

Based on the minimum energy principle of magnetization stable state, this paper presents a definition of three-dimensional (3-D) hysteresis operator by analyzing the magnetization state of the magnetic particles. The direction determination rule of magnetization is given in detail. The critical surface equation of hysteresis operator is described in the space sphere coordinate system according to the modeling principle of the hybrid vector hysteresis model. The magnetization process of the hysteresis operator is studied under the alternating and rotational magnetic field, respectively. Moreover, a 3-D isotropic hysteresis model is initiatively established, and the hysteresis properties of SMC material are simulated and compared with the experimental results. This study lays a foundation to build the 3-D hybrid vector hysteresis model and simulate the 3-D magnetic properties.


I. INTRODUCTION
Modeling of magnetic properties is a hot topic in the field of electrical engineering. Preisach model, Jiles-Atherton (J-A) model, Stoner-Wohlfarth (S-W) model and Enokizono and Soda (E&S) model et al. are typical modeling methods for magnetic properties study at present. [1]- [8] Each model has its own characteristics and applicable occasions in different ways.
After years of development, the research on these models can be divided into three categories: scalar hysteresis model, vector hysteresis model and hybrid vector hysteresis model. The scalar model is essentially a special case of the vector model. The hybrid hysteresis model combining two or more basic models has attracted the attention of scholars recently. It combines the advantages of the basic models, while avoids disadvantages, which can simulate the hysteresis properties more accurately. The idea of hybrid vector hysteresis model was proposed earlier by G. Friedman and I. D. Mayergoyz. [9] It was derived from a general vector Preisach model by extending the S-W model and Mayergoyz's vector Preisach model. C. Michelakis et al. proposed a hybrid vector hysteresis model for single domain single axis particles by combining the classical Preisach model and S-W model. [10], [11] E. Della Torre and E. Cardelli presented the general vector hysteresis model (GHM, Generalized hysteresis model) and DPC (Della Torre, Pinzaglia and Cardelli) model. The GHM model calculates the irreversible magnetization components based on the closed critical surface of each hysteresis operator. [12]- [14] While the DPC model can analyze the magnetization characteristics of unit single domain particles. [15]- [17] Dingsheng Lin et al. presented an anisotropic vector hysteresis model based on the improved isotropic vector play operators. [18] The author's research group has been keen on the study in the field of hysteresis characteristic simulation. The definition method of vector hysteresis operator was presented and a 2-D vector hysteresis mathematical model was established to simulate the hysteresis characteristic for soft magnetic composites in previous publications. [19]- [21] On the basis of the results of previous research, a definition of hysteresis operator in 3-D case is proposed and expanded in detail in this paper. The magnetization state of magnetic particles is analyzed based on the principle of minimum energy of magnetized stable state. According to the astroid properties, the 3-D hysteresis operators are established for anisotropic and isotropic materials in the spherical coordinate system, respectively. Meanwhile, the influence of the energy parameter on the critical surface of the hysteresis operator is analyzed. And then, the magnetization process of hysteresis operator is studied under alternating and rotational magnetic field. Finally, a 3-D isotropic hysteresis model is preliminarily established based on the 3-D hysteresis operator. The rotational hysteresis properties of SMC material are simulated and compared with the experimental data to verify the availability of the 3-D model.

II. ANALYSIS OF MAGNETIZATION STATE OF MAGNETIC PARTICLES
According to the minimum energy principle of magnetization stable state, the magnetization state of magnetic particles is analyzed. In the case of single domain single axis anisotropic particles, the direction of the magnetization is consistent with the easy axis direction under zero magnetic field. After applying the field, the direction of the magnetization is rotated off the easy axis towards the direction of the applied magnetic field. The vector magnetization M, applied magnetic field Ha and their projections on the xoy plane are shown in Fig. 1.
In Fig. 1, z axis is assumed to be the easy axis, M is the angle between M and z axis, H is the angle between Ha and z axis, and  is the angle between Mxoy or Ha-xoy (the projection of M or Ha on the xoy plane) and the x axis.
The total free energy E of the single domain single axis anisotropic particles with volume V is equal to the sum of uniaxial anisotropic energy and field energy, see (1).
According to the principle of minimum energy, the energy parameter e0 must satisfy (4).
Equation (6) can also be directly derived from Fig. 1.
The relationship among hx, hy and hz can be derived by (5) When the energy is taken to the extreme value, the critical surface can be described as (9), and the corresponding critical surface named 3-D astroid surface is shown in Fig. 2.
where e is the energy parameter.
As shown in Fig. 1 and (10), when M = −/2, 0, /2 and , the surface formed by the equation can not form a closed surface, none of them satisfies the requirements of establishing 3-D hysteresis operator. Therefore, neither sinM nor cosM is 0, which also verifies the assumption.
According to (8) and (10), it is known that all points on the equipotential surface are in a stable state when e < −0.5.
And when the effect of anisotropic constant is considered, (11) can be obtained by using e0=E/2V, ha=0MsHa/2 and dividing ha into three axes to simplify (2). The 3-D astroid surface can also be obtained when the energy is taken to the extreme value, see (12) and Fig. 3.
It is also known that all points on the equipotential surface are in a stable state when e < −(K/2). If the anisotropic constant is not considered, it is the same situation as taking the anisotropic coefficient K to 1, this paper takes K = 1.
In the same way, the magnetization state of isotropic magnetic particlesis is analyzed, too. The total free energy of the single domain single axis isotropic particles can be equal to the field energy, see (14).
where m is the unit vector magnetization and its modulus is 1. The coordinate algorithm of vector dot product is used to transform (14) to (15 By using e0=E/2, ha=0MsHa/2 and dividing ha into three axes, (15) can be simplified. 0 sin cos sin sin cos Samely, the energy parameter e0 must satisfies (4). It can be also assumed that neither sinM nor cosM is 0, so (17) can be be deduced from ∂e0 / ∂M = 0. cos sin 0 sin cos Equation (18) (18) The relationship among hx, hy , and hz can be derived by (17) and (18), see (19).
sin cos cos sin sin cos Equation (20) can be be deduced from Substituting hx, hy and hz into (16) respectively, and the equation of equipotential surface can be got from de0=0, see (21). sin cos sin sin As can be seen from Fig. 1 and (21), when M = −/2, 0, /2 and , the surface formed by the equation is not closed, so they cannot satisfy the requirements of establishing 3-D hysteresis operator. Therefore, neither sinM nor cosM is 0, it also verifies the assumption. According to (20) and (21), all points on the equipotential surface are testified in a stable state when e < 0.

III. DEFINITION OF HYSTERESIS OPERATOR IN 3-D CASE
The definition of hysteresis operator plays a very important role in the vector hysteresis modeling and can be concluded from the previous study. The magnetization direction of the points on the equipotential surface in the stable state satisfies the astroid properties of the S-W model. [17], [19]- [21] The direction of spontaneous magnetization should satisfy the principle of minimum energy in S-W model. The principle of hysteresis operator definition is to make the magnetization direction easy to be determined and the numerical realization of the whole model simple and feasible.
Here, the equation of equipotential surface in stable state is selected as the critical surface of hysteresis operator. The reason why hysteresis operator is defined in this way and the determination of magnetization direction are analyzed as follows. The determination of the magnetization direction is shown in Fig. 4.
As can be seen from Fig. 4, for any point P in H space, a line PB tangents to the astroid curve, and the intersection point of the tangent line and the Hxoy axis is P0(sinM, 0). The tangent line PB intersects with the astroid curve at A(sin 3 M, −cos 3 M). When 0 < M < , it is obvious that PA segment is in a stable state and AB segment is in an unstable state. Similarly, when − < M < 0, PC segment is in an unstable state and CD segment is in a stable state. The stable part of the tangent line is perpendicular to the equipotential curve of the stable state. Since de0 = −0mdha = 0 when the field changes along the isopotential curve, the magnetization direction is perpendicular to the equipotential curve.
It is easier to determine the direction of magnetization by selecting the equipotential surface equation of stable magnetization state of magnetic particles as the critical surface of the operator than by selecting S-W astroid surface.
From what has been discussed above, the equipotential surface equations under magnetized stable state (when e < −0.5 in (10)) are selected as the critical surface formulation for anisotropic hysteresis operators. Similarly, the critical surface of isotropic hysteresis operators is formed by the equations (when e < 0 in the (21)). Hysteresis operator is defined as a closed domain in H space, which is convenient for numerical analysis and has certain physical significance. The relationship between hzmax (the maximum of hz along z axis) with rmax (the maximum radius of cross profile paralleled to xoy plane) is obtained from the analysis of Fig.  5(a). When e takes different values, the relationship of rmax and hzmax satisfy the function hzmax = rmax − 0.5. Whatever e value decreases, the longitudinal cross section of anisotropic hysteresis operator is not circular but will be infinitely close to circle. While in Fig. 5(b), the relationship between hzmax and rmax satisfies the function rmax = hzmax when e takes different values. That is, no matter how the value of e decreases, the longitudinal cross section of the isotropic hysteresis operator is always circular.
Considering the principle of Preisach model, the entry point is the opening threshold, and the exit point is the closing threshold. When the applied magnetic field passes through the equipotential surface from different directions, only when e < −0.5 and passing through z axis is the shortest distance, which satisfies the assumption that z axis is the easy magnetization axis.
As shown in Fig. 4, for any point P in H space, the magnetization direction is perpendicular to the equipotential curve. According to these results and astroid properties in S-W model, one conclusion is obtained that the magnetization direction is depending on the locations of the applied field. When the applied field is outside of the operator's critical surface, the magnetization directionis perpendicular to the critical surface. When the magnetic field is inside the operator, the magnetization direction is fixed with the same direction as the entry point where the applied field enters the operator's critical surface. There will occur a Barkhausen jump when the applied field exits the operator's critical surface, the magnetization direction mutates to the normal direction of the tangent plane of the exit point.
The defined 3-D hysteresis operator is also in accordance with the second principle of thermodynamics, and with the loss property, congruency property and deletion property.

A. Alternating Magnetization
The magnetization process of anisotropic and isotropic hysteresis operators is analyzed under alternating magnetic field, respectively. The change of applied field path needs to reflect the process of entering and exiting the operator, so the setting of the maximum value of the applied field should be larger than the size of the operator. When the applied field entering the operator, the magnetization direction is fixed with the same direction as the entry point. And when the applied field exiting the operator, the magnetization direction mutates to the normal direction of the tangent plane of the exit point. Hence, when the change of the applied field from the positive maximum to the negative maximum along a certain axis and back to the positive maximum corresponds to a period of magnetization, and there will be a corresponding loop between M and H.
The magnetization paths of a single anisotropic hysteresis operator and the corresponding hysteresis loops are shown in Fig. 6 when e = −1. And Fig. 7 shows the same case of the isotropic hysteresis operator. It is assumed that the paths of alternating magnetic field are (1.5, −1.5, −1.5) -(−1.5, 1.5, 1.5) -(1.5, −1.5, −1.5) and the hysteresis loops are given from the direction of x, y and z. The simulation steps of the operator are similar to the hysteresis modeling, which will be explained in section V.
As shown in Fig. 6, anisotropic hysteresis operators have the same hysteresis loop in the x and y directions, while that of the z direction is bigger. This is because the magnetization path inside the operator has the same projection in the x and y directions, while the z direction is slightly bigger. When the applied magnetic field crosses the critical surface of the operator, the magnetization M increases in the x and y directions, while decreases in the z direction. The reason is that as the applied magnetic field increases, the corresponding unit magnetization direction should be perpendicular to the operator's critical surface. So the component of unit magnetization increases in the x and y directions and decreases in the z direction.
As shown in Fig. 7, the hysteresis loops for isotropic operators in the x, y and z directions can be plotted similarly. It is deduced that the magnetization M remains unchanged along with the three directions when the applied magnetic field crosses the critical surface of the operator.
In order to keep the study more general possible, more 3-D hysteresis operators are taken into account. The magnetization paths of three concentric anisotropic hysteresis operators and the corresponding hysteresis loops are shown in Fig. 8 when e = −0.6, −1.0 and −1.4. And Fig.  9 shows the same case of the isotropic hysteresis operator. It is assumed that the paths of alternating magnetic field are (1.5, −1.5, −1.5) -(−1.5, 1.5, 1.5) -(1.5, −1.5, −1.5).
Comparing the hysteresis loops derived from a single hysteresis operator and three ones, it can be seen that a single operator has only one jump during magnetization, while the three concentric anisotropic hysteresis operators have three jumps through the critical surface during the magnetization. The results show that the energy parameter e can influence the hysteresis loop.
When considering the effect of the interaction field Hi, the total magnetic field is equal to the vector sum of Ha and Hi. The equation of equipotential surface of the hysteresis operators can be transformed to (24) Three anisotropic and isotropic hysteresis operators with the effect of interaction field are considered, respectively. The alternating magnetization paths of three non-concentric anisotropic hysteresis operators and the corresponding hysteresis loops are shown in Fig. 10 when e = −1. And Fig.11 shows the same case of three isotropic hysteresis operators. It is assumed that the paths of alternating magnetic field are (4, −4, −1) -(−4, 4, 1) -(4, −4, −1), and Hi =2.5.
As illustrated in Fig. 10 and 11, the interaction field affects the distribution position of hysteresis operator and thereby affects the hysteresis loop. As the number of operators increases and the distribution tends to be more uniform, the hysteresis loop will become more and more smooth. Therefore, once the number of operators and the distribution function is determined, the 3-D hybrid hysteresis model can be initially established. The results show the interaction field Hi can affect the hysteresis loop.
Analysis of the magnetization process shows that the factors influencing the hysteresis loop of the operator are energy parameters (related to coercive field), interaction field (affect the distribution of operators) and the number of operators.

B. Rotational Magnetization
The magnetization process of anisotropic and isotropic hysteresis operators is analyzed under rotational magnetic field. The excitation field is applied in xoy, yoz and xoz planes respectively. Three different sizes rotational magnetic field are applied to hysteresis operators, whose maximum magnetic field are set as 1.6, 3.2 and 4.8. The rotating direction of the rotational magnetic field is anticlockwise and Hi = 3.5.
The rotational magnetization paths of five anisotropic hysteresis operators and the corresponding hysteresis loops are shown in Fig. 12. And Fig. 13 shows the same case of five isotropic hysteresis operators.
It can be seen from the Fig. 12 and 13 that there has hysteresis phenomenon in the magnetization curves. When the applied field enters the critical surface of operators, the magnetization direction is frozen until it occurs barkhausen jump when the applied field exits the operator's critical surface. Conclusively, the rotational magnetization of hysteresis operator satisfies the modeling principle of hybrid vector hysteresis model.

V. PRELIMINARY STUDY AND EXPERIMENTAL VERIFICATION OF 3-D HYSTERESIS MODEL
The 3-D isotropic operator is firstly used in the 3-D modeling study due to its 3-D isotropy property. In order to simplify the calculation, the 3-D hysteresis model is regarded as the 2-D magnetization of its projection on the xoy, yoz and xoz planes. The previous numerical simulation method in [22] was improved by adjusting the relevant parameters and setting the cycle stop condition, and the closed magnetic characteristic curves could be obtained by the improved method.
To verify the validity of the 3-D model, the hysteresis properties of SMC material are simulated and compared with the experimental ones. The SMC material used in the experiment is SOMALOY TM 500, which maximum magnetic flux density is 2.1T, and the sample is a 22mm×22mm×22mm cube. The measuring device used in the experiment is the 3-D magnetic properties tester. [23], [24] The physical quantities measured in the experiment are the magnetic flux density and magnetic field, while what is used in the model are the magnetization and magnetic field, so the magnetic flux density B should be converted to magnetization M using (26).
The magnetic characteristics of the SMC material were measured by using the 3-D magnetic tester under the condition of Bm = 0.5 T, 0.7 T, 1 T, 1.4 T and 1.6 T.
By using the 3-D hysteresis model, 3-D magnetic properties of the SMC material are simulated. The main simulation steps are as follows: (1) A three-dimensional interaction field space is defined and discretized into nx  ny  nz points, the space boundary is HB, and the coordinate of the center point of each hysteresis operator is (Hix, Hiy, Hiz). VOLUME XX, 2017 6 (2) The shortest distance from the center point of the hysteresis operator (Hix, Hiy, Hiz) to the space boundary HB is set as R.
(3) The radius r considered can be obtained by discretizing the parameter R, from which the distribution law of the hysteresis operator can be obtained. Therefore, the number of hysteresis operators can be calculated by the equation (27). Among them, for each center point (Hix, Hiy,  Hiz), the value of nr depends on the size of r and HB. Therefore, when different sizes of external field are applied, the total number of operators is not a constant value.
(4) The applied external magnetic field is discretized, and the magnetization paths Hx, Hy and Hz are obtained. When calculating the magnetization, the external field is discretized as an input.
(5) Calculating the magnetization (mx, my, mz) of each (Hx, Hy, Hz) with respect to N hysteresis operators, then the final M is the formula (28).
where Ms is the saturation magnetization of the material, and  is the correction coefficient.
According to the above steps, after converting the input external magnetic field into the output magnetization, the simulated magnetization curve can be obtained.
The alternating simulation results are compared with the experimental results, as shown in Fig. 14. The simulation results in the z direction is not better than that in the x and y directions. The reason is that the z direction is the rolling direction, it will also bring the influence of some stress factors although the material is isotropic in theory.
The rotational simulation results are compared with the experimental results, as shown in Fig. 15. The black curves in the figure are the projection of the experimental measurement results on three planes, and the color curves are the projection of the simulation results on three planes. The parameters of the model need to be determined by several groups of 3-D B-H data. In order to simulate hysteresis more accurately, the least square method was used to fit the correction coefficients of each group of simulation results, which can simulate the data that has not been used before. In the future, intelligent algorithms will be considered for parameter identification.
In Fig. 15(a), the simulation results are more accurate when the magnetization is low, and the simulation accuracy decreases slightly with the increase of the magnetization. As shown in Fig. 15(b), with the increase of the magnetic field intensity, the simulation accuracy decreases obviously when Bm is greater than 1 T. The reason is that the relevant parameters of operators cannot accurately reflect the real distribution of magnetic domains of the material. Further study should be carried out to optimize the parameters by means of the experimental data.

VI. CONCLUSION
This paper proposed a method for the definition of hysteresis operator in 3-D case. According to the minimum energy principle of magnetization stable state and the modeling principle of hybrid vector hysteresis model, the critical surface equations of anisotropic and isotropic hysteresis operator were given in the space sphere coordinate system, respectively. Hysteresis operator is defined as a closed surface in H space. The magnetization direction of hysteresis operator is depending on the locations of the applied field. When the applied field is outside, the magnetization direction is perpendicular to the critical surface of operators. When the applied field is inside, the magnetization direction is frozen on the direction of the entry point until it occurs barkhausen jump to the direction of the exit point.
The magnetization of anisotropic and isotropic hysteresis operators was studied under alternating and rotational magnetic field, respectively. The results show that the factors influencing the hysteresis loop of the operator are energy parameters, interaction field and the number of operators. A 3-D hysteresis model was preliminarily established based on the 3-D isotropic hysteresis operator, and the experimental verification of SMC material was carried out under different magnetic field, which verified the validity of the 3-D hysteresis operator.