A mixed multiscale FEM for the eddy current problem with T,Φ-Φ and vector hysteresis Open Access

The mixed multiscale finite element method (MMSFEM) has already been successfully introduced for linear and nonlinear eddy current problems (ECPs), see for example Hollaus (2019) and Hollaus and Schöbinger (2020). In addition, the method has been applied to scalar hysteresis as it occurs in ferromagnetism (Schöbinger et al., 2019). Most notably, it has been demonstrated that the hysteresis phenomenon cannot be neglected when comparing simulation results to measurement data.

A 2 D/1D approach considering hysteresis with a comparison of simulation results with measurement data can be found in Bottauscio and Chiampi (2002). A b-conform and a h-conform homogenisation techniques for the ECP in laminated cores is presented in Dular (2008). A two-step technique has been proposed in Bíró et al. (2005). In the first step, the laminated medium is assumed to have an anisotropic electric conductivity and in the second step, the eddy currents are computed individually in each sheet.

The aim of this paper is to efficiently simulate eddy currents in laminated iron cores considering vector hysteresis. Thanks to the mixed T-Φ, Φ formulation, with a current vector potential (CVP) T and a magnetic scalar potential (MSP) Φ, the vector Preisach model (VPM) can be used in the forward mode. Using a magnetic vector potential (MVP) A would require an inverse mode of the VPM, which is computationally expensive.

The developed technique is evaluated by a transformer excited by known currents in coils, see also Hollaus (2019). The excitation is considered by the corresponding Biot-Savart field.

First, a brief explanation of the computationally optimised scalar Preisach model (SPM) and the VPM and their integration into the finite element method (FEM) are given in Sec. II. Then, the standard formulation for T-Φ, Φ to solve the nonlinear ECP by the time stepping method with the FEM is presented in Sec. III. Next, the MMSFEM is introduced to substantially reduce the overall computational costs of the considered ECP in Sec. IV. Simulation results in Sec. V demonstrate that the results obtained by the MMSFEM agree very well with the results obtained by the reference solution.

The Preisach model describes a hysteresis phenomenon (Mayergoyz, 1991).

In ferromagnetic materials, hysteresis occurs between the magnetic field strength H and the magnetic flux density B. The original version of the Preisach model considered scalar hysteresis only. The fundamental idea for the SPM consists of describing the hysteresis effect through an infinite number of weighted two-state operators γ_αβ[H(t)]:  →{0,1}, where α and β denote the upper and lower threshold for switching the current state. The operators are weighted by the Preisach function μ(α,β), which uniquely defines a specific material. The integration of these weighted operators over the Preisach plane determines the magnetic flux density

The Preisach plane T_max := T(H_max, − H_max) is defined as the triangle

for a maximal value of the magnetic field strength H_max, which sets the limits of the Preisach model. Possible larger values must be handled separately in a feasible numerical scheme.The integration over a triangle T(α,β) ⊆ T_max yields the Everett function

Using the Everett function in (1) yields an increase in performance. Moreover, the approach of perfect demagnetisation reduces the average computational costs tremendously (Tousignant et al., 2017). For the perfect demagnetisation approach the maximal absolute input value over time of the magnetic field strength

is stored. Finally, the time-variant magnetic flux density is

The values M_k and m_k denote the essential maxima and minima in the input sequence, respectively.

An additional improvement of the performance can be achieved by storing the subtotals of the sum in (5). This approach avoids recalculation of previously calculated results and is particularly useful when the input varies in a limited range without wiping out previous extrema.

For describing the material, the Lorentzian Preisach function is used (Schöbinger et al., 2019). Its Everett function is given by

The parameters in Table 1 are obtained by solving an inverse problem using measurement data at a frequency of 50 Hz. The remaining integral in (6) does not have an analytic representation and has to be calculated numerically. A major hysteresis loop along with the initial magnetisation curve obtained by these parameters is shown in Figure 1. The discretisation of the Preisach plane is done in 701 steps along the α- and the β-axis, respectively. A critical aspect of magnetic hysteresis is that the magnetic permeability $μ=BH$ is not defined for the magnetic remanence (H = 0) leading to problems in numerical schemes. However, the differential permeability

where ΔH describes a sufficiently small local discrete step size, is greater than zero in any case.

The VPM is a superposition of an infinite number of SPMs. In the three-dimensional case, the SPMs are distributed on the surface of a unit sphere. Hence, each SPM has a position vector e_R which affects the scalar magnetic flux density $BR=Γ^[H·eR]$ of the corresponding SPM. The integration over the full unit sphere surface defines the vectorial magnetic flux density

To uniformly distribute a finite number of SPMs on the unit sphere, Lebedev coordinates with the direction vectors e_R,i and the associated weights w_i have been chosen for the numerical integration of (8) as shown in Figure 2 (Lebedev, 1976). Experiments have shown that Lebedev coordinates are well suited for the distribution of SPM on the surface of a unit sphere. The number of SPMs has been selected with N = 73 as a good compromise between accuracy and computation costs. The numerical integration of (8) yields the approximation for the vectorial flux density

It is worth mentioning that the Everett function (3) of a SPM has to be adapted when it is used in the VPM (Mayergoyz, 1991).

The tensor-valued differential permeability

is a superposition of the scalar differential permeabilities $μR,iΔ$ in (7) weighted by w_i according to the Lebedev coordinates and the outer product of the direction vectors e_R,i pointing in the direction of each SPM.

The ECP couples a static magnetic field in an electrically non-conducting domain Ω₀ with a quasi-static magnetic field in an electrically conducting domain Ω_c. Since the VPM is implemented for the forward mode only, the T,Φ-Φ formulation is used (Bíró, 1991).

Considering Ampere’s law ∇ × H = J + J_BS, with an impressed current density J_BS, the divergence of the current density

is identically zero. Therefore, the approach

with the CVPs T and T_BS can be applied. The Biot-Savart field

is used for J_BS, where r’ and r denote source and field point, respectively. Thus, the magnetic field strength is

with a MSP Φ. Since no eddy currents occur in Ω₀, the CVP T vanishes.

Therefore, the boundary value problem (BVP) with the mixed T,Φ-Φ formulation for the quasi-static magnetic field in Ω_c is described by

Moreover, the static magnetic field in Ω₀ is described by

Finally, the interface conditions are

where Γ_0c is the interface between Ω₀ and Ω_c. The indices E, H or B mean that the tangential components of E or H or the normal component of B are prescribed. The indices 0 or c denote the non-conducting or the conducting domain, respectively.

The nonlinear system is split into a linear and a nonlinear part (Bottauscio and Chiampi, 2002). The fixed-point permeability

is set to be the mean value of the smallest and largest differential permeability. Further, the nonlinear parts are shifted to the right-hand side, leading to constant matrices on the left-hand side. Since the differential permeability of the VPM $μ¯¯Δ$ is tensor-valued, the fixed-point differential permeability $μ¯¯FPΔ$ has to be tensor-valued as well. Therefore, the fixed-point differential permeability is a 3 × 3 diagonal matrix with each diagonal element being equivalent to $μFPΔ$ of the SPM.

The differential approach

with H^Δ= H⁽ⁿ⁺¹⁾ − H⁽ⁿ⁾ has been used in the weak formulation in order to exploit the differential permeability $μ¯¯Δ$⁠. The full weak formulation for the standard finite element method (SFEM) is shown in App. I.

To avoid the necessity to resolve the individual sheets of the iron core, the mixed multiscale approach

with the micro-shape function ϕ₂, which is shown in Figure 3, is selected. The second order Gauss-Lobatto polynomial

is chosen as micro-shape function (Hollaus, 2019). The average MSP Φ₀ takes into account the solution on the large scale while the CVP T₂ along with the periodic micro-shape function ϕ₂ considers the oscillating variation on the small scale (Hollaus and Schöbinger, 2020).

In the mixed-multiscale approach the insulated sheets are considered as bulk material. To obtain the specific material parameters in a global integration point, an additional integration with local integration points is needed, see Figure 4. In each local integration point an independent VPM is set up and is updated in every time instant.

The weak formulation for MMSFEM can be found in App. II.

The ECP of the laminated iron core shown in Figures 5 and 6 is investigated. The thickness of an iron sheet is d_Fe = 0.5 mm and the width of the gap d₀ = 0.01 mm, yielding a fill-factor $kf=dFedFe+d0=0.9804$ and the core is composed of 184 sheets. The electric conductivity of iron is selected with σ = 2 · 10⁶S/m. The excitation of the problem is considered by the Biot-Savart field of four symmetric cylindrical coils with 60 turns each. The frequency is selected with f = 50Hz. The arrangement of the core with the coils exhibits three planes of symmetry. Handmade structured hexahedral finite element meshes have been used. The same discretisation by finite elements in the x, y-plane have been used in the models for the SFEM and MSFEM to ensure a fair comparison. One period of time is discretised in 600 steps. The material of the iron sheets is considered to be demagnetised initially. The simulations have been done by Netgen/NGSolve (Schöberl, 2021).

The reference solution with SFEM has been computed to verify the results obtained by the MMSFEM. Simulations with different impressed currents I₀ and different number of iron sheets have been carried out. For every time instant the eddy current losses

have been calculated. Eddy current losses using the BH-curve (initial magnetisation curve) are shown in Figures 7 and 8 and those considering hysteresis are presented in Figures 9 and 10. The additional peaks in Figures 7 and 8 are due to the convex-concave nature of the BH-curve. The eddy current losses are obtained for simulations with I₀ = 1 A and I₀ = 3 A and with 184 sheets. For the sake of visibility, only every third data point is shown.

Further, the mean value

has been calculated in the steady state. With the aid of P the error

is calculated, which is not very sensitive. For a stricter criterion the error

using the time instants of the eddy current losses p(t), has been introduced. A comparison of the local solutions are shown in Figures 11 and 12 representing the magnetic flux density B and the current density J in the same layer. Further simulation results without hysteresis are summarised in Table 2 and with hysteresis in Table 3. Compared to the nonlinear simulations, the number of iterations increases slightly for the simulations with hysteresis, see Table 4. The reduction of degrees of freedom N_DOF increases with the number of sheets in the laminated iron core for both MMSFEM and SFEM. However, MMSFEM requires essentially less N_DOF. The same holds for the number of VPMs N_VPM, as well as for the required computation time. The number of needed VPMs for small problems with just a few sheets is relatively high, which results in long simulation times t_sim also for MMSFEM.

The MMSFEM introduced in the frequency domain for a linear ECP (Hollaus and Schöbinger, 2020) has been extended for materials with vector hysteresis in the time domain. Simulation results of SFEM with resolved sheets and MMSFEM with a coarse finite element mesh show a very high agreement. The efficiency of MMSFEM compared to SFEM substantially grows with the number of iron sheets.

This work was supported by the Austrian Science Fund (FWF) under Project P 31926.

The weak formulation for SFEM reads as:

Find (T⁽ⁿ⁺¹⁾,Φ⁽ⁿ⁺¹⁾) ∈V_D := {T⁽ⁿ⁺¹⁾,Φ⁽ⁿ⁺¹⁾} :T ∈ u,Φ ∈ v and T⁽ⁿ⁺¹⁾ × n = 0 on Γ_0c ∪ Γ_Hc,Φ⁽ⁿ⁺¹⁾ = 0 on Γ_Hc ∪ Γ_H0}, such that

and

for all (T’,Φ’) ∈ V₀, where u and v are finite element subspaces of H(curl,Ω_c) and H¹(Ω), respectively (Schöberl and Zaglmayr, 2005). The fixed-point method in the weak formulation requires Φ^Δ,(n) and T^Δ,(n), cf. (24), of the former time instant.

The weak formulation for MMSFEM reads as:

Find $(T2(n+1), Φ0(n+1))∈VD:={(T2(n+1), Φ0(n+1)):T2(n+1)∈U, Φ0(n+1)∈V, T2(n+1)×n=0 on Γ0c∖ΓT2∪ΓHc, Φ0(n+1)=0 on ΓHc∪ΓH0}$⁠, such that

and

for all (T’₂,Φ’₀) ∈ V₀, where u ⊂ H(curl,Ω_c), v ⊂ H¹(Ω) and $ϕ2∈Hper1(Ωc)$ have been selected. The interface $ΓT2$ is the part of Γ_0c which represents the smooth surface of the laminated iron core, see Figure 5. The averaged coefficients are denoted by the bar and are calculated as described in (Hollaus and Schöberl, 2018). The fixed-point method is exploited to solve the nonlinear problem, see Sec. III-B. The VPM is applied to describe the ferromagnetic material, see Sec. II-B. To overcome the singularity of the tensor-valued permeability in the magnetic remanence, a differential approach is chosen, see (7).