# Application of the Natural Element Method to Model Moving Electromagnetic Devices

Lounès Illoul, Julien Yvonnet, Francisco Chinesta, Stéphane Clénet

► To cite this version:

Lounès Illoul, Julien Yvonnet, Francisco Chinesta, Stéphane Clénet. Application of the Natural Element Method to Model Moving Electromagnetic Devices. IEEE Transactions on Communications, 2006, 42 (4), pp.727-730. 10.1109/TMAG.2006.871587. hal-01005258

HAL Id: hal-01005258

https://hal.science/hal-01005258v1

Submitted on 26 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 CC0 1.0 - Universal - International License
# Application of the Natural-Element Method to Model Moving Electromagnetic Devices

L. Illoul¹, J. Yvonet¹, F. Chinesta¹, and S. Clénet²

¹LMSP, ENSAM, CER Paris, 75013 Paris, France

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

The natural-element method, which belongs to the family of meshless methods, is applied in the context of two-dimensional magnetostatics with moving parts. The method is reviewed and its interest for handling discontinuities in electromagnet devices with moving parts is illustrated through a numerical example.

**Index Terms**—Magnetostatics, movement, natural-element method (NEM).

## I. INTRODUCTION

EFFICIENT methods have been already developed in two dimensions to handle the motion of electromagnetic devices with rotating parts in finite-element model (FEM) [1], [2]. When the motion becomes complex (combination of rotation and translation), it leads to large mesh deformations. A remeshing step is then required, which increases the complexity of the method and can also induce noise due to the gap between the field projections from the old to the new meshes. Meshless methods provide new tools to take motion of interfaces into account. Nevertheless, in most mesh-free methods, imposition of essential boundary conditions is a difficult task. The natural-element method (NEM), which is a great interest in the domain of mechanics, enables us to easily overcome the aforementioned issues [3], [4]. In this paper, we propose the use of the NEM to model an electromagnetic device involving moving parts. First, the magnetostatics problem is summarized. Then, the NEM is reviewed. Finally, the interest of the NEM to handle moving discontinuities is evaluated in the context of a simple electromagnetic problem.

## II. MAGNETOSTATICS PROBLEM

Let denote $\mathbb{D}$ the domain studied and $\Gamma$ its surface. In magnetostatics, the magnetic field $\mathbf{H}$ and the magnetic flux density $\mathbf{B}$ verify

$$ \mathbf{curl} \mathbf{H} = \mathbf{J} \quad \mathrm{div} \mathbf{B} = 0 \qquad (1) $$

where $\mathbf{J}$ is the current density whose distribution is assumed known in $\mathbb{D}$. Both fields are related by the constitutive equation

$$ \mathbf{B} = \mu \mathbf{H} \qquad (2) $$

with $\mu$ the permeability. On the surface $\Gamma$, the following boundary conditions are prescribed:

$$ \mathbf{n} \cdot \mathbf{B} = 0 \text{ on } \Gamma_b \quad \mathbf{n} \times \mathbf{H} = 0 \text{ on } \Gamma_h \qquad (3) $$

with $\Gamma_b$ and $\Gamma_h$, two complementary parts of $\Gamma$, and $\mathbf{n}$, the outward normal unity vector defined on $\Gamma$. In magnetostatics, the

problem is generally solved using the vector potential $\mathbf{A}$ which satisfies

$$ \mathbf{B} = \mathbf{curl} \mathbf{A} \text{ with } \mathbf{n} \times \mathbf{A} = \mathbf{0} \text{ on } \Gamma_b \qquad (4) $$

The vector potential $\mathbf{A}$ has only one component in two dimensions, $\mathbf{A}(\mathbf{x})$ with $\mathbf{x}$ the coordinates of a point $\mathbf{x}$ in the domain $\mathbb{D}$. The equation to be solved then reads

$$ \frac{\partial}{\partial x} \left( \mu \frac{\partial \mathbf{A}}{\partial x} \right) + \frac{\partial}{\partial y} \left( \mu \frac{\partial \mathbf{A}}{\partial y} \right) = -\mathbf{J} \quad \text{and} \quad \mathbf{A} = \mathbf{0} \text{ on } \Gamma_b. \qquad (5) $$

## III. DISCRETISATION

A set of $N$ nodes $\{n_1, \dots, n_N\}$, with coordinates $\mathbf{x}_1, \dots, \mathbf{x}_N$, respectively, is distributed on the whole domain $\mathbb{D}$. A scalar function $\phi_i(\mathbf{x})$ is associated with each node in $\mathbb{D}$. The vector potential $\mathbf{A}$ is then approximated in the usual form

$$ \mathbf{A}(\mathbf{x}) = \sum_{i=1}^{N} A_i \phi_i(\mathbf{x}) \qquad (6) $$

with $A_i$ representing the nodal degrees of freedom which are obtained solving the linear system obtained by applying the Galerkin method to the weak formulation related to (5). In the FEM context, the nodal degree of freedom $A_i$ corresponds with the vector potential at node $n_i$. The interpolation function $\phi_i(\mathbf{x})$ verifies the Kronecker's delta property, i.e. $\phi_i(\mathbf{x}_j) = \delta_{ij}$. Thus, the essential boundary conditions can be imposed by prescribing the vector potential at the nodes located on $\Gamma_b$ (the shape functions related to internal nodes vanish on the domain boundary). In most mesh-free methods, it is not possible to directly impose essential boundary conditions because the shape functions do not satisfy the Kronecker property, and the shape functions related to the interior nodes do not vanish on the domain boundary. In the next section, we show how the NEM enables us to overcome this issue and allows handling discontinuities in the domain in a very simple manner.

In recent years, mesh-free methods for the solution of partial differential equations have significantly matured and have been used in various fields of engineering science. One of the reasons for this development is the fact that mesh-free discretizations and particle models are often better suited to
cope with geometrical changes of the domain of interest than
mesh-based discretization techniques, such as finite differences,
finite elements, or finite volumes techniques. Furthermore, the
computational costs associated with mesh generation are al-
leviated in mesh-free approaches since they are based only
on a cloud of nodes without any geometrical restriction with
respect to their distribution, in contrast with the finite-element
method, where elements that are too distorted lead to poor
accuracy. The so-called meshless or mesh-free methods have
been investigated and used by many researchers for treating a
large variety of engineering problems, usually involving large
displacements as encountered, for example, in forming process
simulations (free surface problems, moving boundaries prob-
lems, moving interfaces, and cracks propagation, among many
others). In these problems, accurate finite-element solutions
require significant computational efforts in remeshing steps.
In contrast, meshless methods require only nodal data without
explicit connectivity between nodes. The first generation of
meshless methods, the smooth particle hydrodynamics (SPH)
method, originally proposed by Lucy [5], introduces kernel
approximations and imposes the conservation laws using a
collocation technique. Recently, other mesh-free methods have
been proposed, such as the element-free Galerkin (EFG) [6],
the reproducing kernel particle method (RKPM) [7], and the
NEM [8], among many others, and they have been success-
fully used to discretize the Galerkin weak form associated
with the different models. Several methods for the analysis of
discontinuities or interfaces have been developed in both mesh
and meshless techniques. The most common approach lies in
explicitly tracking the interface motion. Within the interface
tracking approach, two main alternatives exist: the moving
mesh methods and the mixed Eulerian–Lagrangian methods.
Moving finite-element mesh methods conform element bound-
aries to the interface as it evolves. Although these methods
are very accurate, they are limited by severe mesh distortion.
Thus, frequent remeshing is needed, with the associated field
projections between successive meshes. Moreover, remeshing
is, even today, a delicate task in three dimensions. To alleviate
remeshing efforts a number of Eulerian–Lagrangian methods
have been developed recently that track the interface while
solving the equations on a fixed grid. Many of these methods
effectively smear the discontinuity over a few grid cells, and,
therefore, are not capable of accurately representing the solution
across the interface. A new approach for representing localized
behaviors has recently emerged in the field of the finite-element
method, known as extended finite-element method (X-FEM),
in which the interface evolution can be properly represented
on a fixed background mesh, just by adding an appropriate
enrichment in the functional approximation in the elements
that are intersected by the moving discontinuity. Alternative
to the aforementioned methods, the constrained NEM [4] has
been proposed to handle discontinuities and moving interfaces
in mesh-free methods in a very simple manner. The technique
uses the features of the NEM [8], [3] in tandem with a visibility
criterion. Unlike in most mesh-free methods, the constructed
approximation is strictly conforming and continuous across
interfaces (with discontinuous normal derivative), which avoids
the use of discontinuous enrichment techniques.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Construction of the Sibson shape functions. On left, the point x is located inside the domain. On right, the point x is located on the boundary.</figcaption></figure>

IV. NEM

In this section, after a brief review of the Voronoï-based in-
terpolants, the utility of the technique for handling moving in-
terfaces is presented. We briefly touch upon the foundation of
Sibson’s natural neighbor coordinates (shape functions) that are
used in the NEM. For a more in-depth discussion on the Sibson
interpolant and its application for solving second-order partial
differential equations, the interested reader can refer to [8] and
the references therein. The NEM interpolant is constructed on
the basis of the Voronoï diagram. The Delaunay tessellation is
the topological dual of the Voronoï diagram.

Consider a set of nodes S = {n₁, n₂, ..., nN} ∈ ℝ² [for the sake of simplicity, we will consider only two-dimensional (2-D) models, but all the results can be directly extended to the three-dimensional (3-D) case]. The Voronoï diagram is the subdivision of ℝ² into regions Tᵢ (Voronoï cells) defined by

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

where d(.) denotes the Euclidean distance. Thus, points in Tᵢ
are closer to the node nᵢ than to any other node in the cloud. Let
consider a point x inside the domain studied (Fig. 1). To calcu-
late the value of the Sibson coordinates φₖᵢ(x) of x with respect
to a natural neighbor nᵢ (see Fig. 1), the set Sₓ of neighbor nodes
is determined. Nodes nᵢ, nⱼ, nₖ are natural neighbors of x if the
point x is within the circumcircle of the triangle (nᵢ, nⱼ, nₖ)
coming from the Delaunay tessellation. In Fig. 1, only natural
neighbors of x are represented. The second-order Voronoï cell
Tᵢⱼ is given by

$$
T_{ij} = \{ \mathbf{y} \in \mathbb{R}^2 : d(\mathbf{y}, \mathbf{n}_i) < d(\mathbf{y}, \mathbf{n}_j) \\
\qquad < d(\mathbf{y}, \mathbf{n}_k), \forall k \neq i, \forall k \neq j \}. \quad (8)
$$

In Fig. 1, the second-order cell T<sub>xn1</sub> related to the point x
and the node n<sub>1</sub> is the polygon (afghe). If we denote K<sub>ni</sub>(x), the
surface of the cell T<sub>xni</sub>, the Sibson coordinate φ<sub>ni</sub>(x) is given
by

$$
\phi_{\text{ni}}(\mathbf{x}) = \frac{K_{\text{ni}}(\mathbf{x})}{\sum_{\mathbf{n}_{\text{j}} \in S_x} K_{\text{jn}}(\mathbf{x})}. \quad (9)
$$

For example, in the case of the Fig. 1, $\phi_{\text{ni}}(\mathbf{x})$ is given by

$$
\phi_{n1}(\mathbf{x}) = \frac{\text{Area}(\text{afghe})}{\text{Area}(abcde)}. \quad (10)
$$

If the point x coincides with the node nᵢ, i.e., x = xᵢ, then
φᵢ(xᵢ) = 1, and all other shape functions are zero, i.e., φᵢ(xⱼ) =
$\delta_{ij}(\delta_{ij}$ being the Kronecker delta). If a node $n_i$ is not a neighbor of the point $\mathbf{x}$, the function $\phi_i(\mathbf{x})$ vanishes. The support of the function $\phi_i(\mathbf{x}_i)$ is bounded.

The properties of positivity, interpolation, and partition of unity are then verified [8]

$$0 \le \phi_i(\mathbf{x}) \le 1 \text{ and } \phi_i(\mathbf{x}_j) = \delta_{ij} \text{ and } \sum_{i=1}^{n} \phi_i(\mathbf{x}) = 1. \quad (9)$$

The natural neighbor shape functions also satisfy the local coordinate property, namely

$$\mathbf{x} = \sum_{i=1}^{n} \phi_i(\mathbf{x})x_i \quad (10)$$

which combined with (9), implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness). Sibson natural neighbor shape functions are $C^1$ at any point except at the nodes, where they are only $C^0$. The $C^1$ continuity can be improved by using special classes of natural neighbor shape functions. Another important property of this interpolant is its strict linearity over the boundary of convex domains [8]. An illustration is depicted in Fig. 1 (right): As the areas of the Voronoï cells associated with the nodes located on the domain boundary become infinite, the contribution of internal nodes (with a finite area) vanish in the limit when the point approaches the convex boundary, and the shape functions associated with nodes $n_1$ and $n_2$ become linear on the segment $n_1 - n_2$. This is not true in the case of nonconvex boundaries, but the use of a special treatment such as a visibility criterion in the context of the NEM generalizes this property to any kind of domain boundary. The resulting technique (which involves a visibility criterion to limit the number of neighbor nodes) is called the constrained natural-element method (C-NEM) [4]. To handle fixed or moving interfaces the C-NEM proceeds by defining the different interfaces as a collection of linear segments (or triangular facets in three dimensions), which can be viewed as opaque boundaries. Now, the C-NEM functional interpolation is defined by considering at point $\mathbf{x}$ the contribution of all neighbor and visible nodes (the interfaces being considered as opaque). Thus, we can write

$$A(\mathbf{x}) = \sum_{i=1}^{V} A_i \phi_i^c(\mathbf{x}) \quad (12)$$

where $V$ is the number of natural neighbors visible from point $\mathbf{x}$. The computation of the constrained natural neighbor shape functions $\phi_i^c(\mathbf{x})$ is similar to the natural neighbor shape function, when one proceeds using the constrained Voronoï diagram [4]. Due to the inherent meshless character, these interfaces can be added in arbitrary clouds of nodes and the nodes defining the interface can move freely on the background cloud without any geometrical restriction. In this context, the constructed approximation is strictly continuous across any interface, being its normal derivative discontinuous. To illustrate this behavior, we consider the situation depicted in Fig. 3, where the point $\mathbf{x}$ moves from $\Omega_1$ to $\Omega_2$. If $\mathbf{x}$ is in $\Omega_1$ then the interpolated field is constructed using the neighbor visible nodes from point $\mathbf{x}$ ($\Gamma_I$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Sibson shape function related to a node surrounded by eight other nodes.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Interface tracking: reproducing transmission conditions.</figcaption></figure>

is assumed opaque). If $\mathbf{x}$ is on $\Gamma_I$, according to the previous discussion, the interpolated field is strictly linear because it only depends on the two neighbor nodes located on $\Gamma_I$. Finally, when $\mathbf{x}$ is in $\Omega_2$, the interpolated field is defined using the neighbor and visible nodes from point $\mathbf{x}$ (assumed opaque). The continuity of the interpolated field is then guaranteed, but a discontinuity appears in the normal derivative, because of a sudden change in the neighbor nodes across the interface.

## V. MOVEMENT

In the framework of the NEM, updating, addition or suppression of nodes is quite easy. So, to move a part $D_P$ in $D$, we proceed as follows. First, nodes are distributed on the whole domain $D$ without considering the moving part $D_P$. Then, the moving part is added to the domain with all its associated nodes. Then, whatever the position of $D_P$ in $D$, all the nodes of $D$ shadowing by $D_P$ are suppressed. The problem is solved with the remaining nodes of $D$ and the nodes of $D_P$. This technique is
<figure><img src="image_5.png" /><figcaption>Fig. 4. Electromagnet under study and the distribution of nodes used for the discretization with the NEM.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the flux versus the angle.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of the attractive torque versus the angle.</figcaption></figure>

simple to carry out whatever the movement (translation, rotation, or a combination of both).

## VI. NUMERICAL EXAMPLE

The electromagnet depicted in Fig. 4, which involves a moving part, is analyzed. The static part is composed of a ferromagnetic core with two coils in series. The moving part is made up ferromagnetic material. The relative permeability of the ferromagnetic material is constant and equal to 1000. The nodal distribution used for the discretization of the electromagnet is illustrated in Fig. 4. When the coil is supplied and the moving part rotates around the point 0 (see Fig. 4), we have calculated the flux $\Phi$ (Fig. 5) and the attractive torque (Fig. 6). The torque $T$ is calculated using the Maxwell stress tensor. The results obtained are very coherent because the torque $T$ is close to the torque $T'$ calculated using another method based on the

derivative of the coenergy (i.e., by derivating the flux versus the angle)

$$T'_{n} = \frac{1}{2} i \frac{\Phi_{n+1} - \Phi_{n-1}}{2\Delta\theta} \quad (12)$$

where $T'_n$ is the torque, $\Phi_n$ is the flux, both calculated at the n'th angular position, and $\Delta\theta$ is the angle step. Concerning the computational effort, it is not possible at present to compare of both the FEM and the NEM because due to the novelty of the last one (i.e., the numerical algorithms are not yet optimized). The construction of the NEM shape functions are more expensive that those of the FEM; however, the accuracy of the NEM is significantly higher and no remeshing is required when one uses the NEM (remeshing procedures are expensive in two dimensions and very delicate in three dimensions). In a former work [9], we have proposed a new fast algorithm for computing the NEM shape functions that is not much more expensive that the one used in the FEM context.

## VII. CONCLUSION

In this work, we have proved the ability of a novel meshless method, the C-NEM, for simulating electromagnetic devices involving moving interfaces in two dimensions. All the constructions described in this paper can be directly extended to the 3-D case for nodal interpolation [10]. We are conscious that the extension to the 3-D treatment (i.e., edge elements) of electromagnetic models requires more in deep developments.

## REFERENCES

[1] B. Davat, Z. Ren, and M. Lajoie-Mazenc, "The movement in field modeling," *IEEE Trans. Magn.*, vol. MAG-21, no. 6, pp. 2296-2298, Nov. 1985.

[2] B. Boualem and F. Piriou, "Numerical models for rotor cage induction machines using finite element method," *IEEE Trans. Magn.*, vol. 34, no. 5, pp. 3202-3205, Sep. 1998.

[3] E. Cueto, N. Sukumar, B. Calvo, J. Cegoñino, and M. Doblaré, "Overview and recent developments in natural neighbor galerkin methods," *Arch. Comput. Meth. Eng.*, vol. 10/4, pp. 307-384, 2003.

[4] J. Yvonnet, D. Ryckelynck, P. Lorong, and F. Chinesta, "A New extension of the natural element method for nonconvex and discontinuous problems: The constrained natural element method (C-NEM)," *Int. J. Numer. Meth. Eng.*, vol. 60/8, pp. 1451-1474, 2004.

[5] L. B. Lucy, "A Numerical approach to the testing of fusion process," *Astron. J.*, vol. 88, pp. 1013-1024, 1997.

[6] T. Belytschko, Y. Lu, and L. Gu, "Element-free galerkin methods," *Int. J. Numer. Meth. Eng.*, vol. 37, pp. 229-256, 1994.

[7] W. K. Liu, S. Jun, and Y. F. Zhang, "Reproducing kernel particle methods," *Int. J. Numer. Meth. Fluids*, vol. 21, pp. 1081-1106, 1995.

[8] N. Sukumar, B. Moran, and T. Belytschko, "The natural element method in solid mechanics," *Int. J. Numer. Meth. Eng.*, vol. 43, pp. 839-887, 1998.

[9] J. Yvonnet and F. Chinesta, "An Hybrid element free galerkin and natural element meshfree method for direct imposition of boundary conditions and faster three-dimensional computations," presented at the 3rd MIT Conf. Computational Fluid and Solid Mechanics, Cambridge, MA, 2005.

[10] I. Alfaro, D. Bel, E. Cueto, M. Doblaré, and F. Chinesta, "Three-dimensional simulation of aluminum extrusion by the alpha-shape based natural element method," in *Proc. Comput. Meth. Appl. Mech. Eng.*, to be published.