An Equilibrated Error Estimator for the 2D/1D MSFEM T-Formulation of the Eddy Current Problem

The 2D/1D multiscale finite element method (MSFEM) is an efficient way to simulate rotating machines in which each iron sheet is exposed to the same field. It allows the reduction of the three dimensional sheet to a two dimensional cross-section by resolving the dependence along the thickness of the sheet with a polynomial expansion. This work presents an equilibrated error estimator based on flux equilibration and the theorem of Prager and Synge for the T-formulation of the eddy current problem in a 2D/1D MSFEM setting. The estimator is shown to give both a good approximation of the total error and to allow for adaptive mesh refinement by correctly estimating the local error distribution.


I. INTRODUCTION
T HE simulation of eddy currents in electrical machines consisting of many steel sheets with the finite element method quickly leads to infeasibly large equation systems.In many machines each sheet in the active zone is exposed to the same field, which allows for a great reduction in computational effort by simulating only a single sheet.However, this reduced problem is still far from trivial.One method to further simplify the problem while maintaining a good approximation of the solution is by spacial decomposition.
The thickness of one sheet is less than a millimeter while the length and width are in the range of meters.A method to treat the two dimensional (2D) cross section and the one dimensional (1D) thickness of the sheet as two coupled problems has been presented in [1].It solves the two problems iteratively until convergence is reached.The nature of this coupling has been analyzed in more detail in [2].
In [3] and [4] different approaches have been presented which isolate the one dimensional problem as a pre-processing step in order to obtain parameters for the two dimensional one.
The 2D/1D multiscale finite element method (MSFEM) presented in [5] uses ideas from the multiscale finite element method to use classic finite element functions for the two dimensional problem while approximating the dependence on the third axis with pre-defined polynomial shape functions, similar to the method presented in [6] which is based on trigonometrical shape functions.This enables the solution of the problem within a single iteration while requiring only a mesh for the two dimensional cross section of the sheet.It is also able to include the insulation layers between sheets and correctly treat the edge effect [7].
This paper presents an error estimator for the T-formulation of the 2D/1D MSFEM.It is based on flux equilibration and based on the same theory as the error estimator for the Tformulation for the MSFEM presented in [8].In order to fit within the 2D/1D MSFEM framework it has been restructured so both the construction and the evaluation of the estimator require only the two dimensional mesh while being valid in the complete three dimensional domain.
A numerical example shows that the estimator gives a good approximation of the error in both a global and a local sense.This latter property is used to implement adaptive mesh refinement which allows for a high accuracy of the 2D/1D MSFEM solution while requiring significantly less degrees of freedom than uniform mesh refinement.

II. THE T − Φ-FORMULATION
We use the T − Φ formulation for the reference solution of the eddy current problem as described in [9].The problem domain Ω is split into the conducting domain Ω c , consisting of the steel sheet, and the non-conducting domain Ω 0 , consisting of the air regions and the insulation layers.The sheet is assumed to be axis-aligned with the cross-section in the x − y plane and the thickness aligned with the z axis.The total thickness of Ω is d = d F e + d 0 with the thickness of the sheet d F e and the thickness of the insulation layer d 0 .
The magnetic field strength H is written as with the current vector potential T ∈ H(curl, Ω c ) fulfilling curl T = J, the magnetic scalar potential Φ ∈ H 1 (Ω) and a prescribed Biot-Savart field H BS .The strong formulation of the eddy current problem in the frequency domain is given as where ρ = σ −1 is the electric resistivity with the electric conductivity σ, µ is the magnetic permeability, ω = 2πf with the frequency f and i is the imaginary unit.Multiplication with a test function and integration by parts, together with the auxiliary condition div B = 0, lead to the weak formulation: Find T ∈ H(curl, Ω c ) and Φ ∈ H 1 (Ω) so that arXiv:2302.01601v1 [math.NA] 3 Feb 2023 for all v ∈ H(curl, Ω c ) and all q ∈ H 1 (Ω).

III. THE 2D/1D MSFEM T-FORMULATION
This paper uses the 2D/1D MSFEM approach for the Tformulation which has been described in detail in [5] and [10].The three dimensional unknown components T − ∇Φ are approximated by where Here and in the following, coordinates x, y or z in the index denote the individual components of a vectorvalued function.The shape functions φ 0 and φ 2 are predefined piecewise polynomial of order 0 and 2, respectively, see also Appendix A. The two dimensional rotation operator curl 2D of a two dimensional vector function V = (V x (x, y), V y (x, y)) T is defined as The discretization of the space H(curl 2D ) is discussed in detail in [11].
The full magnetic field strength is then given by For later reference, the (three dimensional) rotation of To obtain the weak formulation, ( 4) is used in (3) for both the trial function and the test function.Note that H 2D/1D only depends on z via the shape functions φ 0 and φ 2 , which are known a-priori.Therefore integration over z can be carried out analytically.This yields the weak 2D/1D MSFEM formulation: for all V 2 ∈ H(curl 2D , Ω 2D,c ) and q ∈ H 1 (Ω 2D ) where a bar denotes that the respective function has been integrated with respect to z.

IV. ERROR ESTIMATION
The proposed error estimator is based on the theorem of Prager and Synge and the theory presented in [12], which can be adapted to obtain the following identity which is the basis for all further calculations: where T is the strong solution of the eddy current problem (2) and γ an equilibrated flux fulfilling the condition The energy norm .ρ can be interpreted as a measurement for the eddy current losses, i.e. for the current density J there holds where the asterisk denotes the complex conjugate.A variant of ( 9) for the two dimensional scalar Tformulation has been proven in [8] and for the vector-valued magnetostatic case in [12].The proof of ( 9) is analogous.
Note that the first term on the left hand side of ( 9) is the error of the 2D/1D MSFEM solution measured in the norm of the eddy current losses.Assuming a suitable γ is known, the right hand side of ( 9) can be calculated.Given that all terms on the left hand side are guaranteed to be non-negative, the right hand side provides an upper bound for the error.
The main problem is the construction of a suitable γ which needs to fulfill (10) on Ω while at the same time being able to be constructed using only Ω 2D .If the error estimator required the full three dimensional domain Ω, it would be much more computationally expensive than the calculation of T 2D/1D and nullify the advantages of using a 2D/1D MSFEM.Similarly, the evaluation of the estimator, as defined by the three dimensional integral on the right hand side of (9), needs to be doable using only Ω 2D .This is achieved by using a 2D/1D MSFEM approach for γ as well.More specifically, we set with φ1 := φ 0 dz and φ3 := φ 2 dz (13) and the unknowns γ 0 , γ 2 ∈ H 1 (Ω 2D ) and γ 1 , γ 3 ∈ H(curl 2D , Ω 2D ) to be determined.Note that the estimator on the right hand side of (9) consists of σγ, which is equal to zero in the insulation because of the conductivity, and curl T 2D/1D , which has only components containing the shape function φ 2 , see (7), which is also zero in the insulation.Therefore it suffices to consider only the domain of the conducting material for the construction γ, i.e. γ 0 , γ 2 ∈ H 1 (Ω 2D,c ) and The rotation of γ is given by Writing out the condition (10) using both ( 4) and ( 14) and comparing the coefficients with respect to the shape functions yields the equations From ( 19) and (20) it follows that γ 1 = ∇Φ 1 and ).With this the remaining equations can be rewritten as Note that (21) and ( 22) do not uniquely define all components of γ.Every solution yields a valid error estimator, but the overestimation (given by the second term on the left hand side of ( 9)) may become arbitrarily large.As can be seen from ( 9), because the error is independent of γ, minimizing the overestimation is equivalent to minimizing the estimator.For this purpose additional conditions are imposed.
Because the estimator is small if σγ is a good approximation of curl T 2D/1D , a comparison of ( 7) and ( 12) suggests that should hold.
In this paper we solve (21) under the constraint and ( 22) under the constraint While this does not yield the optimal minimizer of the estimator because the interdependence of the components is neglected, the numerical example shows that this suffices to achieve an acceptable amount of overestimation.The main advantage of this approach is that instead of one big minimization problem one only has to solve two smaller ones, which is both faster in itself and can even be done in parallel.
The weak formulation for the problem (21) and ( 27) reads as: Find γ 0 , Φ 1 ∈ H 1 (Ω 2D,c ) and a Lagrange multiplier ), where, according to the de Rham complex, the Lagrange multiplier space H λ (Ω 2D,c ) is given as the space H(div, Ω 2D,c ) restricted to divergence-free functions.
Similarly, the weak formulation for the problem ( 22) and (28) reads as: Find Once all components are calculated, the total estimator can be evaluated on the two dimensional mesh as The integrand can also be used locally to identify the finite elements with the highest contribution to the total error.

V. NUMERICAL EXAMPLE
Consider the machine shown in Fig. 1.Using rotational symmetries, only one twelfth of the entire machine has to be simulated.For the steel sheet a magnetic permeability of µ = 1000µ 0 and an electric conductivity of σ = 2.08MS is prescribed.The frequency is 50Hz.The sources are not resolved in the finite element mesh and only included via their Biot-Savart fields.All calculations were done using the opensource software Netgen/NGSolve [13].The calculations start with the coarsest possible mesh for the given geometry, see Fig. 2. In each iteration, the estimator is evaluated for each individual finite element.Then, all elements where this evaluation yields at least half of the maximum encountered estimator, are refined.In this process some adjacent elements might get refined as well in order to avoid hanging nodes.
As can be seen, the refinements are concentrated at the inner edges (where the currents turn around due to the edge effect), at the corners (where the fields peak, see also Fig. 3) and at the inner and outer boundaries (where the boundary conditions need to be resolved correctly).Note also that almost no refinement happens along the vertical symmetry line where the fields are perfectly parallel and easy to resolve.
A qualitative evaluation of the estimator is shown in Fig. 3 where both the error (compared to a high order three dimensional reference solution) and the estimator are depicted after two mesh refinements.It can be seen that the estimator correctly identifies the regions where the error is concentrated, further justifying the refinements.
Finally, as a quantitative evaluation Fig. 4 shows the total error of the 2D/1D MSFEM solution compared to the required degrees of freedom (nDoF) in the finite element problem for both adaptive refinement and uniform refinement.As can be  seen, the adaptive refinement leads to a great increase in the rate of convergence.Furthermore, the estimator gives a good approximation of the behavior of the error with only a small overestimation.Numerical examples show that it gives reliable estimates of the error in both a global and a local sense.This makes it an efficient tool for adaptive mesh refinement to increase the rate of convergence of the 2D/1D MSFEM solution.

Fig. 1 .
Fig. 1.One twelfth of a fictitious machine, consisting of steel sheets (grey) and air domains (left blank).Positive and negative sources are drawn in red and blue, respectively.All measurements are in mm.The thickness of the sheet is d = 0.5mm with a fill factor of 0.95.

Fig. 3 .
Fig.3.The absolute value of J drawn for the reference solution (left), the error with respect to the reference solution (middle) and the error estimator (right) after two refinement steps, drawn on one half of the domain.

Fig. 4 .
Fig.4.The total error, measured in the norm of the eddy current losses, for both adaptive and uniform mesh refinement as well as the total estimated error for the adaptive approach.