Abstract
This paper presents in detail the extension of the - formulation for eddy currents based on higher-order hierarchical basis functions so that it can automatically deal with conductors of arbitrary topology. To this aim, we supplement the classical hierarchical basis functions with nonlocal basis functions spanning the first de Rham cohomology group of the insulating region. Such nonlocal basis functions may be efficiently and automatically found in negligible time with the recently introduced Dłotko–Specogna (DS) algorithm. The approach presented in this paper merges techniques together which to some extent already existed in literature but they were never grouped together and tested as a single unit.
1. Introduction
This paper considers the problem of computing the so-called eddy currents [1] that come up in an electrically conducting body when it is immersed in a slowly varying magnetic field . Electromagnetic phenomena are governed by Maxwell’s equations [2] and material constitutive relations. For slowly time-varying fields whose change in magnetic field energy is dominant and electromagnetic wave propagation can be ignored, it is wise to neglect the displacement current density in the Ampère–Maxwell’s equation [2]. Thus, the equations to be solved in the time-harmonic case [1] arewhere is the electric field, is the electric conductivity (that vanishes inside the insulators), is the magnetic permeability, is the angular frequency, and is the imaginary unit.
In this paper we focus on the widely used - Finite Element formulation [1] for solving eddy currents problems. This formulation is hardly the only choice for the solution of this particular physical problem; nevertheless it is one of the most popular and computationally attractive solutions given that it reduces the number of unknowns by exploiting the fact that the magnetic field is irrotational in the insulating region (where vanishes).
Let us assume that the conductors form a combinatorial three-dimensional manifold [3] and that , where represents the insulating region and the whole computational domain. To make the presentation easier, we assume that is simply connected. This is the case in the great majority of practically relevant problems. Furthermore this assumption does not affect the generality of the results of this paper, given that it may be removed as described in [4].
If one wants to solve the eddy current problem with the magnetic scalar potential in configurations containing topologically nontrivial conductors (e.g., conductors with “handles” like a torus), representatives of first cohomology group generators of insulator are the most efficient mean to make the problem well defined. They are typically addressed as cuts in the community.
A common misconception is that the problem of finding cuts has been solved many years ago. Yet, there seems to be not a single proposed solution available in literature which is fully automatic (in the sense that the user needs not to provide any homological information as input), provably correct, and acceptably fast. To back up this statement let us review the approaches available in the literature so far.
One pioneering attempt to solve this issue for the lowest order - Finite Element formulation was proposed in [5, 6]. In these papers Kotiuga proposed various algorithms to produce a set of cuts, whose elements are representatives of generators of a second relative homology group basis [3] of realized as discrete surfaces in . Yet, those are restricted to non-self-intersecting surfaces, which is a computationally expensive choice. As a result, the algorithms, although a considerable progress, are not to be considered practical for industrial size engineering problems. Incidentally, the drawbacks of these approaches is one of the main reasons why the authors chose to use cohomology [3] theory instead.
In principle, one can use rigorous methods from standard algebraic topology, like the ones presented in [7], to compute representatives of a cohomology basis. These approaches are based on reduction of the complex followed by a Smith Normal Form (SNF) computation on the boundary matrix of the reduced complex. Unfortunately, this route to finding cohomology generators easily leads to the solution of the topological problem to be the bottleneck of the whole simulation (see [8] for a recent comparison); therefore such methods to solve this issue are not used in electromagnetic solvers.
On the other hand, there is a large body of literature in which the problem of dealing with nonsimply connected conductors is claimed to be solved by using heuristic approaches. Most of these papers present a general formulation for nonsimply connected conductors, but the required topological preprocessing is not actually completely dealt with. For example in [9] the unknown source fields which have to be precomputed to use the formulation must be still provided by hand. In fact, in the vast majority of papers, an algorithm to produce these unknown source fields is omitted, and even a precise mathematical definition of what these fields are is totally missing. Some papers like [10, 11] propose a partial solution, in the sense that the computation of the unknown source fields is addressed. Yet, they still expect to have some topological information as input. In particular, the algorithm of both papers assumes to have one current distribution as input for each unknown source field to be computed. It is easily seen from nuclear fusion engineering applications, in which the conductor is a -fold torus with being , that finding of such current distributions is a very hard problem by itself which is not addressed in any paper. Addressing this issue is a computationally nontrivial task: since we assume that the meshes of the conductor and of the insulator are such that the sum of them is topologically trivial, the so-called exact sequence [12] of the pair of homology groups of insulating and conducting complexes yields to a one-to-one correspondence between the topology of the conductor and the topology of the insulator. In other words, finding independent currents to compute the unknown source fields is to some extent as difficult as finding the unknown source fields. In simple academic cases (like a simple coil, i.e., a torus) those independent currents can be found by hand or by solving a direct-current (DC) problem, but for more complicated cases (as the aforementioned -fold torus) a dedicated algorithm is needed.
In this paper, we propose to perform the topological preprocessing with the DS algorithm (see [4, 13]) which is, to the best of the authors knowledge, the first and so far the only fully automatic, provably correct, and fast way of obtaining cohomology generators. Yet, the set of representatives provided by the DS algorithm is a lazy cohomology basis [4, 13]: 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 [4] or the lean cohomology computation described in [14]) given a lazy one. However, it has been verified that these techniques which produce a standard cohomology basis do 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. We remark that the linear system matrix is singular when a set of lazy generators is used and therefore a linear solver able to deal with singular systems is required. We used the Intel MKL Pardiso which is able to solve general complex and symmetric systems (by setting the matrix type input parameter to 6).
The DS algorithm was originally developed in the framework of the Discrete Geometric Approach (DGA) described in [15]. Yet, both the DGA and the standard Finite Element counterpart based on Whitney edge elements [16] 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. There are various possibilities of obtaining a high order of convergence. The hierarchical basis functions introduced in [17, 18] are particularly appealing. They allow having 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 -adaptivity, automatically by a mesh refinement scheme. Furthermore, these hierarchical higher-order basis functions have been applied to the solution of eddy currents problems using the - Finite Element formulation in [19]. Yet, the authors of [19] assume that the conducting region is simply connected.
In this contribution we will supplement the basis functions introduced in [19] with nonlocal basis functions spanning the first de Rham cohomology group [3] of the insulating region.
We were able to find in literature the closest attempt to solve the problem reported in [20]. In particular, the main contribution of [20] and our paper is exactly to combine the use of cohomology theory [3] to the well-known - Finite Element formulation using hierarchical higher-order basis functions. Yet, in Section 2.4, we point out several shortcomings for the solution proposed in [20], which is based on the algorithm introduced in [16] and also used in [21], to demonstrate that the solution of the problem addressed in this paper is not yet documented in literature. Different from [20], in this paper, we propose to perform the topological preprocessing with the DS algorithm [4], whose typical computational complexity is linear and its correctness is proved for any possible input.
The paper, which details the material presented in [22], is organized as follows. In Section 2 we introduce the novel - formulation employing a lazy cohomology basis together with higher-order hierarchical basis functions to solve eddy current problems with conductors that have arbitrary topology. In the same section, we also survey the algorithm used to extract the lazy cohomology basis. Section 3 presents the results obtained when solving two Testing Electromagnetic Analysis Methods (TEAM) problem benchmarks. Finally, in Section 4, the conclusions are drawn.
2. - Formulation with Higher-Order Hierarchical Basis Functions and Nonlocal Lazy Cohomology Basis Functions
2.1. De Rham Cohomology
If is simply connected (or topologically trivial), as in Figure 1(a), a curl-free vector field, like the magnetic field in an insulator, can be expressed as the gradient of the magnetic scalar potential .
(a)
(b)
(c)
Topology starts playing a role in the - formulation when is not simply connected. Let us consider an example where is the complement of a conductive torus with respect to a box as in Figure 1(b). The magnetic field is curl-free in but its circulation on the loop of Figure 1(b), because of Ampère’s law, has to match the electric current that flows around . The consequence is that the gradient of the magnetic scalar potential is not enough to span the space where the magnetic field lives.
In the following we use concepts of algebraic topology that due to the limited space cannot be reproduced here. Please consult [23] for a formal introduction or [4, 6, 15, 24] for an informal one.
In algebraic topology, the first de Rham cohomology group of , is exactly the space of curl-free vector fields in that are not gradients [3]. Therefore, by its very definition, is the space we have to add to the space generated by the well-known hierarchical basis functions described in [25]. It is known that the dimension of this space is equal to the first Betti number of [3].
Focusing again on the example in Figure 1(b), , thus the basis for the space is composed by just one element called generator. Let us construct the generator as a curl-free field that has a circulation on a loop linking the conductor as in Figure 1(b). If one now considers the product , by varying the independent current one is able to span the whole space. Therefore, 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 [4, 26] for a more formal explanation. Here we informally show what happens in a more complicated example, that is, let us consider a twofold solid torus as ; see Figure 1(c). Here, (the conductor has two handles); thus there are two independent currents and . The space is generated by , where the two generators are such that where is the tangent vector of the loop and is the Kronecker delta.
In general, the magnetic field is represented as where the space of gradients of the magnetic scalar potential has been already defined in [25] 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.
2.2. Construction of Generators
Let us try to replicate the properties of in the example of Figure 1(b) at the discrete level. Let us call the discrete counterpart of . First, the circulation of on a loop made of edges that links the conductor like has to be as for . Therefore, it is natural to define as a discrete field having integer coefficients assigned to mesh edges. Second, we want to be curl-free as its continuous equivalent. This translates in the discrete case to the condition , where is the face-edge incidence matrix of the mesh restricted to . That is, the discrete circulation of 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 [4, 26]. Then, to retrieve , we just need to interpolate with the standard Whitney edge elements basis functions [1]. Thus, it is evident that there is no need to use higher-order basis functions to span the de Rham space.
are realized using the thick cut technique [28]; see Figure 2.
(a)
(b)
The advantage in switching to cohomology over the integers is that the generators of 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 related to can be accounted for when using higher order by adding it to the scalar first-order dofs on the edges of . In what follows we give a detailed description of an algorithm that performs this task.
2.3. Efficient Computation of Generators
We find the required nonlocal basis functions spanning the first de Rham cohomology group space by exploiting the recently introduced idea of lazy cohomology generators [4]. In what follows we recall the DS algorithm that finds the lazy generators of the first cohomology group :(1)First the discrete surface (see Figure 3(a)) is extracted and its first cohomology group generators are computed with a combinatorial algorithm with linear worst-case complexity [4]. In Figure 3(a) the supports of the dual in of two possible representatives and of the cohomology generators are shown.(2)Thinned currents are found by premultiplying the generators of by the incidence matrix between face and edge pairs [4] restricted to .(3)Finally, a vectorialized version of the ESTT algorithm [4] is run on the whole complex for all thinned currents at once. The ESTT algorithm is a general version of the Webb–Forghani iterative algorithm [25] 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 forms the required lazy cohomology generators.
(a)
(b)
The output lazy generators span the cohomology group, but they do not form a base given that they are dependent. In the example, the generator obtained by processing span , whereas the generator obtained with is not useful given that the dual of bounds a discrete surface inside (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 in Figure 3(b). Thus, in general there are no easier ways to do it than the change of basis described in the appendix of [4].
The novel idea of lazy cohomology generators [4] 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 it provides great speedup in the topological preprocessing part. Concerning the quality of the solution, the results with a standard versus lazy basis in terms of eddy currents are the same up to solver tolerance.
This technique is appealing first of all because the topological preprocessing 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.
2.4. Critical Survey on an Alternative Solution
In this section we explain why the technique proposed by He et al. [20] lacks a number of decisive details. Without them, it is impossible to assess its correctness nor to implement and use it. The algorithm presented in [20] consists of two steps, Step and Step , summarized below.
In Step , the algorithm of [20] finds a cut in the two-dimensional boundary of by searching for a maximal collection of 2-simplices such that does not have any nontrivial cycle. One can prove that by doing so the sum of supports of a first cohomology group basis of is the edges in the difference . However, this technique does not specify how to isolate each cut from their union.
Let us now suppose that we are able in some way to separate the each cut. In this case, let us assume that Step of [20] algorithm starts from the two canonical cohomology generators of a torus, like and of Figure 3(a). Then, the tetrahedra from are iteratively added to the surface . One can think about this as process of thickening the surface in the direction of the insulator.
The problem appears when this thickened surface collides with itself. A tetrahedron located in the spot of the collision will not be added to the thickened given that it would close a loop. This produces an additional sprue to the output of Step that has to be eliminated to obtain a correct output; see also [29, Figure 5]. No technique to remove it is discussed in [20]. Let us assume, however, that there is a way to remove all sprues. What we can get in this case is a sum of cut surfaces. Once again it is not discussed how the individual cuts may be obtained from their union. Note that obtaining separated cuts is essential for the electromagnetic computations. Even if there exists a way to separate those cuts, for some surfaces , one simply cannot get the generators by using this procedure. Let us consider an example of a Sudanese Möbius band cut surface represented in Figure 4. It can be obtained when the complement of the surface in the boundary of the conductor contains a noncanonical generator of the torus, as in Figure 3(b). In this case, the “cut” one may obtain from it is a Möbius band; it is therefore nonorientable. Since it is not orientable, it cannot be used to fulfill the global Ampère’s law implicitly. In this case, again, the method as presented in [20] fails.
3. Numerical Results
The DS algorithm together with the - formulation with hierarchical basis functions has been implemented inside the EMS solver (https://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 [30]. It consists of a racetrack shaped coil driven by a time-harmonic current (amplitude AT, frequency Hz) over a square aluminum plate. Figure 5 represents the eddy currents resulting from the solution with first and second order. Table 1 contains the comparison of the computed heat losses , where is the current density vector field and is the resistivity.
(a) 1st order
(b) 2nd order
As a second benchmark, we consider the TEAM Workshop problem , see [31]. Figure 6 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.
(a)
(b)
4. Conclusions
Lazy cohomology generators technology has been integrated inside a commercial 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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.