An Alternative Method for Calculating the Eddy Current Loss in the Sleeve of a Sealless Pump

A sealless pump, also known as a wet rotor pump or a canned pump, requires a stationary sleeve in the air gap to protect the stator from a medium that flows around the rotor and the pump impeller. Since the sleeve is typically made from a non-magnetic electrically conductive material, the time-varying magnetic flux density in the air gap creates an eddy current loss in the sleeve. Precise assessment of this loss is crucial for the design of the pump. This paper presents a method for calculating the eddy current loss in such sleeves by using only a two-dimensional (2D) finite element method (FEM) solver. The basic idea is to use the similar structure of Ampère’s circuital law and Faraday’s law of induction to solve eddy current problems with a magnetostatic solver. The theoretical background behind the proposed method is explained and applied to the sleeve of a sealless pump. Finally, the results obtained by a 2D FEM approach are verified by three-dimensional FEM transient simulations.


Introduction
An electrical machine, either synchronous or induction, with an internal rotor and a stationary sleeve fixed to the inner side of the stator is typically used as a sealless pump [1]. Optionally, an additional sleeve is mounted on the rotor to reduce the friction caused by rotor cavities. Sealless pumps are used in circulating systems that require a leak-proof enclosure, where glands and seals may not be reliable enough. Although such pumps have clear advantages in hazardous applications, they face the challenge of finding a suitable non-magnetic material for the sleeves. Since electrically non-conductive materials such as carbon graphite may not be suitable in certain applications, materials with a low electrical conductivity such as titanium, silicone steel, SUS304, SUS316L, Inconel 718, Hastelloy C, or stainless steel 1.4571 might be used [2][3][4]. Regardless of the type of the machine, the time-varying magnetic flux density in the air gap leads to substantial eddy current losses in electrically conductive sleeves. These losses can significantly outweigh the copper and iron losses, thereby reducing the efficiency of the machine. Thus, minimizing the eddy current loss in each such sleeve should be the main focus of the design and optimization of every sealless pump [4,5].
An assessment of the eddy current loss in thin electrically conductive sleeves can be approached purely theoretically as in [6]. In this work the application of the field theory leads to a classical boundary-value problem that accounts for the dimensions of the sleeve and the length of the overhang. The eddy current loss problem is commonly formulated numerically and adopted for use with boundary and finite element methods [7,8]. On the contrary, in some cases a volume integral formulation using facet elements is translated into an equivalent lumped element network as demonstrated in [9].
Analytical models are typically fast to evaluate, but often lack accuracy due to some geometrical assumptions and a limited number of spatial harmonics of the magnetic flux density considered in the calculation of the eddy current loss. On the other hand, the twoor three-dimensional (3D or 2D) finite element method (FEM) provides a variable accuracy that depends on the number of finite elements employed in the calculation. As such, it can provide a high accuracy at the expense of the computation time and complexity. This paper tries to make a compromise between speed and accuracy by using a magnetostatic 2D FEM solver while taking advantage of the similarity between the differential forms of Faraday's law of induction and Ampère's circuital law for magnetostatics. The functionality of the proposed method is demonstrated for a permanent magnet synchronous machine with embedded rotor magnets. The main advantage of the proposed method is its easy integration in the simulation workflow of genetic algorithms for multi-objective optimizations, such as those in the System Model Space (SyMSpace) [10].

Structural Similarity of Ampère's Circuital Law and Faraday's Law of Induction
Ampère's circuital law defines the relationship between the magnetic flux density (B), the current density (J), and the electric field (E). Its differential form can be written as where µ 0 is the magnetic permeability of free space and 0 is the electrical permittivity of free space. The term 0 ∂E ∂t is also referred to as the displacement current. In electrical machines this displacement current can be neglected if conductors with high electrical conductivity σ are used, such that σ ω (2) applies, where ω denotes the angular speed. Therefore, the electromagnetic problem is reduced to a magnetostatic one which can be easily handled by a finite element solver [11]. For a more efficient and smart simulation workflow it is beneficial to be able to solve eddy current problems by using a magnetostatic solver. Solutions of eddy current problems are obtained by solving Faraday's law of induction, which can be written in its differential form as Equation (4) describes a time-varying problem which cannot be directly solved by a magnetostatic solver. Nevertheless, the structure of (3) and (4) can be generalized in the form of ∇ × x = y.
Therefore, the solution of eddy current problems (4) can be evaluated using a solver for magnetostatic problems (3). For this approach, the field quantities must be reassigned according to After substituting J by the right-hand side of (4), the numerical value of the solution of (3) refers to E instead of B.

Calculation of Eddy Current Loss
The magnetic energy (W m ) can be obtained from the magnetic energy density (w m ) in the volume of the sleeve (V) as In general, w m can be obtained from the magnitude of the magnetic field strength (H) as In magnetically non-conductive materials used for the aforementioned sleeves, the relationship between H and the magnitude of B (B) can be taken as B = µ 0 H, so that (8) can be written in terms of B as Thus, if the solution of (3) obtained by a magnetostatic solver refers to E, then (7) according to (9) can be expressed in terms of the magnitude of E (E) as and solved by an FEM solver such as FEMM [11]. The eddy current loss (P ec ) in V is defined by Joule's law based on E and the eddy current density (J ec ) or either of them via the electrical conductivity of the sleeve (σ sl ) as Assuming constant σ sl the magnetostatic solution of (10) can be combined with (11), yielding P ec = 2µ 0 σ sl W m .
The attenuation of the changing excitation field caused by the magnetic field created by the induced eddy currents that opposes the changing excitation field is neglected, which in turn limits the application of this method to a certain frequency range. This limitation is demonstrated in the simulation results.

Simulation Methodology
The proposed method is implemented in SyMSpace for the geometry shown in Figure 1. Although in Figure 1 it can be seen that the machine has 2 pole pairs (p z ) and 36 slots, the remaining machine parameters including the sleeve parameters are listed in Tables 1 and 2.  Since the axial length of the lamination stack (l Fe ) determines the nominal power of the machine (P n ), the three variants listed in Table 2 were analyzed to determine the influence of the axial length of the sleeve (l sl ) on P ec .
To calculate P ec , evaluation points shown in Figure 1 are defined in the SyMSpace project as coordinates at the mid radius of the sleeve (r sl = d slo + d sli /4) along its circumference. Each evaluation point (p k ) is defined based onr sl , the total number of evaluation points along the circumference of the sleeve (K) that needs to be specified, and the number of each evaluation point (k) as After being defined, the evaluation points are evaluated in SyMSpace by means of FEMM. From the obtained results for B in the x and the y direction (B x and B y ), the radial component of B (B rad ) at each evaluation point is calculated as while the tangential component of B (B tan ) is calculated as Based on (14) and (15), variations in B rad and B tan alongr sl at the specified evaluation points can be visualized, which for different load points is shown in Figure 2, where the direct and the quadrature components of the stator current (i sd and i sq ) are expressed in terms of the nominal stator phase current (i sn ). For l sl = 40 mm i sn amounts to 6.65 A, while for l sl = 80 mm it amounts to 13.34 A, and for l sl = 120 mm to 20.46 A. The slotting effects can be seen as the peaks in both B rad and B tan at the beginning and the end of each pole shoe. Harmonic spectra of the waveforms of B rad and B tan in Figure 2 are shown in Figure 3, where h o denotes the harmonic order.
Due to the negligible radial thickness of the sleeve, B rad is the component mainly responsible for P ec . Hence, the partial derivative of B with respect to time in (6) refers to B rad , which needs to be observed as a function of the electrical angular position of the rotor (ϕ e ), as shown in Figure 4.     Based on ϕ e and the electrical angular speed (ω e ), which corresponds to the electrical frequency of the rotating magnetic flux density in the air gap, the partial derivative of B rad with respect to time can be calculated as With the total number of evaluation steps per electric period (M) the partial derivative of B rad with respect to the rotor angle can be approximated by For each p k , (16) can be written in terms of B rad , p z , the mechanical speed (n) and the number of each evaluation step (m) as After obtaining the partial derivative of B rad with respect to time, a new 2D magnetostatic problem in FEMM, in which the sleeve is unrolled as shown in Figure 5, is used for the calculation of P ec .  The sleeve is split up into small surface areas, each with the circumferential width of 2πr sl /K and the axial length of l sl , as it can be seen in Figure 5. The sleeve can be axially positioned, and l sl and l Fe do not have to be equal. The eddy current excitation at p k (i ec,k ) for the calculation of P ec in the new 2D magnetostatic problem in FEMM is obtained based on the corresponding magnitude of J ec (J ec,k ) in the surface area (A) between two neighboring evaluation points according to (6) and (18) as Finally, the average value of P ec that represents the average 2D magnetostatic solution of FEMM (P ec,2D ) is obtained from the value of P ec at each evaluation step during a full electric period (P ec,m ) asP where each P ec,m is obtained according to (12).

Simulation Results
Simulations were performed for K = 180 (evaluation points) and M = 72 (evaluation steps), which in the presented case seemed like a good tradeoff between the accuracy and the calculation speed. The results ofP ec,2D for l sl of 40 mm, 80 mm, and 120 mm, at the same load points for which B rad is presented in Figures 2 and 4 are listed in Table 3. In Table 3 it can be noticed that using only i sq increases P ec , while negative values of i sd reduce P ec . The latter implies that topologies with embedded magnets, which are suitable for field weakening to use the reluctance torque, might be better suited for these applications.

Verification of 2D Simulation Results by 3D Simulations
For a validation of the 2D simulation results obtained by the proposed method, the 2D model from SyMSpace was implemented as a 3D model in Ansys Electronics Desktop 2019 R3 (AED) [12] and used for 3D FEM transient simulations to obtain reference values of P ec . Since it is sufficient to observe the no-load case, the AED model was simulated without the stator windings by using both tangential and axial symmetry. That means that a single pole over a half of the axial length of the machine was simulated, as shown for l sl = 40 mm in Figure 6, to minimize the computation time. The distribution of J ec in the sleeve (colored in teal in Figure 6) is shown for different values of l sl in Figure 7. The transient 3D FEM simulation results of P ec (P ec,3D ) obtained by AED at 4000 rpm for l sl of 40 mm, 80 mm, and 120 mm, are shown in Figure 8. To avoid the influence of the transient character of P ec,3D in Figure 8, the mean average value of P ec,3D (P ec,3D ) over the last half of the electrical period of each characteristic was taken as the representative value of each 3D simulation. A comparison of the no-load results ofP ec,2D listed in Table 3 withP ec,3D is given in Table 4, where the average magnitude of J ec in V (J ecV ) was obtained from AED, while the relative 2D percent error in P ec (ε 2D ) was calculated as ε 2D =P ec,2D −P ec,3D P ec,3D · 100%.
(21) l sl /2 l sl /2r sl r sl Figure 6. The 3D AED model used for the calculation ofP ec,3D in the sleeve.
The values ofP ec,3D listed in Table 4 are smaller thanP ec,2D because the transient solver in AED takes into account the influence of the magnetic field created by the induced eddy currents on the excitation field. From Table 4 can also be seen that by doubling l sl ,J ecV increases by about 50%, which in this case explains the quadratic behavior, since P ec is defined by (11) where according to the microscopic form of Ohm's law E = J ec /σ sl .

Valid Frequency Range of the Proposed Method
The application of the proposed method is limited due to the neglect of the attenuation of the changing excitation field caused by the magnetic field created by the induced eddy currents that opposes the changing excitation field. To demonstrate the frequency limitation P ec,2D andP ec,3D were calculated for l sl of 40 mm at different frequencies in the range from 10 Hz to 10 kHz and presented in Figure 9.
In Figure 9 it can be seen that the magnetic field created by the induced eddy currents in the sleeve begins to attenuate the changing excitation field more significantly as the frequency increases above 2 kHz, which in turn reducesP ec,3D .  The frequency range presented in Figure 9 corresponds in this case to the range from 300 rpm to 300 krpm, which is significantly above n n specified in Table 2 as well as the upper limit of typical industrial applications. Hence, for a meaningful comparison,P ec,2D andP ec,3D together with their difference (∆P ec =P ec,2D −P ec,3D ) and ε 2D are shown in Figure 10 in the range from 300 rpm to n n for l sl of 40 mm, 80 mm, and 120 mm.
From the presented results it can be seen that the discrepancies betweenP ec,2D and P ec,3D increase with a decrease in l sl . The main cause of that is presumably the inability of the 2D FEM solver to calculate the axial stray magnetic field.
The preference for axially long sleeves in 2D calculations of P ec is also implied by the axial-to-radial ratio ( h ) presented per harmonic (h = h o /p z ) in [13] as which is a part of a relatively accurate analytical description of P ec (P ec,h ) for all harmonics of B rad presented in the form of where the overhang correction factor (k ξ,h ) is defined in terms of h and the overhang factor (ξ) as in which ξ defines the relationship between l sl and l Fe as From (23) it can be seen that P ec increases with the square of ω e and B rad that often cannot be reduced. More importantly, P ec increases with the cube ofr sl , what should be considered in the design.
An additional comparison of no-load simulation results ofP ec,2D ,P ec,3D , and P ec,h obtained by (23) from the harmonic spectra of B rad shown in Figure 3 is shown in Figure 11 for l sl in the range from 40 mm to 200 mm at 4000 rpm together with ε 2D and the relative analytical percent error in P ec (ε h ). The lack of smoothness in ε 2D and ε h in Figure 11 is primarily caused by the coarseness of the mesh of the 3D FEM solver that creates noise in the simulation results, which can also be seen in Figure 8. Nevertheless, it is clearly visible that the relative error decreases with increasing sleeve length due to the decreasing impact of the neglected end effects.
From the tendency of ε 2D and ε h it can be seen that the discrepancies between P ec,h and P ec,3D for values of l sl smaller than approximately 70 mm are greater than those between P ec,2D andP ec,3D in the same range. For l sl greater than approximately 160 mm, the share of the axial stray magnetic field in the total excitation field is relatively small, which results in the absolute value of ε 2D below 0.25%.  l sl / mm Figure 11. A no-load comparison ofP ec,2D ,P ec,3D , and P ec,h for l sl in the range from 40 mm to 200 mm with ε 2D and ε h at 4000 rpm.

Conclusions
This paper has successfully demonstrated that a 2D FEA solver for magnetostatic problems can be used for eddy current loss calculation. The method is based on the similarity between the differential forms of Faraday's law of induction and Ampère's circuital law for magnetostatics. By reassigning the field quantities, the numerical solution of B refers to E. Therefore P ec can be evaluated efficiently for a thin sleeve situated in the air gap of a sealless pump.
According to the simulation results listed in Table 3, using only i sq increases P ec in the sleeve primarily due to the increase in the fundamental harmonic but also because of the increased harmonic content of B rad introduced by the stator magnetic flux density, as it can be seen in Figure 3. In contrast, field weakening reduces P ec in the sleeve due to the drop in the fundamental harmonic of B rad mainly produced by the rotor magnets, which can as well be seen in Figure 3. This also suggests that rotor topologies with embedded magnets that exhibit reluctance torque might be better suited for these applications when a control strategy that uses field weakening is being employed. However, increasing i sd and i sq increases the harmonic content of B tan , which can be seen in Figure 3.
For the frequency range of typical industrial applications in which our proposed method is intended to be applied, the reaction field can be neglected. Furthermore, it is derived in (12) and (23) that the relationship between eddy current loss and sleeve conductivity σ sl is linear. Combined with the fact that the reaction field is negligible it can be concluded that the eddy current loss increases linear with the sleeve conductivity.
The comparison of 2D and 3D no-load simulation results listed in Table 4 shows that ε 2D increases as l sl decreases. That is mainly caused by the inability of the 2D FEM solver to correctly calculate the share of the axial stray magnetic field of the sleeve in the total excitation field. For smaller values of l sl these end effects have a higher fraction relative to the actual total excitation field, thereby making P ec seemingly higher. From (23) it can also be seen that P ec increases with the square of ω e and B rad and with the cube ofr sl .
Since ω e and B rad may not be easily reduced,r sl is an important parameter for an efficient pump design.