Dynamic Modeling of a Generator With Anisotropic Nonlinear Permanent Magnets in Finite Element Method Software

In this article, a method of dynamic modeling of nonlinear permanent magnets (PMs) with recoil lines in 2-D finite element analysis (FEA) software was presented. COMSOL Multiphysics 6.0 FEA software was used in this study. The method is implemented through the variable utilities. The simulation results of a spoke-type synchronous generator for a wind turbine with anisotropic aluminum–nickel–cobalt (Alnico) 5, 8B, and 9 grades were used to exemplify the model and compared. The proposed methodology can be used for the simulation of nonlinear PMs with recoil lines and includes reversible and irreversible losses of magnetization of nonlinear PMs. The effect of the magnetic field from the stator winding on nonlinear PMs during normal operation and short circuits was studied. The modeling results were compared to the model without any demagnetization and a previous study with recoil lines and averaged minimum magnetic flux points. The no-load (NL) voltages were compared before and after a demagnetization. The dynamic model showed considerable demagnetization of Alnico magnets during normal operational and three-phase short circuits. Alnico 5 and 9 showed higher sensitivity to short-circuit currents and the short-circuit currents caused remagnetization of the upper part of the magnet in the opposite direction. The anisotropy of the PM implemented in the model improved the magnetic field simulation inside the magnet and partially protected the magnet from demagnetization by inclined fields. At last, the method was experimentally verified by tests on an iron core.


I. INTRODUCTION
T HE current increase in the demand for energy [1] and the dependence of some countries on the import of energy resources from abroad push the development of new, and improvement of existing renewable technologies [2], [3].
Currently, wind energy and solar energy play the main role in the replacement of fossil fuel-based electricity production [4].Wind energy is a more effective and cheaper option for countries with low solar irradiance and close location to the sea or ocean.Together with the development of the design of wind turbines and electrical generators, the improvement of magnetic materials [5] contributes to the higher efficiency of the machines.
Nowadays neodymium (NdFeB) magnets are considered the best magnets for electrical machines due to the high remanence, energy product, and low and fairly constant relative permeability in the second quadrant before the knee point.The disadvantages of NdFeB magnets are the environmental impact during mining and production of magnets, low Curie temperature, monopolization of the market, and impact on the national security [6].
Discovering new magnets, improving existing magnets, and improvement of the design of electrical machines [7] can help replace NdFeB magnets with more environmentally friendly and more accessible magnets [8], [9].Some of these alternative permanent magnets (PMs) have nonlinear B-H curve Manuscript received 24  in the operating region of electrical machines under normal conditions and it makes them prone to demagnetization due to an external magnetic field [10], [11].One of the PMs with a nonlinear B-H curve is aluminum-nickel-cobalt (Alnico) magnets.The development of Alnico magnets began in the 1930s.Alnico magnets have higher Curie temperature, higher remanence, higher resistivity to corrosion, and lower cost in comparison to NdFeB magnets.Alnico magnets are sometimes used in motors since they have the advantage of easier control of the magnetic field in the airgap [12], [13].Disadvantages of Alnico magnets are easy demagnetization, the nonlinearity of the B-H curve in the operating region of electrical machines under normal conditions, lower energy product, and lower coercivity.There are two main solutions to handle the demagnetization of Alnico magnets in electrical machines: protection from demagnetization [14] or fast remagnetization of the magnets as in memory machines [15].There are many ways to model PMs [16].Nonlinear PMs can be modeled with a main curve interpolation, a linear recoil curve, or a recoil loop.The main B-H loop model which is the default option in COMSOL Multiphysics 1 6.0 finite element analysis (FEA) software for nonlinear isotropic PMs utilizes only the interpolated first and second quadrants of the main loop of the B-H curve.It does not take into account demagnetization and the operating point moves along the main loop of the B-H curve.
Most of the models with demagnetization of PMs focus on demagnetization due to temperature [17], [18], demagnetization of the magnets with constant permeability of the main loop [19], [20], calculated minimum B point at the beginning 1 Registered trademark.This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/ of the simulation [21], [22] etc.The models can be dynamic where the minimum magnetic flux density is calculated at each step of the simulation or in some models it is calculated only at the beginning of the simulation.In some studies, the model details are not presented at all and in other studies, nonlinear PMs were modeled with constant permeability at the operating region and the permeability was chosen to be equal to the recoil line permeability [23], [24].This type of model of Alnico magnets can be found in the COMSOL library.
The recoil loop model requires more computational power.Some magnets, such as Alnico 8, have recoil loops mostly in the third quadrant and the recoil lines in the second quadrant of the B-H curve [25].The PMs of electrical machines most of the time operate in the second quadrant.
This article describes the methodology of the implementation of a nonlinear PM model with recoil lines in finite element method (FEM) software.The generator was simulated at NL, with a normal load, and under short-circuit conditions.The demagnetization of the magnets was observed.The model was tested on a 2-D model of a synchronous generator with Alnico magnets.It was experimentally verified with a 3-D model of an iron core.

II. THEORY A. Permanent Magnets
The magnetic flux density of PMs can be formulated as where B r is the remanent magnetic flux density, µ 0 is the permeability of free space, µ r is the relative permeability, and H is the internal magnetic field.It can alternatively be rewritten as where M is magnetization.Conventional PMs, such as NdFeB and ferrite, have a constant relative permeability µ r and thereby a linear behavior if the opening point does not reach the knee point.
For the magnets with a sharp knee of the B-H curve, such as NdFeB and ferrites, the knee point is the point at the second and fourth quadrants of the B-H curve where permeability changes drastically.Nonlinear PMs, such as Alnico and MnAl, do not have a drastic change in the permeability along the B-H curve.The knee point of nonlinear PMs is the point where local permeability is larger than recoil permeability and the remanence of the partially demagnetized magnet cannot be higher than the remanence of the fully magnetized magnet.
Parameters of Alnico magnets from Arnold Magnetic Technologies are shown in Table I.The relative permeability of the hard direction in the table was approximated from [26].

B. Linear Recoil Curve
The relative permeability of nonlinear PM changes during operation and depends on the applied magnetic field and the history of the magnet.Nonlinear PMs have recoil loops of the B-H curve that can be simplified to recoil lines with a constant slope (see Fig. 1) due to a comparably small difference between a minor loop and recoil line values of  magnetic flux density in the second quadrant, according to the measurement results [25].The operating point follows the recoil line trajectory only if the magnetic field density B is lower than the point on the B-H curve where µ pm > µ recoil and higher than the minimum B field that the magnet experienced after being fully magnetized.Otherwise, it will be on the main B-H curve.

III. MAGNET MODEL
The methodology of modeling nonlinear PMs in FEM-based software presented in this article is using the linear recoil curve model implemented through the expression operator and state variables that describe the behavior of Alnico magnet according to Fig. 2 and (2).
The synchronous generator with Alnico magnets was modeled in the 2-D rotating machinery component of COMSOL Multiphysics 6.0.The demagnetization is calculated during the simulation for each point of the mesh.The Alnico magnets are modeled as anisotropic PMs.The easy (preferred) direction is implemented as a lookup table with interpolation for the main curve and a constant permeability for recoil lines.The nonlinear PM from the standard COMSOL library with modifications was used for the preferred direction and data for the lookup table was taken from the Arnold Magnetics website by digitalizing the second quadrant of the B-H curve plot from the datasheet (see Fig. 3).The B-H curve of Flowchart of the PM 2-D model with recoil lines.x refers to the magnetization direction and y is transverse (perpendicular to the easy direction) to the magnetization direction, and min refers to the minimum value.the hard (transverse) direction was modeled with a constant permeability (see Table I).It is assumed that the PM will not be magnetized in the y-direction since most of the magnet will operate in the second quadrant of the B-H curve in the easy (preferred) x-direction.The values for the permeability in the hard (transverse) direction were approximated from [26] since there is no hard (transverse) direction B-H curve for given grades of Alnico magnets by the supplier.
Modeling starts with the input of initial values, parameters, and constants (see Fig. 2).The constitutive B-H relation is implemented as a magnetization model with the relation from ( 2).The magnetic field is calculated with state values of magnetic flux density B x and B y from the model.

IV. GENERATOR
A spoke-type generator with the magnetization of PMs in a tangential direction was modeled since that topology has better demagnetization protection in comparison to a surfacemounted PMs generator [7].The geometry and mesh of the generator are presented in Fig. 4. The silicon steel NGO 50PN400 was used for the stator and rotor and its electrical conductivity was set to 0 S/m due to lamination.Setting the electrical conductivity to the value from the datasheet will give inaccurate results for 2-D models with high eddy currents and additional B fields.
The main parameters of the synchronous generator designed for a wind turbine are shown in Table II.The nominal rotational speed is 500 r/min, and the nominal electrical frequency is 50 Hz, but both can vary during the operation since rotational speed depends on the wind speed and load.Also, the generator is connected to the turbine through a gearbox and the grid through a power converter, so the output voltage can have a wide range of voltage amplitude and electrical frequency levels.

V. SIMULATION SETUP
A 2-D model of the generator was used since most of the rotating machines have plane symmetry and the field distribution is the same in parallel cross sections in the axial direction.It simplifies the modeling and decreases the modeling time but does not take into account the end effect.Only one pole pitch of the machine or (1/12)th was modeled to reduce the simulation time.
The no-load (NL) case was simulated with 3 k and the three-phase short-circuit case was simulated with 0.1 µ resistance per phase.The transition of the load from normal load to the short circuit is 0.1 µs.The other short-circuit types have not been simulated since the three-phase short circuit has been shown to have the highest demagnetizing effect on the magnet [27] and the focus of this article is to study the PM model.The load angle in the simulation of the normal load case was set to 25 • for stability reasons and was confirmed through the analytical calculation and graphical simulation results.The simulation starts with 0 rotational speed and ramps to the nominal 500 r/min.This ramping of the rotational speed helps with the convergence of the model at the beginning of the simulation.The stator of the generator was rotated instead of the rotor in order to simplify the modeling of the magnets.It should not be done if a mechanical analysis is included in the modeling results.The maximum magnetic flux density B in the airgap is measured in the middle of the airgap and the average B in the magnet was measured to observe the demagnetization.
Two cases were modeled: 1) NL-normal-NL (from 0.1 to 0.2 s NL, from 0.2 to 0.3 s normal load, and from 0.3 to 0.4 s NL) and 2) NL-SC-NL-normal (from 0.1 to 0.2 s NL, from 0.2 to 0.3 s short circuit, from 0.3 to 0.35 s NL, and from 0.35 to 0.4 s normal load).The time-dependent study was simulated with the help of the fully coupled attribute node (damped Newton's method) with constant damping.The fully coupled attribute node with the automatic highly nonlinear damping requires less computational time than the constant damping but might have problems with convergence.The relative tolerance of the model is 0.01 as suggested for FEM modeling.
The anisotropic main B-H loop model without recoil lines (model B) was developed together with the anisotropic main B-H loop model with recoil lines (model A) to compare the simulation results of those two models.The model without recoil lines is an improved version of the built-in model (model C).The built-in model (model C) of the isotropic magnet utilizes only the first and second quadrants of the B-H curve and the rest of the B-H curve is extrapolated.

VI. SIMULATION RESULTS AND DISCUSSION
The simulation results for the synchronous generator with Alnico magnets are presented below in Tables III and IV.There is a large difference between the anisotropic main B-H loop model without recoil lines (model B) and the anisotropic main B-H loop model with recoil lines (model A) even after a normal load operation.Both models A and B were developed from the default isotropic model of the nonlinear PM without recoil lines from the COMSOL library (model C).The NL voltage, normal load terminal voltage, and armature current of those two models also slightly differ.It can be explained by small variations of the magnetic field during rotation due to  a change the reluctance of the magnetic circuit.The results for the recoil line model are slightly lower, as expected.
The NL voltages in models B and C go back to their initial values after normal load and short circuit since the built-in nonlinear PM model in COMSOL does not model demagnetization correctly (model C).It does not experience permanent demagnetization.There is a tutorial on COMSOL's website on the development of a demagnetization model of the PM but it is not dynamic.The demagnetization in that tutorial is calculated at the beginning of the time-dependent study and the magnet itself is modeled with the recoil permeability only.
The analytically calculated NL voltage of the generator with Alnico 8B magnets using the generator equation is 715.2V.The simulation result for the NL case after the normal load is 672.2V which is 6% lower.It can be assumed that this difference is due to magnetic flux leakage, geometrical limitations of the generator, and nonsinusoidal magnetic flux density.
The generator was simulated with PMs to check if the simulated short-circuit current is reasonable since the subtransient inductance was higher than expected.For a machine with exactly the same geometry but linear PMs and with a relative permeability of 1, the subtransient inductance was 14.4%.The subtransient reactance depends on the B-H curve shape (relative permeability) of the magnet and operating point.
The model was tested with different transition times from the normal load to the short circuit.No significant change in the subtransient current amplitude was noticed when the transition time was changed from 10 to 0.1 µs.
Alnico 5 and 9 demagnetized more after a short circuit in comparison to Alnico 8B.The demagnetization of the magnet depends on the permeability of the magnet to the left of the operating point on the B-H curve and its coercivity.
The demagnetization of Alnico 8B and 9 magnets can be seen and confirmed by the simulation results in Figs.5-8.Also, the simulation results confirm that Alnico magnets experience larger demagnetization after a short circuit than after normal operation.
The average magnetic flux density of the Alnico magnet after the short circuit is higher than the minimum magnetic density during the short circuit.It can be assumed that the average operating point of the magnet is located on the recoil line, as expected.The magnet demagnetizes closer to the upper corners and clamp edges.The demagnetization is not symmetrical after an operation of the generator under normal load conditions due to the load angle.It should be noticed that the magnetic flux direction inside the magnet is different from the isotropic main B-H loop with averaged recoil lines model (model D) from the previous study [22] (see Figs. [9][10][11].It can be explained by the larger difference between neighboring parts of the magnets in [22] which changes the direction of the magnetic field.The model could be more accurate if the magnet in [22] is divided into smaller parts.The anisotropy of the PM protects the magnet from demagnetization in the easy direction.It also could be a reason for a more uniform direction of the magnetic field inside the magnet in comparison to [22].
After a short circuit, the upper part of the magnet experiences larger demagnetization than the rest of the magnet because of the close location to the stator winding.The upper part of Alnico 5 and 9 magnets was remagnetized in the opposite direction due to a combination of high remanence and low coercivity.Also, it can be noticed that the y-direction component of the B field in the remagnetized part is almost 0 after the short circuit even if it was significant during the short circuit.The explanation is that the hard direction of the PM was modeled as a soft magnetic material.

A. Experimental Setup
The demagnetization of nine pieces of Alnico 8 LNGT40 magnets was tested in an iron core with a coil with 290 turns (see Figs. 12 and 13).The iron core is made of laminated sheets of M700-100A.The size of each magnet is 10 × 10 × 10 mm.The magnetic flux density was measured in the middle of the airgap (4 mm) with the help of an FW Bell 5180 Gauss/Tesla meter.The specified accuracy for values over 30 mT is ±(0.6% + 3 counts), where one count is 0.1 mT below 300 and 1 mT for larger values.Below 30 mT specified  accuracy is ±(0.8% + 0.04 mT).The B-H curve and recoil lines of Alnico 8 LNGT40 dimensioned 3 × 3 × 3 mm were measured in vibrating sample magnetometer (VSM) [25].be significant in this case.This model was developed in the magnetic fields (MF) module of COMSOL Multiphysics 6.1.It is suggested to use 6.1 and a higher version of COMSOL for the time-dependent simulation models with coils and recoil lines.

B. Simulation and Experimental Results
The measurement and simulation results (model A) of demagnetization of the Alnico 8 LNGT40 magnets in an  iron core can be seen in Fig. 14.The airgap in the model was adjusted to 3.4 mm and the magnetic flux density was decreased to 2% in the B-H curve of the silicone steel of the iron core in order to make B field in the airgap match at 0 A. The difference between the simulation and measurement results can be explained by the lamination of the iron core steel sheets, the difference of the B-H curves of the rolling and transverse directions, and a small remanence of the iron core.Also cutting the iron core might have affected its edges.
The load line of the magnet is in the second quadrant of the B-H curve before, during, and after the experiment since the magnetic flux density does not reach negative values.The magnetic flux density after an application of a demagnetizing field will not reach the initial value when the applied field is removed (see Fig. 14).This confirms the partial demagnetization of the magnets.The operating point of the magnet follows the recoil loop instead of the main B-H loop.The slope of the recoil line of the simulation and measurement results slightly differ.It can be assumed that the operating point in the simulation is located higher on the B-H curves of PMs and steel in comparison to reality.As a result, the permeability in the simulation is lower than in the measurements.Another difference of the current model is an implementation of the magnetic anisotropy that helps the magnet to maintain the magnetization in the easy direction.Model A was verified through experimental tests on an iron core.The model of the steel needs to be improved in order to get perfectly matching experimental and simulation results.
The simulation results confirm that drastic demagnetization of the permanent Alnico magnets occurs during short circuits and it depends on the subtransient currents.It can be concluded that the B-H curve shape of the magnet and load impedance change rate, affect the demagnetization of the magnets.The B-H curve of the magnet affects the maximum current during a short circuit or, in other words, affects the subtransient reactance of the generator.
Alnico 5 is sensitive to short circuits and loses most of the magnetization due to subtransient currents.A synchronous generator with Alnico 9 has a higher terminal voltage but lower rated power than a PM synchronous generator (PMSG) with Alnico 8B magnets.
In conclusion, for Alnico magnets to be considered for use in generators, the magnetic circuit needs to be carefully designed to decrease the demagnetizing effects of short-circuit currents presented here.

Fig. 1 .
Fig. 1.Hysteresis curves of a nonlinear PM Alnico 8B with a recoil line.

Fig. 2 .
Fig. 2.Flowchart of the PM 2-D model with recoil lines.x refers to the magnetization direction and y is transverse (perpendicular to the easy direction) to the magnetization direction, and min refers to the minimum value.

Fig. 4 .
Fig. 4. Model of the spoke-type synchronous generator with Alnico magnets.(a) Geometry of the machine: white-air; light gray-laminated silicon steel NGO 50PN400; dark gray-nonmagnetic steel; purple-PM; and red, yellow, and blue-copper conductors (three phases).(b) Mesh-free triangular.

Fig. 9 .
Fig. 9. Alnico 5 x-direction component of the magnetic flux density B x in T. (a) NL R = 3 k and t = 0.13 s.(b) NL after normal load 35 and t = 0.38 s.(c) Short circuit 0.1 µ and t = 0.2095 s.(d) NL R = 3 k after short circuit and t = 0.43 s.Arrows show magnetic flux density B-direction.

Fig. 10 .
Fig. 10.Alnico 8B x-direction component of the magnetic flux density B x in T. (a) NL R = 3 k and t = 0.13 s.(b) NL after normal load 20 and t = 0.38 s.(c) Short circuit µ and t = 0.2095 s.(d) NL R = 3 k after short circuit and t = 0.43 s.Arrows show magnetic flux density B-direction.

Fig. 11 .
Fig. 11.Alnico 9 x-direction component of the magnetic flux density B x in T. (a) NL R = 3 k and t = 0.13 s.(b) NL after normal load 25 and t = 0.38 s.(c) Short circuit 0.1 µ and t = 0.2091 s.(d) NL R = 3 k after short circuit and t = 0.43 s.Arrows show magnetic flux density B-direction.

Fig. 12 .
Fig. 12. Iron core geometry.The coil is installed on the left side of the iron core.The depth is 30 mm everywhere.

Fig. 14 .
Fig. 14.Experimental and simulation results (model A) of demagnetization of the Alnico 8 LNGT40 magnets in the iron core.Magnetic flux density B in the airgap was measured with different coil currents.
VIII.CONCLUSIONThe presented novel dynamic linear recoil line model of demagnetization of anisotropic nonlinear magnets (model A) provides more accurate simulation results in comparison to the other three models discussed in this article (models B, C, and D).Model B refers to the anisotropic main B-H loop model without recoil lines, model C is the default model in COMSOL, and model D is the model with averaged recoil line model from the previous study.The difference can be explained by the dynamic tracking of demagnetization and as a result lower subtransient and transient short-circuit current amplitudes.

TABLE I PARAMETERS
OF ALNICO MAGNETS