- Email: [email protected]
A challenging hysteresis operator for the simulation of Goss-textured magnetic materials
Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Contents lists available at ScienceDirect
Journal of Magnetism and Magnetic Materials journal homepage: www.elsevier.com/locate/jmmm
A challenging hysteresis operator for the simulation of Goss-textured magnetic materials Ermanno Cardelli a,b, Antonio Faba b,c, Antonino Laudani d, Michele Pompei a, Simone Quondam Antonio a,⇑, Francesco Riganti Fulginei d, Alessandro Salvini d a
Dipartimento di Ingegneria, Università di Perugia, Via G. Duranti, 67, 06125 Perugia, Italy Centre for Electric and Magnetic Applied Research, Perugia, Italy Polo Didattico Scientifico di Terni, Strada Pentima Bassa n. 4, 05100 Terni, Italy d Università Roma tre, Via Ostiense, 159, 00154 Roma, Italy b c
a r t i c l e
i n f o
Article history: Received 6 June 2016 Accepted 23 January 2017 Available online 25 January 2017 Keywords: Goss-texture Hysteresis operators Magnetocrystalline anisotropy Vector hysteresis modeling
a b s t r a c t A new hysteresis operator for the simulation of Goss-textured ferromagnets is here defined. The operator is derived from the classic Stoner–Wohlfarth model, where the anisotropy energy is assumed to be cubic instead of uniaxial, in order to reproduce the magnetic behavior of Goss textured ferromagnetic materials, such as grain-oriented Fe–Si alloys, Ni–Fe alloys, and Ni–Co alloys. A vector hysteresis model based on a single hysteresis operator is then implemented and used for the prediction of the rotational magnetizations that have been measured in a sample of grain-oriented electrical steel. This is especially promising for FEM based calculations, where the magnetization state in each point must be recalculated at each time step. Finally, the computed loops, as well as the magnetic losses, are compared to the measured data. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Vector hysteresis modeling for engineering applications is being developed by means of phenomenological approaches [1,2], due to the limitation of the computational cost required to solve micromagnetic equations into a macroscopic domain [3]. Typical phenomenological approaches are either based on probabilistic or deterministic distributions of hysteresis operators (hysterons) [4– 12]. An assembly of hysteresis operators in general can’t model with accuracy the magnetic behavior of ferromagnetic materials at macroscopic scale-length, unless a suitable magnetization dependent distribution (moving function) of the hysteresis operators is introduced [6,13–15]. This fact, however, leads at least to a couple of critical questions. First of all, high memory occupation and computation time, since the number of hysterons to use to have a good accuracy in results is almost high and the magnetization status of all the hysterons must be kept in memory at each time-step [16]. Furthermore, the definition of the distribution function and its dependence on the magnetization (moving modeling) is an intrinsically twice ill-posed problem, and we have to face ⇑ Corresponding author. E-mail address: [email protected] (S. Quondam Antonio). http://dx.doi.org/10.1016/j.jmmm.2017.01.068 0304-8853/Ó 2017 Elsevier B.V. All rights reserved.
in general a multi-solution problem. Although some papers have discussed possible strategies to avoid moving modeling, using a fixed assembly of hysterons [10], or conditioned moving modeling techniques [13], the complete solution of the problem is actually far from our knowledge. We are confident that an improvement could be achieved if the applied hysteresis operator is made closer to the physics of the magnetization. The need for complicated moving modeling and consequent ill posed identification of the model parameters is especially due to the use of hysterons with null of limited intrinsic property to describe individually magnetic anisotropy, such in case of hysterons with spherical or ellipsoidal critical surface [10,11]. More complex hysterons, based on the Stoner–Wohlfarth (SW) theory, seem to have more efficient capability in reproducing at macromagnetic scale-length the experimental anisotropic behavior of several ferromagnetic materials [17–25]. In this paper Goss-textured ferromagnetic materials, such as for the grain oriented iron silicon steels, and some Ni–Fe and Ni–Co alloys, are taken into account. The planar structure of these materials is realized by aligning their grains along the [0 0 1] direction with only one minimum energy axis in-plane, while the other two easy axes are oriented at 45 degrees respect to the plane of the structure. Consequently, a new hysteron is mathematically
15
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
defined in this paper, assuming to model circular magnetization processes in Goss-textured planar structures using a 2D operator, where the cubic magnetic anisotropy is described in the (1 1 0) plane, instead of the uniaxial anisotropy of the classical SW theory. The approach here presented seems especially promising because we may think to use a single hysteron to reproduce the experimental behavior of the material so that the computational time is strongly reduced. In addition, the moving strategy can be easily defined and the model parameters can be deduced in univocal way. These are important advantages if the model of Gosstextured ferromagnetic materials is matched with a FEM solver, as in the case of design of electrical machines, because the distribution of the magnetization state must be computed at each time step, and the use of an assembly of hysterons makes the calculation times too long and the memory allocation too large for a practical usage. In this paper the performance of the new hysteresis operator is checked and the reconstruction of rotational magnetizations measured in a sample of grain-oriented electrical steel (GOES) is presented. The experimental analysis takes into account also the estimation of the static magnetic losses and the comparison between measured and computed data will be presented. 2. Vector hysteresis operator The Stoner–Wohlfarth approximation is assumed here to define the normalized Gibb’s free energy per unit of volume in case of cubic magnetic anisotropy. Let R(O,x,y,z) be the fundamental frame of reference, and m the unitary magnetization vector measured respect to R. The general expansion of the cubic anisotropy energy as function of the unitary magnetization m in cartesian coordinates is given by
~ Ecubic AN ðmÞ ¼ K 1 ðm2x m2y þ m2y m2z þ m2z m2x Þ þ K 2 m2x m2y m2z V
ð1Þ
where V is the volume of the particle, while K1 and K2 the anisotropy constants. The values for the iron are K1 = 4.8 104 J/m and K2 = 0.5 104 J/m, while for the typical Fe–Si alloys with few percentage points of Si in Fe, both the values are slightly lower. It can be shown that the maximum error in the computation of the anisotropy energy that is committed if the second term of (1) is neglected is around 0.1%. For this reason we can put K2 = 0. Now, let R0 (O0 ,x0 ,y0 ,z0 ) be a second frame of reference rotated clockwise respect to R by an amount of 45° around the z axis. If the anisotropy energy is written as function of the magnetization m0 expressed respect to R0 , then the x0 -z0 plane coincides with one of the (1 1 0) equivalent planes. Refer to Fig. 1 for more clarity.
Fig. 1. (1 1 0) Miller plane in a cubic cell. In Goss-textured materials the rolling direction coincides with one of the h0 0 1i equivalent directions, while the lamination plane coincides with one of the {1 1 0} equivalent planes.
The normalized anisotropy energy can be written in the following form:
~0 Ecubic 1 1 4 2 AN ðm Þ ¼ sin h þ sin h cos2 h 2K 1 V 8 2
ð2Þ
where h is the angle of the unitary magnetization vector respect to the positive half-axis aligned with the unique in-plane easy direction. The obtained anisotropy energy presents two equivalent absolute minima for h = 0° and h = 180°, when the magnetization vector lie along the easy axis, but two additional relative energy minima can be found for h = 90° and h = 270°, due to the off-plane orientation of the other two easy axes. Four other equilibrium points that are energy maxima can be found for h ffi 55°, 125°, 235°, 305°. The normalized anisotropy energy as function of the magnetization angle h is pictured in Fig. 2. The normalized Gibb’s free energy for a given normalized external magnetic field hext = Hext/HAN is
g L ðh; hÞ ¼ h cosðhH hÞ þ
1 1 4 2 sin h þ sin h cos2 h 8 2
ð3Þ
where the anisotropy field is defined as HAN = 2K1/l0MS. The equilibrium state for the unitary magnetization vector is found by minimizing the expression (3) respect to h for a given external magnetic field hext that is aligned with the angle hH respect to the easy-axis direction. As for the classic SW operator [17], the number of stable states for the magnetization depends on the value of the applied magnetic field vector. In particular the system is characterized by a bifurcation set that is more complex than the classic SW astroid, and it is defined as the solution of the coupled equations gL0 = 0, gL00 = 0:
8 8 @g L UðhÞ VðhÞ tan h > > > > ¼ 0 > > < < hext;x ¼ @h 1 þ tan2 h ) 2 > > VðhÞ þ UðhÞ tan h > @ gL > > > : 2 ¼0 : hext;y ¼ 1 þ tan2 h @h 2
4
ð4Þ
3
h where UðhÞ ¼ 92 sin h cos h 2sin cos3 h and VðhÞ ¼ 12 sin h þ cos h
sin h cos2 h: The geometrical representation of the bifurcation curve, given in parametrical form by (4), is the closed curve on the hx hy plane pictured in Fig. 3. In the outer region A the external magnetic field is dominant respect to the anisotropy and there is only one stable state for the magnetization; in this case the calculation of the
Fig. 2. Normalized anisotropy energy for (1 1 0) in-plane anisotropy.
16
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Fig. 3. Bifurcation curve for the SW derived operator in case of cubic anisotropy cut along the (1 1 0) plane The curve divides the plane into the four regions A, B, C and D. The number of equilibrium states for the magnetization vector depends on the region the applied field belongs to.
magnetization angle is univocal. Conversely, for smaller fields the number of stable states is greater and the computation of the magnetization angle requires an opportune initialization of the previous magnetization state. In particular, if the normalized applied field hext belongs to the region B two stable states exists, in C three and finally in the inner D region four possible stable states are equivalent points of minimum energy for m. Anyway, for a given external field, the stable states for the magnetization are oriented along the angles between the tangents to the bifurcation curve passing through the point hext and the positive hx half-axis. One should conclude that, if an external field outside to the region A is applied and the previous magnetization state is unknown, then the actual stable state can’t be found. Therefore the starting point can either be near the saturation region, otherwise the magnetization angle must be initialized first. As previously mentioned each equilibrium state for the magnetization is given by the tangent to the bifurcation curve, that passes through the extrema point of the applied field vector on the hx hy plane. The number of the stable states is always equal to the number of the unstable states and, as already explained, it depends on the region where the applied field lie. If hext = 0 (gL = EAN/2K1V) the normalized field lies in D and the equilibrium states are given by those lines tangent to the bifurcation curve from the origin of the h-plane. It can be shown that in this case there are eight tangents corresponding to the eight equilibrium states of the anisotropy energy reported in Fig. 2. Suppose to simulate a quasi-stationary process characterized by a finite sequence of magnetic field vectors H1, H2, H3, . . ., Hn. Then the corresponding M1, M2, M3, . . ., Mn can be found when either H12A or the previous magnetization is defined. If at least one of these requirements is given, the magnetization can be computed uniquely for each value of the applied magnetic field. When the bifurcation set curve is crossed from the outer region to one of the B, C or D region, then the number of stable states becomes greater than one; we implemented a selection criterion based on the evaluation of the derivative gL0 (h) = dgL/dh on an opportune interval centered at the previous stable state of the magnetization [hp e, hp + e]. For e ? 0 the function gL0 can be either negative or positive, because the magnetic field has changed respect to the previous applied field and the new stable state for the magnetization is selected among the computed ones by choosing the closest
Fig. 4. Scalar hysteresis cycles for the SW derived operator obtained by applying a scalar normalized magnetic field hext = Hext/HAN along the directions hH = 0°, 30°, 45°, 60°, 90°: tangential component of the unitary magnetization (panel a) and orthogonal component (panel b).
angle that is found when the magnetization vector rotates toward lower energy values. In order to better illustrate the response of the hysteresis operator here accounted to an external applied field, some scalar cycles for hH = 0°, 30°, 45°, 60°, 90° have been computed, and in Fig. 4 both the tangential and orthogonal component of the magnetization are shown. It must underlined that this operator exhibits hysteresis also along the hard direction, contrarily to the traditional SW theory, because the hard axis is not symmetric respect the two adjacent easy axis. Furthermore, the presence of a secondary minimum anisotropy energy for hH = 90° is clearly represented by the associated squared loop pictured in the Fig. 4, but in this case the coercive field is less than the anisotropy field because the minimum values associated are higher than the principal one. 3. Single hysteron model and identification procedure The single hysteron model here presented is based on the hysteresis operator discussed in the previous section. The unitary magnetization is calculated according to the minimization of the normalized Gibb’s free energy (3). It must be noticed that the hysteron discussed is related to a single particle magnetic system and only the Zeeman energy and the anisotropy energy are considered in this case. The other physical processes that take place in a magnetic material at macroscopic scale-length, such as the domains formation and the Bloch walls motion, are not intrinsically included in the model and must be somehow modeled. In addition, the macroscopic behavior of Goss-textured ferromagnets is strongly affected by the quality of the grain boundaries, as well as the inclusions, whose effects is difficult, or almost impossible,
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
to account for, according to the physical theory. In our approach the operator, whose physical origin, even if strongly approximated is undoubted, requires to be mixed with a phenomenological quantity in order to overcome the necessity to model the mean effects due to the other energy terms. The moving technique developed here consists in the introduction of a phenomenological quantity, that according to the common usage in literature we can call ‘‘interaction field vector” Hi, that is added to the magnetic field H in order to evaluate the equivalent magnetic field, ‘‘seen” by the operator. The interaction field is obtained by measured data and whose practical identification with a particular physical situation is rather hard to get, being phenomenologically derived. It can be viewed as the mean field that results from the other energy terms, such as the magneto-elastic energy, the exchange energy, the demagnetizing energy, as well as the effects of the temperature, inclusions, impurities, crystallographic discontinuities such as grain boundaries. But, of course its relationship with each one of these additional items, some of them have to be stochastically treated, is not clear. It must be outlined that the vector H is the magnetic field that is measured at the surface of the material; theoretically it can be expressed as the sum of the field measured in vacuum Hext and the demagnetizing field Hdem which represents the response of the material. Therefore we have Heq = Hext + Hdem + Hi. Clearly, a more accurate modeling of such physical phenomena is far from the purpose to build a vector hysteresis model for bulk materials, because the introduction of a more accurate modeling brings to not acceptable calculation times and memory allocations, for real devices, such as electrical machines. In our viewpoint, the most important requirement is to short the CPU time and memory required for the calculation of the mean magnetization field of a ferromagnetic material at macroscopic scale length, in order to insert the model in FEM software for the design of electrical machines. The identification of the single-hysteron model consists, first of all, in the evaluation of the interaction field for a given magnetic field paths and relative magnetizations. In addition, the variation of the magnetization vector amplitude must be also modeled, since the single hysteron is a single domain particle characterized by a constant magnetic moment. We recognized that the most significant closed paths for H and M vectors, which can be used to identify the model, are those relative to circular magnetizations. For this reason we have done a set of measurements using a Round Rotational Single Sheet Tester (RRSST) with an opportune feedback technique which imposes the magnetization value to be constant during each cycle. The effect of the anisotropy and the demagnetizing effect are therefore accounted by the magnetic field that must be applied in order to obtain pure circular magnetizations. The identification procedure starts from the family of the measured magnetic field closed paths. The constant value of the magnetization vector is associated to each closed path and then the interaction field is computed. This is deducible from the measurements as the magnetic field, due to internal phenomena, that must be added to the measured field in order to make the unitary magnetization, computed by the anisotropic operator, aligned with the measured one. Finally, the amplitude correction is realized in the model by multiplying the unitary magnetization with the modulus that is associated to the corresponding H path. The computation of the interaction field for each measured field is not trivial, since, according to the theory, there are infinite magnetic fields that produce the same magnetization vector, that are all the fields tangent to the bifurcation curve associated to the given magnetization. It is reasonable that the interaction field must be as small as possible, that is to assume that the equivalent field must be as close as possible to the measured field. We have applied this criterion to
17
Fig. 5. Minimum normalized interaction field that must be added to the measured normalized field h in order to obtain the computed magnetization aligned along h = 75°.
evaluate the interaction field. An example for that is drawn in Fig. 5 for a magnetization angle of 75°. Once known both the components Hx, Hy, and the angle h of the magnetization, the normalized interaction field vector can be calculated as follows:
8 hy þ cot h hx VðhÞ > > hx < hi;x ¼ tan h þ cot h > h þ cot h hx VðhÞ > : hi;y ¼ y hy þ VðhÞ 1 þ cot2 h
ð5Þ
The identification algorithm has been written in MatlabÓ and it consists in a function that accepts as arguments the components of the measured field Hx, Hy and the components of the measured magnetization Mx, My. The function computes the angle h = a tan (My/Mx) and afterward the interaction field components using the Eq. (5). We reported the developed code in Appendix I. Once that the identification procedure is carried out for a known series of rotational magnetization paths, the singlehysteron model should be able to reconstruct the magnetic response of the magnetic material at macroscopic scale length when a generic sequence of quasi-static field vectors is applied. Let H(n) = Hncos hHn1x + Hnsin hHn1y be the n-th sample of the measured field vector sequence. If H(n) belongs to one of the closed paths that have been used for the identification of the material, both the interaction field vector Hi(n) and the magnetization modulus M(n) are immediately determined. Otherwise, the magnetic field lies between two consecutive closed paths, that we call here lower path and upper path, Hlow(hH), Hup(hH). A constant value of the magnetization amplitude and a hH-dependent interaction field are associated to each closed path. In this case the computation of both the interaction field vector and the magnetization magnitude requires a suitable interpolation technique, that is described below. Consider the Fig. 6. The line Hy = tan hH0Hx crosses the lower and the upper closed path in the points A and B, and the extrema point of the magnetic field H(n) divides the segment AB in two parts having length dlow and dup respectively. Let D = dlow + dup, we can define the relative distances clow = dlow/D, cup = dup/D that are the constants we used for the interpolation. Since the identification procedure gives the interaction field components for each closed path Hixlow(hH), Hixlow(hH), that are both functions of the magnetic field angle hH, we can now interpolate the value for the current field H.
18
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Fig. 7. Single-hysteron model block diagram: the equivalent field Heq = H + Hi is applied to the vector hysteresis operator (V.H.O.) and the computed unitary magnetization is multiplied by the amplitude M. Both M and Hi are given by the identification procedure.
multiplied by the numerical amplitude M. A block diagram of the single-hysteron model is reported in the figure below (see Fig. 7). 4. Simulations and measurements
Fig. 6. Interpolation parameters for a magnetic field with magnitude Hn along the direction specified by the angle hHn. The identified closed paths that, for hH = hHn are immediately lower and higher respect to the field, are indicated with the red dashed line and the blue line, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(
Hix ¼ 12 ðclow Hix;low ðhH0 Þ þ cup Hix;up ðhH0 ÞÞ Hiy ¼ 12 ðclow Hiy;low ðhH0 Þ þ cup Hiy;up ðhH0 ÞÞ
ð6Þ
The interpolation technique for the magnitude is being realized by using the same constants.
M¼
1 ðclow M low þ cup M up Þ 2
ð7Þ
We must point out that, theoretically, the rotational hysteresis can be simulated directly by means of a data interpolation between the measured cycles, but, in practice, this numerical approach is not applicable for some reasons. First of all, it is seen that a very much larger set of rotational cycles should be used, and the development of a finer set of measurements is a very difficult and onerous task, due to the high anisotropy of the material. In fact, if for two adjacent loops the magnetization amplitudes are very close, this is absolutely not true for the respective amplitudes of the magnetic field. In addition, the loop shape strongly changes when the materials is approaching to saturation. Finally, the complete and accurate sampling in the H and M plane, necessary for a good interpolation method, leads to a very high number of data, making the computation very slow. Taking into account such considerations, we have recognized the numerical advantages of a model, which is able to take into account the mean anisotropy of the material, respect to a pure numerical tool, like an interpolator. The single-hysteron model has been implemented in MatlabÓ; the numerical code it is composed by two functions and a principal M file. One function is used to calculate the equivalent field and to extrapolate the magnetization vector amplitude from the components of the magnetic field. The other one calculates the unitary magnetization according to the developed vector hysteresis operator; this last function is the core of the developed hysteron and we reported it in Appendix II. Finally, the unitary magnetization is
The single-hysteron model is now used for the simulation of rotational magnetization processes in a macroscopic sample of GOES sheet. The material under test is a disk having radius r = 3.2 cm, thickness t = 0.3 mm, silicon content around 3%, average grain size D = 1 cm. The measurement region is the central squared region with area A = 2 cm 2 cm. We can determine the values of both the magneto-crystalline anisotropy constant and the saturation polarization, from the silicon in iron content [26], and we obtain K1 = 3 105 J/m3 and JS = l0MS = 2.0 T. The measurements have been carried out by means of a RRSST system, as we have already anticipated, and they consists in a series of thirteen measured magnetic field paths Hi and the respective circular magnetic inductions Bi for i = 1, 2, . . ., 13. While the magnetic field locii are rather distorted due to the magnetic anisotropy, the feedback technique was able to make the magnetic induction perfectly circular, whit magnitude B = 0.20, 0.40, 0.60, 0.80, 0.90, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35 T. The odd measured loops have been used for the identification of the material, while the even loops have been used for the test. First of all the identification procedure is applied to the seven odd loops, obtained for i = 1, 3, 5, 7, 9, 11, 13. In Fig. 8 both the computed equivalent field and the measured field are displayed in the case of four identified loops. It is interesting to observe that along both the easy and hard directions we have Heq = H, that means Hi = 0. This interesting feature is verified not only at large magnetic induction, when the cubic (1 1 0) magnetic anisotropy of the Goss-texture is manifest, but also for lower magnetic induction, when the material shows a magnetic anisotropy that is quite uniaxial. The physical derivation of the hysteron discussed here, that is able to take into account the magnetic anisotropy, contributes to make the interaction field change very smooth and regular. This feature is no longer achievable by means of other hysterons, such as the Preisach type operators, since they have no connection with the physical mechanisms of the magnetization phenomena. In order to verify the congruency of the identification procedure we first present the simulations of the GOES starting from the identified (odd) rotational loops. When the measured magnetic field belongs to one of the closed paths used for the identification, the exact punctual values of the interaction field vector and the magnetization magnitude has been found, as expected, and one can observe the complete agreement between simulated and computed data in the left panel of Fig. 9. Then, the simulations have been carried out for the series of the even rotational loops, obtained for i = 2, 4, 6, 8, 10, 12, because the model must be able
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
19
Fig. 8. Measured magnetic field and computed equivalent field for the odd rotational loops B = 0.20 T (panel a), B = 0.90 T (panel b), B = 1.15 T (panel c) and B = 1.35 T (panel d).
Fig. 9. Reconstruction of circular magnetizations: comparison between simulated and measured data for both the rotational loops used for the identification (left panel) and those ones used for the test (right panel).
to reconstruct the magnetization processes mainly when the magnetic field path is different from the ones that have been used for the identification. In the right panel of Fig. 9 the computed and measured magnetizations are compared. For the test loops, a reliable reconstruction of the measured cycles is achieved, both in
amplitude and angle of the magnetization vector. These quantities are both of primary importance for the computation of the static magnetic losses, since the energy loss per cycle is proportional to the area of the hysteresis cycles (Mx–Hx and My–Hy curves); therefore, the model should be able not only to predict the shape of the
20
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Fig. 10. Simulated and measured lag angle Dh(hH) relative to three even rotational loops obtained by applying rotating magnetic inductions having amplitude 0.40 T (upper panel), 1.00 T (center panel) and 1.30 T (lower panel).
closed paths of the magnetization for a set of applied fields, but also the magnetization angle or the phase displacement between M and H. Such quantity, called lag angle, is usually evaluated as function of the magnetic field angle hH. In Fig. 10 the lag angle relative to three even loops (B = 0.40 T, B = 1.00 T, B = 1.30 T) is plotted versus the corresponding magnetic field angle hH. The good agreement between the simulated and the measured values is clearly visible. It is interesting to note that the behavior of the lag angle versus the magnetic field angle shows an interesting feature of this material, which presents an in-plane uniaxial anisotropy for low inductions, although the Goss-texture should be characterized by a cubic anisotropy. In the upper panel of Fig. 10 the maximum delay of the magnetization is found for hH ? 90° when the global magnetization tends to remain aligned with the positive easy semi-axis. For hH immediately greater than 90° the magnetization vector jumps toward the negative easy semi-axis and remains in advance respect to the field. In the interval [180°, +180°] the lag angle curve presents two minima and two maxima points. When the magnetic induction increases, the lag angle tends to be distorted (center panel of Fig. 10) but only for B > 1.10 T the real cubic anisotropy is manifest, as in the loop B = 1.30 T shown in the lower panel of Fig. 10. In this case the maximum delay of M is found approximately around 45°, followed by a slight anticipation due to the secondary energy minima that is found along the y axis for the (1 1 0) in-plane cubic anisotropy. The number of minima and maxima point of the lag angle are doubled. Our analysis is here concluded with a discussions about the magnetic losses. The energy lost per cycle is computed using the equation:
Eloss ¼ Eloss;x þ Eloss;y ¼
1 cT
Z dBx dBy Hx þ Hy dt dt dt T
ð8Þ
Fig. 11. Energy loss per cycle as function of the magnetic induction magnitude calculated by means of both simulated and measured data.
The integral (8) can be solved by evaluating numerically the area of the hysteresis cycles along x and y directions. The rotational losses present the typical behavior shown in Fig. 11, where the comparative analysis between the computed and measured data is reported. For low inductions, the measured field is also very low and the area contribution, that depends on the average lag angle during one cycle, is clearly small. If B < 1.05 T the magnetic anisotropy tends to be similar to the uniaxial one. In this interval the energy loss is an increasing function of the magnetic induction: the major contribution to the energy loss is due to the increment of the fields, while the average lag angle does not vary considerably from one cycle to the other, as
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
one can see in Fig. 10. If B > 1.05 T the cubic anisotropy behavior is manifest and the lag angle tends to be almost constant. Although the maximum misalignment between the magnetic field and the magnetic induction decreases slightly, the average phase displacement during one cycle tends to zero. In this region the lag angle is dominant over the absolute values of the vector fields and the magnetic losses decrease. For this reason it is evident that an accurate reconstruction of the lag angle is necessary in the correct prediction of the magnetic losses. Of course the computed losses associated to the odd rotational loops, used to identify the model, are closer to the measured data, respect to the ones associated to the even rotational loops. However the maximum displacement defined as e = |Eloss,sim Eloss,meas| is below 4 mJ/kg. As final technical comment, let us underline that the simulation of the thirteen rotational loops were carried out in the same computer, equipped with a CPU Intel Core i7 @ 2.20 GHz. The single hysteron model required only few seconds to reconstruct the entire series of rotational magnetization processes.
21
The approach is derived by the classical SW theory applied to the case of cubic magnetic anisotropy defined in the (1 1 0) plane. Vector hysteresis in Goss-textured ferromagnetic material can be accurately modeled at macro-magnetic scale length by means of a single hysteron. The moving modeling strategy defined in this case leads to an univocal calculation of the model parameters, deduced by vector measurements. The comparison with experimental data of GOES indicates that the new hysteron can predict with accuracy the magnetization trajectories, the lag angle and the static losses for circular magnetizations. The hysteron presented seems especially promising for a matching of the Goss-textured ferromagnetic materials modeling with FEM analysis in time domain, because puts together numerical accuracy and strong computation time and memory allocation saving. Scalar magnetization modeling deserves further studies. We are considering the possibility to model scalar hysteresis by means of the reconstruction of elliptic loops with increasing values of their aspect ratio, but this is a future work. Appendix I
5. Conclusions and further Research In this paper we have presented a new approach for the phenomenological modeling of Goss-textured ferromagnetic materials.
In this section we list the MatlabÓ code that we have written to implement the identification function of the single hysteron model.
22
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Appendix II In this section we list the MatlabÓ code that we have written to implement the magnetization calculation according to the single hysteron model.
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
References [1] C. Appino, C. Ragusa, F. Fiorillo, Can rotational magnetization be theoretically assessed ?, Int J. Appl. Electromagnet. Mech. 44 (2014) 355–370. [2] T. Matsuo, Anisotropic vector hysteresis model using an isotropic vector play model, IEEE Trans. Magn. 46 (8) (2010) 3041–3044 5512912. [3] E. Cardelli, A. Faba, G. Finocchio, B. Azzerboni, Mathematical modelling of magnetic hysteresis in exchange-bias spin valves, IEEE Trans. Magn. 48 (2012) 3367–3370. [4] F. Preisach, Über die magnetische Nachwirkung, Z. Angew. Phys. 94 (1935) 277–302. [5] A. Visintin, Differential models of hysteresis, in: Applied Mathematical Sciences, 111, Springer, Berlin, 1994. [6] E. Della Torre, Magnetic Hysteresis, IEEE Press, Piscataway, NJ, 1999. [7] E. Della Torre, E. Pinzaglia, E. Cardelli, Vector modeling – part I: generalized hysteresis model, Phys. B 372 (2006) 111–114. [8] E. Della Torre, E. Pinzaglia, E. Cardelli, Vector modeling: Part II, ellipsoidal vector hysteresis model – numerical application to a 2-d case, Phys. B 372 (2006) 115–119. [9] L. Dupré, G. Bàn, M. Von Rauch, J. Melkebeek, Relation between the microstructural properties of electrical steels and the Preisach modelling, J. Magn. Magn. Mater. 195 (1999) 233–249. [10] E. Cardelli, E. Della Torre, A. Faba, A general vector hysteresis operator: extension to the 3-D case, IEEE Trans. Magn. 46 (12) (December 2010) 3990– 4000. [11] E. Cardelli, A general hysteresis operator for the modeling of vector fields, IEEE Trans. Magn. 47 (8) (2011) 2056–2067. [12] Z. Szabo, J. Füzi, Implementation and identification of Preisach type hysteresis models with Everett Function in closed form, J. Magn. Magn. Mater. 406 (2016) 251–258. [13] E. Cardelli, A. Faba, Genetic algorithm identification of a H-moving vector hysteresis model, Phys. B 435 (February 2014) 11–15.
23
[14] E. Cardelli, A. Faba, A. Laudani, F. Riganti Fulginei, A. Salvini, A moving approach to magnetic modeling of electrical steels in 2-d, Int. J. Appl. Electromagnet. Mech. 48 (2,3) (2015) 263–270. [15] E. Cardelli, A. Faba, A. Laudani, S. Quondam Antonio, F. Riganti Fulginei, A. Salvini, A moving approach for the Vector Hysteron Model, Phys. B 486 (1) (April 2016) 92–96. [16] E. Cardelli, E. Della Torre, A. Faba, Numerical implementation of the DPC model, IEEE Trans. Magn. 45 (3) (2009) 1186–1189. [17] E.C. Stoner, E.P. Wohlfarth, Philos. Trans. R. Soc. London, Ser. A, Math. Phys. Sci. 240 (826) (1948) 599–642. [18] M. de Campos, F. Sampaio da Silva, E. Perigo, J. de Castro, Stoner–Wohlfarth model for the anisotropic case, J. Magn. Magn. Mater. 345 (2013) 147–152. [19] A. Kuncser, V. Kuncser, Magnetization reversal via a Stoner–Wohlfarth model with bi-dimensional angular distribution of easy axis, J. Magn. Magn. Mater. 395 (2015) 34–40. [20] C. Chang, Micromagnetic studies of coherent rotation with quartic crystalline anisotropy, J. Appl. Phys. 69 (1991) 2431. [21] A. Thiaville, Extensions of the geometric solution of the two dimensional coherent magnetization rotation model, J. Magn. Magn. Mater. 182 (1998) 5– 18. [22] A. Thiaville, Coherent rotation of magnetization in three dimensions. A geometrical approach, Phys. Rev. B 61 (2000) 12221–12232. [23] A. Stancu, L. Spinu, Transverse susceptibility for single-domain particle with cubic anisotropy, J. Magn. Magn. Mater. 266 (2003) 200–206. [24] J. Petrila, A. Stancu, Analytical vector generalization of the classical Stoner– Wohlfarth hysteron, J. Phys.: Condens. Matter 23 (2011) 076002. [25] D. Cimpoesu, L. Stoleriu, A. Stancu, Generalized Stoner–Wohlfarth model accurately describing the switching processes in pseudo-single ferromagnetic particles, J. Appl. Phys. 114 (2013) 223901. [26] Bozorth, Ferromagnetism, IEEE Press, 1993.
Contents lists available at ScienceDirect
Journal of Magnetism and Magnetic Materials journal homepage: www.elsevier.com/locate/jmmm
A challenging hysteresis operator for the simulation of Goss-textured magnetic materials Ermanno Cardelli a,b, Antonio Faba b,c, Antonino Laudani d, Michele Pompei a, Simone Quondam Antonio a,⇑, Francesco Riganti Fulginei d, Alessandro Salvini d a
Dipartimento di Ingegneria, Università di Perugia, Via G. Duranti, 67, 06125 Perugia, Italy Centre for Electric and Magnetic Applied Research, Perugia, Italy Polo Didattico Scientifico di Terni, Strada Pentima Bassa n. 4, 05100 Terni, Italy d Università Roma tre, Via Ostiense, 159, 00154 Roma, Italy b c
a r t i c l e
i n f o
Article history: Received 6 June 2016 Accepted 23 January 2017 Available online 25 January 2017 Keywords: Goss-texture Hysteresis operators Magnetocrystalline anisotropy Vector hysteresis modeling
a b s t r a c t A new hysteresis operator for the simulation of Goss-textured ferromagnets is here defined. The operator is derived from the classic Stoner–Wohlfarth model, where the anisotropy energy is assumed to be cubic instead of uniaxial, in order to reproduce the magnetic behavior of Goss textured ferromagnetic materials, such as grain-oriented Fe–Si alloys, Ni–Fe alloys, and Ni–Co alloys. A vector hysteresis model based on a single hysteresis operator is then implemented and used for the prediction of the rotational magnetizations that have been measured in a sample of grain-oriented electrical steel. This is especially promising for FEM based calculations, where the magnetization state in each point must be recalculated at each time step. Finally, the computed loops, as well as the magnetic losses, are compared to the measured data. Ó 2017 Elsevier B.V. All rights reserved.
1. Introduction Vector hysteresis modeling for engineering applications is being developed by means of phenomenological approaches [1,2], due to the limitation of the computational cost required to solve micromagnetic equations into a macroscopic domain [3]. Typical phenomenological approaches are either based on probabilistic or deterministic distributions of hysteresis operators (hysterons) [4– 12]. An assembly of hysteresis operators in general can’t model with accuracy the magnetic behavior of ferromagnetic materials at macroscopic scale-length, unless a suitable magnetization dependent distribution (moving function) of the hysteresis operators is introduced [6,13–15]. This fact, however, leads at least to a couple of critical questions. First of all, high memory occupation and computation time, since the number of hysterons to use to have a good accuracy in results is almost high and the magnetization status of all the hysterons must be kept in memory at each time-step [16]. Furthermore, the definition of the distribution function and its dependence on the magnetization (moving modeling) is an intrinsically twice ill-posed problem, and we have to face ⇑ Corresponding author. E-mail address: [email protected] (S. Quondam Antonio). http://dx.doi.org/10.1016/j.jmmm.2017.01.068 0304-8853/Ó 2017 Elsevier B.V. All rights reserved.
in general a multi-solution problem. Although some papers have discussed possible strategies to avoid moving modeling, using a fixed assembly of hysterons [10], or conditioned moving modeling techniques [13], the complete solution of the problem is actually far from our knowledge. We are confident that an improvement could be achieved if the applied hysteresis operator is made closer to the physics of the magnetization. The need for complicated moving modeling and consequent ill posed identification of the model parameters is especially due to the use of hysterons with null of limited intrinsic property to describe individually magnetic anisotropy, such in case of hysterons with spherical or ellipsoidal critical surface [10,11]. More complex hysterons, based on the Stoner–Wohlfarth (SW) theory, seem to have more efficient capability in reproducing at macromagnetic scale-length the experimental anisotropic behavior of several ferromagnetic materials [17–25]. In this paper Goss-textured ferromagnetic materials, such as for the grain oriented iron silicon steels, and some Ni–Fe and Ni–Co alloys, are taken into account. The planar structure of these materials is realized by aligning their grains along the [0 0 1] direction with only one minimum energy axis in-plane, while the other two easy axes are oriented at 45 degrees respect to the plane of the structure. Consequently, a new hysteron is mathematically
15
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
defined in this paper, assuming to model circular magnetization processes in Goss-textured planar structures using a 2D operator, where the cubic magnetic anisotropy is described in the (1 1 0) plane, instead of the uniaxial anisotropy of the classical SW theory. The approach here presented seems especially promising because we may think to use a single hysteron to reproduce the experimental behavior of the material so that the computational time is strongly reduced. In addition, the moving strategy can be easily defined and the model parameters can be deduced in univocal way. These are important advantages if the model of Gosstextured ferromagnetic materials is matched with a FEM solver, as in the case of design of electrical machines, because the distribution of the magnetization state must be computed at each time step, and the use of an assembly of hysterons makes the calculation times too long and the memory allocation too large for a practical usage. In this paper the performance of the new hysteresis operator is checked and the reconstruction of rotational magnetizations measured in a sample of grain-oriented electrical steel (GOES) is presented. The experimental analysis takes into account also the estimation of the static magnetic losses and the comparison between measured and computed data will be presented. 2. Vector hysteresis operator The Stoner–Wohlfarth approximation is assumed here to define the normalized Gibb’s free energy per unit of volume in case of cubic magnetic anisotropy. Let R(O,x,y,z) be the fundamental frame of reference, and m the unitary magnetization vector measured respect to R. The general expansion of the cubic anisotropy energy as function of the unitary magnetization m in cartesian coordinates is given by
~ Ecubic AN ðmÞ ¼ K 1 ðm2x m2y þ m2y m2z þ m2z m2x Þ þ K 2 m2x m2y m2z V
ð1Þ
where V is the volume of the particle, while K1 and K2 the anisotropy constants. The values for the iron are K1 = 4.8 104 J/m and K2 = 0.5 104 J/m, while for the typical Fe–Si alloys with few percentage points of Si in Fe, both the values are slightly lower. It can be shown that the maximum error in the computation of the anisotropy energy that is committed if the second term of (1) is neglected is around 0.1%. For this reason we can put K2 = 0. Now, let R0 (O0 ,x0 ,y0 ,z0 ) be a second frame of reference rotated clockwise respect to R by an amount of 45° around the z axis. If the anisotropy energy is written as function of the magnetization m0 expressed respect to R0 , then the x0 -z0 plane coincides with one of the (1 1 0) equivalent planes. Refer to Fig. 1 for more clarity.
Fig. 1. (1 1 0) Miller plane in a cubic cell. In Goss-textured materials the rolling direction coincides with one of the h0 0 1i equivalent directions, while the lamination plane coincides with one of the {1 1 0} equivalent planes.
The normalized anisotropy energy can be written in the following form:
~0 Ecubic 1 1 4 2 AN ðm Þ ¼ sin h þ sin h cos2 h 2K 1 V 8 2
ð2Þ
where h is the angle of the unitary magnetization vector respect to the positive half-axis aligned with the unique in-plane easy direction. The obtained anisotropy energy presents two equivalent absolute minima for h = 0° and h = 180°, when the magnetization vector lie along the easy axis, but two additional relative energy minima can be found for h = 90° and h = 270°, due to the off-plane orientation of the other two easy axes. Four other equilibrium points that are energy maxima can be found for h ffi 55°, 125°, 235°, 305°. The normalized anisotropy energy as function of the magnetization angle h is pictured in Fig. 2. The normalized Gibb’s free energy for a given normalized external magnetic field hext = Hext/HAN is
g L ðh; hÞ ¼ h cosðhH hÞ þ
1 1 4 2 sin h þ sin h cos2 h 8 2
ð3Þ
where the anisotropy field is defined as HAN = 2K1/l0MS. The equilibrium state for the unitary magnetization vector is found by minimizing the expression (3) respect to h for a given external magnetic field hext that is aligned with the angle hH respect to the easy-axis direction. As for the classic SW operator [17], the number of stable states for the magnetization depends on the value of the applied magnetic field vector. In particular the system is characterized by a bifurcation set that is more complex than the classic SW astroid, and it is defined as the solution of the coupled equations gL0 = 0, gL00 = 0:
8 8 @g L UðhÞ VðhÞ tan h > > > > ¼ 0 > > < < hext;x ¼ @h 1 þ tan2 h ) 2 > > VðhÞ þ UðhÞ tan h > @ gL > > > : 2 ¼0 : hext;y ¼ 1 þ tan2 h @h 2
4
ð4Þ
3
h where UðhÞ ¼ 92 sin h cos h 2sin cos3 h and VðhÞ ¼ 12 sin h þ cos h
sin h cos2 h: The geometrical representation of the bifurcation curve, given in parametrical form by (4), is the closed curve on the hx hy plane pictured in Fig. 3. In the outer region A the external magnetic field is dominant respect to the anisotropy and there is only one stable state for the magnetization; in this case the calculation of the
Fig. 2. Normalized anisotropy energy for (1 1 0) in-plane anisotropy.
16
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Fig. 3. Bifurcation curve for the SW derived operator in case of cubic anisotropy cut along the (1 1 0) plane The curve divides the plane into the four regions A, B, C and D. The number of equilibrium states for the magnetization vector depends on the region the applied field belongs to.
magnetization angle is univocal. Conversely, for smaller fields the number of stable states is greater and the computation of the magnetization angle requires an opportune initialization of the previous magnetization state. In particular, if the normalized applied field hext belongs to the region B two stable states exists, in C three and finally in the inner D region four possible stable states are equivalent points of minimum energy for m. Anyway, for a given external field, the stable states for the magnetization are oriented along the angles between the tangents to the bifurcation curve passing through the point hext and the positive hx half-axis. One should conclude that, if an external field outside to the region A is applied and the previous magnetization state is unknown, then the actual stable state can’t be found. Therefore the starting point can either be near the saturation region, otherwise the magnetization angle must be initialized first. As previously mentioned each equilibrium state for the magnetization is given by the tangent to the bifurcation curve, that passes through the extrema point of the applied field vector on the hx hy plane. The number of the stable states is always equal to the number of the unstable states and, as already explained, it depends on the region where the applied field lie. If hext = 0 (gL = EAN/2K1V) the normalized field lies in D and the equilibrium states are given by those lines tangent to the bifurcation curve from the origin of the h-plane. It can be shown that in this case there are eight tangents corresponding to the eight equilibrium states of the anisotropy energy reported in Fig. 2. Suppose to simulate a quasi-stationary process characterized by a finite sequence of magnetic field vectors H1, H2, H3, . . ., Hn. Then the corresponding M1, M2, M3, . . ., Mn can be found when either H12A or the previous magnetization is defined. If at least one of these requirements is given, the magnetization can be computed uniquely for each value of the applied magnetic field. When the bifurcation set curve is crossed from the outer region to one of the B, C or D region, then the number of stable states becomes greater than one; we implemented a selection criterion based on the evaluation of the derivative gL0 (h) = dgL/dh on an opportune interval centered at the previous stable state of the magnetization [hp e, hp + e]. For e ? 0 the function gL0 can be either negative or positive, because the magnetic field has changed respect to the previous applied field and the new stable state for the magnetization is selected among the computed ones by choosing the closest
Fig. 4. Scalar hysteresis cycles for the SW derived operator obtained by applying a scalar normalized magnetic field hext = Hext/HAN along the directions hH = 0°, 30°, 45°, 60°, 90°: tangential component of the unitary magnetization (panel a) and orthogonal component (panel b).
angle that is found when the magnetization vector rotates toward lower energy values. In order to better illustrate the response of the hysteresis operator here accounted to an external applied field, some scalar cycles for hH = 0°, 30°, 45°, 60°, 90° have been computed, and in Fig. 4 both the tangential and orthogonal component of the magnetization are shown. It must underlined that this operator exhibits hysteresis also along the hard direction, contrarily to the traditional SW theory, because the hard axis is not symmetric respect the two adjacent easy axis. Furthermore, the presence of a secondary minimum anisotropy energy for hH = 90° is clearly represented by the associated squared loop pictured in the Fig. 4, but in this case the coercive field is less than the anisotropy field because the minimum values associated are higher than the principal one. 3. Single hysteron model and identification procedure The single hysteron model here presented is based on the hysteresis operator discussed in the previous section. The unitary magnetization is calculated according to the minimization of the normalized Gibb’s free energy (3). It must be noticed that the hysteron discussed is related to a single particle magnetic system and only the Zeeman energy and the anisotropy energy are considered in this case. The other physical processes that take place in a magnetic material at macroscopic scale-length, such as the domains formation and the Bloch walls motion, are not intrinsically included in the model and must be somehow modeled. In addition, the macroscopic behavior of Goss-textured ferromagnets is strongly affected by the quality of the grain boundaries, as well as the inclusions, whose effects is difficult, or almost impossible,
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
to account for, according to the physical theory. In our approach the operator, whose physical origin, even if strongly approximated is undoubted, requires to be mixed with a phenomenological quantity in order to overcome the necessity to model the mean effects due to the other energy terms. The moving technique developed here consists in the introduction of a phenomenological quantity, that according to the common usage in literature we can call ‘‘interaction field vector” Hi, that is added to the magnetic field H in order to evaluate the equivalent magnetic field, ‘‘seen” by the operator. The interaction field is obtained by measured data and whose practical identification with a particular physical situation is rather hard to get, being phenomenologically derived. It can be viewed as the mean field that results from the other energy terms, such as the magneto-elastic energy, the exchange energy, the demagnetizing energy, as well as the effects of the temperature, inclusions, impurities, crystallographic discontinuities such as grain boundaries. But, of course its relationship with each one of these additional items, some of them have to be stochastically treated, is not clear. It must be outlined that the vector H is the magnetic field that is measured at the surface of the material; theoretically it can be expressed as the sum of the field measured in vacuum Hext and the demagnetizing field Hdem which represents the response of the material. Therefore we have Heq = Hext + Hdem + Hi. Clearly, a more accurate modeling of such physical phenomena is far from the purpose to build a vector hysteresis model for bulk materials, because the introduction of a more accurate modeling brings to not acceptable calculation times and memory allocations, for real devices, such as electrical machines. In our viewpoint, the most important requirement is to short the CPU time and memory required for the calculation of the mean magnetization field of a ferromagnetic material at macroscopic scale length, in order to insert the model in FEM software for the design of electrical machines. The identification of the single-hysteron model consists, first of all, in the evaluation of the interaction field for a given magnetic field paths and relative magnetizations. In addition, the variation of the magnetization vector amplitude must be also modeled, since the single hysteron is a single domain particle characterized by a constant magnetic moment. We recognized that the most significant closed paths for H and M vectors, which can be used to identify the model, are those relative to circular magnetizations. For this reason we have done a set of measurements using a Round Rotational Single Sheet Tester (RRSST) with an opportune feedback technique which imposes the magnetization value to be constant during each cycle. The effect of the anisotropy and the demagnetizing effect are therefore accounted by the magnetic field that must be applied in order to obtain pure circular magnetizations. The identification procedure starts from the family of the measured magnetic field closed paths. The constant value of the magnetization vector is associated to each closed path and then the interaction field is computed. This is deducible from the measurements as the magnetic field, due to internal phenomena, that must be added to the measured field in order to make the unitary magnetization, computed by the anisotropic operator, aligned with the measured one. Finally, the amplitude correction is realized in the model by multiplying the unitary magnetization with the modulus that is associated to the corresponding H path. The computation of the interaction field for each measured field is not trivial, since, according to the theory, there are infinite magnetic fields that produce the same magnetization vector, that are all the fields tangent to the bifurcation curve associated to the given magnetization. It is reasonable that the interaction field must be as small as possible, that is to assume that the equivalent field must be as close as possible to the measured field. We have applied this criterion to
17
Fig. 5. Minimum normalized interaction field that must be added to the measured normalized field h in order to obtain the computed magnetization aligned along h = 75°.
evaluate the interaction field. An example for that is drawn in Fig. 5 for a magnetization angle of 75°. Once known both the components Hx, Hy, and the angle h of the magnetization, the normalized interaction field vector can be calculated as follows:
8 hy þ cot h hx VðhÞ > > hx < hi;x ¼ tan h þ cot h > h þ cot h hx VðhÞ > : hi;y ¼ y hy þ VðhÞ 1 þ cot2 h
ð5Þ
The identification algorithm has been written in MatlabÓ and it consists in a function that accepts as arguments the components of the measured field Hx, Hy and the components of the measured magnetization Mx, My. The function computes the angle h = a tan (My/Mx) and afterward the interaction field components using the Eq. (5). We reported the developed code in Appendix I. Once that the identification procedure is carried out for a known series of rotational magnetization paths, the singlehysteron model should be able to reconstruct the magnetic response of the magnetic material at macroscopic scale length when a generic sequence of quasi-static field vectors is applied. Let H(n) = Hncos hHn1x + Hnsin hHn1y be the n-th sample of the measured field vector sequence. If H(n) belongs to one of the closed paths that have been used for the identification of the material, both the interaction field vector Hi(n) and the magnetization modulus M(n) are immediately determined. Otherwise, the magnetic field lies between two consecutive closed paths, that we call here lower path and upper path, Hlow(hH), Hup(hH). A constant value of the magnetization amplitude and a hH-dependent interaction field are associated to each closed path. In this case the computation of both the interaction field vector and the magnetization magnitude requires a suitable interpolation technique, that is described below. Consider the Fig. 6. The line Hy = tan hH0Hx crosses the lower and the upper closed path in the points A and B, and the extrema point of the magnetic field H(n) divides the segment AB in two parts having length dlow and dup respectively. Let D = dlow + dup, we can define the relative distances clow = dlow/D, cup = dup/D that are the constants we used for the interpolation. Since the identification procedure gives the interaction field components for each closed path Hixlow(hH), Hixlow(hH), that are both functions of the magnetic field angle hH, we can now interpolate the value for the current field H.
18
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Fig. 7. Single-hysteron model block diagram: the equivalent field Heq = H + Hi is applied to the vector hysteresis operator (V.H.O.) and the computed unitary magnetization is multiplied by the amplitude M. Both M and Hi are given by the identification procedure.
multiplied by the numerical amplitude M. A block diagram of the single-hysteron model is reported in the figure below (see Fig. 7). 4. Simulations and measurements
Fig. 6. Interpolation parameters for a magnetic field with magnitude Hn along the direction specified by the angle hHn. The identified closed paths that, for hH = hHn are immediately lower and higher respect to the field, are indicated with the red dashed line and the blue line, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(
Hix ¼ 12 ðclow Hix;low ðhH0 Þ þ cup Hix;up ðhH0 ÞÞ Hiy ¼ 12 ðclow Hiy;low ðhH0 Þ þ cup Hiy;up ðhH0 ÞÞ
ð6Þ
The interpolation technique for the magnitude is being realized by using the same constants.
M¼
1 ðclow M low þ cup M up Þ 2
ð7Þ
We must point out that, theoretically, the rotational hysteresis can be simulated directly by means of a data interpolation between the measured cycles, but, in practice, this numerical approach is not applicable for some reasons. First of all, it is seen that a very much larger set of rotational cycles should be used, and the development of a finer set of measurements is a very difficult and onerous task, due to the high anisotropy of the material. In fact, if for two adjacent loops the magnetization amplitudes are very close, this is absolutely not true for the respective amplitudes of the magnetic field. In addition, the loop shape strongly changes when the materials is approaching to saturation. Finally, the complete and accurate sampling in the H and M plane, necessary for a good interpolation method, leads to a very high number of data, making the computation very slow. Taking into account such considerations, we have recognized the numerical advantages of a model, which is able to take into account the mean anisotropy of the material, respect to a pure numerical tool, like an interpolator. The single-hysteron model has been implemented in MatlabÓ; the numerical code it is composed by two functions and a principal M file. One function is used to calculate the equivalent field and to extrapolate the magnetization vector amplitude from the components of the magnetic field. The other one calculates the unitary magnetization according to the developed vector hysteresis operator; this last function is the core of the developed hysteron and we reported it in Appendix II. Finally, the unitary magnetization is
The single-hysteron model is now used for the simulation of rotational magnetization processes in a macroscopic sample of GOES sheet. The material under test is a disk having radius r = 3.2 cm, thickness t = 0.3 mm, silicon content around 3%, average grain size D = 1 cm. The measurement region is the central squared region with area A = 2 cm 2 cm. We can determine the values of both the magneto-crystalline anisotropy constant and the saturation polarization, from the silicon in iron content [26], and we obtain K1 = 3 105 J/m3 and JS = l0MS = 2.0 T. The measurements have been carried out by means of a RRSST system, as we have already anticipated, and they consists in a series of thirteen measured magnetic field paths Hi and the respective circular magnetic inductions Bi for i = 1, 2, . . ., 13. While the magnetic field locii are rather distorted due to the magnetic anisotropy, the feedback technique was able to make the magnetic induction perfectly circular, whit magnitude B = 0.20, 0.40, 0.60, 0.80, 0.90, 1.00, 1.05, 1.10, 1.15, 1.20, 1.25, 1.30, 1.35 T. The odd measured loops have been used for the identification of the material, while the even loops have been used for the test. First of all the identification procedure is applied to the seven odd loops, obtained for i = 1, 3, 5, 7, 9, 11, 13. In Fig. 8 both the computed equivalent field and the measured field are displayed in the case of four identified loops. It is interesting to observe that along both the easy and hard directions we have Heq = H, that means Hi = 0. This interesting feature is verified not only at large magnetic induction, when the cubic (1 1 0) magnetic anisotropy of the Goss-texture is manifest, but also for lower magnetic induction, when the material shows a magnetic anisotropy that is quite uniaxial. The physical derivation of the hysteron discussed here, that is able to take into account the magnetic anisotropy, contributes to make the interaction field change very smooth and regular. This feature is no longer achievable by means of other hysterons, such as the Preisach type operators, since they have no connection with the physical mechanisms of the magnetization phenomena. In order to verify the congruency of the identification procedure we first present the simulations of the GOES starting from the identified (odd) rotational loops. When the measured magnetic field belongs to one of the closed paths used for the identification, the exact punctual values of the interaction field vector and the magnetization magnitude has been found, as expected, and one can observe the complete agreement between simulated and computed data in the left panel of Fig. 9. Then, the simulations have been carried out for the series of the even rotational loops, obtained for i = 2, 4, 6, 8, 10, 12, because the model must be able
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
19
Fig. 8. Measured magnetic field and computed equivalent field for the odd rotational loops B = 0.20 T (panel a), B = 0.90 T (panel b), B = 1.15 T (panel c) and B = 1.35 T (panel d).
Fig. 9. Reconstruction of circular magnetizations: comparison between simulated and measured data for both the rotational loops used for the identification (left panel) and those ones used for the test (right panel).
to reconstruct the magnetization processes mainly when the magnetic field path is different from the ones that have been used for the identification. In the right panel of Fig. 9 the computed and measured magnetizations are compared. For the test loops, a reliable reconstruction of the measured cycles is achieved, both in
amplitude and angle of the magnetization vector. These quantities are both of primary importance for the computation of the static magnetic losses, since the energy loss per cycle is proportional to the area of the hysteresis cycles (Mx–Hx and My–Hy curves); therefore, the model should be able not only to predict the shape of the
20
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Fig. 10. Simulated and measured lag angle Dh(hH) relative to three even rotational loops obtained by applying rotating magnetic inductions having amplitude 0.40 T (upper panel), 1.00 T (center panel) and 1.30 T (lower panel).
closed paths of the magnetization for a set of applied fields, but also the magnetization angle or the phase displacement between M and H. Such quantity, called lag angle, is usually evaluated as function of the magnetic field angle hH. In Fig. 10 the lag angle relative to three even loops (B = 0.40 T, B = 1.00 T, B = 1.30 T) is plotted versus the corresponding magnetic field angle hH. The good agreement between the simulated and the measured values is clearly visible. It is interesting to note that the behavior of the lag angle versus the magnetic field angle shows an interesting feature of this material, which presents an in-plane uniaxial anisotropy for low inductions, although the Goss-texture should be characterized by a cubic anisotropy. In the upper panel of Fig. 10 the maximum delay of the magnetization is found for hH ? 90° when the global magnetization tends to remain aligned with the positive easy semi-axis. For hH immediately greater than 90° the magnetization vector jumps toward the negative easy semi-axis and remains in advance respect to the field. In the interval [180°, +180°] the lag angle curve presents two minima and two maxima points. When the magnetic induction increases, the lag angle tends to be distorted (center panel of Fig. 10) but only for B > 1.10 T the real cubic anisotropy is manifest, as in the loop B = 1.30 T shown in the lower panel of Fig. 10. In this case the maximum delay of M is found approximately around 45°, followed by a slight anticipation due to the secondary energy minima that is found along the y axis for the (1 1 0) in-plane cubic anisotropy. The number of minima and maxima point of the lag angle are doubled. Our analysis is here concluded with a discussions about the magnetic losses. The energy lost per cycle is computed using the equation:
Eloss ¼ Eloss;x þ Eloss;y ¼
1 cT
Z dBx dBy Hx þ Hy dt dt dt T
ð8Þ
Fig. 11. Energy loss per cycle as function of the magnetic induction magnitude calculated by means of both simulated and measured data.
The integral (8) can be solved by evaluating numerically the area of the hysteresis cycles along x and y directions. The rotational losses present the typical behavior shown in Fig. 11, where the comparative analysis between the computed and measured data is reported. For low inductions, the measured field is also very low and the area contribution, that depends on the average lag angle during one cycle, is clearly small. If B < 1.05 T the magnetic anisotropy tends to be similar to the uniaxial one. In this interval the energy loss is an increasing function of the magnetic induction: the major contribution to the energy loss is due to the increment of the fields, while the average lag angle does not vary considerably from one cycle to the other, as
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
one can see in Fig. 10. If B > 1.05 T the cubic anisotropy behavior is manifest and the lag angle tends to be almost constant. Although the maximum misalignment between the magnetic field and the magnetic induction decreases slightly, the average phase displacement during one cycle tends to zero. In this region the lag angle is dominant over the absolute values of the vector fields and the magnetic losses decrease. For this reason it is evident that an accurate reconstruction of the lag angle is necessary in the correct prediction of the magnetic losses. Of course the computed losses associated to the odd rotational loops, used to identify the model, are closer to the measured data, respect to the ones associated to the even rotational loops. However the maximum displacement defined as e = |Eloss,sim Eloss,meas| is below 4 mJ/kg. As final technical comment, let us underline that the simulation of the thirteen rotational loops were carried out in the same computer, equipped with a CPU Intel Core i7 @ 2.20 GHz. The single hysteron model required only few seconds to reconstruct the entire series of rotational magnetization processes.
21
The approach is derived by the classical SW theory applied to the case of cubic magnetic anisotropy defined in the (1 1 0) plane. Vector hysteresis in Goss-textured ferromagnetic material can be accurately modeled at macro-magnetic scale length by means of a single hysteron. The moving modeling strategy defined in this case leads to an univocal calculation of the model parameters, deduced by vector measurements. The comparison with experimental data of GOES indicates that the new hysteron can predict with accuracy the magnetization trajectories, the lag angle and the static losses for circular magnetizations. The hysteron presented seems especially promising for a matching of the Goss-textured ferromagnetic materials modeling with FEM analysis in time domain, because puts together numerical accuracy and strong computation time and memory allocation saving. Scalar magnetization modeling deserves further studies. We are considering the possibility to model scalar hysteresis by means of the reconstruction of elliptic loops with increasing values of their aspect ratio, but this is a future work. Appendix I
5. Conclusions and further Research In this paper we have presented a new approach for the phenomenological modeling of Goss-textured ferromagnetic materials.
In this section we list the MatlabÓ code that we have written to implement the identification function of the single hysteron model.
22
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
Appendix II In this section we list the MatlabÓ code that we have written to implement the magnetization calculation according to the single hysteron model.
E. Cardelli et al. / Journal of Magnetism and Magnetic Materials 432 (2017) 14–23
References [1] C. Appino, C. Ragusa, F. Fiorillo, Can rotational magnetization be theoretically assessed ?, Int J. Appl. Electromagnet. Mech. 44 (2014) 355–370. [2] T. Matsuo, Anisotropic vector hysteresis model using an isotropic vector play model, IEEE Trans. Magn. 46 (8) (2010) 3041–3044 5512912. [3] E. Cardelli, A. Faba, G. Finocchio, B. Azzerboni, Mathematical modelling of magnetic hysteresis in exchange-bias spin valves, IEEE Trans. Magn. 48 (2012) 3367–3370. [4] F. Preisach, Über die magnetische Nachwirkung, Z. Angew. Phys. 94 (1935) 277–302. [5] A. Visintin, Differential models of hysteresis, in: Applied Mathematical Sciences, 111, Springer, Berlin, 1994. [6] E. Della Torre, Magnetic Hysteresis, IEEE Press, Piscataway, NJ, 1999. [7] E. Della Torre, E. Pinzaglia, E. Cardelli, Vector modeling – part I: generalized hysteresis model, Phys. B 372 (2006) 111–114. [8] E. Della Torre, E. Pinzaglia, E. Cardelli, Vector modeling: Part II, ellipsoidal vector hysteresis model – numerical application to a 2-d case, Phys. B 372 (2006) 115–119. [9] L. Dupré, G. Bàn, M. Von Rauch, J. Melkebeek, Relation between the microstructural properties of electrical steels and the Preisach modelling, J. Magn. Magn. Mater. 195 (1999) 233–249. [10] E. Cardelli, E. Della Torre, A. Faba, A general vector hysteresis operator: extension to the 3-D case, IEEE Trans. Magn. 46 (12) (December 2010) 3990– 4000. [11] E. Cardelli, A general hysteresis operator for the modeling of vector fields, IEEE Trans. Magn. 47 (8) (2011) 2056–2067. [12] Z. Szabo, J. Füzi, Implementation and identification of Preisach type hysteresis models with Everett Function in closed form, J. Magn. Magn. Mater. 406 (2016) 251–258. [13] E. Cardelli, A. Faba, Genetic algorithm identification of a H-moving vector hysteresis model, Phys. B 435 (February 2014) 11–15.
23
[14] E. Cardelli, A. Faba, A. Laudani, F. Riganti Fulginei, A. Salvini, A moving approach to magnetic modeling of electrical steels in 2-d, Int. J. Appl. Electromagnet. Mech. 48 (2,3) (2015) 263–270. [15] E. Cardelli, A. Faba, A. Laudani, S. Quondam Antonio, F. Riganti Fulginei, A. Salvini, A moving approach for the Vector Hysteron Model, Phys. B 486 (1) (April 2016) 92–96. [16] E. Cardelli, E. Della Torre, A. Faba, Numerical implementation of the DPC model, IEEE Trans. Magn. 45 (3) (2009) 1186–1189. [17] E.C. Stoner, E.P. Wohlfarth, Philos. Trans. R. Soc. London, Ser. A, Math. Phys. Sci. 240 (826) (1948) 599–642. [18] M. de Campos, F. Sampaio da Silva, E. Perigo, J. de Castro, Stoner–Wohlfarth model for the anisotropic case, J. Magn. Magn. Mater. 345 (2013) 147–152. [19] A. Kuncser, V. Kuncser, Magnetization reversal via a Stoner–Wohlfarth model with bi-dimensional angular distribution of easy axis, J. Magn. Magn. Mater. 395 (2015) 34–40. [20] C. Chang, Micromagnetic studies of coherent rotation with quartic crystalline anisotropy, J. Appl. Phys. 69 (1991) 2431. [21] A. Thiaville, Extensions of the geometric solution of the two dimensional coherent magnetization rotation model, J. Magn. Magn. Mater. 182 (1998) 5– 18. [22] A. Thiaville, Coherent rotation of magnetization in three dimensions. A geometrical approach, Phys. Rev. B 61 (2000) 12221–12232. [23] A. Stancu, L. Spinu, Transverse susceptibility for single-domain particle with cubic anisotropy, J. Magn. Magn. Mater. 266 (2003) 200–206. [24] J. Petrila, A. Stancu, Analytical vector generalization of the classical Stoner– Wohlfarth hysteron, J. Phys.: Condens. Matter 23 (2011) 076002. [25] D. Cimpoesu, L. Stoleriu, A. Stancu, Generalized Stoner–Wohlfarth model accurately describing the switching processes in pseudo-single ferromagnetic particles, J. Appl. Phys. 114 (2013) 223901. [26] Bozorth, Ferromagnetism, IEEE Press, 1993.