Modeling of history‐dependent magnetization in the finite element method on the example of a postassembly rotor magnetizer

Demagnetization characteristics of permanent magnets play a central role for the design and optimization process of permanent magnet electrical machines and magnetic actuators. Currently, practical issues, such as loss of performance of electromagnetic energy conversion devices, due to demagnetization effects in permanent magnets or due to reduced servo‐dynamics in magnetic actuators, cannot be addressed by using conventional material models or are investigated by using rough approximations valid in restricted operation regions, if at all. This paper strives to integrate demagnetizing characteristics dependent on previous magnetizing field strength in the finite element method. The discussed examples are a magnetizer of a single permanent magnet and a spoke‐type permanent magnet rotor.


| Classical approaches
The relative permeability of fully magnetized rare-earth permanent magnets is close to one in the first and second quadrant of B(H) characteristic, which is why there occurs almost no irreversibility. 4 The simplest models correspond to an analytical linear or nonlinear equation. A single hyperbolic tangent function already forms a good approximation of the scalar magnetic polarization in preferred direction over the magnetic field strength. In Zhou et al, 5 the scalar demagnetization curve in preferred direction, J in T, of fully magnetized permanent magnets is illustrated by two hyperbolic tangent functions as follows: The coefficients b 0,1 describe the magnetic remanence and h 0,1 the occurring knee field strength. Due to Ampere's law, a coercive field strength can be recalculated into an equivalent source current, which is used in standard finite element method (FEM) approaches to model nonlinear permanent magnet behavior in third and fourth quadrants as follows: Taking this additional source into account, on the excitation side of Equation (4), a remanent magnetic polarization is represented. The equivalent current source moves the magnetization curve inside the permanent magnet to the origin, which produces a symmetrical scalar magnetization curve in the first and third quadrants, equivalent to the modeling of soft magnetic materials. 6 This single characteristic is not variable during the simulation, and therefore, only represents a reversible demagnetization curve without dependency from previous magnetization. In addition, various hysteresis models have been transferred to permanent magnets or developed especially for them, such as the ferromagnetic positive-feedback theory. 7 Simple mathematical models are pice-wise-linear 8 and transplantation type, such as Tellinen 9 or Zirka. 10 Popular models are the Jiles-Atherton 11 or Preisach 12 and Play. 13 Unfortunately, most hysteresis models only produce symmetric curves, which are not capable of representing inner virgin magnetization loops, as the ones shown in Figure 1A. Therefore, the focus of this paper is on the nonlinear magnetization, field-dependent demagnetization characteristic, which later still can be integrated into a mathematical hysteresis model, such as the Tellinen, which is related to Equation 1.

| Continuous history-dependent (de-)magnetization surface
Permanent magnets are operated as a field source in the second quadrant of B(H) characteristic. 4 For that, the virgin magnetization curve must be run through to saturation first. Depending on the maximum field strength experienced of the permanent magnet, different demagnetization characteristics occur, as illustrated in Figure 1A. Based on a common virgin curve in the first quadrant, depending on the maximum field strength seen by the hard magnetic material during magnetization process, the demagnetization curves split with nonlinear courses in the second quadrant. While for high magnetizing field strength, the demagnetization curves are almost linear up to the knee point; the demagnetization curves at low magnetizing field strength have saddled and turning points, Figure 1A. According to the previous elaborations, it can be concluded that the previously presented analytical model with only one demagnetization curve is not suitable for the description of a magnetization process of permanent magnet material.
In this paper, a history-dependent permanent magnet material model is proposed, which stores the magnetization field strength occurring in each element during the transient magnetization process. First, while increasing impulse magnetization, the virgin magnetization curve is utilized. Second, as soon as the magnetization process reaches its peak, a demagnetizing curve is selected in each element, which is related to the maximum magnetization field strength. Subsequently, the material equation is solved as follows: For simplicity, it is assumed that a remanent magnetic polarization only appears in the preferred magnetization direction. All other directions have a linear characteristic with relative permeability close to one. The necessary material data for calculation, as shown in Figure 1B, is derived by interpolation of the measured material data from the Vacuumschmelze catalog, as shown in Figure 1A. The interpolation surface is directly derived from measurements, which may contain measurement deviations, causing nonphysical behavior. Thus, a regularization algorithm has to be applied to avoid numerical instability. A regularization method can be the application of an analytical equation, such as Equation (1) to the scalar magnetization in preferred direction J(H) and a sigmoid function to the corresponding J(H mag ). In this paper, a two-dimensional smoothing spline is directly applied to the measurements, which result in the surface shown in Figure 1B. The flow chart for the determination of the magnetic polarization is depicted in Figure 1C.

| FINITE ELEMENT METHOD
In the considered case of a rotor magnetizer, the general Maxwell equations can be simplified and merged into a quasi-static magnetic vector potential formulation to calculate the transient electromagnetic fields, as in Equation 4. The equation is derived from the combination of Ampere's law with the neglection of displacement current density 15 and the magnetic material law. The reluctivity ν can be a scalar, a scalar function, or a tensor of second rank. The time derivative of the magnetic vector potential is calculated by an implicit Euler scheme. The previously defined equation must be fulfilled in the entire solution area in compliance with all constraints, such as sub-areas with different materials, boundary conditions, and impressed excitations. Complex considered geometries have to be solved either in two or three dimensions. The magnetic vector potential A is a function of space and time and is discretized over geometry by means of finite elements. The magnetic vector potential A ! represents an approximation to the true Solution A ! * . The absolute error is called residual and is calculated by Equation (5).
The material behavior depicted in the last chapter needs to be considered for a sufficient solution, which is why a suitable iteration method for nonlinearity needs to be applied. 16 Because of its fast convergence speed, the Newton method is applied. 17 The Newton method calculates the correction Δ A ! by linearizing the material behavior in the current working point, represented by the Jacobian matrix P. The error in Equation (5) is minimized by the following equation: The derivative of the resudial in terms of the magnetic vector potential ∂R ∂A is called Jacobi matrix. 18 The evaluation of the Jacobi matrix, which is referred to as P, requires the correct evaluation of the derivatives 19 The entries within P k are defined according to the different material regions. With a full consideration of the material properties, especially in the case of anisotropic materials, all entries must be evaluated. In this work, the assumption is made that the nonpreferred direction of permanent magnet regions behaves like air and is therefore not coupled to the preferred direction, which leads to a diagonal Jacobian matrix. Due to these restrictions, the Newton method does not evaluate the exact Jacobi matrix, but an approximation, which is why this modified method is a quasi Newton.

| Impulse magnetizer for single magnet
The geometry of an impulse magnetizer, taken from Przybylski et al, 20 for a single magnet is shown in Figure 2A  The magnetizing circuit, as shown in Figure 3A, consists of an interconnection of a capacity C1 and a magnetizing coil L1. The capacity is initially charged and then shorted over the coil via a thyristor Thy1. A corresponding pulse current is shown in Figure 3B. The current amplitude is 1200 A multiplied by 440 winding turns. The simulation time is divided into 500 equidistant time steps. The capacity must be big enough to provide the necessary magnetization energy and the additional (Ohmic) losses. The requirement for the magnetizing coil is to provide the required field distribution in the most homogeneous way possible, despite inhomogeneities due to geometry and eddy currents.
The 2-D model in Cartesian coordinates is an approximation of the axial symmetric geometry and causes deviations in local field and eddy currents. Nevertheless, the 2-D and 3-D simulations are in good agreement, see Figures 4 and 5. The major effect is the magnetizing field-dependent demagnetization in addition to a minor effect due to local eddy currents. The influence of the magnetizing field strength-dependent demagnetization is demonstrated by simulations at scaled magnetizing current I peak from 20% up to 100%. The values for the magnetic flux density and magnetic field strength are extracted at the points K (1, center; 3, edge; and 2, in-between) inside the permanent magnet in preferred magnetization direction; the points P are the corresponding ones at the surface of the permanent magnet. As an example, the magnetic flux for the impulse current factor of 0.7 is shown at three points inside the permanent magnet over the full simulation time in Figure 4B and for the peak, as well as the end of impulse current over the geometry in Figure 4C and 4D, respectively. The measurements in Przybylski et al 20 correspond to the simulation in this paper, with the difference of another NdFeB material (VACODYM 362/140). A validation of the method concerning the same material is given in Bavendiek et al. 21 At the first glance of Figure 4A, it is contradictory that although the outer field acts only in one direction, the remanent magnetization of the magnet changes its polarity from center (P1) to edge (P3). The reason for that is the inhomogeneous magnetizing field strength during transient impulse magnetization, increased by local eddy currents, causing different demagnetization resistances. That results in self-degaussing of the strong magnetized edges through the weak magnetized center of the permanent magnet. This inhomogeneity of magnetic flux density is particularly visible in the flux distributions of the 3-D simulations in Figure 5. The permanent magnet magnetizes with its virgin curve at increasing pulse current. The magnetic flux density reaches its maximum almost at the peak value of the pulse current, Figure 5A and 5B. The corresponding maximum magnetization field strength determines the demagnetization curve at descending pulse current. At the end of the pulse current, only the remanent magnetization remains, Figure 5C and 5D. However, the inhomogeneity then increases due to different demagnetization resistance. This effect leads to the fact that the local remanence considerably deviates from the measured mean magnetization measured by a flux meter or hysterograph. The goal is to ensure the maximum demagnetizing field strength in all regions of the magnet with minimal energy.

| Postassembly magnetizer for spoke-type rotor
This section deals with the magnetizing process of a spoke-type permanent magnet rotor. To reduce the self-degaussing and to simplify the installation of permanent magnets in the rotors, the magnetization is postassembly. This simplifies the installation of the magnets as they are inserted into pockets in the rotor or glued to the surface. In the literature, 22,23 a magnetic circuit model is applied for modeling the magnetization process. In comparison, in this paper, the approach described above is applied to a transient FEM simulation.
The sample topology of a spoke-type rotor of a permanent magnet synchronous machine is shown in Figure 6. Despite the geometry is not corresponding to any real machine in terms of its dimensions, all the fundamental effects can be discussed here. Due to symmetry, the simulation is carried out by using only a quarter model with antiperiodic boundary conditions assigned to the limiting nodes that lie on the margin lines at ±45ř. This links the nodes so that the flow on those lines has a reverse sign. This corresponds to the opposing field of the next poles of the geometry at ±90ř. The inner rotor radius is 40 mm, and the outer one 100 mm. The two magnetization coils have a radius of 10 mm and This leads to a varying thickness of the surrounding electrical steel bridges, which are 2.5 mm at the horizontal symmetry axis. The permanent magnet material is VACODYM 362/140 and the electrical steel is an arbitrary one with 2.0 T saturation polarization. In order to magnetize the magnets, the coils are excited with the current illustrated in Figure 3 b. The current was produced in 20 by an impulse magnetizer to magnetize NdFeB magnets in an air-cored coil. The amplitude of the impulse current at 100% is 1200 A and is varied by multiplication factors mag mul to achieve different magnetization sates. The excitation and therefore the transient simulation is divided into 500 equidistant time steps to be sufficiently accurate. It is assumed that the laminated electrical steel of the rotor is nonconductive. The electrical conductivity of the magnet is assumed to 70 000 S/m. In Figure 7A-F, the local flux density distribution of the spoke-type rotor are depicted for different amplitudes of impulse magnetization currents, I p eak * mag mul . At the beginning of the transient simulation, the magnetic flux concentrates at the small steel bridges surrounding the magnet, because of their lower magnetic resistance. This is due to the higher magnetic resistance of the magnetizing curve of the magnet as well as induced eddy currents inside the magnet, which counteract the penetrating field. As a result, the magnet begins to magnetize only after the rotor bridges are saturated. Since the eddy currents strongly counteract the penetration of the field, the magnetic flux has to pass the saturated rotor bridges. Because of the smaller cross-section, the right-hand bar first goes into saturation and the magnetic flux density enters from the right side into the magnet. This results in an inhomogeneous field distribution inside the magnet. Some time steps later, also the left bridge is saturated and the magnetic flux density enters from the left side of the magnet. The different change of the magnetic flux penetration at both sides creates eddy currents that are distributed unsymmetrically to the center of the magnet. To sum up the geometry, determining the reluctance network has a strong influence on the flux density and eddy-current distribution. The flux caused by the inner coil encloses a smaller area compared with the outer coil, whereby the flux density of the inner rotor surface is higher than on the outer rotor surface, which spreads back to the magnets due to the curvature of the rotor. For the lowest current amplitude, there is a small remanence on the right side of the magnet, but the magnet is self-degaussing over the two steel bridges. In all simulations, the center of the magnetization is slightly shifted to the right, which is recognizable by the distribution of the equipotential lines of the magnetic vector potential. The distance between the equipotential lines decreases to the edges of the magnet, corresponding to an increasing flux density at the edges. The shift of the center of magnetization from the center of the magnet decreases by increasing the amplitude of impulse magnetization current. A homogeneous magnetization is almost reached at the highest current amplitude. The magnetized magnet locally demagnetizes at the edges over the steel bridges. The impact of eddy currents is visualized in Figure 8, which sets the simulations with and without eddy currents against each other. It can be clearly seen that the eddy currents reduce the magnetic field in the center and increase it at the edges of the permanent magnet. A purely static consideration of the maximum magnetization field can therefore lead to an erroneous assessment of the magnetization requirement. The mentioned rule of thumb, concerning at least three or up to five times of the coercivity for full magnetization, can now be tailored exactly to the material and the geometry.

| CONCLUSIONS
The magnetizing field strength-dependent demagnetization behavior of permanent magnets in the first and second quadrants is modeled by the history-dependent model. For this purpose, the maximum experienced field strength is stored in each element. Based on an anisotropic behavior, only a magnetic polarization is considered in the preferred direction and is neglected in the other directions. To handle the nonlinear material behavior, a quasi Newton method for the FEM is introduced. Eddy currents in the material are considered and their influence on the local distribution of magnetization is investigated. In the magnetization process of the spoke-type rotor of a permanent magnet synchronous machine, eddy currents lead to inhomogeneous magnetization as long as the magnetizing field does not reach a necessary saturation threshold. An inhomogeneous magnetization by eddy currents during magnetization is amplified at demagnetization by the different demagnetization characteristics. An extension of the model around return lines and remagnetization is conceivable. Based on this, investigations on the magnetization of rotor magnets by multiple magnetization can be done. 3 How to cite this article: Bavendiek G, Müller F, Steentjes S, Hameyer K. Modeling of history-dependent magnetization in the finite element method on the example of a postassembly rotor magnetizer. Int J Numer Model. 2020;33:e2674. https://doi.org/10.1002/jnm.2674