# Comparison Between NEM and FEM in 2D Magnetostatics Using an Error Estimator

Antoine Bruyère, Lounès Illoul, Francisco Chinesta, Stéphane Clénet

► To cite this version:

Antoine Bruyère, Lounès Illoul, Francisco Chinesta, Stéphane Clénet. Comparison Between NEM and FEM in 2D Magnetostatics Using an Error Estimator. IEEE Transactions on Magnetics, 2008, 44 (6), pp.1342-1345. 10.1109/TMAG.2007.915951. hal-01005866

HAL Id: hal-01005866

https://hal.science/hal-01005866

Submitted on 2 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comparison Between NEM and FEM in 2-D Magnetostatics
Using an Error Estimator

A. Bruyère¹, L. Illoul², F. Chinesta², and S. Clénet¹

¹L2EP, ENSAM, CER Lille, 59046 Lille Cedex, France

²LMSP, ENSAM, CER Paris, 75013 Paris, France

This communication deals with a comparison between two methods of discretization: the well known finite element method and the natural element method that is a meshless method. An error estimator, based on the nonverification of the constitutive law, is used. This estimation has been applied to two examples: a device with permanent magnets and a variable reluctance machine.

*Index Terms*—Error estimation, finite element method, magnetostatics, natural element method.

## I. INTRODUCTION

NOWADAYS, the finite element method (FEM) is well established and allows solving many problems in applied physics. However, difficulties can appear with the FEM when the system includes moving parts. Actually, a remeshing step can be required to avoid a significant distortion of elements or its overlapping. Recently, novel methods have been proposed to circumvent, or at least alleviate, the main difficulties related to the mesh constraints. The most usual mesh constraints are related to the need of efficient remeshing procedures, which in the 3-D case constitutes a difficulty nowadays under active investigation, as well as to the field projection from the old mesh to the new mesh that induces an inevitable numerical diffusion whose unfavorable impact on the solution after numerous resmeshing steps is widely accepted. It was in this context that the so-called “meshless methods” were introduced one decade ago. Their applications in several domains experienced a spectacular progression in the last decade, after proving its ability for avoiding, or at least for alleviating significantly, the remeshing requirements. The natural element method [1] (from now on referred as NEM) belongs to this new and vast family of meshless methods and its interest compared with the vast majority of meshless techniques concern its interpolation property that allows enforcing essential boundary conditions in an easy way, as in the finite element method. We recall that the imposition of essential boundary conditions is the main issue in the other meshless methods. On the other hand, it has been proved that the solution accuracy does not depend on the quality of the subjacent triangulation (“meshless character”). Thus, we could affirm, crudely speaking, that the NEM is a kind of FEM in which the triangles could be distorted without a significant impact on the solution accuracy.

The NEM, which was firstly applied in the framework of solid and fluid mechanics [2], [3], has already been successfully tested in electromagnetic systems for accounting, in a simple manner, the existence of moving parts [4]. However, in the electromagnetic framework the NEM and FEM solution

accuracy comparison has never been addressed. This is a question needing a clear response and consequently justifying or disappointing the use of meshless approaches in this field. In this communication, we compare both methods NEM and FEM by using an error estimator based on the non verification of the constitutive law on two examples [5], [6]. The first example concerns two permanent magnets in opposition and the second example a Variable Reluctance Machine.

## II. PROBLEM DESCRIPTION AND ERROR ESTIMATOR

The magnetostatic problem, defined in the domain D, can be written as

$$ \{ \operatorname{div} \mathbf{B} = 0 \text{ and } \mathbf{B} \cdot \mathbf{n} = 0 \text{ on } S_b \quad (1) $$

$$ \{ \operatorname{curl} \mathbf{H} = \mathbf{J} \text{ and } \mathbf{H} \times \mathbf{n} = 0 \text{ on } S_h. \quad (2) $$

$\mathbf{B}$ is the magnetic flux density, $\mathbf{H}$ is the magnetic field, and $\mathbf{J}$ is the current density that is assumed to be known in the domain $D$. $D$ is bounded by the surface $S = S_h \cup S_b$ where $\mathbf{n}$ represents the outward normal unity vector. Moreover, $\mathbf{B}$ and $\mathbf{H}$ verify the constitutive relationship given by:

$$ \mathbf{B} = \mu(\mathbf{H} + \mathbf{H}_c) \quad (3) $$

with $\mathbf{H}_c$ the coercitive magnetic field and $\mu$ the permeability. In the sequel, the problem is assumed to be linear. Let's call “admissible” fields $\mathbf{B}_{ad}$ and $\mathbf{H}_{ad}$ that verify (1) and (2) respectively. It is obvious that an admissible field couple $(\mathbf{B}_{ad}, \mathbf{H}_{ad})$ verifying the constitutive relationship (3) is the exact solution denoted $(\mathbf{B}_{ex}, \mathbf{H}_{ex})$. We consider now an admissible couple and we define a constitutive law error ‘e’:

$$ e = \| \mathbf{B}_{ad} - \mu(\mathbf{H}_{ad} + \mathbf{H}_c) \|_{D} \quad (4) $$

$$ \text{with } \|\mathbf{X}\|_{D}^{2} = \int_{D} \mu^{-1} \mathbf{X} \cdot \mathbf{X} dD $$

It can be proven [5] that this error e is equal to the distance between an average field $\mathbf{B}_m$ and the exact solution $\mathbf{B}_{ex}$:

$$ e^2 = 4 \|\mathbf{B}_{ex} - \mathbf{B}_m\|_{D}^{2} \text{ with } \mathbf{B}_m = \frac{\mu(\mathbf{H}_{ad} + \mathbf{H}_c) + \mathbf{B}_{ad}}{2}. \quad (5) $$

Therefore, the scalar e is a numerical error estimator for the admissible solution. To obtain the admissible field couple, several methods have been already proposed. The more accurate method consists in solving both potential formulations. The vector potential $\mathbf{A}$ formulation $(\mathbf{B} = \operatorname{curl}\mathbf{A}$ and $\mathbf{A}\mathbf{x}\mathbf{n} = 0$ on $S_b$) yields an
<figure><img src="image_2.png" /><figcaption>Fig. 1. Construction of the Voronoï cell related to node 5.</figcaption></figure>

admissible induction field $\mathbf{B}_{ad}$. The scalar potential $\varphi$ formulation ($\mathbf{H} = \mathbf{H}_s - \text{grad}\varphi$ and $\varphi = \text{C}^{\text{ste}}$ on $S_h$ with $\mathbf{H}_s$ such that $\text{curl} \mathbf{H}_s = \mathbf{J}$ and $\mathbf{H}_s \times \mathbf{n} = 0$ on $S_h$) yields an admissible magnetic field $\mathbf{H}_{ad}$. Then, using (4), an estimation of the numerical error can be calculated. To evaluate the global solution quality, let us also define a global relative error estimator

$$ \hat{\epsilon}^2 = \frac{\|\mathbf{B}_{ad} - \mu(\mathbf{H}_{ad} + \mathbf{H}_c)\|_D^2}{\|\mathbf{B}_{ad} + \mu(\mathbf{H}_{ad} + \mathbf{H}_c)\|_D^2}. \quad (6) $$

### III. DISCRETIZATION

In 2-D, both potential formulations derived from the model defined in (1), (2) and (3), can be solved by discretizing the associated Galerkin's weak form. For this purpose we need to define the approximation to be considered for both the trial and the test functions. This approximation is usually built from the values that the unknown fields have at $N$ nodes distributed on the domain $D$. A shape function $\phi_i(\mathbf{x})$ ($\mathbf{x}$ is the position vector) is associated to each node $i$ ($1 \le i \le N$). Thus, both potentials can be approximated from a linear combination of such shape functions

$$ \mathbf{A} = \sum_{i=1}^{N} \mathbf{A}_i \phi_i(\mathbf{x}) \quad \varphi = \sum_{i=1}^{N} \varphi_i \phi_i(\mathbf{x}). \quad (7) $$

The only difference between FEM and NEM lies on the expression of the shape functions that in both methods are defined in a different manner. The ones related to the FEM are well known and therefore they do not need additional explanations. On the contrary, as the ones related to the NEM being less known, a brief introduction seems to be suitable. The NEM shape functions $\phi_i(\mathbf{x})$ have been introduced by Sibson and are calculated using the Voronoï diagram [1]. We are illustrating the concepts related to the construction of such shape functions by considering the simple distribution of five nodes depicted in Fig. 1. Some basic constructions, as the Delaunay triangulation and the associated Voronoi diagram, that serve to define the NEM shape functions in 2-D will be illustrated.

#### A. Voronoï Diagram

We define the Voronoï diagram associated to a set of $N$ nodes $S = \{n_1, n_2, n_3, \dots, n_N\}$ as the geometric subdivision of the

<figure><img src="image_3.png" /><figcaption>Fig. 2. Construction of the Voronoï diagram.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Shape function calculation related to the node 5, at 2 different positions $\mathbf{x}, T_{x5}$ being grey area.</figcaption></figure>

plane into regions $T_i$, called first order Voronoï cells, each one related to a distinct node $n_i$, such that

$$ T_i = \{\mathbf{x} \in \mathbb{R}^2 : d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j), \forall j \neq i \}, \forall i \quad (8) $$

with $\mathbf{x}_i$ the coordinates of the node $n_i$ and $d(\mathbf{x}, \mathbf{x}_j)$ the distance between node $n_j$ and point $\mathbf{x}$. Thus, $T_i$ is the region of the plane that contains the points $\mathbf{x}$ closest to node $n_i$ than to any other node in $S$. In Fig. 1 the nodes are numbered from "1" to "5" and the Voronoï cell related to node 5 is depicted. The construction of that Voronoï cell is quite simple. We start identifying the natural neighbors [1] of node 5. Now, we consider the segments joining node 5 to each neighbor node, and a straight line, normal to each one of these segments, traced at the central point of those segments as depicted in Fig. 1. These straight lines defined the Voronoï cell. Fig. 2 illustrates the construction of the Voronoï cells related to the other nodes. For additional details the interested reader can refer to [1] and the references therein.

#### B. Sibson Shape Functions

Once the Voronoï diagram is constructed, we can calculate the value of the NEM shape function $\phi_i(\mathbf{x})$, related to the node $i$ at the position $\mathbf{x}$. For this purpose we consider, in the scenario depicted in Fig. 3, the insertion of a point $\mathbf{x}$ and the computation of the shape function related to node 5 at that position, that is, $\phi_5(\mathbf{x})$. We consider two possible placements of such
<figure><img src="image_5.png" /><figcaption>Fig. 4. (a) Voronoï diagram (heavy lines), Delaunay triangulation (thin lines) and the associated Delaunay circumcircles. (b) Sibson shape function.</figcaption></figure>

point [Fig. 3(a) and Fig. 3(b)]. After introduction of point $x$, the Voronoï diagram is modified because the consideration of the cell related to the point just introduced that in our example consists of the polygon “abcd”. The cell just inserted, the one related to point $x$, intersects the former cells related to the original nodes. In Fig. 3 we emphasize in grey color the intersection between this cell and the one initially related to node 5. Now, the shape function related to node 5 at position $\mathbf{x}$ is computed from the ratio between the intersected area (the grey area in Fig. 3) and the one related to new cell (polygon “abcd”), that is

$$ \phi_5(\mathbf{x}) = \frac{\text{Area}(\text{"ebcf"})}{\text{Area}(\text{"abcd"})} \qquad (9) $$

It is easy to prove [1] that the support of the shape function related to a node $n_i$ is given by the union of the circles passing through this node and any other two neighbor nodes of $n_i$ such that no other node falls within the circle. Fig. 4(a) depicts those circles and the corresponding support for node 5. By computing $\phi_5(\mathbf{x})$, by using the procedure just described at different positions $\mathbf{x}$, the shape function can be represented. Fig. 4(b) depicts the shape function related to node 5. One of the reasons of the robustness of the Sibson interpolation against the distortion of the background mesh defining the Delaunay triangles and the associated Voronoï diagram, is that even if a triangle becomes very distorted, the support of the shape functions is not restricted to this triangle (as it is the case in the finite element method). Nowadays there is not a mathematical proof of the higher robustness of this kind of interpolation; however, there are numerous examples showing the low sensitivity of NEM to the mesh distortion [7]. Several simulations in solid and fluid mechanics [2], [3] have been successfully performed using the same cloud of nodes throughout the whole simulation.

### C. Nodal Sibson Shape Function Properties

The NEM shape functions $\phi_i(\mathbf{x})$ are $C^\infty$ except at the positions of nodes where they are $C_0$ and on the boundary of the Delaunay circles where they are $C^1$. The shape functions have a compact support delimited approximately by the neighbour nodes, leading to a sparse system. The band width of the resulting linear systems is a little bit higher than the bandwidth associated to the FEM. That is precisely one of the reason of the

<figure><img src="image_6.png" /><figcaption>Fig. 5. Studied system: two magnets surrounded by an air box.</figcaption></figure>

natural element interpolation robustness. The shape functions have other interesting properties

$$ \phi_i(\mathbf{x}_j) = \delta_{ij} \sum_{i=1}^{N} \phi_i(\mathbf{x}) = 1 \qquad \sum_{i=1}^{N} \phi_i(\mathbf{x})\mathbf{x}_i = \mathbf{x} \quad (10) $$

where $\mathbf{x}_j$ is the position of node $j$ and $\delta_{ij}$ the Kronecker's delta. The so-called delta property (the first one in (10)) enables to enforce essential boundary conditions as in FEM. The other two properties in (10) prove the linear completeness, that is, the Sibson interpolation can represent exactly linear functions. Moreover, the partition of unity, second relation in (10) allows functional enrichments in the context of the partition of unity paradigm.

The computing time related to the construction of the Sibson shape functions is higher than the one associated with the definition of such functions in the finite element framework. However, nowadays, there are fast algorithms coming from the computational geometry community that are robust, fast even in 3-D, and some of them freely distributed. If we take the CPU time as comparison criterion and consider simulations involving massively remeshing operations (that define the natural context of application of meshless methods), finite elements simulations are nowadays faster than natural element method, because the existence of efficient and robust remeshing procedures in 2-D. If one consider convergence and accuracy criteria, the convergence rate for linear FEM and NEM is the same (first order in the energy norm) but as we proved in some of our former works in mechanics [8] (and we will prove also in magnetostatics later in this paper) the NEM accuracy is higher than the one related to the finite elements. Moreover, in models requiring field transfer after each remeshing process the superiority of the NEM is more relevant as argued previously.

## IV. APPLICATIONS

First, two magnets face to face are considered (Fig. 5). Ten meshes with an increasing number of nodes were used. Fig. 5 gives the Voronoï diagram associated to the clouds consisting of 107, 800 and 2200 nodes. For the FEM simulation, linear interpolation on the Delaunay triangles was used.

For each mesh, both potential formulations have been solved using both the NEM and the FEM. Fig. 6 depicts the evolution of energy with the number of nodes. We can notice that the energies calculated with both formulations bound the exact energy. However, the energies computed from NEM seem to be more accurate than the ones coming from FEM calculations.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of the energy with the number of nodes.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Global relative error estimator.</figcaption></figure>

By representing the evolution of the numerical error $\hat{e}^2$ computed from (6) versus the number of nodes in a log-log scale (see Fig. 7) for both the finite element and the natural element solutions, we notice that both evolutions are linear and that the convergence rate, slope of such evolutions, is the same and approximately equal to one. We can also notice the higher accuracy of the NEM solution.

A more complex model in now considered. It concerns the variable reluctance machine (VRM) depicted in Fig. 8. The modelling of VRM is a very trick problem because of regions with sharp angles where the concentration of the magnetic flux density is very dramatic. Some simple error indicators are available [9] allowing the nodal refinement. The relative permeability of the rotor and stator yokes is equal to 1000. Different nodal distributions have been tested (the coarsest one containing 6000 nodes and the finest one more than one million) and both potential formulations solved by using the NEM and the FEM. The errors $e_{NEM}$ and $e_{FEM}$ have been computed from (5). The error ratio ($e_{FEM}/e_{NEM}$) evolves from a value close to 2 related to the coarsest nodal distribution, to 2.5 for the finest one.

## V. CONCLUSION

In this paper, the accuracy of meshless natural element method has been compared to the FEM one. For the comparison purposes we proposed the use of an error estimator

<figure><img src="image_9.png" /><figcaption>Fig. 8. Variable Reluctance Machine. The windings that have been supplied for the test is emphasized.</figcaption></figure>

based on the non verification of constitutive equation. For the considered examples, superiority in terms of accuracy of the NEM solutions has been noticed.

## REFERENCES

[1] N. Sukumar, “The natural element method in solid mechanics,” Ph.D. dissertation, Northwestern Univ., Evanston, IL, Jun. 1998.

[2] M. A. Martinez, E. Cueto, I. Alfaro, M. Doblaré, and F. Chinesta, “Up-dated Lagrangian free surface flow simulations with the natural neighbour Galerkin methods,” *Int. J. Numer. Meth. Eng.*, vol. 60, no. 13, pp. 2105–2129, 2004.

[3] I. Alfaro, D. Gonzalez, D. Bel, E. Cueto, M. Doblare, and F. Chinesta, “Recent advances in the meshless simulation of extrusion and other related forming processes,” *Archives of Computational Methods in Engineering, State of the Art Reviews*, vol. 13, no. 1, pp. 3–44, 2006.

[4] L. Illoul, J. Yvonnet, F. Chinesta, and S. Clénet, “Application of the natural element method to model moving electromagnetic device,” *IEEE Trans. Magn.*, vol. 42, no. 2, pp. 727–730, 2006.

[5] J. L. Synge, “The method of the hypercircle in function-space for boundary-value problems,” *Proc. Roy. Soc. London. Series A, Math. Phys. Sci.*, vol. 191, no. 1027, pp. 447–467, Dec. 3, 1947.

[6] J. Rikabi, C. F. Bryant, and F. Freeman, “An error based approach to complementary formulations of static field solutions,” *Int. J. Numer. Meth. Eng.*, vol. 26, pp. 1963–1987, 1988.

[7] J. Yvonnet, D. Ryckelynck, P. Lorong, and F. Chinesta, “The meshless constrained natural element method (C-NEM) for treating thermal models involving moving interfaces,” *Int. J. Thermal Sci.*, vol. 44, pp. 559–569, 2005.

[8] D. Misoum-Benziane, D. Ryckelynck, and F. Chinesta, “A new fully coupled two-scales modelling for mechanical problems involving microstructure,” *Comput. Meth. Appl. Mech. Eng.*, vol. 196, pp. 2325–2337, 2007.

[9] J. Yvonnet, D. Ryckelynck, G. Coffignal, P. Lorong, and F. Chinesta, “A simple error indicator for meshfree methods based on natural neighbors,” *Comput. Struct.*, vol. 84, no. 21, pp. 1301–1312, 2006.