© EDP Sciences, 2018

1 Introduction

One of the key physical quantities while designing an electrical motor [1,2] is the flux linkage of the inductors since it allows to compute the back electromotive force [3,4]. To do so, the finite elements method (FEM) is widely used at an industrial level (see pre-design softwares such as Flux Motor issued by Altair motorering). The FEM yields fairly good results but the computation time can sometimes be a bit long in a pre-design objective [5]. While the FEM needs the whole geometry to be meshed (air, active materials, sources, etc…), the volume integral method (VIM) only requires the meshing of active materials. For a few years, the VIM has been widely studied [6–8] and is know today as a good alternative to the FEM for solving magnetostatic problems in the presence of non-linear materials thanks to the robust acceleration method for its resolution. The magnetic field in the magnetic materials will thus be computed thanks to a 2D VIM presented in [9]. Consequently, the information regarding the magnetic field outside of the active materials is unavailable without an expensive post computation. It is thus necessary to use a flux computation method that requires only the information in the active materials. The method for the flux computation is inspired from [10] and adapted in 2D. This method will be tested with the calculation of the flux linkage of the first phase of a three-phased synchronous permanent magnet electrical motor. The results are compared to those obtained with the FEM.

2 Flux computation method

The flux computation method proposed is complementary to the VIM since it only requires the magnetic field inside the active materials. The vector potential can be split into A = A₀ + A_m where A₀ is the vector potential created by the coils in the vacuum and A_m is the vector potential in the magnetic regions. The magnetic flux linkage of a coil k associated to a domain Ω_0k can be expressed as [10]: $${\mathrm{\Phi }}_k={\int }_{{\mathrm{\Omega }}_{0k}}{\mathbf{\text{j}}}_{0k}\times \mathbf{\text{A}}\text{d}{\mathrm{\Omega }}_{0k}={\mathrm{\Phi }}_{0k}+{\mathrm{\Phi }}_{mk}$$(1) where j_0k is the current density vector that creates a current of 1 A in the coil k. Φ_0k is the contribution of all coils and Φ_mk is the contribution of the magnetized part. Let us calculate the contribution of the flux coming from A₀. One can write A₀ as: $${\text{A}}_0=\frac{{\mu }_0}{2\pi }{\displaystyle \sum _{l=0}^{N_{coil}}}\left({\text{i}}_l{\int }_{{\mathrm{\Omega }}_{0l}}{\text{j}}_{0l}G(\text{r})\text{d}{\mathrm{\Omega }}_{0l}\right)$$(2)

with G(r) the 2D Green kernel: $G(r)=\log (1/\parallel \text{r}\parallel )$, with r distance between the source point and the target point of integration. The contribution from the coils to the flux linkage of the coil k is: $$\begin{array}{ll}{\mathrm{\Phi }}_{0k}\hfill & =\frac{{\mu }_0}{2\pi }{\int }_{{\mathrm{\Omega }}_{0k}}{\text{j}}_{0k}\left\{{\displaystyle \sum _{l=0}^{N_{coil}}}\left(i_l{\int }_{{\mathrm{\Omega }}_{0l}}{\text{j}}_{0l}G(\text{r})\text{d}{\mathrm{\Omega }}_{0l}\right)\right\}\text{d}{\mathrm{\Omega }}_{0k}\hfill \\ \mathrm{}\hfill & ={\displaystyle \sum _{l=0}^{N_{coil}}}{L}_{kl}{i}_l\hfill \end{array}$$(3)

where i_l is the current through the coil l and: $$\begin{array}{ll}{L}_{kl}\hfill & =\frac{{\mu }_0}{2\pi }{\int }_{{\mathrm{\Omega }}_{0k}}{\text{j}}_{0k}\cdot {\int }_{{\mathrm{\Omega }}_{0l}}{\text{j}}_{0l}G(\text{r})\text{d}{\mathrm{\Omega }}_{0l}\text{d}{\mathrm{\Omega }}_{0k}\hfill \\ \mathrm{}\hfill & ={\int }_{{\mathrm{\Omega }}_{0k}}{\text{j}}_{0k}{\text{a}}_{0l}\text{d}{\mathrm{\Omega }}_{0k}\hfill \end{array}$$(4)

where a_0l represents the vector potential generated by the coil l in the vacuum with a 1 A current and L_kl the self (for k = l) and mutual inductance of the coils in the vacuum.

The contribution from the magnetic parts can be written as: $${\mathrm{\Phi }}_{mk}={\int }_{{\mathrm{\Omega }}_{0k}}{\text{j}}_{0k}\cdot {\text{A}}_m\text{d}{\mathrm{\Omega }}_{0k}$$(5) that can be modified by the divergence theorem: $${\mathrm{\Phi }}_{mk}={\int }_{{\mathrm{\Omega }}_{0k}}{\text{H}}_{0k}\cdot \nabla \times {\text{A}}_m\text{d}{\mathrm{\Omega }}_{0k}+{\int }_{{\mathrm{\Gamma }}_{0k}}\text{n}\cdot ({\text{A}}_m\times {\text{H}}_{0k})\text{d}{\mathrm{\Gamma }}_{0k}$$(6)where H_0k is the magnetic field created by the coil k when crossed by a 1 A current, Γ_0k is the boundary of the domain Ω_0k. The second part of equation (6) is equal to zero since the potential and magnetic field are conserved through Γ_0k [10]. We can thus obtain: $${\mathrm{\Phi }}_{mk}={\int }_{{\mathrm{\Omega }}_{0k}}{\text{H}}_{0k}\cdot \nabla \times {\text{A}}_m\text{d}{\mathrm{\Omega }}_{0k}$$(7)that can be transformed into [10]: $${\mathrm{\Phi }}_{mk}={\int }_{{\mathrm{\Omega }}_m}{\mu }_0{\text{H}}_{0k}\cdot M\text{d}{\mathrm{\Omega }}_m$$(8)where Ω_m represents the active materials domain and M is the magnetisation. We finally obtain a general formula for the flux linkage of a coil k: $${\mathrm{\Phi }}_k={\displaystyle \sum _{l=0}^{N_{bob}}}{\int }_{{\mathrm{\Omega }}_{0k}}{\text{j}}_{0k}\cdot {\text{a}}_{0l}\text{d}{\mathrm{\Omega }}_{0k}{i}_l+{\int }_{{\mathrm{\Omega }}_m}{\mu }_0{\text{H}}_{0k}\cdot M\text{d}{\mathrm{\Omega }}_m$$(9)

One can notice that the formula necessitates only the information of the magnetic field in the active materials (Ω_m). This makes this formula highly compatible with the VIM. Nevertheless, the vector potential and the magnetic field created by the coils are required in the formula (9). The second part of this paper introduces an analytical method to obtain those values.

3 Magnetic field and vector potential induced by a 2D meshed coil

Since the shape of the coils within an electrical rotating machine are not basic shapes (e.g. squares, circles), computing the magnetic field and vector potential generated by those coils is not easy. We have adopted a semi-analytical approach that consists on meshing the inductors with triangles and computing analytical magnetic field and vector potential per triangle Figure 1. Contrary to the 3D situation where Urankar's formulas are well known [11–14], There are no analytical formulas for the field and potential created by a triangle in 2D. Let T be a triangle of the mesh, we can shift and rotate this triangle so that it is in a reference frame (longest edge along the O_x axis) and cut it in two rectangle triangles (see Fig. 2). We can thus compute analytically the field and vector potential created by both triangles T₁ and T₂ (see Figs. 2 and 3). We could have chosen to work with a normalised reference triangle and a geometrical Jacobian transformation [15] but our approach is more straightforward and easier to implement, that is why we decided to present it here. Moreover, the use of a jacobian transformation is less effective regarding the computation time since more mathematical operations are required for the computation of the field and potential created by each triangle. We use the complex potential created by each triangle defined as [16]: $$P(z)=-\frac{J}{2\pi }{\int }_{\mathrm{\Omega }}\log (z-z_i)\text{d}\mathrm{\Omega }$$(10) where Ω is the triangle surface crossed by a current density J and z the complex coordinate of the target point. From this potential can be extracted the magnetic field and vector potential: $$\frac{\text{d}P(z)}{\text{d}z}=-i{H}^{*}=-H_y-i{H}_x$$(11) $$A=\Re e(P(z))$$(12)where i is the imaginary number.

Let us compute the complex potential created by T₁ and T₂ described in Figure 3. The formula for the complex potential yields:

$P{(z)}_{T1}=-\frac{J}{2\pi }{\int }_0^a{\int }_0^{b\xi /a}\log (z-(\xi +i\eta ))\text{d}\eta \text{d}\xi$(13)

$$P{(z)}_{T2}=-\frac{J}{2\pi }{\int }_0^a{\int }_0^{b-b\xi /a}\log (z-(\xi +i\eta ))\text{d}\eta \text{d}\xi$$(14)

Let us focus on the potential created by T₁. We use the fact that a primitive of log(u) is u[log(u) − 1]: $$\begin{array}{ll}P{(z)}_{T1}\hfill & =\frac{J}{2\pi }{\int }_0^a{[(z-\xi -i\eta )(\log (z-\xi -i\eta )-1)]}_0^{b\xi /a}\text{d}\xi \hfill \\ \hfill & =\frac{J}{2\pi }{\int }_0^a\left\{-\xi \frac{ib}{a}-[z-\xi ]\log (z-\xi )+\left[z-\xi \left(1+\frac{ib}{a}\right)\right]\log \left(z-\xi \left(1+\frac{ib}{a}\right)\right)\right\}\text{d}\xi \hfill \end{array}$$(15) which leads to, with z₀ = 0, z₁ = a, z₂ = a + ib: $$P{(z)}_{T1}=\frac{J}{2\pi }\left[i{C}^{te}+\frac{{(z-z_1)}^2}{2}\left(\log (z-z_1)-\frac{1}{2}\right)-\frac{{(z-z_0)}^2}{2}\left(\log (z-z_0)-\frac{1}{2}\right)-\frac{1}{1+\frac{ib}{a}}\left(\frac{{(z-z_2)}^2}{2}\left(\log (z-z_2)-\frac{1}{2}\right)-\frac{{(z-z_0)}^2}{2}\left(\log (z-z_0)-\frac{1}{2}\right)\right)\right]$$(16)

We then use the equations given in equation (11) to get the vector potential: $$\frac{2\pi A}{J}=D_0\frac{{\alpha }_0}{2}+P_0\left(\frac{1}{2}-\log (r_0)\right)-D_1\frac{{\alpha }_1}{2}+P_1\left(\log (r_1)-\frac{1}{2}\right)+\frac{1}{1+\frac{b^2}{a^2}}\left(D_2\frac{{\alpha }_2}{2}+P_2\left(\frac{1}{2}-\log (r_2)\right)-D_0\frac{{\alpha }_0}{2}+P_0\left(\log (r_0)-\frac{1}{2}\right)\right)+\frac{\frac{b}{a}}{1+\frac{b^2}{a^2}}\left(D_2\left(\frac{1}{4}-\frac{\log (r_2)}{2}\right)-P_2{\alpha }_2+D_0\left(\frac{\log (r_0)}{2}-\frac{1}{4}\right)+P_0{\alpha }_0\right)$$(17) where D_i = (y − y_i) ² − (x − x_i) ², and P_i = (x − x_i) (y − y_i). For the magnetic field: $$\frac{-2\pi H_x}{J}=x\log \left(\frac{r_0}{r_1}\right)+y({\alpha }_1-{\alpha }_0)+a\log (r_1)+\frac{x\log \left(\frac{r_2}{r_0}\right)+y({\alpha }_0-{\alpha }_2)-a\log (r_2)+b{\alpha }_2}{1+\frac{b^2}{a^2}}+\frac{b}{a}\left(\frac{y\log \left(\frac{r_2}{r_0}\right)+x({\alpha }_2-{\alpha }_0)-b\log (r_2)-a{\alpha }_2}{1+\frac{b^2}{a^2}}\right)$$(18) $$\frac{2\pi H_y}{J}=y\log \left(\frac{r_0}{r_1}\right)+x({\alpha }_0-{\alpha }_1)+a{\alpha }_1+\frac{y\log \left(\frac{r_2}{r_0}\right)+x({\alpha }_2-{\alpha }_0)-a{\alpha }_2}{1+\frac{b^2}{a^2}}+\frac{b}{a}\left(\frac{x\log \left(\frac{r_0}{r_2}\right)+y({\alpha }_2-{\alpha }_0)-b{\alpha }_2}{1+\frac{b^2}{a^2}}\right)$$(19)

The analytical formulas for the vector potential and the field created by T₂ are obtained the same way: $$\frac{2\pi A}{J}=D_0\frac{{\alpha }_0}{2}+P_0\left(\frac{1}{2}-\log (r_0)\right)-D_1\frac{{\alpha }_1}{2}+P_1\left(\log (r_1)-\frac{1}{2}\right)+\frac{1}{1+\frac{b^2}{a^2}}\left(D_1\frac{{\alpha }_1}{2}+P_1\left(\frac{1}{2}-\log (r_2)\right)-D_2\frac{{\alpha }_2}{2}+P_2\left(\log (r_2)-\frac{1}{2}\right)\right)-\frac{\frac{b}{a}}{1+\frac{b^2}{a^2}}\left(D_1\left(\frac{1}{4}-\frac{\log (r_1)}{2}\right)-P_1{\alpha }_1+D_2\left(\frac{\log (r_2)}{2}-\frac{1}{4}\right)+P_2{\alpha }_2\right)$$(20) $$\frac{-2\pi H_x}{J}=x\log \left(\frac{r_0}{r_1}\right)+y({\alpha }_1-{\alpha }_0)+a\log (r1)+\frac{x\log \left(\frac{r_1}{r_2}\right)+y({\alpha }_2-{\alpha }_1)-a\log (r1)-b{\alpha }_2}{1+\frac{b^2}{a^2}}+\frac{b}{a}\left(\frac{y\log \left(\frac{r_2}{r_1}\right)+x({\alpha }_2-{\alpha }_1)-b\log (r2)+a{\alpha }_1}{1+\frac{b^2}{a^2}}\right)$$(21) $$\frac{2\pi H_y}{J}=y\log \left(\frac{r_0}{r_1}\right)+x({\alpha }_0-{\alpha }_1)+a{\alpha }_1+\frac{y\log \left(\frac{r_1}{r_2}\right)+x({\alpha }_1-{\alpha }_2)-a{\alpha }_1+b\log (r_2)}{1+\frac{b^2}{a^2}}+\frac{b}{a}\left(\frac{x\log \left(\frac{r_1}{r_2}\right)+y({\alpha }_2-{\alpha }_1)-b{\alpha }_2-a\log (r_1)}{1+\frac{b^2}{a^2}}\right)$$(22)

Thanks to those formulas, it is possible to compute the magnetic field and vector potential created by any unspecified triangle and thus by any 2D triangularly meshed coil. Consequently, those formulas allow us to compute the flux with the semi-analytical equation (9).

Fig. 1 Example of the meshed the coil of an electrical motor.

Fig. 2 Shifting and rotation of the triangle T.

Fig. 3 Parametrization of the two rectangle triangles.

4 Validation of the analytical formulas

We validate our analytical formulas with the test case presented Figure 4. It consists of three triangularly shaped coils with a current of 1000 A and we computed the field and vector potential along the curve C. We compared the values obtained with our analytical approach to an FEM solution: We plotted the mean;1; difference between the analytical solution and the solution yielded by the FEM for several mesh qualities: $$\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}{H}_x=\frac{1}{N}{\displaystyle \sum _{i=1}^N}\left|\frac{{\text{H}}_{\text{anal}}^i\cdot {\text{e}}_x-{\text{H}}_{\text{FEM}}^i\cdot {\text{e}}_x}{{\displaystyle \underset{i=\mathrm{1,..},N}{{\displaystyle \max }}}({\text{H}}_{\text{FEM}}^i\cdot {\text{e}}_x)}\right|$$ $$\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}{H}_y=\frac{1}{N}{\displaystyle \sum _{i=1}^N}\left|\frac{{\text{H}}_{\text{anal}}^i\cdot {\text{e}}_y-{\text{H}}_{\text{FEM}}^i\cdot {\text{e}}_y}{{\displaystyle \underset{i=\mathrm{1,..},N}{{\displaystyle \max }}}({\text{H}}_{\text{FEM}}^i\cdot {\text{e}}_y)}\right|$$ $$\mathrm{M}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{E}\mathrm{r}\mathrm{r}\mathrm{o}\mathrm{r}A=\frac{1}{N}{\displaystyle \sum _{i=1}^N}\left|\frac{A_{\text{anal}}^i-A_{\text{FEM}}^i}{{\displaystyle \underset{i=\mathrm{1,..},N}{{\displaystyle \max }}}(A_{\text{FEM}}^i)}\right|$$(23) where N is the number of points along the curve C (we chosed N = 360). The results are displayed Figure 5. The graph shows that the better the FEM problem is meshed, the closer its solution is to the analytical solution, which validates our analytical formulas. The reader can notice that the error on the potential vector is an order of magnitude less than on the magnetic field. This can be explained by the relative “smootheness” of the potential vector due to the fact that it is a primitive to the magnetic field.

Fig. 4 Description of test case.

Fig. 5 Results of the test case.

5 Flux calculation in the vacuum

Since the integration of Φ_0k is semi-analytical (a_0l is computed analytically with the formula described above), we decide to use an adaptative number of Gauss points to ensure the convergence of the calculation, which mean we stop increasing the number of Gauss points when the variation of the flux is less than a certain value (in this paper this value is 0.001%). This allows us to use the coarsest mesh for the coils and still get an accurate value for mutual and self inductances L_kl in equation (4).

We use the algorithm shown Figure 6 to do the calculation, where $L_{kl}^i$ is the value of the mutual inductance between coil k and coil l computed with a number of Gauss points i on the Ω_0k domain. We used a value of 10⁻⁵ for ε. One has to keep in mind that the further the two coils are from each other, the less Gauss points the calculation needs.

To validate our flux in vacuum computation method, we compute the flux in the coil B of the test case Figure 7 with our adaptative semi-analytical method. To validate our approach we use an FEM reference solution, obtained with a sufficiently fine mesh (around 300 000 second order elements). The FEM yielded a value of the flux of 1.792688 × 10⁻⁵ Wb while our adaptative computation yielded 1.792684 × 10⁻⁵ Wb (less than 3 × 10⁻⁵% error). We plot Figure 8 the relative error and the computation time of our adaptative approach for several numbers of Gauss points. We can see that we used 25 Gauss points, computing the result in less than 0.05 s. To obtain this level of precision with the FEM, at least 10 s and a very refined mesh (200 000 elements and 3 gauss points per elements, almost as refined as our “reference solution”) are required while we used only 18 triangles Figure 8.

Fig. 6 Algorithm for the computation of the inductances.

Fig. 7 Test case for the flux computation in the vacuum.

Fig. 8 Error and computation time for several Gauss points number used for integration.

6 Flux computation through the phase of an electrical motor

In order to validate our flux computation method, we used a 24 slot − 8 poles permanent magnet electrical motor (see Fig. 9). We calculated the flux in the phase A of the motor with no current in the coils for several positions of the rotor (from 0 to π/4). Three qualities of mesh have been used for the resolution of the magnetostatic problem via the VIM: 70, 205 and 371 unknowns for one anti-periodicity of the motor. A comparison between the 70 and 371 unknowns mesh can be found Figure 10.

We compared the solution obtained by our flux computation method with one solution obtained via the FEM with a very high quality mesh (around 400 000 elements). The plot of those four curves (one for the FEM reference and three for our method at different qualities of mesh) can be found Figures 11 and 12 while the relative differences between the reference and our method can be found Figure 13. The plot of the computation times required for the FEM and the VIM and the analytical method approach can be found Figure 14.

One can notice that even with a very low number of unknowns, our semi-analytical method yields very accurate results for a very reduced amount of calculation: indeed, while the mean difference between the FEM reference and the solution obtained with 70 unknowns is 2.28% with a maximum difference of 3.81%, the computation time is less than 0.3 s per position of the rotor. This low computation time makes this method very effective in a pre-design context (which was the aim of this study). In opposite, the FEM takes at least 1.5 s per position for the same level of precision (we used a 5000 elements mesh) compared to the reference.

Fig. 9 Test case for the flux computation.

Fig. 10 A total of 70 and 371 unknowns mesh for the used for the flux computation.

Fig. 11 Flux linkage of the phase A of the test motor for several positions of the rotor.

Fig. 12 Flux linkage of the phase A of the test motor for several positions of the rotor (zoomed).

Fig. 13 Difference between the reference and the semi-analytic solution for several positions of the rotor.

Fig. 14 Error and computation time of FEM and VIM + the proposed semi-analytical approach (per position for the computation time).

7 Conclusions

Original analytical formulas for the computation of the magnetic field and the potential created by a 2D unspecified triangle have been introduced. Those formulas have been rigorously validated with a Finite Elements model. They allow us to calculate semi-analytically the flux linkage of a 2D meshed coil of any shape. This semi-analytical method of flux computation is very well suited for the post processing of a solution yielded by a VIM. We thus used the coupling of the VIM and the proposed method to compute the flux linkage of the phase of a permanent magnet synchronous electrical motor. This coupling allows to gather good quality results rather quickly compared to the FEM. We can see that in addition to being as precise as the FEM for normally meshed problems [7,8,17] is very efficient in a pre-design context, which was the aim of this study.

Author contribution statement

The vast majority of the work has been carried out by Quentin Debray who was then a PhD student under the direction of Dr. Gerard Meunier (Univ. Grenoble Alpes, CNRS, Grenoble INP) and supervised by Dr. Olivier Chadebec (Univ. Grenoble Alpes, CNRS, Grenoble INP), Dr. Jean Louis Coulomb (Univ. Grenoble Alpes, CNRS, Grenoble INP) and Dr. Anthony Carpentier (Altair Engineering). Nevertheless, Drs. Meunier, Chadebec, Coulomb and Carpentier played an essential role of counseling during both the research and redaction phase.

References

All Figures

	Fig. 1 Example of the meshed the coil of an electrical motor.
In the text

	Fig. 2 Shifting and rotation of the triangle T.
In the text

	Fig. 3 Parametrization of the two rectangle triangles.
In the text

	Fig. 4 Description of test case.
In the text

	Fig. 5 Results of the test case.
In the text

	Fig. 6 Algorithm for the computation of the inductances.
In the text

	Fig. 7 Test case for the flux computation in the vacuum.
In the text

	Fig. 8 Error and computation time for several Gauss points number used for integration.
In the text

	Fig. 9 Test case for the flux computation.
In the text

	Fig. 10 A total of 70 and 371 unknowns mesh for the used for the flux computation.
In the text

	Fig. 11 Flux linkage of the phase A of the test motor for several positions of the rotor.
In the text

	Fig. 12 Flux linkage of the phase A of the test motor for several positions of the rotor (zoomed).
In the text

	Fig. 13 Difference between the reference and the semi-analytic solution for several positions of the rotor.
In the text

	Fig. 14 Error and computation time of FEM and VIM + the proposed semi-analytical approach (per position for the computation time).
In the text