$T$-$\Omega$ formulation with higher order hierarchical basis functions for non simply connected conductors

TL;DR

This paper enhances the $T$-$ Omega$ formulation for eddy current problems by incorporating higher order hierarchical basis functions and non-local basis functions to handle conductors with complex topologies.

Contribution

It introduces a novel extension of the $T$-$ Omega$ formulation using non-local basis functions derived from the first de Rham cohomology group, enabling analysis of arbitrary conductor topologies.

Findings

Effective handling of arbitrary topology conductors

Efficient computation of non-local basis functions with DS algorithm

Extension of $T$-$ Omega$ formulation for complex geometries

Abstract

This paper extends the $T$-$\Omega$ formulation for eddy currents based on higher order hierarchical basis functions so that it can deal with conductors of arbitrary topology. To this aim we supplement the classical hierarchical basis functions with non-local basis functions spanning the first de Rham cohomology group of the insulating region. Such non-local basis functions may be efficiently found in negligible time with the recently introduced D{\l}otko--Specogna (DS) algorithm.

Figures

Tables

Table 1. Table 1: TEAM problem 7 benchmark

Number of tetrahedra	Losses 1st order [W]	Losses 2nd order [W]
52 539	7.11	8.02
126 615	7.59	8.83
225 478	9.02	9.37
423 841	9.32	9.60
755 333	9.66	-

Table 2. Table 2: Loss in Watt for the TEAM problem 21 benchmark

Problem	Tetrahedra	Loss 1st ord.	Loss 2nd ord.	Measured
21a-1	45 177	3.19	3.40	3.40
21a-2	29 577	1.58	1.66	1.68
21a-3	199 781	0.95	1.04	1.25

Equations

$${\int }_{c_i}{g}_j\cdot \hat{t}dl={\delta }_{ij},$$

$$h=\nabla \Omega +\sum _{k=1}^{{\beta }_1({\mathcal{K}}_a)}{i}_{c_k}{g}_k,$$

Full text

keywords: finite element method, eddy currents, high order, hierarchical basis functions, cuts, cohomology.

$T$-$\Omega$ Finite Element formulation for solving eddy currents problems is computationally very attractive given that it reduces the number of unknowns by exploiting the fact that the magnetic field is irrotational in the insulating region. Yet, the standard Finite Element formulation based on Whitney edge elements [1] and the Discrete Geometric Approach (DGA) counterpart described in [2] provide a current density which is only uniform inside each mesh element. Higher order basis functions are therefore very attractive since they yield greater accuracy for a given computational cost and smoother current density vector field.

Among the various possibilities to obtain a high order of convergence, the hierarchical basis functions introduced in [3] and [4] are particularly appealing. They allow to have a good control over the distribution of degrees of freedom (dofs) given that different orders can coexist on the same mesh. Moreover, these different orders may be set by hand or, as in $p$-adaptivity, automatically by a mesh refinement scheme.

The $T$-$\Omega$ Finite Element formulation for solving eddy currents problems using hierarchical high order basis functions has been introduced in [5]. Yet, the authors of [5] assume that the conducting region is simply connected. Since in the engineering practice in most cases this assumption is overly restrictive, the aim of this paper is to fill this gap by proposing a formulation valid for manifold conductors of arbitrary topology. A previous attempt to solve this issue is reported in [6], even though necessary details for such a solution are missing, see Sec. 1.3.

Let us assume that the conductors form a combinatorial three-manifold $\mathcal{K}_{c}$ and that $\mathcal{K}=\mathcal{K}_{c}\cup\mathcal{K}_{a}$, where $\mathcal{K}_{a}$ represents the insulating region and $\mathcal{K}$ the whole computational domain. We assume that $\mathcal{K}$ is simply connected since this is the case in the great majority of problems while also making the presentation easier. However, this assumption does not affect the generality of the results of this paper, given that it may be removed as described in [9].

If one wants to solve the eddy current problem with the magnetic scalar potential in configurations containing topologically nontrivial conductors (for example, conductors with “handles” like a torus), representatives of first cohomology group generators of insulator are the only objects that make the problem well defined. Thus, the idea proposed in this contribution is to supplement the basis functions introduced in [5] with non-local basis functions spanning the first de Rham cohomology group $H^{1}_{dR}(\mathcal{K}_{a})$ [7] of the insulating region.

To achieve good overall performances, it is mandatory to address this topological pre-processing efficiently. With classical algebraic methods like [8] based on reduction of the complex $\mathcal{K}_{a}$ followed by a Smith Normal Form (SNF) computation on the reduced complex, finding such non-local basis functions easily becomes the bottleneck of the whole simulation, see [11]. This is the main reason why such rigorous methods to solve this issue did not take off in electromagnetic solvers.

In this paper we propose to perform the topological pre-processing with the DS algorithm [9], whose typical computational complexity is linear. Yet, the set of representatives provided by the DS algorithm is a lazy cohomology basis [9], [10]: the provided set of representatives span the needed cohomology group, but contains additional, dependent elements. The size of the lazy basis is no more than twice the size of a standard cohomology basis and with moderate effort one may produce a standard cohomology basis (see [9]) given a lazy one. However, it has been verified that this technique does not provide any speedup in the solution of the electromagnetic problem while giving exactly the same solution in terms of induced currents up to linear solver tolerance.

The paper is organized as follows. In Section II we introduce the novel $T$-$\Omega$ formulation to solve eddy current problems with high order hierarchical basis functions which works also when conductors have arbitrary topology. Moreover, we survey the algorithm used to extract the lazy cohomology basis. Section III presents the results obtained when solving two TEAM problem benchmarks. Finally, in Section IV, the conclusions are drawn.

1 Novel $T$-$\Omega$ formulation

1.1 de Rham cohomology

If $\mathcal{K}_{a}$ is simply-connected (or topologically trivial), as in Fig. 1a, a curl-free vector field, like the magnetic field $h$, can be expressed as the gradient of the magnetic scalar potential $\Omega$.

Topology starts playing a role in the $T$-$\Omega$ formulation when $\mathcal{K}_{a}$ is not simply-connected. Let us consider an example where $\mathcal{K}_{a}$ is the complement of a conductive torus $\mathcal{K}_{c}$ with respect to a box as in Fig. 1b. The magnetic field $h$ is curl-free in $\mathcal{K}_{a}$ but its circulation on the loop $c$ of Fig. 1b, because of Ampère’s law, has to match the electric current $i_{c}$ that flows around $\mathcal{K}_{c}$. The consequence is that the gradient of the magnetic scalar potential is not enough to span the space where the magnetic field $h$ lives.

In the following we use concepts of algebraic topology that due to the limited space cannot be reproduced here. Please consult [13] for a formal introduction or [12], [2], [9] for an informal one.

In algebraic topology, the first de Rham cohomology group $H^{1}_{dR}(\mathcal{K}_{a})$ of $\mathcal{K}_{a}$, is exactly the space of curl-free vector fields in $\mathcal{K}_{a}$ that are not gradients [7]. Therefore, by its very definition, $H^{1}_{dR}(\mathcal{K}_{a})$ is the space we have to add to the space generated by the well known hierarchical basis functions described in [17]. It is known that the dimension of this space is equal to the first Betti number $\beta_{1}(\mathcal{K}_{a})$ of $\mathcal{K}_{a}$ [7].

Focusing again on the example in Fig. 1b, $\beta_{1}(\mathcal{K}_{a})=1$, thus the basis for the $H^{1}_{dR}(\mathcal{K}_{a})$ space is composed by just one element called generator. Let us construct the generator $g\in H^{1}_{dR}(\mathcal{K}_{a})$ as a curl-free field that has a circulation $1$ on a loop linking the conductor as $c$ in Fig. 1b. If one now considers the product $i_{c}\,g\in H^{1}_{dR}(\mathcal{K}_{a})$, by varying the independent current $i_{c}$ one is able to span the whole $H^{1}_{dR}(\mathcal{K}_{a})$ space. Therefore, $i_{c}$ is a new degree of freedom that has to be added as an additional unknown of the eddy current problem.

This property is generalizable to more complicated topologies, see [14], [9] for a more formal explanation. Here we informally show what happens in a more complicated example, i.e. let us consider a two-fold solid torus as $\mathcal{K}_{c}$, see Fig. 1c. Here, $\beta_{1}(\mathcal{K}_{a})=2$ (the conductor has two handles), thus there are two independent currents $i_{c_{1}}$ and $i_{c_{2}}$. The $H^{1}_{dR}(\mathcal{K}_{a})$ space is generated by $i_{c_{1}}\,g_{1}+i_{c_{2}}\,g_{2}$, where the two generators $g_{1},g_{2}$ are such that

[TABLE]

$\hat{t}$ is the tangent vector of the loop $c_{i}$ and $\delta_{ij}$ is the Kronecker delta.

In general, the magnetic field $h$ is represented as

[TABLE]

where the space of gradients of the magnetic scalar potential $\Omega$ has been already defined in [17] and each independent current has to be included as an additional degree of freedom to span the whole space of the eddy current problem solution.

1.2 Construction of $H^{1}_{dR}(\mathcal{K}_{a})$ generators

Let us try to replicate the properties of $g$ in the example of Fig. 1b at the discrete level. Let us call $\mathbf{g}$ the discrete counterpart of $g$. First, the circulation of $\mathbf{g}$ on a loop $\mathbf{c}$ made of edges that links the conductor like $c$ has to be $1$ as for $g$. Therefore, it is natural to define $\mathbf{g}$ as a discrete field having integer coefficients assigned to mesh edges. Second, we want $\mathbf{g}$ to be curl-free as its continuous equivalent. This translates in the discrete case to the condition $\mathbf{C_{a}}\mathbf{g}=\mathbf{0}$, where $\mathbf{C_{a}}$ is the face-edge incidence matrix of the mesh restricted to $\mathcal{K}_{a}$. That is, the discrete circulation of $\mathbf{g}$ over the edges in the boundary of each face of the mesh must be zero. It is possible to prove that such edge discrete fields may be interpreted as the representatives of the generators of the first cohomology group over integers $H^{1}(\mathcal{K}_{a},\mathbb{Z})$ [14], [9]. Then, to retrieve $g$, we just need to interpolate $\mathbf{g}$ with the standard Whitney edge elements basis functions [15]. Thus, it is evident that there is no need to use higher order basis functions to span the de Rham space.

The $H^{1}(\mathcal{K}_{a},\mathbb{Z})$ are realized using the thick cut technique [16], see Fig. 2.

The advantage in switching to cohomology over the integers is that the generators of $H^{1}(\mathcal{K}_{a},\mathbb{Z})$ can be rigorously and efficiently computed. In this respect, it is important to note that the property of the basis functions of being hierarchical is of fundamental importance: thus, the term of $h$ related to $H^{1}(\mathcal{K}_{a},\mathbb{Z})$, can be accounted for when using higher order by adding it to the scalar first order dofs on the edges of $\mathcal{K}_{a}$. In what follows we give a detailed description of an algorithm that performs this task.

1.3 Efficient computation of $H^{1}(\mathcal{K}_{a},\mathbb{Z})$ generators

We find the required non-local basis functions spanning the first de Rham cohomology group space by exploiting the recently introduced idea of lazy cohomology generators [9]. In what follows we recall the DS algorithm that finds the lazy generators of the first cohomology group $H^{1}(\mathcal{K}_{a},\mathbb{Z})$:

First the discrete surface $\mathcal{S}=\mathcal{K}_{c}\cap\mathcal{K}_{a}$ (see Fig. 3a) is extracted and its first cohomology group generators $H^{1}(\mathcal{S},\mathbb{Z})$ are computed with a combinatorial algorithm with linear worst-case complexity [9]. In Fig. 3a the supports of the dual in $\mathcal{S}$ of two possible representatives $\mathbf{c}^{1}$ and $\mathbf{c}^{2}$ of the cohomology generators are shown. 2. 2.

Thinned currents are found by pre-multiplying the generators of $\mathcal{S}$ by the incidence matrix $\mathbf{C}_{c}$ between face and edge pairs [9] restricted to $\mathcal{K}_{c}$. 3. 3.

Finally, a vectorialized version of the ESTT algorithm [9] is run on the whole complex $\mathcal{K}$ for all thinned currents at once. The ESTT algorithm is a general version of the Webb–Forghani iterative algorithm [17] to obtain a discrete field whose discrete curl has to match the thinned current vectors obtained at the previous step. Its typical complexity is linear, even though the worst-case complexity is cubical. The output of the ESTT restricted to $\mathcal{K}_{a}$ form the required lazy cohomology generators.

The output lazy generators span the $H^{1}(\mathcal{K}_{a},\mathbb{Z})$ cohomology group, but they do not form a base given that they are dependent. In the example, the generator obtained by processing $\mathbf{c}^{1}$ span $H^{1}(\mathcal{K}_{a},\mathbb{Z})$, whereas the generator obtained with $\mathbf{c}^{2}$ is not useful given that the dual of $\mathbf{c}^{2}$ bounds a discrete surface inside $\mathcal{K}_{c}$ (therefore it is not topologically interesting but it is used anyway in the final system of equations).

Disentangling boundary generators is far from being trivial especially because it is in general not enough to pick half of the generators to produce the suitable basis. This is because the generators of the basis may be “mixed”, meaning that one of them is a linear combination of the two as $\mathbf{c}^{3}$ in Fig. 3c. Thus, in general there are no easier ways to do it than the change of basis described in the appendix of [9].

The novel idea of lazy cohomology generators [9] is exactly that all boundary generators produce a cut, even if half of them may be eliminated. This approach does not provide a sensible slowdown in the system solver time, whereas greatly speedup the topological pre-processing part. Concerning the quality of the solution, the results with a standard vs lazy basis in terms of eddy currents are the same up to solver tolerance.

This technique is appealing first of all because the topological pre-processing requires mere seconds even on very complicated problems. Moreover, the DS algorithm is proved to be general, therefore it works for every possible input, no matter how complicated the geometry or the topology of the insulating region is.

Finally, we note that the technique proposed by Ren [1] and used later by He et al [6] lack a number of decisive details. The description of the algorithm does not specify if surface cuts (i.e. boundary generators) are disentangled or not (and, if they are then how). It is also not clear whether the surface cuts (i.e. cohomology generators) are extended in $\mathcal{K}_{a}$ or in the whole $\mathcal{K}$. The idea of extending topologically trivial sets that seems to be at the root of the algorithm in [1] suffers from very serious drawbacks that may prevent successful termination of the algorithm. However, given the level of generality of the presentation in [1], it is not possible to assess its correctness nor to implement and use it.

2 Numerical results

The DS algorithm together with the $T$-$\Omega$ formulation with hierarchical basis functions have been implemented inside the EMS solver (www.emworks.com). We validated the software and assessed the performances of first and second order approximation in solving various eddy current problems.

As a first benchmark, we use the TEAM Workshop problem 7, see [18]. It consists of a racetrack shaped coil driven by a time harmonic current (amplitude $2742$AT, frequency $f=$200Hz) over a square aluminum plate. Fig. 4 represents the eddy currents resulting from the solution with first and second order. Table 1 contains the comparison of the computed heat losses $P=\int_{\mathcal{K}}\rho j^{2}$, where $j$ is the current density vector field and $\rho$ is the resistivity.

As a second benchmark, we consider the TEAM Workshop problem 21, see [19]. Fig 5 contains the eddy currents resulting from the solution with first and second order in case of benchmark 21a-3. Table 2 contains the comparison of the computed heat losses for the cases 21a-1, 21a-2 and 21a-3.

3 Conclusions

Lazy cohomology generators technology has been integrated inside an eddy current solver employing second order hierarchical basis functions. As expected, the second order formulation indeed provides more accurate results than its first order counterpart. The inclusion of automatic mesh adaptivity to render the second order more efficient is left for further studies.