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

Abstract

Purpose

This work introduces an efficient and accurate technique to solve the eddy current problem in laminated iron cores considering vector hysteresis.

Design/methodology/approach

The mixed multiscale finite element method based on the based on the T,Φ-Φ formulation, with the current vector potential T and the magnetic scalar potential Φ allows the laminated core to be modelled as a single homogeneous block. This means that the individual sheets do not have to be resolved, which saves a lot of computing time and reduces the demands on the computer system enormously.

Findings

As a representative numerical example, a single-phase transformer with 4, 20 and 184 sheets is simulated with great success. The eddy current losses of the simulation using the standard finite element method and the simulation using the mixed multiscale finite element method agree very well and the required simulation time is tremendously reduced.

Originality/value

The vector Preisach model is used to account for vector hysteresis and is integrated into the mixed multiscale finite element method for the first time.

Citation

Hanser, V., Schöbinger, M. and Hollaus, K. (2022), "A mixed multiscale FEM for the eddy current problem with T,Φ-Φ and vector hysteresis", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. 41 No. 3, pp. 852-866. https://doi.org/10.1108/COMPEL-02-2021-0053

1. Introduction

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.

2. Preisach model

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

2.1 Scalar Preisach model

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

(1) B(t)=Γ^[H(t)]=∬Tmaxμ(α,β)γαβ[H(t)] dαdβ.

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

(2) T(α′,β′):={(α,β):α≥β,α≤α′,β≥β′}

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

(3) E(α,β)=2∬T(α,β)μ(α,β) dαdβ.

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

(4) Hmax,in(t)=max⁡t′∈[0,t]|H(t′)|

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

(5) B(t)=−12E(Hmax,in,−Hmax,in)+∑k=0N−1(E(Mk,mk−1)−E(Mk,mk)).

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

(6) E(α,β)=K1b2atan(β+ab)(atan(β−ab)−atan(α−ab))+K2e(atan(βe)−atan(αe))+f(β−α)−∫αβK1b3atan(ξ+ab)ξ2−2aξ+a2+b2 dξ.

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

(7) μΔ=Γ^[H(t)]−Γ^[H(t)±ΔH]ΔH,

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

2.2 Vector preisach model

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

(8) B(t)=∫θ=0π∫φ=02πeRΓ^[H(t)·eR]dφdθ.

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

(9) B≈∑i=0N−1wiΓ^[eR,i·H]eR,i.

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

(10) μ¯¯Δ=∑i=0N−1wi μR,iΔ eR,ieR,iT,

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.

3. Standard formulation

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).

3.1 T,Φ-Φ formulation

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

(11) ∇·(J+JBS)=∇·∇×H≡0

is identically zero. Therefore, the approach

(12) J+JBS=∇×(T+TBS)

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

(13) TBS(r)=14π∫ΩSJBS(r′)×r−r′|r−r′|3 dΩ

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

(14) H={T+TBS−∇Φ in ΩcTBS−∇Φin Ω0

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

(15) ∇×(ρ∇×T)+∂tμ¯¯T−∂tμ¯¯∇Φ=−∇×(ρ∇×TBS)−∂tμ¯¯TBS in  Ωc

(16) −∂t∇·μ¯¯∇Φ+∂t∇·μ¯¯T=−∂t∇·μ¯¯TBS   in  Ωc

(17)       T×n=0,       Φ=0     on  ΓHc

(18) (ρ∇×T)×n=−(ρ∇×TBS)×n,μ¯¯(T−∇Φ)·n=−μ¯¯TBS·n on  ΓE.

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

(19) −∂t∇·(μ0∇Φ)=−∂t∇·(μ0TBS)  in Ω0

(20)    −μ0∇Φ·n=−μ0TBS·n  on  ΓB

(21)           Φ=0       on  ΓH0.

Finally, the interface conditions are

(22) μ¯¯(T−∇Φ)·nc−μ0∇Φ·n0=−μ¯¯TBS·nc−μ0TBS·n0T×nc=0 on Γ0c,

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.

3.2 Fixed-point method

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

(23) μFPΔ=μmaxΔ+μminΔ2

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.

3.3 Weak formulation

The differential approach

(24) ∂tB=∂B∂H∂H∂t ≈ μ¯¯Δ1ΔtHΔ

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.

4. Mixed multiscale formulation

4.1 Mixed multiscale approach

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

(25) H˜={TBS+T2ϕ2−∇Φ0 in ΩcTBS−∇Φ0in Ω0

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

(26) ϕ2:z∈[−d/2,d/2]→ℝz↦{0  for z∉[−dFe/2,dFe/2]38(4dFe2z2−1)                else

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).

4.2 Material parameters

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.

4.3 Weak formulation

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

5. Numerical example

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

(27) p(t)=∫ΩcρJ(t)·J(t) dΩ

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

(28) P=1T∫t′t′+Tp(t) dt

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

(29) ε1=|PSFEM−PMMSFEMPSFEM|100%

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

(30) ε2=100%T∫t′t′+T|pSFEM−pMMSFEMpSFEM| dt

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.

6. Conclusion

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.

Figures

Figure 1.

Initial magnetisation curve and the hysteresis loop with an Everett function based on the Lorentzian function.

Figure 2.

Lebedev distribution for points on the surface of a unit sphere

Figure 3.

Second order Gauss-Lobatto polynomial as micro-shape function ϕ₂ in a sheet with an air gap of width d₀

Figure 4.

Global (left) and local (right) integration points in the bulk core

Figure 5.

One eighth of the laminated iron core (gray) with the coils (red), not drawn to scale

Figure 6.

Top (left) and front (right) view of the geometry, not drawn to scale, dimensions are in mm, planes of symmetry are x = 0, y = 0 and z = 0

Figure 7.

Eddy current losses for I₀ = 1 A peak without hysteresis

Figure 8.

Eddy current losses for I₀ = 3 A peak without hysteresis

Figure 9.

Eddy current losses for I₀ = 1 A peak with hysteresis

Figure 10.

Eddy current losses for I₀ = 3 A peak with hysteresis

Figure 11.

Magnetic flux density B for I₀ = 3 A peak-to-peak at z ≈ 45.9 mm and t = 25ms, reference solution with SFEM on the left and MMSFEM on the right

Figure 12.

Current Density J for I₀ = 3 A peak-to-peak at z ≈ 45.9 mm and t = 50ms, reference solution with SFEM on the left and MMSFEM on the right

Table 1.

Parameters for the Lorentzian Everett function for 50 Hz (M400-50A)

a	–8.18773707 · 101
b	4.13538892 · 101
K1	2.00442033 · 10– 2
K2	2.49345353 · 10– 1
e	1.33306276 · 102
f	5.57513398 · 10– 3

Table 2.

Numerical data for various simulations without hysteresis

No. Sheets		4	20	184	184
I 0 in A		3	3	1	3
P in W	SFEM	0.336	1.753	3.956	16.18
	MMSFEM	0.339	1.771	4.009	16.383
ε1 in %		0.881	1.034	−1.337	1.25
ε2 in %		3.51	8.21	1.39	7.28
N DOF	SFEM	58,868	195,940	1,600,928	1,600,928
	MMSFEM	49,859	74,021	122,345	122,345
NVPM	SFEM	14,336	71,680	659,456	659,456
	MMSFEM	100,352	200,704	401,408	401,408
tsim in h	SFEM	6.0	22.2	208.5	210.2
	MMSFEM	7.3	12.9	24.7	23.8

Table 3.

Numerical data for various simulations with hysteresis

No. Sheets		4	20	184	184
I 0 in A		3	3	1	3
P in W	SFEM	0.408	2.086	5.008	19.206
	MMSFEM	0.407	2.079	5.004	19.216
ε1 in %		0.22	0.33	0.1	0.05
ε2 in %		6.37	11.62	1.37	11.72
N DOF	SFEM	58,868	195,940	1,600,928	1,600,928
	MMSFEM	49,859	74,021	122,345	122,345
NVPM	SFEM	14,336	71,680	659,456	659,456
	MMSFEM	100,352	200,704	401,408	401,408
tsim in h	SFEM	10.7	34.1	295.5	294.9
	MMSFEM	12.9	18.3	29.5	30.8

Table 4.

Number of nonlinear iterations

		SFEM		MMSFEM
No. Sheets	I 0	mean	max.	mean	max.
Nonlinear
4	3	1.71	4	2.04	4
20	3	1.83	4	2.07	4
184	1	1.82	3	2.05	3
184	3	1.86	3	2.08	3
Hysteresis
4	3	1.68	8	2.05	10
20	3	1.77	8	2.05	10
184	1	1.80	8	2.1	9
184	3	1.81	3	2.08	3

Appendix I. Weak formulation SFEM

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

(31) Δt∫Ωcρ∇×T(n+1)·∇×T′ dΩ+∫ΩcμFP(TΔ−∇ΦΔ)·T′ dΩ  in Ωc =−∫Ωcμ¯¯ΔTBSΔ·T′ dΩ+∫Ωc(μ¯¯FP−μ¯¯Δ)(TΔ,(n)−∇ΦΔ,(n))·T′ dΩ

and

(32a) ∫Ω0μ0∇ΦΔ·∇Φ′ dΩ=∫Ω0μ0TBSΔ·∇Φ′ dΩ in Ω0

(32b) ∫ΩcμFP(∇ΦΔ−TΔ)·∇Φ′ dΩ=∫Ωcμ¯¯ΔTBSΔ·∇Φ′ dΩ    in Ωc +∫Ωc(μ¯¯FP−μ¯¯Δ)(∇ΦΔ,(n)−TΔ,(n))·∇Φ′ dΩ

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.

Appendix II. Weak formulation MMSFEM

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

(33) Δt∫Ωcρϕ2ϕ2¯ ∇×T2(n+1)·∇×T′2 dΩ+Δt∫Ωcρ∂zϕ2∂zϕ2¯ (T2,x(n+1)T′2,x+T2,y(n+1)T′2,y) dΩ+∫ΩcμFPΔϕ2ϕ2¯ T2Δ·T′2 dΩ−∫ΩcμFPΔϕ2¯ ∇Φ0Δ·T′2 dΩ =∫Ωc(μ¯¯FPΔ−μ¯¯Δ)ϕ2ϕ2¯ T2Δ,(n)·T′2 dΩ  −∫Ωc(μ¯¯FPΔ−μ¯¯Δ)ϕ2¯ ∇Φ0Δ,(n)·T′2 dΩ  −∫Ωcμ¯¯Δϕ2¯ TBSΔ·T2′ dΩ

and

(34a) ∫Ω0μ0 ∇Φ0Δ·∇Φ0′ dΩ  =∫Ω0μ0TBSΔ·∇Φ0′ dΩ in Ω0

(34b) ∫ΩcμFPΔ¯ ∇Φ0Δ·∇Φ0′ dΩ−∫ΩcμFPΔϕ2¯ T2Δ·∇Φ0′ dΩ=∫Ωc(μ¯¯FPΔ−μ¯¯Δ)¯ ∇Φ0Δ,(n)·∇Φ0′ dΩin Ωc−∫Ωc(μ¯¯FPΔ−μ¯¯Δ)ϕ2¯ T2Δ,(n)·∇Φ0′ dΩ  +∫Ωcμ¯¯Δ¯ TBSΔ·∇Φ0′ dΩ

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).