# A mixed finite element/meshless natural element method for simulating rotative electromagnetic machines

Lounes Illoul, Yvonnick Le-Menach, Stéphane Clénet, Francisco Chinesta

► To cite this version:

Lounes Illoul, Yvonnick Le-Menach, Stéphane Clénet, Francisco Chinesta. A mixed finite element/meshless natural element method for simulating rotative electromagnetic machines. European Physical Journal: Applied Physics, 2008, 43, pp.197-208. 10.1051/epjap:2008102. hal-01007369

HAL Id: hal-01007369

https://hal.science/hal-01007369

Submitted on 5 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
# A mixed finite element/meshless natural element method for simulating rotative electromagnetic machines

L.A. Illoul<fnref n="1" />, Y. Le Menach<fnref n="2" />, S. Clenet<fnref n="3" />,ᵃ, and F. Chinesta<fnref n="1" />

**Abstract.** Efficient numerical models have been already developed in 2D to take into account the movement of electromagnetic devices with rotating parts in the framework of the Finite Element Method (FEM). When the movement becomes complex it leads to large mesh distortions. A remeshing step is then required which increases the computational complexity and can also lead, in some cases, to numerical ripples on forces and torques due to the field projections between old and new meshes. Moreover, remeshing procedures in 3D remain an open topic. Meshless methods seem an appealing choice for alleviating the mesh constraints. The Natural Element Method (NEM) which, has known a growing interest in the domain of mechanics, allows to proceed in the meshless framework, avoiding one of the main drawbacks related to the vast majority of meshless techniques, as is the imposition of essential boundary conditions. In this paper, a variant of the NEM, known as constrained natural element method (C-NEM) is applied for simulating electromagnetic machines involving rotating parts. A new mixed strategy combining the finite element and the constrained natural element methods is proposed and then tested by using an appropriate error estimator.

Magnetostatics; magnetic shielding, magnetic induction, boundary-value problems, Finite-element and Galerkin methods

# 1 Introduction

## 1.1 Numerical modelling of electromagnetic machines

To accurately model electrical machines, one needs to solve the Maxwell equations taking into account the movement of mobile parts. Generally, an electrical rotating machine is made up of two parts. The moving part (the rotor) and the static one (the stator) are separated by a thin clearance (the air gap) to enable the rotor movement. In fact, most of the energy of the machine is concentrated in the air-gap. Thus, it is necessary to take into account very accurately the relative movement between both parts.

At present, the Finite Element Method (FEM) is widely used to solve numerically the Maxwell equations. To take the movement into account, several methods have already been proposed: the macro element [18], the moving band [5], the slip surface [16], the use of Lagrange multipliers [8], the overlapping technique [27], the nodal interpolation strategy [15], the mortar element [17], among many others. Meshless techniques could be envisaged for alleviating the computational complexity related to these approaches.

The main aim of the present work is the analysis of meshless approaches as well as the proposal and evaluation of a new mixed strategy combining the simplicity of finite element descriptions with the robustness of meshless descriptions in the regions involving large mesh deformations.

## 1.2 Meshless approaches

In the past years meshfree methods for the solution of partial differential equations have significantly matured and are used in various fields of applications. One of the reasons for this development is the fact that meshfree 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. Furthermore, the computational cost associated with mesh generation is alleviated in meshfree approaches since they are based only on a cloud of nodes, without any geometrical restriction concerning their relative positioning.

Among the main meshless approaches we highlight the Smooth Particle Hydrodynamics (SPH) [10], the Diffuse Element Method (DEM) [13], the Element Free

<footnote n="1">¹ LMSP, UMR 8106 CNRS ENSAM-ESEM, 151 boulevard de l'Hôpital, 75013 Paris, France</footnote>

<footnote n="2">² L2EP, USTL, Bât. P2, 59655 Villeneuve d'Ascq Cedex, France</footnote>

<footnote n="3">³ L2EP, ENSAM CER Lille, 8 boulevard Louis XIV, 59046 Lille Cedex, France</footnote>
Galerkin (EFG) [1], the Reproducing Kernel Particle Method(RKPM) [9], the h-p Cloud Method [6] and the Natural Element Method (NEM) [25].

Some of these meshless techniques are still nowadays confronted to some difficulties related to: (i) the quality of the approximation as well as the conditioning of the global system is pathologically dependent on the size of the support of the shape functions; (ii) imposing of essential boundary conditions needs special treatments; (iii) integration is not accurate enough, and (iv) physical discontinuities across interfaces must be accurately enforced.

Several techniques for the analysis of problems involving discontinuities or interfaces have been proposed in the finite element [2,12] and in the meshless [14,28] frameworks. The extended finite element method (X-FEM) [12] makes use of the Partition of Unity paradigm [11,12] for approximating discontinuities. The goal of the extended finite element method is to expand the finite element approximation functions for including some features of the solution that are known.

Alternative to the aforementioned methods, the constrained natural element method (C-NEM) [29] has been proposed to handle discontinuities and moving interfaces in meshfree methods in a very simple manner [30]. The technique uses the features of the natural element method in tandem with a visibility criterion. In opposition to the vast majority of meshless strategies, the C-NEM allows to enforce exactly essential boundary conditions on the domain boundary, and then to account for transmission conditions across fixed or moving interfaces.

In this paper, the constrained natural element method will be applied to model rotating electrical machines. For this purpose, we first summarize the constrained natural element method (C-NEM). Then, the equations to be solved are given and the method to take the movement into account is described. Finally, a Variable Reluctance Machine (VRM) is modelled using the constrained natural element method and results are compared with results obtained by using the FEM. A mixed formulation combining the versatility of finite elements and the robustness of the meshless approaches to address large mesh deformations is then proposed and evaluated.

For comparison purposes we use the error estimator introduced in [3] in the natural element method framework, which is based on the solution of both complementary potential formulations. It will be applied to quantify the accuracy of the different solutions as well as to compare C-NEM and FEM solutions. A mixed approach combining the versatility of finite elements and the robustness of the meshless approaches to address large mesh deformations will be also proposed and its accuracy checked. This comparison is one of the main original contributions of this paper, an one of the first exhaustive comparative analysis of both approximations to take into account movement, at least in our knowledge.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Voronoi diagram.</figcaption></figure>

## 2 The meshfree constrained natural element method (C-NEM) revisited

In this section, the ability of the C-NEM for describing moving interfaces and discontinuities in a fixed cloud of nodes is discussed. After a brief review of the Voronoi-based interpolants, we introduce the constrained Voronoi diagram which is used to compute the shape functions in any domain. To avoid duplication with some of our former publications, different references to our former works will be addressed.

### 2.1 Natural neighbor interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. 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 Sambridge and Braun [19], and Sukumar et al. [25]. The NEM interpolant is constructed on the basis of the Voronoi diagram (see Fig. 1). The Delaunay tessellation is the topological dual of the Voronoi diagram.

From now on, we focus only in the 2D case, the 3D case being a direct extension. Consider a set of nodes $S = \{n_1, n_2, ..., n_N\}$ in $\mathbb{R}^2$. The Voronoi diagram is the partition of $\mathbb{R}^2$ into regions $T_i$ (Voronoi cells) defined $\forall i \in S$ by:

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

where $d(\cdot)$ denotes a distance. The Voronoi cells related to neighbor nodes have a common edge.

In order to define the natural neighbour coordinates it is necessary to introduce the second-order Voronoi diagram of the cloud defined for each couple of nodes $(i, j) \in S^2$ with $i \neq j$ as

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

Sibson [24] defined the natural neighbor coordinates of a point $\mathbf{x}$ with respect to one of its neighbors $n_i$ as the ratio of the cell $T_i$ that is transferred to $T_x$ when adding $\mathbf{x}$ to the
initial cloud of points, to the total volume of $T_x$. In other words, if $\kappa(\mathbf{x})$ and $\kappa_i(\mathbf{x})$ are the Lebesgue measures of $T_x$ and $T_{xi}$ respectively, the natural neighbor coordinates of $\mathbf{x}$ with respect to the node $n_i$ is defined as

$$ \phi_i(\mathbf{x}) = \frac{\kappa_i(\mathbf{x})}{\kappa(\mathbf{x})}. \qquad (3) $$

Figure 2 illustrates the construction of $\phi_1(\mathbf{x})$, that in this case is given by:

$$ \phi_1(\mathbf{x}) = \frac{\text{Area}(afghe)}{\text{Area}(abcde)}. \qquad (4) $$

If the point $\mathbf{x}$ coincides with the node $n_i$, i.e. ($\mathbf{x} = \mathbf{x}_i$), $\phi_i(\mathbf{x}_i) = 1$, and all other shape functions are zero, i.e. $\phi_j(\mathbf{x}_i) = \delta_{ij}$ ($\delta_{ij}$ being the Kronecker's delta). The properties of positivity, interpolation, and partition of unity are then verified [25]:

$$ \left\{ \begin{array}{l} 0 \le \phi_i(\mathbf{x}) \le 1 \\ \phi_i(\mathbf{x}_j) = \delta_{ij} \\ \sum_{i=1}^{N} \phi_i(\mathbf{x}) = 1. \end{array} \right. \qquad (5) $$

The natural neighbor shape functions also satisfy the local coordinate property [24], namely:

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

which combined with equation (5), 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 continuity can be enhanced by using special classes of natural neighbor shape functions [7].

Another important property of this interpolant is its strict linearity over the boundary of convex domains. The proof can be found in Sukumar et al. [25]. An illustration is depicted in Figure 2b: as the areas associated to points on the boundary become infinite, the contribution of internal nodes 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 joining both nodes ($S_{n_1 \to n_2}$). This result is essential to guarantee strict continuity of the approximation across material interfaces. The lack of this property is an important issue in most meshfree methods which require special numerical strategies to circumvent this drawback [14]. Nevertheless, the property of linearity of the interpolant is only satisfied along convex boundaries [25]. The non-convex case needs further developments that will be addressed in the next section.

Consider an interpolation scheme for a scalar function $u(\mathbf{x}) : \Omega \subset \mathbb{R}^2 \to \mathbb{R}$, in the form:

$$ u^h(\mathbf{x}) = \sum_{i \in \mathcal{I}} \phi_i(\mathbf{x}) u_i \qquad (7) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Construction of the Sibson shape functions.</figcaption></figure>

where $u_i$ are the nodal values of the field $u$ at the $n(\mathbf{x})$ neighbor nodes of point $\mathbf{x}$, all them grouped in the set $\mathcal{I}$, and $\phi_i(\mathbf{x})$ are the shape functions at that point associated with each neighbor node. It is noted that equation (7) defines a local interpolation scheme. Thus, the trial and test functions used in the discretization of the variational formulation related to the electromagnetic model treated in this paper could be approximated by equation (7).

## 2.2 The constrained natural element method

### 2.2.1 Constrained Voronoi diagram

It was proved in [25,29] and [4] that spurious influences between “non-visible” nodes and loss of linearity in the interpolation along boundaries of non-convex domains appear in the framework of the NEM. In order to avoid this drawback and to recover all properties of the method for any geometry (including non-convex domains containing fixed or moving interfaces), a visibility criterion was introduced in order to restrict influent nodes among natural neighbors [29]. The computation of the shape functions is then done on the basis of the so-called constrained (or extended) Voronoi diagram (CVD), introduced by Seidel in [22] which is the strict dual of the constrained Delaunay triangulation.

To illustrate this behavior, we consider the situation depicted in Figure 3, where the point $\mathbf{x}$ moves from $\Omega_1$ to $\Omega_2$. If $\mathbf{x}$ is in $\Omega_1$, the interpolated field is constructed from equation (7) using the neighbor visible nodes from point $\mathbf{x}$ ($\Gamma_I$ is assumed opaque). If $\mathbf{x}$ is on $\Gamma_I$ 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 nodes, visible from point $\mathbf{x}$ ($\Gamma_I$ being opaque). The continuity of the interpolated field is then guaranteed, but a
<figure><img src="image_4.png" /><figcaption>Fig. 3. Reproducing discontinuous derivatives using the constrained Voronoi diagram.</figcaption></figure>

discontinuity appears in the field derivatives, because of a sudden change in the neighbor nodes across the interface. Thus, we can then reproduce a continuous field with discontinuous normal derivative across the interface.

The intersection between the constrained Voronoi diagram and the domain closure is composed of cells $T_i^C$, one for each node $n_i$, such that any point $\mathbf{x}$ inside $T_i^C$ is closer to $n_i$ than to any other node $n_j$ visible from $\mathbf{x}$. The constrained Voronoi cells are defined formally $\forall i \in S$ by:

$$T_i^C = \{ \mathbf{x} \in \Re^n : d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j);$$

$$\forall j \in S, j \neq i; S_{x \to n_i} \cap \Gamma = \emptyset; S_{x \to n_j} \cap \Gamma = \emptyset \} \quad (8)$$

where $\Gamma$ is the domain boundary and $S_{a \to b}$ denotes the segment between the points $a$ and $b$.

All these constructions can be generalized to 3D [20,21].

### 2.2.2 The constrained natural element approximation

In order to solve partial differential equations defined on non-convex domains, we consider the following approximation for both the trial and the test functions:

$$\begin{align}
u^h(\mathbf{x}) &= \sum_{i \in \mathcal{I}_V} \phi_i^C(\mathbf{x}) u_i \nonumber \\
\delta u^h(\mathbf{x}) &= \sum_{i \in \mathcal{I}_V} \phi_i^C(\mathbf{x}) \delta u_i \tag{9}
\end{align}$$

where the set $\mathcal{I}_V$ contains the $V(\mathbf{x})$ natural neighbors nodes visible from point $\mathbf{x}$, and $\phi_i^C$ represents the constrained natural neighbor shape function related to the $i$th node at point $\mathbf{x}$. The computation of the C-n-n (constrained natural neighbor) shape functions is similar to the natural neighbor ones, when one proceeds using the constrained Voronoi diagram introduced previously. It was proved in [29] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation, allowing the extension of the linearity of the shape functions on the convex domain boundaries, to any geometry, convex or not.

The ability of the C-NEM to treat problems involving cracks was illustrated in [29] and for moving interfaces

in thermal problems in [30]. In the present paper, the domain is partitioned in some regions with different material properties. Each subdomain is discretized using a cloud of nodes and the interfaces between the different regions are described by a polygonal curve defined by a set of nodes. Then, a constrained Voronoi diagram is defined in each subdomain with respect to the domain boundary and the interfaces. The attractive feature of the present technique is the possibility to move the interfaces without special care for the shape of the underlying Delaunay triangles because the interpolation accuracy does not depend on the geometrical quality of the Delaunay triangles, in contrast to the FEM. In this manner, the continuity of the approximation is guaranteed by the strict linearity of the interpolation across the interfaces, that are defined by a set of interface nodes.

## 2.3 Taking into account movement

It is well known in the context of the FEM that frequent remeshing procedures are required for accounting accurately the rotation. To model the rotation of the electrical machine different methods are available. Remeshing can be very simple when using a regular mesh. Thus, one must shift the unknowns located on the boundary of the moving part. Nevertheless, the rotation step must be constant, that can be a drawback if one wants to model a machine working at variable speed.

To overcome this task, other methods are applied on a slip line using an interpolation technique to obtain a relation between the static and mobile unknowns. Thus, the nodal values on the slip line of the mobile part are expressed from the ones on the fixed part and then no additional unknowns are introduced. The interpolation can be linear, expressing a moving unknown from two static unknowns [15]. It can be also quadratic [23] where three static nodes are necessary to define the interpolation. Another more general approach is given by the “mortar” method [17] which allows expressing the moving unknowns as a function of all the static ones. The solution defined on the moving surface is projected on the fixed one using the Galerkin technique. It is obvious that the “mortar” method gives the best results but the size of the linear system bandwidth increases significantly, in contrast with the two first techniques which do not increase the linear system bandwidth but lead to less accurate
solutions, in particular the linear interpolation. Still on the slip line, Lagrange Multipliers can be considered [8]. This technique consists in introducing test functions for minimizing the gap in the nodal values (related to moving and static nodes) computed on the slip line. In this case all the nodal unknowns must be computed, and consequently the size of the linear system increases significantly and its conditioning is degraded, like when interpolation methods are used.

Other techniques use a special mesh in the air gap such as the moving band [5]. In this method, the mesh is distorted, due to the movement, until an unacceptable degree of deformation. Then, the surface is re-meshed and the nodal connectivity updated as in the locked step approach (regular mesh in the air gap). Unfortunately some numerical ripples can be observed on the torque when remeshing is used. For the macro element technique, an analytic solution is introduced in the air-gap [18]. However, the associated linear system bandwidth increases too and the analytic solution can be difficult to obtain for some kind of machines. Finally the overlapping method uses two meshes in the air-gap [27]. One mesh is associated to the static part and the whole air-gap region. A Dirichlet condition is imposed at the nodes of the mesh that belong to the mobile part. The same operation is carried out with a second mesh (that contains the mobile part and the air-gap) prescribing the field at the nodes located in the static domain. Thus, two different meshes cover the air gap. This method leads to very accurate results but the mesh computation in the airgap requires significant efforts.

In the C-NEM context, accounting for movement is quite natural and very simple to implement. One could just move the nodes defining the interface because the geometrical quality of the resulting Delaunay triangles does not affect significantly the accuracy of the resulting approximation. In the simulations reported in this work we decided to rotate all the nodes within the inner airgap circle. Then the constrained Voronoi diagram is updated and the interpolation shape functions computed again.

## 3 Simulating rotating electrical machines

To test the capability of the C-NEM to model rotating electrical machines, we consider a 6-4 Variable Reluctance Machine (VRM). The calculation of the torque for each relative position between stator and rotor, that is parameterized using the angle $\theta$, is a very tricky problem because for computing the torque (as we explain later) one must proceed in the air gap region where sharp angles due to the presence of rotor and stator teeth are encountered. If the mesh is not regular and fine enough, parasitic ripples can appear in the computed torque (as illustrated later) when one proceeds in the finite element framework.

In Figure 4, a 3D view of the machine is given, as well as the 2D model that is solved. The domain under study is a cross section of the machine defined as $\Omega = \{(x, y) : x^2 + y^2 < R_{ext}^2\}$, being $R_{ext}$ the stator external radius. The stator has 6 teeth and holds three phase windings. Each one is composed of two coils that are wound around

<figure><img src="image_5.png" /><figcaption>Fig. 4. (Color online) Scheme of a rotative electromagnetic machine: 3D and 2D views (units in metric system). In the bottom figure we can see the regions where the current density is non zero (in red and blue). Only one phase among the three phases is supplied.</figcaption></figure>

one tooth. In our study only one phase is supplied by a DC current. Two regions are then defined the first is flown by a current density $+J$ and the other by a current density $-J$ ($|J| = 1 \text{ A mm}^{-2}$) (see Fig. 4). Elsewhere the current density is equal to zero. The eddy currents are also neglected.

The electromagnetic field distribution can be calculated by using either the vector potential or the scalar potential formulations. To introduce the notations we consider the basic electromagnetic equations:

$$ \left\{ \begin{array}{l} \nabla \mathbf{B} = 0 \\ \nabla \times \mathbf{H} = \mathbf{J} \end{array} \right. \qquad (10) $$

where $\mathbf{B}$ is the magnetic flux density, $\mathbf{H}$ the magnetic field and $\mathbf{J}$ the current density. There is a constitutive equation that links $\mathbf{B}$ and $\mathbf{H}$ that writes

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

where $\mu > 0$ represents the magnetic permeability.
On the boundary $\partial\Omega$ of $\Omega$ we assume that the normal component of $\mathbf{B}$ vanishes, that is, $\mathbf{B}\mathbf{n} = 0$, where $\mathbf{n}$ is the unit outwards vector defined on $\partial\Omega$.

A field $\mathbf{B}$ is said admissible, $\mathbf{B}_{adm}$, if it verifies the first relationship in (10). Analogously, a field $\mathbf{H}$ is said admissible, $\mathbf{H}_{adm}$, if it verifies the second relationship in (10). Obviously a couple $(\mathbf{B}_{adm}, \mathbf{H}_{adm})$ that verifies the constitutive equation (11) is the exact solution of the electromagnetic problem.

In the 2D problem that we are analyzing we can assume that the magnetic flux density has a zero out-of-plane component, that is $B_z = 0$, that implies that it exists a field $\mathbf{A}$, with components $\mathbf{A}^T = (0, 0, A_z) = (0, 0, A)$ such that $\mathbf{B} = \nabla \times \mathbf{A}$. Thus the magnetic field can be defined everywhere using the constitutive equation (11): $\mathbf{H} = \nabla \times \mathbf{A}/\mu$, that introduced in the second relation of (10) leads to the 2D vector potential formulation:

$$ \left\{ \begin{array}{ll} -\nabla \left( \frac{1}{\mu} \nabla A \right) = J & \text{in } \Omega \\ A = 0 & \text{on } \partial \Omega \end{array} \right. \qquad (12) $$

where the vector current density $\mathbf{J}$ has only a z-component, i.e. $\mathbf{J}^T = (0, 0, J_z) = (0, 0, J)$, with both $A$ and $J$ independent on the z-coordinate.

Obviously, by construction, its solution allows to define an admissible magnetic flux density $\mathbf{B}_{adm}$ that will be referred as $\mathbf{B}_A$. The associated magnetic field writes $\mathbf{H}_A = \mathbf{B}_A/\mu$.

In the same way, another complementary potential formulation can be introduced. Let $\mathbf{H}_s$ be a field such that $\nabla \times \mathbf{H}_s = \mathbf{J}$. Since $\mathbf{H} - \mathbf{H}_s$ is curl free, we can introduce a scalar potential $\omega$ such that:

$$ \mathbf{H} = \mathbf{H}_s - \nabla \omega. \qquad (13) $$

Then, the equation to be solved can be written:

$$ \left\{ \begin{array}{ll} -\nabla(\mu \nabla\omega) = -\nabla(\mu \mathbf{H}_s) & \text{in } \Omega \\ \mu (\mathbf{H}_s - \nabla\omega) \mathbf{n} = 0 & \text{on } \partial\Omega. \end{array} \right. \qquad (14) $$

This formulation leads to the admissible magnetic field denoted by $\mathbf{H}_\omega$, but the associated magnetic flux density $\mathbf{B} = \mu\mathbf{H}_\omega$ is not admissible anymore because it does not verify $\nabla\mathbf{B} = 0$. It will be denoted in what follows by $\mathbf{B}_\omega = \mu\mathbf{H}_\omega$.

In the electromagnetic system here considered the permeability of the ferromagnetic part (stator and rotor) is assumed to be constant $\mu_s = \mu_r = 4\pi \times 10^{-4} \text{ H m}^{-1}$ (saturation and hysteresis effects are neglected) and the permeabilities of the air and of the coil are $\mu_a = \mu_c = 4\pi \times 10^{-7} \text{ H m}^{-1}$.

The weak formulation related to equation (12) results: Find $A \in H_0^1(\Omega)$ such that

$$ \int_{\Omega} \frac{1}{\mu} \nabla A \nabla A \, d\Omega = \int_{\Omega} J \, \delta A \, d\Omega, \forall \delta A \in H_0^1(\Omega) \quad (15) $$

where $H_0^1(\Omega)$ is the usual Sobolev functional space.

In the same way the weak formulation related to equation (14) writes:

$$ \int_{\Omega} \mu \nabla\omega \nabla\delta\omega \, d\Omega = \int_{\Omega} \mu \mathbf{H}_s \nabla\delta\omega \, d\Omega, \quad \forall \delta\omega \in H^1(\Omega). \qquad (16) $$

In Figures 5–7 we can observe the cloud of nodes considered for the discretization of the electromagnetic equations as well as the associated Voronoi cells. The interfaces between the different regions are described by a polygonal curve defined by a set of nodes (thick lines in Fig. 7). In order to appreciate the nodal refinement in the air gap between stator and rotor we depict in Figures 6 and 7 two zooms of this area. We can observe that a constrained Voronoi diagram has been associated with each subdomain $\Omega_s$ (stator), $\Omega_c$ (coil), $\Omega_a$ (air) and $\Omega_r$ (rotor) with respect to the corresponding interfaces. Thus, the continuity of $A$ and $\omega$ across the interfaces is guaranteed, allowing the discontinuity of the tangential component of the magnetic flux density across the interface, as discussed previously. Let $D_0$ be the disk centered on the shaft of the electrical machine and with a diameter equal to the diameter of the rotor as it was reported at the end of Section 2. For accounting a rotation $\Delta\theta$ of the rotor all nodes within $D_0$ rotates of that angle $\Delta\theta$. At each rotation step, the constrained Voronoi diagram has to be updated, allowing the computation of the approximation shape functions and then the discretization of the variational formulations (Eqs. (15) and (16)).

Thus, substituting the trial and test functions (both approximated in the C-NEM framework) in the above variational formulations and using the arbitrariness of the fields $\delta A$ and $\delta\omega$, the following linear system of equations is obtained for the vector potential formulation:

$$ \mathbf{K}_A(\theta) \mathbf{A}(\theta) = \mathbf{F}_A(\theta) \qquad (17) $$

and for the scalar potential formulation

$$ \mathbf{K}_\omega(\theta) \omega(\theta) = \mathbf{F}_\omega(\theta) \qquad (18) $$

where $\theta$ refers the relative position between stator and rotor, $\mathbf{A}(\theta)$ the vector containing the nodal values of the vector potential field for that relative position $\theta$, and $\omega(\theta)$ the vector containing the nodal values of the scalar potential at position $\theta$.

The associated torque can be calculated from:

$$ C(\theta) = \int_{0}^{2\pi} \bar{R} F_t(\theta, \bar{R}, \varphi) \, d\varphi \qquad (19) $$

where $\bar{R}$ is the radius of a circle completely defined in the air gap, and $F_t$ the tangential component of the electromagnetic force density $\mathbf{F}$ that is computed from the Maxwell stress tensor. Thus, during a complete turn, the torque varies and due to the localization of the potential gradients significant oscillations appear when the torque is computed using a coarse mesh in the finite element framework if no remeshing is performed, as we illustrate later.

Figure 8 depicts, for 4 relative positions between the rotor and the stator ($\theta_1 = 0^\circ$, $\theta_2 = 15^\circ$, $\theta_3 = 30^\circ$ and
<figure><img src="image_6.png" /><figcaption>Fig. 5. (Color online) 2D scheme of a rotative electromagnetic machine partitioned with constrained Voronoi cells.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. (Color online) 2D scheme of a rotative electromagnetic machine: detail.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. (Color online) 2D scheme of a rotative electromagnetic machine: detail.</figcaption></figure>

$\theta_4 = 45^{\circ}$), the flux path-map given by the vector potential formulation when the phase indicated in Figure 4 (bottom) is supplied with a DC current.

In order to compare the solutions related to different approximation schemes, we are using an appropriate error estimator based on the verification of the constitutive relationship. This error estimator was introduced in [26] and is often used to estimate numerical errors in static electromagnetism. This estimator was recently used to compare finite element and meshless natural element solutions in the case of magnetostatics (without taking into account movement) in [3]. In the present work we are applying this estimator in order to compare both methods when the movement is also taken into account. This estimator

<figure><img src="image_9.png" /><figcaption>Fig. 8. (Color online) Vector potential field $A(\mathbf{x})$ for four different rotor positions ($\theta_1 = 0^{\circ}$, $\theta_2 = 15^{\circ}$, $\theta_3 = 30^{\circ}$ and $\theta_4 = 45^{\circ}$).</figcaption></figure>

is defined by:

$$e_B = \|\mathbf{B}_{adm} - \mu\mathbf{H}_{adm}\| = \|\mathbf{B}_A - \mathbf{B}_{\omega}\| \quad (20)$$

where $\|\mathbf{v}\| = \int_{\Omega}(\mathbf{v}\mathbf{v})/\mu \, d\Omega$.

It was proven in [26] that this error corresponds to the distance between the averaged magnetic flux density $\bar{\mathbf{B}} = 0.5 (\mathbf{B}_{adm} + \mu\mathbf{H}_{adm})$ and the exact field $\mathbf{B}_{ex}$.

Obviously we can define an equivalent error norm operating on the magnetic fields:

$$e_H = \left\| \frac{\mathbf{B}_{adm}}{\mu} - \mathbf{H}_{adm} \right\| = \| \mathbf{H}_A - \mathbf{H}_{\omega} \|. \quad (21)$$

A more appropriate expression for the error consists of normalizing the previous expressions. In that follows a relative error will be computed by using

$$E_H^2 = 2 \frac{\| \mathbf{H}_A - \mathbf{H}_\omega \|^2}{\| \mathbf{H}_A + \mathbf{H}_\omega \|^2} \quad (22)$$

or by using the equivalent norm:

$$\hat{E}_H^2 = 2 \frac{\| \mathbf{H}_A - \mathbf{H}_\omega \|^2}{\max_{\theta} \| \mathbf{H}_A + \mathbf{H}_\omega \|^2}. \quad (23)$$

As described previously the torque can be computed as soon as the fields are available. However, these fields can be computed by using either the vector potential formulation or the scalar potential one. The resulting torques will be denoted by $C_A$ and $C_\omega$ respectively.

Obviously all these fields, $H_A$ and $H_\omega$ and the scalars $C_A$ and $C_\omega$, depend on the relative position between the rotor and stator. To define a torque relative error we must be careful because the torque vanishes for some rotor positions. To overcome this difficulty we could integrate the
<figure><img src="image_10.png" /><figcaption>Fig. 9. (Color online) Error in the magnetic field by using the FEM, the C-NEM and the mixed FEM-NEM discretizations calculated with a cloud of 400 000 nodes.</figcaption></figure>

different norms in the angular coordinate defining the rotor position:

$$ \tilde{E}_C^2 = 2 \frac{\int_0^{\frac{\pi}{4}} |C_A(\theta) - C_\omega(\theta)|^2 d\theta}{\int_0^{\frac{\pi}{4}} |C_A(\theta) + C_\omega(\theta)|^2 d\theta} \quad (24) $$

where the integration domain is reduced to $\theta \in [0, \frac{\pi}{4}]$ by symmetry considerations.

One could expect that the main differences between a finite element and a natural element solutions should be located in the air gap where due to the rotor movement the finite element mesh becomes very distorted. Of course in a precise implementation of the finite element method the mesh is updated locally or globally to avoid too distorted elements. However, in this work we have decided to proceed with the same cloud of nodes to perform both discretizations, the one using the FEM and the other the C-NEM. Thus, a linear finite element interpolation will be applied directly on the Delaunay's triangulation to obtain de finite element shape functions. The same triangulation was also considered to define the Voronoi diagram used to computed the C-NEM interpolation functions.

As the C-NEM becomes piecewise linear on the domain boundary or on any internal interface, as previously illustrated, one could define two interfaces, the first one located on the rotor external radius, and the other one on the internal stator radius. Thus, the region in between both radii defined the air gap, where a maximum mesh distortion is expected. It seems natural that an appealing choice for the interpolation in that region would consist of using the C-NEM approximation. On the other hand, because the CPU time related to the construction of the C-NEM shape functions is much more higher that the one

required to define the finite element approximation, one could prefer to use in the other two remaining regions a linear finite element discretization constructed on the Delaunay triangulation. In this way, both interpolations connect perfectly on the interfaces between the different regions. The mixed discretization strategy uses the C-NEM in the air gap and the FEM in the remaining domain.

Figure 9 depicts the evolution of $\tilde{E}_H(\theta)$ for the three functional approximations: FEM, C-NEM and the mixed FEM/C-NEM, where by symmetry considerations $\theta$ varies in the interval $\theta \in [45^\circ, 90^\circ]$. These solutions were computed using a cloud of 400 000 nodes.

Now, we compute the torque for different rotor positions $\theta \in [45^\circ, 90^\circ]$, different formulations (vector potential and scalar potential), different discretization techniques (C-NEM, FEM and mixed FEM-NEM) and different nodal densities.

Figure 10 compares the computed torques $C(\theta)$ using both potential formulations and the different discretization techniques (the FEM, the C-NEM and the mixed strategy), when a cloud of 6755 nodes was considered. Figures 11–13 compare the different formulations and discretization techniques when a cloud of $N = 27046$, $N = 106804$ and $N = 419116$ nodes are considered respectively.

We can notice that for coarse meshes the FEM generates the expected spurious oscillations whereas the C-NEM solution looks significantly less perturbed. All these solutions come close as the number of nodes increases, illustrating the convergence of all the used formulations. The mixed FEM/C-NEM solution shows less oscillations than the finite element solution, being its
<figure><img src="image_11.png" /><figcaption>Fig. 10. (Color online) Computed torque using the vector and scalar potentials, different discretization techniques and 6755 nodes.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. (Color online) Computed torque using the vector and scalar potentials, different discretization techniques and 27046 nodes.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12. (Color online) Computed torque using the vector and scalar potentials, different discretization techniques and 106 804 nodes.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. (Color online) Computed torque using the vector and scalar potentials, different discretization techniques and 419 116 nodes.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 14. (Color online) Convergence analysis for the magnetic field.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. (Color online) Convergence analysis for the torque.</figcaption></figure>

associated computing cost similar to the one related to the finite element solution.

Finally, we quantify the convergence rates for both the magnetic field and the torque. Concerning the magnetic field we consider a particular position of the rotor, in our simulation $\theta_1 = 0^\circ$, evaluating the error according to equation (22) for different nodal densities. The resulting error evolution is depicted using a log-log scale in Figure 14, where the norm (22) was used. The slope of this representation corresponds to the convergence rate, that for both the C-NEM and the linear FE interpolation here considered is expected to be one for the energy norm here considered. We can notice that despite the same conver-

gence rate that both discretizations exhibit, the C-NEM accuracy is significantly better. The convergence rate of the mixed formulation is also first order and its accuracy is in between the ones related to the FEM and the C-NEM.

Figure 15 depicts the results associated with the torque convergence analysis, where the error is computed according to equation (24). The conclusions are similar to the ones just addressed in the magnetic field convergence analysis. One can notice that the mixed solution is a bit less accurate than the C-NEM solution, but it is much better than the FEM one. However, the computing time related to the mixed solution is very close to the one associated with the FEM.
# 4 Conclusions

In this paper the constrained natural element method has been applied for simulating electromagnetic rotative machines. The use of the aforementioned technique allows to alleviate the difficulties related to remeshing in the finite element context, improving significantly the solution accuracy. A method combining the versatility of the FEM and the robustness of the meshless natural element method to take into account the movement has been proposed and tested. It seems that this combination is a good compromise between computing time and accuracy.

The extension to 3D geometries and the nodal density adaptation using appropriate error estimators are in progress.

# References

1. T. Belytschko, Y.Y. Lu, L. Gu, Int. J. Numer. Meth. Eng. 37, 229 (1994)

2. T. Belytschko, N. Moes, S. Usui, C. Parimi, Int. J. Numer. Meth. Eng. 50, 993 (2001)

3. A. Bruyère, L. Illoul, F. Chinesta, S. Clénet, IEEE Trans. Magn. 46, 1342 (2008)

4. E. Cueto, M. Doblaré, L. Gracia, Int. J. Numer. Meth. Eng. 49, 519 (2000)

5. B. Davat, Z. Ren, M. Lajoie-Mazenc, IEEE Trans. Magn. 21, 2296 (1985)

6. C.A.M. Duarte, J.T. Oden, Comput. Meth. Appl. Mech. Eng. 139, 237 (1996)

7. H. Hiyoshi, K. Sugihara, Computational Geometry 22, 167 (2002)

8. H.C. Lai, D. Rodger, P.J. Leonard, IEEE Trans. Magn. 28, 1732 (1992)

9. W.K. Liu, S. Jun, Y.F. Zhang, Int. J. Numer. Meth. Fluids 21, 1081 (1995)

10. L.B. Lucy, Astronomical J. 88, 1013 (1977)

11. J.M. Melenk, I. Babuška, Comput. Meth. Appl. Mech. Eng. 4, 289 (1996)

12. N. Moës, J. Dolbow, T. Belytschko, Int. J. Numer. Meth. Eng. 39, 131 (1999)

13. B. Nayroles, G. Touzot, P. Villon, Computat. Mech. 10, 307 (1992)

14. D. Organ, M. Fleming, T. Terry, T. Belytschko, Computat. Mech. 18, 1 (1996)

15. R. Perrin-Bit, J.L. Coulomb, IEEE Trans. Magn. 31, 1920 (1995)

16. T.W. Preston, A.B.J. Reece, P.S. Sangha, IEEE Trans. Magn. 24, 471 (1988)

17. F. Rapetti, A. Buffa, F. Bouillaut, Y. Maday, COMPEL 19, 332 (2000)

18. A. Razek, J.L. Coulomb, M. Feliachi, J.C. Sabonnadiere, IEEE Trans. Magn. 17, 3250 (1981)

19. M. Sambridge, J. Braun, M. McQueen, Geophys. J. Int. 122, 837 (1995)

20. J.R. Shewchuck, *Proceedings of the fourteenth annual symposium on computational geometry*, Minneapolis, Minnesota, June 1998, pp. 86–95, association for computing machinery

21. J.R. Shewchuck, *Proceedings of the sixteenth annual symposium on computational geometry*, Hong-Kong, June 2000, pp. 350–359, association for computing machinery

22. R. Seidel, in *1978–1988 Ten Years IIG* (1988), pp. 178–191

23. X. Shi, Y. Le Menach, J.P. Ducreux, F. Piriou, *Proceeding of COMPUMAG*, Shenyang, China, 2005

24. R. Sibson, Math. Proc. Camb. Phil. Soc. 87, 151 (1980)

25. N. Sukumar, B. Moran, T. Belytschko, Int. J. Numer. Meth. Eng. 43, 839 (1998)

26. J.L. Synge, *Proceedings of the Royal Society of London. Series A*, Math. Phys. Sci. 191, 447 (1947)

27. I. Tsukerman, IEEE Trans. Magn. 28, 2247 (1992)

28. G. Ventura, J.X. Xu, T. Belytschko, Int. J. Numer. Meth. Eng. 54, 923 (2002)

29. J. Yvonnet, D. Ryckelynck, P. Lorong, F. Chinesta, Int. J. Numer. Meth. Eng. 60, 1451 (2004)

30. J. Yvonnet, F. Chinesta, P. Lorong, D. Rynckelynck, Int. J. Thermal Sci. 44, 559 (2005)