Mathematical formulation of a new linear surface element for the modeling of axisymmetric magnetic field problems

Abstract

There are currently two main trends in solving magnetic field problems, numerical methods, of which the Finite Element Method (FEM) is undoubtedly the most popular, and analytical methods, of which there is a great diversity of approaches. The main advantage of numerical methods is their flexibility in modeling complicated geometries; however, its calculation times are generally considerable. Analytical methods attain very short calculation times by sacrificing flexibility. This paper presents the mathematical formulation of a new linear surface element for the semi-analytical solution of axisymmetric magnetic field problems. The importance of this formulation is that it is capable of dealing with arbitrary geometries while retaining the flexibility of numerical methods with considerably shorter calculation times. The results of the application of the proposed formulation to a case study are presented in the validation section. The results obtained are compared with those calculated with FEM showing an excellent agreement between them.

Keywords

Magnetic field numerical methods analytical methods surface current

1. Introduction

The precision, speed and flexibility in the calculation of magnetic fields are undoubtedly issues of great relevance in the field of computational electromagnetics [1]. In the last decades an unprecedented improvement has been seen in the calculation capacity of computers; however, the complexity of the problems to be addressed grows with the available computational resources. Therefore, the solution of magnetic field problems involves dealing with a trade-off between these three variables [2].

Despite the great advances in the field of numerical methods in recent years, FEM being the most representative method, it is indisputable that a numerical method capable of solving problems of arbitrarily complicated geometries, with high precision and very short calculation times (on the order of a few seconds) is not yet available. This situation is even worse in the case of 3D transient problems with non-linearities or anisotropies [3–5].

As a consequence, great efforts continue to be made to formulate analytical methods for solving specific magnetic field problems, so that the proposed model is suitable to represent the case of interest with reasonable precision and in relatively short times.

The disadvantage of analytical methods is that their applicability is restricted only to the specific problem for which they were intended. If there is a significant variation in the geometry or the topology of the problem, the calculation methodology must be reformulated. Applications of analytical methods for the solution of the magnetic field in electrical machines are presented in [6–8]. It is important to note that in these works the analytical method is only applicable to the type of device and under the conditions for which it has been formulated.

In a previous paper, the present authors proposed a methodology which can be considered a trade-off between numerical and analytical methods called Semianalytic Integral Method (SAIM) [9]. The main characteristic of this method is that it discretizes the geometry of the problem using elements whose analytical solution is known, and a linear system of equations is proposed to determine the values of the excitations of each element accounting for the boundary conditions [10]. The magnetic axisymmetric problem has been solved using cylindrical [11] and disk-shaped elements [12].

Despite the good performance of SAIM in terms of precision and calculation time in a particular case study, the authors identified that the cylindrical and disk-shaped elements significantly restrict the applicability of SAIM to more general problems, especially those with curved surfaces or surfaces which are not parallel to the axes. This limitation motivated the authors to present the mathematical formulation of a new element with conical geometry in this paper, which allows the modeling of geometries much more complicated than those which can be represented by cylindrical and disk-shaped elements.

The most important contribution of this work is that the mathematical formulation of the new conical element makes it possible to model and solve axisymmetric magnetic field problems with arbitrarily complicated geometries.

A second notable contribution of this work is that a method of integration based on matrices has been proposed to evaluate the interactions between an arbitrary number of elements and evaluation points in a computationally efficient and structured way.

In the final part of the work the proposed conical element is applied along with SAIM to a case study which requires the use of conical elements. In this case study, the magnetic field distribution in the vicinity of a ferrite core under the effect of an external magnetic field generated by a Helmholtz coil is determined. An excellent agreement has been found between the results calculated using the conical element together with SAIM and those calculated by FEM.

1.1. Problem statement

The authors’ research has been focused on the determination of a computationally efficient method to determine the vector potential A _𝜙 and the magnetic field strength in the radial and axial directions H _r and H _z produced by a conical surface current distribution at an arbitrary point in space.

The geometry of the current distribution is shown in Fig. 1. As can be seen, the distribution starts from a ring of radius r ₁ on the axial coordinate z ₁. This ring has a current density K _𝜙1. The current density varies linearly along the surface of the cone until it reaches the value of K _𝜙2 at the coordinates z ₂ and r ₂. The coordinate of the evaluation point P (r, z) is also shown. The coordinates that define the cone are totally arbitrary as long as $\{r,r_{1},r_{2},z,z_{1},z_{2}\}\in \mathbb{R}:r\geq 0,r_{1}>0,r_{2}>0,r_{1}\neq r_{2}$ y $\{K_{{\phi}1},K_{{\phi}2}\}\in \mathbb{C}$ . As it will be seen later, the formulations have a singularity in r ₁ = r ₂ which is inherited from the definition of the slope in the equation of the line. In spite of the above, the validity condition r ₁≠ r ₂ does not represent a disadvantage at all because this particular case was already discussed in a previous work [11]. It is important to mention that the international system of units has been used in this paper.

1.2. Basic mathematical expressions

This section presents the mathematical expressions to be solved to determine the magnetic field produced by the current distribution. The expression for the magnetic vector potential is presented first, and then the expressions for the radial and axial fields are derived from it.

Fig. 1.

Schematic diagram of the finite cone.

1.2.1. Magnetic vector potential

A filamentary circular loop whose axis is the z-axis and located at a distance z ^′ from the x-y plane is examined as first step. The radius of the loop is r ^′ and a current I _𝜙 flows through it. According to the Biot-Savart law the magnetic vector potential produced by the loop at an arbitrary field point P (r, z) can be written as follows [13] $\begin{eqnarray}\displaystyle A_{{\phi}}^{[fc]}=\frac{{\mu}_{0}}{2{\pi}}\int _{0}^{{\pi}}\frac{I_{{\phi}}\cos ({\phi}^{\prime })r^{\prime }\text{d}{\phi}^{\prime }}{\sqrt{r^{2}+r^{\prime 2}+(z-z^{\prime })^{2}-2rr^{\prime }\cos ({\phi}^{\prime })}} & & \displaystyle\end{eqnarray}$ (1) where the primed coordinates are associated with the field source and the non-primed ones are related to the observation point also called field point. The superscript [fc] denotes filamentary current. On the other hand, in the case of a conical current distribution, the current density K _𝜙 is assumed to vary linearly from the coordinate (r ₁, z ₁) to the coordinate (r ₂, z ₂). This linear variation can be written as follows $\begin{eqnarray}\displaystyle K_{{\phi}}=\frac{(K_{{\phi}2}-K_{{\phi}1})(z^{\prime }-z_{1})}{z_{2}-z_{1}}+K_{{\phi}1}. & & \displaystyle\end{eqnarray}$ (2) In the r − z plane, the equation describing the surface of the cone is given by the equation of a straight line as follows $\begin{eqnarray}\displaystyle z^{\prime }=m(r^{\prime }-r_{1})+z_{1} & & \displaystyle\end{eqnarray}$ (3) where $m=\frac{z_{2}-z_{1}}{r_{2}-r_{1}}$ is the slope of the line. From the previous expression of the slope it is evident that the mathematical expressions presented in this work will only be valid when r ₁≠ r ₂, since the equation of the line is not defined when the two radii are equal, i.e., when the cone represents the geometry of a regular cylinder. Fortunately, this particular case has been investigated in detail in [11]. According to [14], the differential length can be written as (see Fig. 1) $\begin{eqnarray}\displaystyle \text{d}l^{\prime }=\left[1+\left(\frac{\text{d}z^{\prime }}{\text{d}r^{\prime }}\right)^{2}\right]^{\frac{1}{2}}\text{d}r^{\prime }=\sqrt{1+m^{2}}\text{d}r^{\prime }={\xi}_{0}\text{d}r^{\prime } & & \displaystyle\end{eqnarray}$ (4) where ${\xi}_{0}=\sqrt{1+m^{2}}$ . The differential current contribution of a circular loop of length dl ^′ is I _𝜙 = K _𝜙 dl ^′, so that the differential current can be written as follows $\begin{eqnarray}\displaystyle I_{{\phi}}={\xi}_{0}\left[\frac{m(K_{{\phi}2}-K_{{\phi}1})(r^{\prime }-r_{1})}{z_{2}-z_{1}}+K_{{\phi}1}\right]\text{d}r^{\prime }. & & \displaystyle\end{eqnarray}$ (5)

To obtain the total magnetic vector potential due to the conical current distribution, the differential current contributions I _𝜙 from r ₁ to r ₂ must be superimposed, which is achieved by integrating (1) between r ₁ and r ₂. $\begin{eqnarray}\displaystyle A_{{\phi}}=\frac{{\mu}_{0}}{2{\pi}}\int _{r_{1}}^{r_{2}}\int _{0}^{{\pi}}\frac{I_{{\phi}}\cos ({\phi}^{\prime })r^{\prime }\text{d}{\phi}^{\prime }}{\sqrt{r^{2}+r^{\prime 2}+(z-z^{\prime })^{2}-2rr^{\prime }\cos ({\phi}^{\prime })}}. & & \displaystyle\end{eqnarray}$ (6)

It should be noted that in (6) only the differential d𝜙^′ appears. This is because the differential dr ^′ is within I _𝜙 as defined in (5). For the sake of simplicity, the integrand of (6) is defined using the variable f _A(𝜙^′, r ^′) as follows $\begin{eqnarray}\displaystyle f_{A}({\phi}^{\prime },r^{\prime })=\frac{{\xi}_{0}\left[\frac{m(K_{{\phi}2}-K_{{\phi}1})(r^{\prime }-r_{1})}{z_{2}-z_{1}}+K_{{\phi}1}\right]\cos ({\phi}^{\prime })r^{\prime }}{\sqrt{r^{2}+r^{\prime 2}+[z-m(r^{\prime }-r_{1})-z_{1}]^{2}-2rr^{\prime }\cos ({\phi}^{\prime })}} & & \displaystyle\end{eqnarray}$ (7) and consequently, $\begin{eqnarray}\displaystyle A_{{\phi}}=\frac{{\mu}_{0}}{2{\pi}}\int _{0}^{{\pi}}\int _{r_{1}}^{r_{2}}f_{A}({\phi}^{\prime },r^{\prime })∼\text{d}r^{\prime }\text{d}{\phi}^{\prime }. & & \displaystyle\end{eqnarray}$ (8)

Note that in (8) the order of integration has been changed because the integral is first evaluated analytically with respect to r ^′. The method to perform the integration with respect to 𝜙^′ will be presented later.

1.2.2. Magnetic field intensity in radial direction

The magnetic flux density $\vec{\boldsymbol{B}}$ is related to the magnetic vector potential by the following relationship $\begin{eqnarray}\displaystyle \vec{\boldsymbol{B}}=\hat{\boldsymbol{r}}B_{r}+\hat{\boldsymbol{{\phi}}}B_{{\phi}}+\hat{\boldsymbol{z}}B_{z}={\nabla}\times \vec{\boldsymbol{A}}. & & \displaystyle\end{eqnarray}$ (9)

Since the problem is axisymmetric, the vector potential only has an azimuthal component 𝜙 ( $\vec{\boldsymbol{A}}=\hat{\boldsymbol{{\phi}}}A_{{\phi}}$ ), then the magnetic flux density in radial direction B _r is obtained from (9) as follows $\begin{eqnarray}\displaystyle B_{r}=-\frac{\partial A_{{\phi}}}{\partial z}. & & \displaystyle\end{eqnarray}$ (10)

By assuming that the current distribution is placed in a homogeneous medium having the permeability of vacuum, the following constitutive relation is valid $\begin{eqnarray}\displaystyle \vec{\boldsymbol{B}}={\mu}_{0}\vec{\boldsymbol{H}}={\mu}_{0}(\hat{\boldsymbol{r}}H_{r}+\hat{\boldsymbol{{\phi}}}H_{{\phi}}+\hat{\boldsymbol{z}}H_{z}). & & \displaystyle\end{eqnarray}$ (11) By substituting (10) in (11) and solving for H _r, the following expression is found $\begin{eqnarray}\displaystyle H_{r}=-\frac{1}{{\mu}_{0}}\frac{\partial A_{{\phi}}}{\partial z}. & & \displaystyle\end{eqnarray}$ (12) The following expression for the magnetic field strength in the radial direction can be written by substituting (8) in (12) and rearranging, $\begin{eqnarray}\displaystyle H_{r}=-\frac{1}{2{\pi}}\int _{0}^{{\pi}}\int _{r_{1}}^{r_{2}}\frac{\partial }{\partial z}f_{A}({\phi}^{\prime },r^{\prime })∼\text{d}r^{\prime }\text{d}{\phi}^{\prime }. & & \displaystyle\end{eqnarray}$ (13)

1.2.3. Magnetic field intensity in axial direction

In this case, (9) gives the magnetic flux density in axial direction B _z $\begin{eqnarray}\displaystyle B_{z}=\frac{1}{r}\frac{\partial (rA_{{\phi}})}{\partial r}. & & \displaystyle\end{eqnarray}$ (14) Hence $\begin{eqnarray}\displaystyle H_{z}=\frac{1}{{\mu}_{0}r}\frac{\partial (rA_{{\phi}})}{\partial r}. & & \displaystyle\end{eqnarray}$ (15) By replacing (8) in (15), the expression of the magnetic field strength in the axial direction is obtained $\begin{eqnarray}\displaystyle H_{z}=\frac{1}{2{\pi}r}\int _{0}^{{\pi}}\int _{r_{1}}^{r_{2}}\frac{\partial }{\partial r}(rf_{A}({\phi}^{\prime },r^{\prime }))∼\text{d}r^{\prime }\text{d}{\phi}^{\prime }. & & \displaystyle\end{eqnarray}$ (16) It should be noted that B _𝜙 = 0, and therefore H _𝜙 = 0 if (9) is applied to the case ${\nabla}\times \hat{\boldsymbol{{\phi}}}A_{{\phi}}$ .

2. Solution methodology

2.1. Analytical integration with respect to r ^′

2.1.1. Magnetic vector potential

It can be shown that the magnetic vector potential can be expressed as the sum of the current densities at the ends of the cone weighed by two geometric constants. These two constants k _𝜙1 and k _𝜙2 are referred to as geometric factors, which only depend on the dimensions of the cone and the coordinates of the field point. According to the above mentioned, the magnetic vector potential can be written as follows $\begin{eqnarray}\displaystyle A_{{\phi}}=k_{{\phi}1}K_{{\phi}1}+k_{{\phi}2}K_{{\phi}2}. & & \displaystyle\end{eqnarray}$ (17)

It is important to note that the form of (17) is very convenient, since in many problems the values of the current densities at the ends of the elements are unknown, so the form of (17) together with the corresponding boundary conditions makes it possible to set up a system of equations to determine the current distribution of the conical elements. The geometric factors can be calculated using the following expressions $\begin{eqnarray}\displaystyle & \displaystyle k_{{\phi}1}=\frac{{\mu}_{0}}{2{\pi}}\int _{0}^{{\pi}}\bigg(s_{a}\bigg|_{r_{i}=r_{2}}\bigg|_{r^{\prime }=r_{1}}^{r^{\prime }=r_{2}}\bigg)∼\text{d}{\phi}^{\prime } & \displaystyle\end{eqnarray}$ (18) $\begin{eqnarray}\displaystyle & \displaystyle k_{{\phi}2}=-\frac{{\mu}_{0}}{2{\pi}}\int _{0}^{{\pi}}\bigg(s_{a}\bigg|_{r_{i}=r_{1}}\bigg|_{r^{\prime }=r_{1}}^{r^{\prime }=r_{2}}\bigg)∼\text{d}{\phi}^{\prime } & \displaystyle\end{eqnarray}$ (19) where r _i is an auxiliary variable which allows using the same expression of s _a to determine both factors k _𝜙1 and k _𝜙2. The only difference between k _𝜙1 and k _𝜙2 is the value of r _i and the signs. In the expressions of the radial and axial field the variable r _i is used in the same way.

The analytical integration with respect to r ^′ and a subsequent simplification yields $\begin{eqnarray}\displaystyle s_{a}={\xi}_{24}\cdot \left(\begin{array}{@{}c@{}}\displaystyle \frac{{\xi}_{3}\log ({\xi}_{13})({\xi}_{12}-4{\xi}_{11}{\xi}_{6}r_{i})}{{{\xi}_{6}}^{5/2}}\\[10.0pt] \displaystyle -\frac{2\sqrt{{\xi}_{10}}{\xi}_{3}(5{\xi}_{11}+{\xi}_{9}-2{\xi}_{6}r_{i})}{{{\xi}_{6}}^{2}}\end{array}\right). & & \displaystyle\end{eqnarray}$ (20)

2.1.2. Magnetic field intensity in radial direction

The procedure followed is the same as in the case of the magnetic vector potential. $\begin{eqnarray}\displaystyle H_{r}=k_{r1}K_{{\phi}1}+k_{r2}K_{{\phi}2}. & & \displaystyle\end{eqnarray}$ (21) Note that the geometric factors for the radial field are now denoted as k _r1 and k _r2, which affect the current densities associated with the cone ends. These geometric factors can be evaluated by the following expressions $\begin{eqnarray}\displaystyle & \displaystyle k_{r1}=-\frac{1}{2{\pi}}\int _{0}^{{\pi}}\bigg(s_{r}\bigg|_{r_{i}=r_{2}}\bigg|_{r^{\prime }=r_{1}}^{r^{\prime }=r_{2}}\bigg)∼\text{d}{\phi}^{\prime } & \displaystyle\end{eqnarray}$ (22) $\begin{eqnarray}\displaystyle & \displaystyle k_{r2}=\frac{1}{2{\pi}}\int _{0}^{{\pi}}\bigg(s_{r}\bigg|_{r_{i}=r_{1}}\bigg|_{r^{\prime }=r_{1}}^{r^{\prime }=r_{2}}\bigg)∼\text{d}{\phi}^{\prime } & \displaystyle\end{eqnarray}$ (23) where $\begin{eqnarray}\displaystyle s_{r}={\xi}_{24}\cdot \left(\begin{array}{@{}l@{}}\displaystyle \frac{{\xi}_{17}}{{{\xi}_{6}}^{2}}-\frac{{\xi}_{15}{\xi}_{3}(5{\xi}_{11}+{\xi}_{9}-2{\xi}_{6}r_{i})}{\sqrt{{\xi}_{10}}{{\xi}_{6}}^{2}}\\[15.0pt] \displaystyle \quad +\frac{{\xi}_{16}({\xi}_{12}-4{\xi}_{11}{\xi}_{6}r_{i})}{{\xi}_{13}{{\xi}_{6}}^{5/2}}\\[10.0pt] \displaystyle \quad +\frac{\log ({\xi}_{13})({\xi}_{14}-4{{\xi}_{3}}^{2}{\xi}_{4}{\xi}_{6}r_{i})}{{{\xi}_{6}}^{5/2}}\end{array}\right). & & \displaystyle\end{eqnarray}$ (24)

2.1.3. Magnetic field intensity in axial direction

Similarly, the following set of expressions for the magnetic field strength in the axial direction can be written. $\begin{eqnarray}\displaystyle & \displaystyle H_{z}=k_{z1}K_{{\phi}1}+k_{z2}K_{{\phi}2} & \displaystyle\end{eqnarray}$ (25) $\begin{eqnarray}\displaystyle & \displaystyle k_{z1}=\frac{1}{2{\pi}r}\int _{0}^{{\pi}}\bigg(s_{z}\bigg|_{r_{i}=r_{2}}\bigg|_{r^{\prime }=r_{1}}^{r^{\prime }=r_{2}}\bigg)∼\text{d}{\phi}^{\prime } & \displaystyle\end{eqnarray}$ (26) $\begin{eqnarray}\displaystyle & \displaystyle k_{z2}=-\frac{1}{2{\pi}r}\int _{0}^{{\pi}}\bigg(s_{z}\bigg|_{r_{i}=r_{1}}\bigg|_{r^{\prime }=r_{1}}^{r^{\prime }=r_{2}}\bigg)∼\text{d}{\phi}^{\prime } & \displaystyle\end{eqnarray}$ (27) where k _z1 and k _z2 are the geometric factors associated with the axial field. The expression of s _z is obtained after integrating with respect to r ^′ and performing a significant amount of simplifications. $\begin{eqnarray}\displaystyle s_{z}={\xi}_{24}\cdot \left(\begin{array}{@{}l@{}}\displaystyle \frac{r{\xi}_{22}}{{{\xi}_{6}}^{2}}+\frac{r\log ({\xi}_{13})({\xi}_{20}-4{{\xi}_{3}}^{3}{\xi}_{6}r_{i}\cos ({\phi}^{\prime }))}{{{\xi}_{6}}^{5/2}}\\[13.0pt] \displaystyle \quad +\frac{r{\xi}_{19}({\xi}_{23}+2{\xi}_{3}{\xi}_{6}r_{i})}{\sqrt{{\xi}_{10}}{{\xi}_{6}}^{2}}+\frac{\log ({\xi}_{13})({\xi}_{18}-4{\xi}_{11}{\xi}_{3}{\xi}_{6}r_{i})}{{{\xi}_{6}}^{5/2}}\\[15.0pt] \displaystyle \quad +\frac{2\sqrt{{\xi}_{10}}({\xi}_{23}+2{\xi}_{3}{\xi}_{6}r_{i})}{{{\xi}_{6}}^{2}}+\frac{r{\xi}_{21}({\xi}_{18}-4{\xi}_{11}{\xi}_{3}{\xi}_{6}r_{i})}{{\xi}_{13}{{\xi}_{6}}^{5/2}}\end{array}\right). & & \displaystyle\end{eqnarray}$ (28) The auxiliary variables 𝜉_i have been used for simplicity, whose expressions are available in Appendix A.2.

2.2. Numerical integration with respect to 𝜙^′

It can be seen in (18), (19), (22), (23), (26) and (27) that the geometric factors are expressed in terms of an integral between 0 to π. This is because the conic configuration produces mathematical expressions that have not known analytic antiderivatives with respect to 𝜙^′.

At this point, we could resort to the use of special functions such as Jacobi, Hankel functions and elliptical integrals. Although these functions can be evaluated using efficient iterative algorithms such as the Arithmetic Geometric Mean (AGM) [11], the fact that they must be evaluated for thousands or millions of interactions compromises the computational performance. That is why a matrix based integration method (MBIM) is proposed in this work, which allows the efficient evaluation of all interactions in a single matrix product. See Appendix A.1 for more details.

3. Results and validation

In order to compare the computational time and the precision of the proposed formulas, a ferrite core has been used as a case study, which is under the influence of the external homogeneous magnetic field produced by a Helmholtz coil [15]. The objective of this section is to show the results of the calculation of the magnetic field distribution in the vicinity of the core by using the proposed formulation combined with SAIM and compare them with those obtained with a commercial tool based on FEM.

3.1. Description of the case study

The case study consists of two ferrite pieces which are separated by a distance g (see Fig. 2). This type of ferrite core configurations can be very useful to shield the cavity from external magnetic fields [16].

Several characteristics of this case study are most suitable to be modeled using conical elements: a. The geometry of the problem is axisymmetric, which is very convenient to represent using conical elements. b. The local field close to the ferrite core is of great interest, especially inside the cavity. The field close to external sources is not of interest. c. It is an open boundary problem, which implies that the boundaries of the problem extend to infinity.

The case study has been depicted in Fig. 2. The relative magnetic permeability of both ferrites has been assumed to be 𝜇_r = 100 and they form a spherical cavity in their interior. Although the magnetic permeability of the core is assumed to be linear and isotropic in this work, the effect of saturation will be introduced in the model in the future.

Fig. 2.

Cutaway drawing of the core of the case study.

Both ferrites are separated by a gap of g = 2 mm, the internal radius of the cavity is 10 mm, the external radius of the core is 20 mm and the total height is 60 mm. The core has a chamfer at 45 degrees at the top and bottom. The core is excited by a Helmholtz coil composed of filamentary conductors located at the coordinates r − z (100, −50) mm and (100, 50) mm respectively. The Helmholtz coil has been modeled by two 10 kA filamentary azimuthal currents.

3.2. Discretization and solution of the case study

Figure 3 shows the final mesh generated by Ansys Maxwell, version 17.2.0. This figure shows only the part of the mesh on a rectangular region near the upper core of ferrite, however it should be noted that the mesh extends beyond the conductors of the Helmholtz coil. Ansys Maxwell needed a total of 8 057 triangular elements to cover the entire solution region creating a system of equations of 16 177 degrees of freedom.

Fig. 3.

Discretization of the problem using Ansys Maxwell (FEM based program).

Fig. 4.

Discretization of the problem using conical elements.

Figure 4 shows the discretization using conical elements. There is a conical element between each pair of consecutive nodes. Note that conical elements are only necessary to be located on the border of the ferrite core and an amount of 57 × 2 elements has been used. The evaluation points associated with a conical element have been represented with x labels in red. The following boundary condition is defined on these field points [9] $\begin{eqnarray}\displaystyle \hat{\mathbf{n}}\times \vec{\boldsymbol{H}}_{b}=\frac{1}{2}\left(\frac{{\mu}_{r}+1}{{\mu}_{r}-1}\right)\vec{\boldsymbol{K}}_{b} & & \displaystyle\end{eqnarray}$ (29) where $\hat{\mathbf{n}}\times \vec{\boldsymbol{H}}_{b}$ is the field intensity tangential to the element and $\vec{\boldsymbol{K}}_{b}$ is the surface current density of the element, while the normal vector $\hat{\mathbf{n}}$ is pointing in all the elements towards the interior of the core.

On the other hand, the field points where the field is to be evaluated have been represented with x labels in blue. These points form a rectangular grid of 20 × 20 evaluation points.

A linear system of equations of the form k ⋅ K = s is set up by means of the boundary condition (29), where the matrix k contains the geometric factors and boundary conditions, the vector s the excitations (currents flowing in the Helmholtz coil) and the vector K the values of the current densities at the ends of each conical element which are initially unknown. The average total calculation time using conical elements was 1.27 s while FEM needed 5.78 s average time using an adaptive error of 0.1%. A 2.8 MHz Core i7 desktop, 8 GB RAM was used for calculating both solutions.

The procedure to set up the system of equations is described in greater detail in section III.B. of [9]. Once the system of equations is solved, the value of the magnetic vector potential in each evaluation point is determined in order to calculate the magnetic flux. Further details about the field matrix calculation from the current density solution values are given in section 3.1 of [10].

3.3. Results using FEM and the proposed formulation

Figure 5 shows the magnetic flux distribution in the vicinity of the ferrite calculated with Ansys Maxwell. It can be seen from the field distribution the shielding effect provided by the ferrite, so that the net flux is less than 50 μWb in the central part of the cavity in the region defined by 0 < r < 5 mm and −5 mm < z < 5 mm. The distortion in the flux lines due to the gap is also visible.

Figure 6 shows the distribution of magnetic flux lines calculated from the superposition of the field of conical elements. The distribution reproduces the result obtained with Ansys Maxwell very closely. It should also be noted that the same number of contour lines (60) have been plotted with both methods so that the result is visually comparable.

Fig. 5.

Distribution of magnetic flux lines calculated using Ansys Maxwell (FEM based).

Fig. 6.

Distribution of magnetic flux lines calculated by superimposing the field produced by conical elements.

Table 1 also shows a numerical comparison between the magnetic flux density calculated with FEM and conical elements in a 5 × 5 grid. Absolute flux density errors in the radial and axial directions are also shown. As can be seen in the table, the absolute value never exceeds 10 mT.

Table 1

Comparison between magnetic flux density calculated with FEM and using conical elements in a 5 × 5 grid

r	z	$B_{r}^{fem}$	$B_{r}^{cone}$	𝜖_Br	$B_{z}^{fem}$	$B_{z}^{cone}$	𝜖_Bz
(mm)	(mm)	(mT)	(mT)	(mT)	(mT)	(mT)	(mT)
0.00	0.00	0.78	−0.30	1.08	56.00	59.58	3.58
0.00	8.75	−0.02	0.01	0.03	25.09	24.96	0.13
0.00	17.50	−0.13	−0.16	0.03	268.71	269.38	0.67
0.00	26.25	0.63	−0.25	0.88	248.82	258.25	9.43
0.00	35.00	0.63	0.44	0.19	187.99	188.24	0.25
6.25	0.00	1.29	0.89	0.40	88.63	89.37	0.74
6.25	8.75	−194.00	−195.45	1.45	237.60	241.50	3.90
6.25	17.50	−12.32	−11.84	0.48	296.70	296.81	0.11
6.25	26.25	28.70	28.33	0.37	288.33	288.01	0.32
6.25	35.00	24.68	23.46	1.22	189.94	189.72	0.22
12.50	0.00	−0.02	0.00	0.02	377.23	375.00	2.23
12.50	8.75	−52.29	−53.38	1.09	390.03	390.12	0.09
12.50	17.50	7.61	7.47	0.14	321.01	321.71	0.70
12.50	26.25	−15.56	−15.01	0.55	322.37	324.85	2.48
12.50	35.00	54.88	55.66	0.78	163.55	163.67	0.12
18.75	0.00	3.70	−1.05	4.75	366.55	368.33	1.78
18.75	8.75	1.81	1.86	0.05	439.54	440.40	0.86
18.75	17.50	82.77	83.28	0.51	337.47	338.75	1.28
18.75	26.25	93.79	94.75	0.96	121.09	120.35	0.74
18.75	35.00	52.58	50.60	1.98	126.22	126.32	0.10
25.00	0.00	−6.50	1.24	7.74	66.55	67.72	1.17
25.00	8.75	10.93	12.77	1.84	39.37	40.36	0.99
25.00	17.50	57.35	58.56	1.21	54.07	54.56	0.49
25.00	26.25	59.40	58.71	0.69	97.23	95.73	1.50
25.00	35.00	39.04	38.37	0.67	109.56	107.25	2.31
			max:	7.74			9.43

The small differences between one method and another can be attributed to the fact that the value of the magnetic flux density calculated with FEM is obtained by interpolation of nodal values, while the value calculated with conical elements corresponds to the magnetic density evaluated at the given field point coordinates.

4. Conclusions

A mathematical formulation of a new conical surface element for the modeling of axisymmetric magnetic field problems has been presented. The corresponding geometric factors due to a linear surface element with conical geometry has been derived. The proposed formulation is used in conjunction with SAIM to determine the field distribution due to a Helmholtz coil in the vicinity of a conical ferrite. The results found with the proposed formulation for this case study have been compared with the results calculated by a commercial program based on FEM and excellent agreement has been found. Additionally, the calculation time using the proposed formulation is about 4 times less than using FEM. The most important contribution of this work is that the proposed formulation allows for the modeling of magnetic field axisymmetric problems with curved geometries. The most important advantage offered by modeling with surface elements is that only the object surface must be discretized, which considerably reduces the size of the problem and computation time compared to that of traditional methods such as FEM. Finally, the Matrix Based Integration Method has been introduced, allowing the determination of all integrals associated with multiple field sources on multiple evaluation points in a structured way by a single matrix product.

Footnotes

Acknowledgements

Special thanks to the Research and Transfer Vice-Rectory of the La Salle University - Bogota Colombia, who supported the publication of this work within the framework of the call for international visibility.

Acknowledgments to the engineers Diego Celis and Nelson Caicedo, who supported the simplification of mathematical expressions during their undergraduate studies.

Appendix

1

vec(M) is equivalent to the M(:) command in Matlab / Octave.

2

mat(M)_m×n is equivalent to the reshape(M,m,n) command in Matlab/Octave.

3

The GaussianQuadratureWeights function of Wolfram Mathematica was used to determine the weights and nodes of the quadrature.

References

Bai

Z.P.

, Numerical Analysis of Electromagnetic Fields, Springer, 1993.

Rylander

Ingelstrom

and Bondeson

, Computational Electromagnetics, Springer Science & Business Media, 2012.

Sun

Chen

Yang

and Li

, Thermal modeling and analysis of bearingless permanent magnet synchronous motors, International Journal of Applied Electromagnetics and Mechanics (1) (2018), 115–130.

Bao

Fang

and Wu

, Analysis of electromagnetic torque of induction motors with two different doubly skewed rotors, International Journal of Applied Electromagnetics and Mechanics (4) (2017), 583–595.

Barglik

Smalcerz

Przylucki

and Doleel

, 3d modeling of induction hardening of gear wheels, Journal of Computational and Applied Mathematics 270 (2014), 231–240, Fourth International Conference on Finite Element Methods in Engineering and Sciences (FEMTEC 2013). [Online]. Available: https://www-sciencedirect-com.web.bisu.edu.cn/science/article/pii/S0377042714000478.

Chen

Gao

Kang

Fang

Cui

Teng

and Lv

, Analytical calculation of no-load magnetic field of external rotor permanent magnet brushless direct current motor used as in-wheel motor of electric vehicle, IEEE Transactions on Magnetics PP(99) (2018), 1–6.

Min

K.C.

Choi

J.Y.

Kim

J.M.

Cho

H.W.

and Jang

S.M.

, Eddy-current loss analysis of noncontact magnetic device with permanent magnets based on analytical field calculations, IEEE Transactions on Magnetics 51(11) (2015), 1–4.

Gillon

Gong

and Bracikowski

, Presentation of a novel transverse-flux permanent magnet linear motor and its magnetic field analysis based on Schwarz Christoffel mapping method, IEEE Transactions on Magnetics 54(3) (2018), 1–4.

Diaz

and Mombello

, Semianalytic integral method for fast solution of current distribution in foil winding transformers, Magnetics, IEEE Transactions on 51(9) (2015), 1–9.

10.

Díaz

Mombello

and Stephan

, Novel technique for the calculation of eddy current losses and lorentz forces in foil winding transformers, International Journal of Applied Electromagnetics and Mechanics (2017), (in press).

11.

Diaz

and Mombello

, New compact and singularity free formulations for the magnetic field produced by a finite cylinder considering linearly varying current density, International Journal of Applied Electromagnetics and Mechanics 50(2) (2016), 483–501.

12.

Diaz

and Mombello

, Magnetic field due to a finite current carrying disk considering a variable current density along its radial dimension, International Journal of Applied Electromagnetics and Mechanics 42(1) (2013), 119–136.

13.

Knoepfel

H.E.

, Magnetic Fields: A Comprehensive Theoretical Treatise for Practical Use, John Wiley & Sons, 2008.

14.

Granville

W.A.

Smith

P.F.

Longley

W.R.

Elements of the Differential and Integral Calculus, Ginn and Co., 1934.

15.

Trout

S.R.

, Use of Helmholtz coils for magnetic measurements, IEEE Transactions on Magnetics 24(4) (1988), 2108–2111.

16.

Park

H.H.

Kwon

J.H.

Kwak

S.I.

and Ahn

, Effect of air-gap between a ferrite plate and metal strips on magnetic shielding, IEEE Transactions on Magnetics 51(11) (2015), 1–4.

17.

Robinson

, Toward an optimal algorithm for matrix multiplication, SIAM News 38(9) (2005), 1–3.

Mathematical formulation of a new linear surface element for the modeling of axisymmetric magnetic field problems

Abstract

Keywords

1. Introduction

1.1. Problem statement

1.2. Basic mathematical expressions

2.1. Analytical integration with respect to r ′

2.1.1. Magnetic vector potential

3. Results and validation

3.1. Description of the case study

Footnotes

Acknowledgements

Appendix

1

2

3

References

2.1. Analytical integration with respect to r ^′