Magnetic transmission gear finite element simulation with iron pole hysteresis

Mattia Filippini; Piergiorgio Alotto; Gregor Glehn; Kay Hameyer

doi:10.1515/phys-2018-0017

Article Open Access

Magnetic transmission gear finite element simulation with iron pole hysteresis

Mattia Filippini , Piergiorgio Alotto , Gregor Glehn and Kay Hameyer

Published/Copyright: April 18, 2018

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 16 Issue 1

Abstract

Ferromagnetic poles in a magnetic transmission gear require particular attention during their design process. Usually, during the numerical simulation of these devices the effects of hysteresis for loss estimation are neglected and considered only during post-processing calculations. Since the literature lacks hysteresis models, this paper adopts a homogenized hysteretic model able to include eddy current and hysteresis losses in 2D laminated materials for iron poles. In this article the results related to the hysteresis in a magnetic gear are presented and compared to the non-hysteretic approach.

Keywords: Magnetic gear; hysteresis models; homogenization

PACS: 75.60.-d; 84.50.+d; 41.20.-q; 02.70.Dh; 85.70.-w

1 Introduction

Magnetic transmission gears are gaining particular attention, since they may offer performances similar to those of conventional gearboxes with the advantages of lower maintenance and losses. Several topologies have been proposed as shown in [1] and depending on the application the best candidate could be different.

In magnetic gears, the magnetic fields produced by the inner and the outer rotor are modulated by the action of ferromagnetic poles. These components play a crucial role on overall device behavior. Several papers in the literature discuss analytical models for magnetic field computation in the gears. For example, in [2] a simplified analytical approach is used within an optimization loop providing an approximated first design of an optimal magnetic gear. To achieve a higher accuracy, the detailed magnetic gear design should be finite element based. In particular, iron poles require an accurate model to assess losses and torque. However, the problem results numerically challenging because of the rotational nature of the magnetic flux density and the nonlinear hysteretic behavior in the ferromagnetic poles.

Li et al. in [3] discuss the role of hysteresis on the torque waveforms of a permanent magnet machine, while in [4] the intrinsic dependence between eddy currents and hysteresis is highlighted: hysteresis should be included for accurate machine modelling.

This paper analyzes the hysteresis effect on a coaxial magnetic gear using a homogenized model embedded inside a finite element model. A simple post processing implementation based on the non-hysteretic FEM results (referred as open loop) is discussed and compared to a direct implementation (referred in the paper as closed loop implementation).

2 FEM implementation

To overcome convergence problems of the simulation due to magnetic material models, a differential reluctivity tensor can be applied to solve the finite element formulation of the magnetic vector potential [5].

Using the 2D A formulation, the magnetic flux density B is defined as B = ∇ × A and the magnetic vector potential A = (0, 0, A_z) is discretized through linear piecewise functions:

Az(x,y,t)=∑j=1Neαj(t)βj(x,y)(1)

where Ne is the number of nodes, α_j(t) is the nodal value of the vector potential z component, β_j(x, y) is the shape function. The vector associated to the shape function is ω_j = (0, 0, β_j) since in the 2D formulation only the z component of the magnetic vector potential is not zero.

Assuming negligible eddy currents in the permanent magnets, the weighted residual approach is applied on Ampere’s law, leading to the weak formulation:

∫ΩH⋅∇×ωi+∮∂τH×ωidτ=∫ΩsJs⋅ωidΩs(2)

where J_s is the source current in a subspace Ω_s of the entire domain Ω. The closed integral on the boundary is equal to zero due to homogeneous Neumann or Dirichlet boundary conditions. To solve (2) in the time domain, a time-stepping technique is applied: this is due to the time dependencies of the hysteretic materials. If A(t_n) is a given state of the magnetic problem the state at the next time instant t_n+1 = t_n + Δt is calculated using iterative Newton-Raphson (NR) method. For each NR iteration A^k = A^k−1 + ΔA^k the increment ΔA^k must be calculated. Therefore Eq. (2) is linearized around A^k−1. This linearization is obtained deriving the equation with respect to α_j, which can be achieved through the differential reluctivity:

dHdαj=dHdB⋅∇×ωj=νd∇×ωj(3)

Substituting (3) in (2), Ampere’s law becomes:

∑j=1NeΔαjk∫Ω(νd⋅∇×ωj)⋅(∇×ωi)dΩ==∫ΩsJ(tn+1)⋅ωidΩ−∫ΩsHk−1⋅∇ωidΩ(4)

In the discrete time-stepping scheme the differential reluctivity can be expressed as νd=ΔHΔB=ΔH⋅ΔBΔB⋅ΔH with ΔH = H^k+1(t_n+1) – H(t_n) and ΔB = B^k+1(t_n+1) – B(t_n).

In all elements where hysteresis is considered, the homogenized parametric algebraic model (PAM) described in [6] is adopted to include eddy currents and hysteresis:

H(B,B˙,pk)=(p0+p1|B|2p2)⋅B+p3B˙+p4B˙p52+|B˙|2(5)

Ḃ is the time derivative of the magnetic flux density and the parameters p₀ – p₅ are material constants that have been found through the identification procedure in [6]. In particular, the parameters p₀, p₁ and p₂ are related to the anhysteretic magnetization curve, p₃ is related to eddy currents in the laminated sheets, p₄ and p₅ are linked to the hysteresis phenomena. Since the term B^2·p2 is not asymptotic, Eq. (5) has been applied below the saturation flux density |B| = B_s, while above saturation the BH curve is assumed to be linear with a slope equal to the vacuum permeability μ₀. B_s is computed as:

Bs=1μ0−p0p1⋅(2p2+1)2p2(6)

Since the waveforms in the magnetic gear are sinusoidal in first approximation, the set of parameters chosen from [6] is the one reported in Table 1.

Table 1

List of parameters for the PAM model retrieved by [6]

Po	P1	P2	P3	P4	P5	B_s
95.9	0.29	11.4	0.041	28.6	8.03	1.67 T

Eq. (4) is computed elementwise for the hysteretic regions providing Helemk=f(Belemk(tn+1),Belem(tn)) which is subsequently substituted in the differential reluctivity expression. In the proposed strategy, (4) is directly implemented in the FEM model, constituting a closed loop integration of the PAM material model in the FEM algorithm. The open loop implementation consists of the simple implementation of (5) as post processing, while the nonlinear BH curve adopted for the FEM simulation is the one extracted by the anhysteretic part in (5).

3 Magnetic gear test case

Figure 1 shows the test case geometry (1/4th of the entire model) and dimensions. The iron poles are modelled using the hysteresis model, while for the yoke area the classical nonlinear BH curve is adopted. The rotational speed of the inner rotor is set to v_in = 955 rpm and to v_in = 3180 rpm, thus the outer rotor speed is v_out = 545.7 rpm and v_out = 1817 rpm respectively. The magnets use a linear model with magnetic remanence Br =1.2 T on both inner and outer rotors and unitary relative permeability. As test case geometry, a magnetic transmission gear with low fractional gear ratio and a high number of inner pole pairs (Pi = 4) is adopted, hence the lowest order of harmonics of the cogging torque on the inner and outer rotor are h₁ = 154 and h₂ = 286 respectively according to [7]. These harmonics are due to the combined interaction between permanent magnets magneto-motive force and reluctances due to iron poles. In this paper, the open and closed loop application of Eq. (5) are compared: the resulting loss magnitude’s order is validated through the dynamic version of the typical loss separation method applied for the steel M330-35HS with density δ = 7650 kg/m³:

Figure 1

Magnetic gear test geom1etry and query points P1, P2, P3.

PSt=khystfBα+keddy(dBdt)2+kexc(dBdt)1.5(7)

where f is the frequency and the material parameters calculated through the fitting are: k_hyst = 0.0194 Wm⁻³T^−αs, k_eddy = 6.78 · 10⁻⁵Wm⁻³T⁻²s², k_eddy = 8.77 · 10⁻⁶Wm⁻³T^−1.5s^1.5, α = 2. Eq. (7), referred as Bertotti’s equation, is only used to compare the results of the PAM model with the most diffused semi-empirical method for loss calculation, but could lead to wrong estimations when applied to rotational loci or with frequencies f > 400 Hz according to [8]. The validity of (5) has been extensively discussed in [6].

4 Results

Figure 2 shows the x and y components of the magnetic flux density in the points P1, P2 and P3 depicted in Figure 1. Since the differences between the waveforms with hysteresis and without hysteresis are hardly distinguishable, two different zooms have been depicted in Figures 3 and 4. In particular Figure 4 shows the multiple inflection points of B_y computed in P1; multiple minor loops are therefore expected in the B_yH_y plane when the the flux density is maximum.

Figure 2

Magnetic flux density waveforms evaluated in the points P1, P2 and P3 of Figure 1. The results for both nonlinear BH curve without and with hysteresis are reported at inner rotational speed v_in = 955 rpm.

Figure 3

Zoom of Figure 2. The difference between the magnetic flux density waveforms is only due to the closed loop implementation of the PAM model. A small inflection point on the B_x component on the point P1 is depicted.

Figure 4

Zoom of Figure 2. Multiple inflection points are shown in the B_y waveform of the point P1.

Figure 5 shows the rotational flux loci computed at the points P1, P2, P3 again in both cases with and without hysteresis.

Figure 5

Magnetic flux density loci at points P1, P2, P3 with and without hysteresis. The rotational behavior of the flux density in the iron poles is clearly visible.

Figure 6 shows the B_x versus H_y waveform at P1. Similar results are obtained for P2 and P3. Thus in this paper only the results relative to P1 has been reported.

Figure 6

B_x versus H_x in point P1 calculated without hysteresis and with the PAM model applied in closed loop and in open loop as post processing. Minor loops are highlighted in the low flux density region.

Figure 7 shows the results for the y component of P1.

Figure 7

B_y versus H_y in P1 calculated without hysteresis and with the PAM model applied in closed loop and in open loop as post processing. Minor loops are highlighted in the high flux density region.

Figure 8 shows the torques developed on inner and outer rotors and on the iron poles when the gear load angle is maximum.

Figure 8

Computed torques on the inner, outer and iron poles rotor when the magnetic gear is at his maximum load capability. The torque ripple is lower than 0.5% because of the high number of pole pairs. The difference due to the hysteresis model is negligible on the torque waveforms in this case. The torque ripple harmonics are in accordance with the predictions in [7].

5 Discussion

As depicted in Figure 2, the field waveforms are composed of the fundamental harmonic with some additional higher order harmonics, due to the interaction between inner and outer magneto motive forces.

The implementation of the PAM model rather than its anhysteretic part during the FEM calculation affects only slightly the magnetic flux density (Figure 3). This small difference justifies the application of hysteresis models as post processing (or in open loop), in an effort to combine the accuracy of the hysteresis models with the efficiency of the nonlinear FEM [9].

The flux loci in Figure 5 is strongly rotational in all the nodes of the iron poles and the influence of the closed loop implementation is clearly visible.

The analysis of Figures 6 and 7 is the key point in this paper that allows us to compare the open loop and closed loop implementation of the hysteresis models. In the case without hysteresis the component wise BH characteristics enclose an area in the first quadrant that is opposite to the one on the third quadrant: this is due to the fact that the nonlinear BH curve is applied at the absolute values of B and H, thus when looking at the x components the y components effects are implicitly included. In the cases with hysteresis the difference between open and closed loop has a remarkable impact on the BH loops: in particular with the closed chain the loop area is bigger than the other case. Some minor loops are also visible, as expected from Figure 3 and 4 where B_x and B_y have an inflection point: minor loops have been observed in both implementations.

The same consideration of Figure 6 applies with the exception of the minor loops that appear at high flux densities again due to the inflection points visible in Figure 2.

According to Figure 8, the torque ripple is lower than ΔR = 0.5% for all the rotors because of the high number of pole pairs of the test case geometry. The shift due to the different model adopted is of the order of ΔT = 0.07%. The shift is directly linked to the power losses obtained in the iron poles when using the closed loop hysteresis implementation. In order to estimate the power loss, the areas of the BH loops have to be calculated:

Ploss=Px+Py=f∮HxdBx+f∮HydBy(8)

In order to improve the integral calculation, the H components are interpolated through piecewise cubic Hermite polynomial functions and then the integration is performed. This procedure allows us to estimate area of arbitrary loops accurately when the number of data points is reduced or crowded in certain regions of the BH space. The minor loop calculation is automatically included.

In Table 2 the results relative to two different gear speeds are reported. Eq. (8) is used using the magnetic field H calculated from FEM in the first column of Table 2 or through (5) using the magnetic flux density B from FEM in the second column. In the case of open loop, P_loss is zero since both the integrals are zero. When adopting the closed loop hysteresis, the post processing calculation of H is performed through the same law adopted in the finite element implementation, thus the calculation will yield the same results: as expected in the open loop case direct integration gives approximatively zero (*) and in the closed loop hysteresis model direct integration and post processing calculation are in agreement (∧ and ∧∧). The most interesting comparison is the one between the losses computed through the open loop and closed loop application of the PAM equation (** and ∧). In Table 2 the losses computed through closed loop are respectively 60% and 50% higher than the ones computed with open loop approximation. This implies that applying the PAM equation as post processing calculation could lead to very misleading results, even if the magnetic flux density B variation due to the closed loop implementation is limited.

Table 2

Comparison between closed and open loop loss calculation. H_p stands for the post processing magnetic field while H is the one adopted in the finite element model. No hyst means that only the anysterecic part of the material model is implemented in FEM, while Hyst means that the FEM material model is the one in (5) with both hysteretic and anhysteretic parts.

W	∫ HdB	∫ H_pdB	P_hyst	P_eddy	P_exc	P_st
No hyst	−0.01^*	3.12^**	3.33	1.66	0.16	5.15
Hyst	4.91^∧	4.95^∧∧	3.31	1.65	0.16	5.12

No hyst	−0.03^*	24.29^**	11.1	18.41	0.99	30.5
Hyst	36.2^∧	36.3^∧∧	11.01	18.1	0.98	30.1

The losses computed through (7) and (8) are in agreement, thus the procedure based on the PAM model provides physically meaningful results and the material coefficients p₀ – p₅ are reliable. In the high speed case, the losses computed through Bertotti's equation are 20% lower than the closed loop implementation while in the low speed case the discrepancy is less noticeable. The mismatch occurs since the loss separation method is a simple procedure normally introduced in the linear material case [8]. Several modifications to the standard equation have been introduced in order to adapt the loss model to the more general cases such as waveforms with minor loops, DC biases and non linearities as shown in [10, 11]. In the magnetic gear case, where the B loci are rotational and minor loops are present, (7) provide a poor estimate of losses. The theoretical rigorous approach for loss computation is the one in (8), where B and H take into account hysteresis, eddy currents and material non linearities. Assuming that the material coefficients p₀ – p₅ are exactly fitted for the case under investigation, (8) should provide a better loss estimation than (7).

The discussion presented in this paper is based on the magnetic gear test case but the results can be extended to the general case: in fact the BH curves only affect the material coefficients p₀ – p₅ while all the other comparisons between the open and closed loop implementations still hold.

6 Conclusions

In this paper a hysteresis model is applied to the magnetic transmission gear iron poles for a more accurate and realistic numerical device simulation. In particular, the open and closed loop results related to the application of a pragmatic homogenized hysteresis model are compared and the resulting losses match the approximated values computed through the well-known Bertotti’s equation. This paper highlights that simple post processing application of PAM model could lead to loss underestimation. Thus, other models should be adopted when hysteresis is not implemented in closed loop in the finite element algorithm. This consideration is valid not only for the case of magnetic gears but is generally true for arbitrary BH waveforms as well.

References

[1] Tlali P. M., Wang R. J., Gerber S., Magnetic gear technologies: A review, 2014, International Conference on Electrical Machines (ICEM), Berlin, 2014, 544-550. 10.1109/ICEL-MACH.2014.6960233Search in Google Scholar

[2] Filippini M., Alotto P., An optimization tool for coaxial magnetic gears, COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, 2017, 36, 1526-1539, 10.1108/COMPEL-02-2017-0068Search in Google Scholar

[3] Li Y. B., Niu S., Ho S. L., Li Y., Fu W. N., Hysteresis Effects of Laminated Steel Materials on Detent Torque in Permanent Magnet Motors, in IEEE Transactions on Magnetics, 2011, 47, 3594-3597, 10.1109/TMAG.2011.2155633Search in Google Scholar

[4] Dlala E., Belahcen A., Pippuri J., Arkkio A., Interdependence of Hysteresis and Eddy-Current Losses in Laminated Magnetic Cores of Electrical Machines, in IEEE Transactions on Magnetics, 2010, 46, 306-309, 10.1109/TMAG.2009.2032930Search in Google Scholar

[5] Krüttgen C., Steentjes S., Glehn G., Hameyer K., Parametric homogenized model for inclusion of eddy currents and hysteresis in 2-D finite element simulation of electrical machines, in IEEE Transactions on Magnetics, 2017, 53, 1-4, 10.1109/TMAG.2017.2660460Search in Google Scholar

[6] Henrotte F., Steentjes S., Hameyer K., Geuzaine C., Pragmatic two-step homogenisation technique for ferromagnetic laminated cores, in IET Science, Measurement and Technology, 2015, 9, 152-159, 10.1049/iet-smt.2014.0201Search in Google Scholar

[7] Niguchi N., Hirata K., Cogging-Torque Analysis of Magnetic Gear, in IEEE Transactions on Industrial Electronics, 2012, 59, 2189-2197, 10.1109/TIE.2011.2159357Search in Google Scholar

[8] Eggers D., Steentjes S., Hameyer K., Advanced Iron-Loss Estimation for Nonlinear Material Behavior, in IEEE Transactions on Magnetics, 2012, 48, 3021-3024, 10.1109/TMAG.2012.2208944Search in Google Scholar

[9] Lin D., Zhou P., Fu W. N., Badics Z., Cendes Z. J., A dynamic core loss model for soft ferromagnetic and power ferrite materials in transient finite element analysis, in IEEE Transactions on Magnetics, 2004, 40, 1318-1321, 10.1109/TMAG.2004.825025Search in Google Scholar

[10] Barbisio E., Fiorillo F., Ragusa C., Predicting loss in magnetic steels under arbitrary induction waveform and with minor hysteresis loops. IEEE Transactions on Magnetics, 2004, 40, 1810-1819, 10.1109/TMAG.2004.830510Search in Google Scholar

[11] Reinlein M., Hubert T., Hoffmann A., Kremser A., Optimization of analytical iron loss approaches for electrical machines, 2013 3rd International Electric Drives Production Conference (EDPC), Nuremberg, 2013, 1-7, 10.1109/EDPC.2013.6689759Search in Google Scholar

Received: 2017-11-01

Accepted: 2017-11-30

Published Online: 2018-04-18

This work is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 License.

Articles in the same Issue

https://doi.org/10.1515/phys-2018-0017

Keywords for this article

Magnetic gear; hysteresis models; homogenization

Creative Commons

BY-NC-ND 4.0