**Enhanced Flexible Algorithm for the Optimization of Slot Filling Factors in Electrical Machines** †

**Armin Dietz 1, Antonino Oscar Di Tommaso 2,\*, Fabrizio Marignetti 3, Rosario Miceli <sup>2</sup> and Claudio Nevoloso <sup>2</sup>**


Received: 14 January 2020; Accepted: 24 February 2020; Published: 26 February 2020

**Abstract:** The continuous development in the field of industrial automation and electric mobility has led to the need for more efficient electrical machines with a high power density. The improvement of electrical machines' slot filling factors is one of the measures to satisfy these requirements. In recent years, this topic has aroused greater interest in the industrial sector, since the evolution of the winding technological manufacturing processes allows an economically sustainable realization of ordered winding arrangements, rather than random ones. Moreover, the manufacture of electrical machines' windings must be preceded by an accurate design phase in which it is possible to evaluate the maximum slot filling factor obtainable for a given wire shape and for its dimensions. For this purpose, this paper presents an algorithmic approach for the evaluation of maximum slot filling factors in electrical machines under an ideal geometric premise. In particular, this algorithm has a greater degree of flexibility with respect to the algorithm approaches found in the literature, since the study has been extended to round, rectangular and hexagonal wire sections. Furthermore, the slot filling factor calculation was carried out both for standard and non-standard slots. The algorithmic approach proposed can be considered as an additional useful tool for the fast design of electrical machine windings.

**Keywords:** electrical motors; sot filling factor; optimization algorithm; windings; magnetic wire; filling factor optimization

#### **1. Introduction**

The development of more and more efficient electrical machines has become a topic of interest for various industrial sectors, such as automation or electric traction. In particular, high efficiency, high power density and cost-effective manufacturing are required in the automotive industry [1–3]. One possible solution to meet these requirements is to optimize the copper filling factor of stator winding [4–11]. In particular, a high copper filling factor involves a more rational and efficient use of copper with economic benefits and improved energy savings. Therefore, the optimization of the slot filling factor is a key focus in winding technology. The improvement of the slot filling factor depends mainly on the winding pattern schemes and the adopted winding manufacturing process.

The simplest type of winding pattern is so-called "random winding". In this case, the random winding process is sustainable for mass production and it is characterized by low manufacturing requirements. The main advantage is represented by its high production speed, while the disadvantages consist, generically, of lower filling factors. The highest possible filling factor is achieved by the "orthocyclic winding pattern" for round wires. This winding pattern presents a high packing density of wire, but the winding process is more complex and, therefore, more costly than the random winding process. In particular, the orthocyclic winding process requires very high manufacturing requirements to obtain an ordered wire positioning within each slot. Another possible winding pattern is the "layer winding", where the wires are uniformly arranged in layers [12]. This winding pattern allows us to obtain higher filling factors than those of the random winding pattern. The choice of typology winding pattern depends on the functions and the design requirements of the electrical machine. Therefore, orthocyclic (or ordered) winding structures are appropriate for high power density applications.

In the past, the realization of commercial solutions that allow the economically sustainable manufacturing of distributed windings with ordered structures was very difficult due to the high economic burdens. The automated winding process technologies available for making distributed windings are: insertion winding technology, flyer winding technology and needle winding technology [13]. In recent years, these winding process technologies have undergone a great technological evolution that has allowed us to reduce manufacturing costs and winding process times with an ordered structure [13,14]. In [12], a new and innovative needle winding method that allows shifted layer winding structures for distributed round wire applications is described, thus significantly increasing the copper filling factor. In this regard, it is of considerable interest to accurately evaluate the slot filling factor obtainable during the design phase of electrical machines. In the previous century, this task was performed by means of manual graphic analysis or by testing the stator of the electrical machine during the pre-production phase. In order to carry out this process, the resources in terms of money, time and technical staff are not indifferent [15–17]. Therefore, a preliminary analysis of the maximum possible value of filling factor is an important step forward for the design of electrical machine windings for a given slot, wire shape and set of dimensions.

In the scientific literature, there are several definitions of slot filling factors and, generally, an electrical slot filling factor *fcu* and a mechanical slot filling factor *fme* are defined. In this study, the electrical slot filling factor *fcu* is given by the ratio between the total wire copper cross-section (*Nw Acu*) and the total slot cross-section *Aslot*. Furthermore, the mechanical slot filling factor *fme* is defined as the ratio between the total wire cross-section (*Nw Aw*) and the effective slot area *Ae*ff, which is defined as the difference between the total slot cross-section *Aslot* and the area occupied by slot insulation. The two filling factors can be described by the following mathematical equations:

$$f\_{cu} = \frac{N\_w \ A\_{cu}}{A\_{slot}}\tag{1}$$

$$f\_{m\epsilon} = \frac{N\_{w\cdot}A\_{w}}{A\_{eff}}\tag{2}$$

where *Nw* is the number of wires contained within a slot, *Acu* is the maximum copper cross-section of a single wire (useful cross-section) and *Aw* is the maximum cross-section of a single wire (including its insulation layer). For the purpose of this work, the electrical slot filling factor values are considered and discussed. This paper proposes an algorithmic approach for the preliminary determination of the slot filling factors. In particular, the goal of the present paper is to extend the work done in [18,19] with reference to various types of wires (this work also considers hexagonal cross-sections) and slots. Moreover, the identification of the best arrangement of wires within the slot for each slot type is suggested, with a technique aimed at the reduction of computation times. In fact, the algorithmic approach described in the previous papers calculates the slot filling factors only for circular and rectangular wire shapes. This algorithmic approach requires the definition of the coordinates of the first wire, that can be useful for the investigation of singular cases where the goal is the evaluation of

the slot filling factors for a specific wire's positioning. In fact, to research the best wire arrangement and, consequently, the highest obtainable values of the slot filling factors, the algorithm requires several executions of the procedure to obtain several values of the first wire position. Since each execution of the algorithm requires a few seconds, the time consumed can be significant.

However, the enhanced algorithm, here described, allows the determination of the maximum slot filling factor for a given slot, wire shape and set of geometrical dimensions. In particular, the algorithm procedure is designed in order to systematically investigate a large number of cases defined for several values of the coordinates of the first wire in order to optimize the computation time and find the best arrangement inside the slot, with an improved accuracy. In more detail, with a given cross-section geometry and given wire dimensions, the algorithm allows the calculation of the maximum number of wires that can be placed inside the slot. In this work, both wire and slot insulation are considered because, together, they cause a reduction of the useful slot area. In detail, with respect to the previous work, three different wire geometries, namely round, rectangular and hexagonal, are taken into consideration and the results are discussed. Since the use of hexagonal wires can present an innovative character in the design of electrical machines, attention has been paid to the comparative study of the filling factors obtainable with wires that have both circular and hexagonal cross-sections, by considering the same cross-section area.

This paper is structured as follows: Section 2 describes the state-of-the-art optimization algorithms of the slot filling factors, Section 3 describes the algorithmic approach proposed and Section 4 illustrates several cases of studies that have been carried out and the relevant obtained results.

#### **2. State-of-the-Art Slot Filling Factor Optimization Algorithms**

At the present time, the scientific literature does not include many works regarding the optimization of the slot filling factors with algorithmic approaches. Furthermore, the algorithmic approaches found in the literature refer to the case where the wire has a round cross-section. The optimal winding pattern is invariant with respect to the axis perpendicular to the stator cross-section. Therefore, the determination of the optimal winding pattern can be carried out with a bi-dimensional approach. A family of widespread algorithms is that of the orthocyclic windings algorithms (OWA). In [20], the results of this kind of algorithm are discussed. A study was carried out on the variation in the value of the mechanical and electrical slot filling factors when the radius of the circular conductor varies. In the first phase, the algorithm sets the coordinates of the center (*xi*; *yi*) of the first wire:

$$w\_i = \begin{pmatrix} x\_i \\ y\_i \end{pmatrix}. \tag{3}$$

Then, the wires of the same layer are plotted and the distance between the centers of two adjacent wires is equal to the diameter of the wire *d*. Once the first layer is finished, the wire coordinates of the upper layer are determined through the following relation:

$$w\_j = \begin{pmatrix} \cos \alpha \\ \sin \alpha \end{pmatrix} d + w\_i \tag{4}$$

with α equal to 60◦ or 120◦ (Figure 1).

The algorithm ends when it is not possible to find a new wire position whose area intercepts the slot contour or its insulation contour; then, the slot is full.

**Figure 1.** Round wires with orthocyclic distribution.

Another algorithm approach, named the filling factor estimation algorithm (FFEA), is described in [15]. The goal of this algorithm is to reduce the design time of the winding by estimating the filling factors. In this case, the first wire is placed centrally at the slot bottom. The next wire is placed by checking all valid positions along a circle around the first wire, with a radius equal to the distance *dw* + *tw*, where *dw* and *tw* are the diameter of the wire and the safety distance or the minimum gap between two wires, respectively. For each new wire positioned, the algorithm checks whether there are overlaps with the contour of the slot, taking into account the presence of the slot insulation. If the new position is valid, the algorithm places the new neighbor wire and continues going around until all positions have been tried. The procedure ends when it is not possible to position wires that do not overlap the slot profile. The result of this algorithm is a random positioning of the wires forming the winding. Since the value of the determined filling factor depends on the initial position of the first wire, the procedure is repeated for several different initial positions.

Another interesting approach is described in [20]. This algorithm is based on the basic rules of the two algorithms described above. This algorithm is called the needle winding simulation algorithm (NWSA) because the objective is to optimize the positioning of each wire in order to simulate the process of needle winding. In this case, for the positioning of wire, a cost function is taken into consideration, accounting for the position and the downward force acting on the wires. In the algorithm, there are constraints to avoid overlapping between wires and between the wire and the slot profile. Optimal positioning is achieved by identifying the global minimum of the cost function. Authors use a genetic algorithm to search the global minimum in order to reduce convergence times.

Beyond the approaches described above, there are the dense packing algorithms, also called the wire inflation algorithms (WIA). The goal of these algorithms is to find the maximum radius of spheres within a given boundary and their corresponding coordinates when the number of spheres is fixed. In [21], an algorithm is described that simulates a system of billiard spheres within a limited space, whose radius is made to grow until it reaches a state whereby it becomes a blockade in the system. This algorithm approach is used for the dense packing of spheres in circular, triangular and hexagonal spaces [22–24]. This algorithm performs an event based on the physical simulation of a billiard system where the coordinates, the number and the speed of the spheres are set in advance as input data. Events are represented by the collision of each sphere with other spheres or with the borders of the system. Depending on the event type (either sphere–obstacle or sphere–sphere interaction) an elastic impact is performed, and the new velocities are calculated. At the increasing radii of the spheres, an eventual jamming occurs, resulting in a dense packing. For the purposes of maximizing the slot filling factor, the objective is to maximize the number of magnetic wires, with given radius, within the slot. In this sense, in [20], the authors, in order to make a comparison with the NWSA, have modified the described approach. In particular, the wires are understood to be charged particles exerting a Coulomb force on each other, leading to the movement of the spheres during the simulation. The results of the comparison are widely discussed.

#### **3. Proposed Enhanced Algorithmic Approach**

The basic rules of the algorithm proposed in this paper are mainly inspired by the approach proposed in [15] and the treatment of the winding patterns with ordered structures in detail [18]. The algorithm is based on a general approach that is valid for different stator slots and for different wire shapes. In detail, for round and hexagonal wires, the algorithm focuses on the windings, with orthocyclic arrangement of the wires within the slot. Moreover, the algorithm has been designed so as to be able to define some critical aspects, such as the possible presence of slot insulation, the positioning of the wire parallel to the flank or to the bottom of the slot, any safety distances between the wires forming the winding and a large number of the possible coordinates of the first wire. In this way, the algorithm presents a high degree of flexibility that makes investigating the maximum slot filling factors obtainable in many cases of study. The algorithm has been implemented in the Matlab environment and it is described in detail below.

#### *3.1. Slot Geometrical Features Definition*

Compared to the works found in the literature, where simplified models of slot have been taken into account, in this work, both a standard slot (STSL) model and the non-standard (NSTSL) ones are used. The slot profiles are defined with the presence of an insulation sheet and, therefore, its thickness *dins* is taken into account. In general, the slot insulation sheet is pre-folded and does not substantially modify the slot cross-section profile available for the positioning of the wires. The geometries and the related contour data of a standard slot, typically used for wires, and a non-standard one are shown in Figure 2a,b [16], respectively.

**Figure 2.** Cross-section of the standard slot (**a**) and non-standard slot (**b**) with contour parameters.

Generally, the slot profiles are supplied with CAD drawings by the manufacturers and, therefore, the geometrical parameters are easily determined. The dimensions of each slot are known and are reported in Tables 1 and 2.


**Table 1.** Standard slot geometric parameters.


**Table 2.** Non-standard slot geometric parameters.

The slot profiles are defined in the *xy* reference frame by means of characteristic points connected by lines and arcs. The implementation of the geometric model of the slot profile is carried out as described in [18]. The cross-sections of slots are calculated in a numeric way by the trapezium rule in the same way as [15]. In particular, the slot areas are divided into different sections whose characteristic points are known (Figure 3a,b). The calculated cross-sections show a maximum deviation of less than one percent from the cross-section given by the specifications.

**Figure 3.** Half cross-section of the standard slot (**a**) and non-standard slot (**b**) divided into sections.

#### *3.2. Magnet Wire Data*

In this study, investigations into round, rectangular and hexagonal wires have been carried out. The geometrical data of the round and rectangular wires used in this paper are obtained from [25,26]. By these standards, the nameplate data, the insulation degree, tolerances and maximum dimensions allowed for each wire shape are reported. For rectangular shapes, the standard defines the shape of the wire with rounding corners. However, there is not an international standard for hexagonal wires. The use of the latter could be an innovative idea, as they can substantially improve the filling factors. Improvements must be compared with the production costs which, compared to the past, have been reduced thanks to the technological evolution of wire manufacturing processes [13,17]. Therefore, the geometric data of the hexagonal wires will be hypothesized. The geometrical reference quantities for each wire shape are defined below for:

	- *dcu* diameter wire without insulation;
	- *dmax* maximum diameter of wire with insulation;
	- *xc* and *yc* coordinates of wire center.
	- *L*1*cu* width of the rectangular wire without insulation;
	- *L*2*cu* height of the rectangular wire without insulation;
	- *L*1*max* maximum width of the rectangular wire with insulation;
	- *L*2*max* maximum height of the rectangular wire with insulation;
	- *rcorner* corner radius;
	- *xc* and *yc* coordinates of the wire center.
	- *rcu* radius of circumference circumscribed to the hexagon without insulation;
	- *rmax* maximum radius of circumference circumscribed to the hexagon with insulation;
	- *xc* and *yc* coordinates of the wire center.
	- *phi* (φ) rotation angle of the hexagon.

Figure 4 shows the wire cross-sections with the related contour data.

**Figure 4.** Wire cross-sections with contour data.

#### *3.3. Constraints for Wire Distribution and Placement*

䭉䈟 ᵚࡠᕅ⭘ⓀDŽ In the algorithm, the slot models have been implemented in a *xy* plane with the bottom or the ground of each slot parallel to the *y*-axis. Each slot profile has been divided into several reference sections Figure 3a,b. These sections are referred to a one-half slot. The distribution and placement conditions of wires are defined for each section and applied in a specular way for each half of the slot. In order to avoid overlapping between the wire and the slot profile and between adjacent wires, several constraints are taken into account. Therefore, for each wire that must be positioned within the slot profile, the following gaps are considered and investigated:


䭉䈟 ᵚࡠᕅ⭘ⓀDŽ In Figure 5, the distances between the slot profile and the round and rectangular wires are shown. In each region of the slot profile, the minimum value of the various distances is evaluated. In particular, for round wires, the distances are defined with respect to the center of the wire and it is imposed that

$$\min\{|d\_1|, |d\_2|, \dots, |d\_n|\} \ge r\_{\max} \tag{5}$$

䭉䈟 ᵚࡠᕅ⭘ⓀDŽ

whereas, in the case of rectangular and hexagonal wires, these distances are evaluated with respect to the vertices of each wire shape (Figure 6) and it is imposed that

$$\min\{|d\_1|, |d\_2|, \dots, |d\_n|\} \ge 0. \tag{6}$$

**Figure 5.** Distances between the wires and the slot profile.

**Figure 6.** Characteristic points of the rectangular and hexagonal wires.

Regarding the distance between wires, it is possible to add an additional safety distance *tw*. The winding positioning pattern adopted in this work is the orthocyclic one that can be obtained in the case of round and hexagonal wires. However, an ordered arrangement was taken into consideration for

䭉䈟 ᵚࡠᕅ⭘ⓀDŽ

䭉䈟 ᵚࡠᕅ⭘ⓀDŽ

rectangular shaped wires. Moreover, it is possible to vary the angle formed between the line joining the centers of the wires and the horizontal one of the *xy* plane (Figure 7). In this case, the new arrangement is obtained by applying the following coordinate transformation to the center of wires for round wires, and also to the vertices for rectangular and hexagonal wires:

**Figure 7.** Rotation of the orthocyclic distribution.

In this sense, it is possible to evaluate the value of the slot filling factors in the case where the wires are arranged parallel to the flank or the bottom of the slot. This study has been performed with particular attention paid to round and rectangular wires since. in the hexagonal wire case, as described above, it is possible to set the rotation angle.

#### *3.4. The Algorithm Procedure*

As described above, the algorithm has been designed to automatically define several positions for the first wire. In this procedure, the position of the first wire is defined by its center (*xw*1, *yw*1), which it is selected in such a way that it will be located in the region of space near the lower wedge of the slot (i.e., in the lower rounded corner):

$$\mathbf{x}\_{w1} = \mathbf{x}\_{init} + d\_{ims} \tag{8}$$

$$y\_{w1} = y\_{init} + d\_{ins} - \frac{w}{2} \tag{9}$$

䭉䈟 ᵚࡠᕅ⭘ⓀDŽ where *dins*, *xint*, *yint* are the slot insulation foil thickness and the initial wire coordinates, respectively, which are chosen in such a way as to avoid overlaps with the slot profile. In order to evaluate the filling factors, several pairs of values of *xint*, *yint* are assigned. For this purpose, it is necessary to evaluate the extremes of the range, within which the initial coordinates must vary. These extremes depend on the slot geometry and the wire geometrical features. In more detail, since the slot geometry is specular with respect to the *x*-axis (see Figure 3) and the first wire is positioned in the region of space near the lower wedge of the slot, the lower extremes are chosen in such way that they are tangent with the lower slot wedge profile, whereas the upper extremes of the range are chosen so that an additional increase in the initial coordinates causes the loss of significant useful wires. Therefore, in the case of circular, hexagonal and rectangular shapes located in the STSL model, the extremes are respectively equal to:

$$\begin{array}{ll} \chi\_{\text{init\\_circ}} = d\_{\text{max}}/2 & \chi\_{\text{init\\_circ}}/hc\_{\text{max}} = \frac{3}{2}d\_{\text{max}}\\ \chi\_{\text{init\\_circ}} = d\_{\text{max}}/2 & y\_{\text{init\\_circ}}/hc\_{\text{max}} = -y\_B + \frac{3}{2}d\_{\text{max}} - d\_{\text{ins}} + \frac{y}{2} \\ \chi\_{\text{init\\_}} = L\_{1\text{max}}/2 & \chi\_{\text{init\\_circ}} = \frac{3}{2}L\_{1\text{max}}\\ y\_{\text{init\\_}} = L\_{2\text{max}}/2 & y\_{\text{init\\_}}rc\_{\text{max}} = -y\_B + \frac{3}{2}L\_{2\text{max}} \end{array}$$

where *yB* is the y-coordinate of the characteristic point B (Figure 3a). On the other hand, for NSTSL, the extremes are equal to, respectively:


The evaluation of the filling factors, when the first wire center position (*xw*1, *yw*1) varies, is carried out by using two loops with indices *k* and *z*, with maximum values *kmax* and *zmax*. They are calculated by the following relationships:

$$k\_{\max} = cei(\frac{\mathbf{x}\_{init\_{\max}} - \mathbf{x}\_{init\_{\min}}}{\delta\_t}) \tag{18}$$

$$z\_{\text{max}} = \text{cell}\left(\frac{y\_{\text{init}\_{\text{max}}} - y\_{\text{init}\_{\text{min}}}}{\delta\_t}\right) \tag{19}$$

where *ceil*(*x*) is a function that rounds *x* to the nearest integer greater or equal to *x* and δ*<sup>t</sup>* is the variation step of the initial coordinates. The last quantity starting value δ<sup>1</sup> is chosen at least 1/50 of the maximum dimension of the considered wire (*dmax* for the circular case, *L*2*max* for the rectangular case and *2rmax* for the hexagonal case). In this way, the first wire center position (*xw*1, *yw*1) varies according the following relationships:

$$\mathbf{x}\_{\rm w1}(1,k) = \mathbf{x}\_{\rm init\_{\rm min}} + d\_{\rm ins} + \delta\_{\rm l} \cdot (k-1) \tag{20}$$

$$y\_{w1}(z,1) = y\_{init\_{min}} + d\_{inv} - \frac{w}{2} + \delta\_{l} \cdot (z-1) \tag{21}$$

Furthermore, in order to investigate the possibility of finding a better wires arrangement when δ*<sup>t</sup>* varies, several iterations with a progressive halving of δ*<sup>t</sup>* are performed until to the maximum value of the number of wires positioned *nmax* presents the same value for three consecutive iterations. Obviously, this choice determines the increasing of computational time, but a more accurate research of the maximum filling factors.

In order to place the next wires, the cross-section of the slot is divided into a grid formed by *i* columns and *j* rows, where both the first column and the first row are defined as starting from the first wire position. The maximum number of columns *imax* and rows *jmax* is calculated according to the slot and wire dimensions. For round, rectangular and hexagonal wires these values are calculated by means of the following relationships, respectively:

$$i\_{\text{max}} = round \left(\frac{h}{\frac{\sqrt{3}}{2}d\_{\text{max}}}\right) \tag{22}$$

$$j\_{\text{max}} = round \left(\frac{w}{d\_{\text{max}}}\right) \tag{23}$$

$$i\_{\text{max}} = ccl \left(\frac{h}{L\_{1\text{max}}}\right) \tag{24}$$

$$j\_{\max} = ccl(\frac{w}{L\_{2\max}}) \tag{25}$$

$$i\_{\max} = cil \left(\frac{h}{\frac{\sqrt{3}}{2}r\_{\max}}\right) \tag{26}$$

$$j\_{\max} = \operatorname{coil}\left(\frac{w}{\frac{3}{2}r\_{\max}}\right) \tag{27}$$

where *round*(*x*) is a function that rounds each element of *x* to the nearest integer. The distribution of the wires occurs through two loops that change the indices *i* and *j* in order to position the wires along with the whole slot profile. In the cases of the round and hexagonal wires, in order to obtain an orthocyclic winding pattern, the coordinates of the successive wire are calculated as previously described in Section 2. The rectangular case is widely described in [13]. The positioning of the wires is considered valid only if it meets the overlapping conditions described above. The algorithm calculates, for each possible pair of coordinates (*xw*1(1,*k*) *yw*1(*z*,1)), the number of positioned wires *n*, the electrical and mechanical slot filling factor values *fcu* and *fme* and the coordinates of the positioned wires (*xwire*, *ywire*). In this way, each of the quantities of interest is represented by an array, in which each element is associated to the coordinates of the first positioned wire. The search for the best wire arrangement for a fixed wire shape, its geometric dimension, type of positioning and slot profile can be achieved by identifying the maximum value of number of wires positioned *nmax* within the array. This procedure, as mentioned above, is repeated until the quantity, *nmax*(δ*t*), presents the same value for three consecutive δ*<sup>t</sup>* values. Therefore, once the coordinates of the first wire (*xw*1, *yw*1) associated to the best wire arrangement is determined, the algorithm provides a graphical representation of the slot containing the wires. Furthermore, the coordinates of the positioned wires (*xwire*, *ywire*) are available and the manufacturer has the ability to know and choose a reference wire for the production phase according to the winding technological process adopted. In short, the algorithm works via the following steps:


In addition, for an accurate description and understanding, the algorithm flowchart is reported in Figure 8.

**Figure 8.** Algorithm flowchart.

#### **4. Cases of Study**

In order to validate the proposed algorithm, several cases of study have been performed. This study has been carried out both for the STSL and for the NSTSL and for each wire shape. From international standards [25,26], it is possible to deduce the nominal dimensions, the insulation grades, the tolerances and the outer maximum dimensions for each typology of wire. In this work, the largest allowed outer dimensions are adopted to consider the worst case for the evaluation of the slot electrical filling factor. An insulation grade (According to [25,26] the grade is defined as "the range of thickness of the insulation wire".) 3 and a grade 2 for the round-shaped wire and for the rectangular-shaped one have been taken into consideration, respectively. The dimensions of the hexagonal wire have been deduced, assuming the cross-section is equal to that of the round wire. In this way, it is possible to compare the slot filling factors and evaluate the benefits derived from the use of hexagonal wires. In this study, six different dimensions have been chosen for each wire shape, respectively. For each dimension, three additional interspaces *tw*, respectively equal to 0, 0.05 and 0.1 mm, were taken into consideration. Furthermore, regarding the standard slot, a study was carried out on the possibility of positioning the wires parallel to the ground or the bottom (PG) and to the flank (PF) of the slot. The dimensions of the wires, used for this study and expressed in mm, are reported in Tables 3 and 4.

**Table 3.** Geometrical features of round and hexagonal wires.



**Table 4.** Geometrical features of rectangular wires.

Figures 9–12, it is possible to notice how the algorithm returns the slot profile with the desired distribution of the wires and also shows the number of inserted wires. In particular, in the case of round wires, the winding pattern presents an orthocyclic structure where the positioning is parallel to the bottom in one case and parallel to the flank of the slot in the other case. In the case of rectangular wire, the winding pattern is ordered. In the case of the hexagonal wire, an orthocyclic winding pattern is obtained, both with *phi* equal to zero and with *phi* equal to π/*2*. Furthermore, it is possible to notice how the structure remains ordered with the use of an additional interspace *tw*.

**Figure 9.** Orthocyclic distribution of the round wires with the parallel ground disposition (**a**) and with the parallel flank disposition in the standard slot (**b**) (*tw* = 0 mm, *dmax* = 1.381 mm).

**Figure 10.** Ordered distribution of the rectangular wires with parallel disposition to the flank in standard slot (**a**) and in non-standard slot (**b**) (*tw* = *0.1* mm, *L*1*max* = 2.17 mm, *L*1*max* = 1.07 mm).

**Figure 11.** Orthocyclic distribution of the hexagonal wires in standard slot (*phi* = 0, *tw* = 0 mm, *rmax* = 0.313 mm) (**a**) and (*phi* = π/2, *tw* = 0 mm, *rmax* = 0.313 mm) (**b**).

**Figure 12.** Orthocyclic distribution of the hexagonal wires in non-standard slot (*phi* = 0, *tw* = 0 mm, *rmax* = 0.618 mm) (**a**) and (*phi* = π/2, *tw* = 0.05 mm, *rmax* = 0.618 mm) (**b**).

In Figure 13, the value of slot electrical filling factors obtained for the case of the round wire, with a parallel bottom disposition, both for the STSL and the NSTSL, are reported. A comparison, between the slot electrical filling factors obtained with a parallel bottom disposition and the one obtained with the parallel flank disposition of the standard slot, are reported in Figure 14. This comparison shows that the slot electrical filling factor is higher for a parallel flank disposition for each additional interspace *tw*. A similar study was performed for rectangular wires. Figure 15 shows the trend of slot electrical filling factors as a function of the cross-section area. On average, the filling factors obtained in the rectangular case are higher than those of the round case for the dimensions chosen. Additionally, in the rectangular case, the slot electrical filling factor is higher in a parallel flank disposition (Figure 16). In this work, the comparison between the slot electrical filling factors obtained in the round and rectangular wire cases, with equal cross-section, has not been reported because it has been widely discussed in [18] and in [19]. Furthermore, the rectangular wires are used for medium and high-power applications, whereas the round ones are mostly used for low-power applications. Figures 17 and 18 show the comparison between the value of electrical slot filling factors obtained for hexagonal wires with *phi* (φ) equal, respectively, to 0 and π/2. This study is performed both for the standard slot and for the non-standard

slot, for each value of the additional interspace *tw*. In the standard slot, the slot electrical filling factors are higher when *phi* = 0 whereas, in the non-standard slot, in some cases, higher electrical slot filling factors are obtained with *phi* = π/2. Particularly interesting is the comparison between the slot filling factors obtained in the case of the round wire and the hexagonal one, with the same cross-section. In Figure 19, the comparison between the slot electrical filling factors of round wires with the parallel flank disposition and the slot electrical filling factors of hexagonal wires with *phi* = 0, is reported. From this comparison, it can be seen that the filling factors obtained in the case of hexagonal wires are always higher than those obtained in the case of round wires. This difference decreases as the additional interspace increases. Therefore, the use of hexagonal wires can provide innovation in the field of electrical machine windings. Obviously, this improvement must be contextualized with any additional costs and the current state-of-the-art of specially shaped wire manufacturing process. Furthermore, with the choice of δ1, equal to 1/50 of the maximum dimension of the considered wire, the worst-case scenario of the study has presented a computation time equal to about 10 min, corresponding to about 8000 different first wire positions. In order to validate the choice of δ<sup>1</sup> equal to 1/50, further investigations have been carried out, both increasing and decreasing the value of δ1.

**Figure 13.** Slot electrical filling factor as a function of wire copper cross-section area for round wires.

**Figure 14.** Comparison between the slot electrical filling factors obtained with round wires distribution parallel to the ground/bottom (PG) and to the flank (PF) of the standard slot.

**Figure 15.** Slot electrical filling factor as a function of wire copper cross-section area for rectangular wires.

**Figure 16.** Comparison between the slot electrical filling factors obtained with rectangular wires distribution parallel to the ground/bottom (PG) and to the flank (PF) of the standard slot.

**Figure 17.** Slot electrical filling factor as a function of wire copper cross-section area for hexagonal wires in a standard slot (STSL).

**Figure 18.** Slot electrical filling factor as a function of wire copper cross-section for hexagonal wires in a non-standard slot (NSTSL).

**Figure 19.** Comparison between the slot electrical filling factors obtained with round wires distributed parallel to the flank of the slot and the one obtained with hexagonal wires as a function of wire copper cross-section area.

In more detail, values of δ<sup>1</sup> smaller than 1/50 only provide longer computational times with the identification of the same best wire arrangements. On the contrary, in almost all analyzed cases, the values of δ*<sup>1</sup>* in the range from 1/50 to 1/25 allow us to identify the same best wire arrangements with shorter computational times compared to those obtained with δ*<sup>1</sup>* being equal to 1/50. Values of δ<sup>1</sup> lower than 1/25 generate longer computational times due the higher number of δ*<sup>i</sup>* iterative variations. Therefore, the choice of δ*<sup>1</sup>* equal to 1/50 of the maximum dimension of the wire under study provides a good compromise between computation time and identification of the optimal wire arrangements.

The case studies shown here demonstrate how the algorithm allows us to accurately estimate the best value of the filling factors and reveal that it is a useful tool for the design of electrical machine windings. The results, like electrical filling factors, are calculated under the assumption of the ideal geometric shape of each wire and slot filling paper. For practical considerations, tolerances of shape and diameter must be taken into account; also, the winding process itself does not give an ideal orthocyclic winding distribution in the slot. This deviations from ideal geometry and winding process can be considered using additional geometric tolerances. With the given algorithm, the influence of non-ideal conditions on the filling factor and maximal turns can be calculated within a few minutes and motor design engineers and process engineers can make decisions based on a reliable calculation concept.

#### **5. Conclusions**

High filling factors make it possible to improve electrical machines' performances, meeting the design requirements of several application fields. This improvement requires an optimal arrangement of the wires inside the slots and a careful evaluation of the maximum slot filling factor obtainable in the design phase. In this paper, an algorithm approach is proposed to perform the calculation of slot filling factors in electrical machines. The algorithm requires, as input data, the geometrical data of the slot, the insulation thickness, the shape of the wire, the dimensions of the wire and the type of disposition inside the slot. From this data, the algorithm determines the maximum slot filling factors, the number of wires positioned, a graphic distribution of the wires inside the slot and the coordinates of the wires positioned. The algorithm has a high degree of flexibility and requires a reasonable computation time (about 10 min in the worst case). The conducted study proves that the algorithm is very simple and can give useful results in the designing processes of winding layouts. Furthermore, the algorithm can be used as investigation tool because it allows us to compare electrical filling factor values when different wire shapes with the same cross-section are employed. From the investigations here presented, it has been shown that the use of hexagonal wires provides higher filling factors than those obtained with the use of round wires. Although the comparison between the slot filling factors obtainable with rectangular wires and hexagonal wires with the same wire cross-section is not presented in this paper, the use of hexagonal wires allows us to obtain higher filling factors than those obtained with the use of rectangular wires. This result is due to the orthocyclic arrangement of the hexagonal wires that optimally occupy the slot area, unlike the ordered arrangement of the rectangular wires. Non-ideal geometry of the magnetic wires and the tolerances of the winding process can be easily considered using additional geometric factors and practical problems can be addressed. These results may be of considerable interest in the optimization of electrical machine windings. Future developments will concern the extension of this study to a greater number of cases in terms of slots and an in-depth study on the producibility of hexagonal wires. In particular, the evaluation of the effects of real cross-section deviation with respect to ideal deviation should be modelled and implemented in future research. Furthermore, experimental investigations will be conducted in order to validate the algorithm results.

**Author Contributions:** Authors contributed equally to this work. Authors of this manuscript jointly conceived the theoretical developments, revised the state-of-the-art algorithms and provided suggestions to obtain a flexible algorithm for the optimization of slot filling factors in the design phase of electrical machines. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** This work was financially supported by MIUR—Ministero dell'Istruzione dell'Università e della Ricerca (Italian Ministry of Education, University and Research) and by SDESLab—Sustainable Development and Energy Saving Laboratory of the University of Palermo.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Non-Salient Brushless Synchronous Generator Main Exciter Design for More Electric Aircraft**

#### **Filip Kutt \*,†, Michał Michna † and Grzegorz Kostro †**

Faculty of Electrical and Control Engineering, Gda ´nsk University of Technology, 80-233 Gda ´nsk, Poland; michal.michna@pg.edu.pl (M.M.); grzegorz.kostro@pg.edu.pl (G.K.)

**\*** Correspondence: filip.kutt@pg.edu.pl; Tel.: +48-58-347-19-39

† These authors contributed equally to this work.

Received: 19 April 2020; Accepted: 26 May 2020; Published: 27 May 2020

**Abstract:** This paper presents a prototype of high speed brushless synchronous generators (BSG) design for the application in autonomous electric power generation systems (e.g., airplane power grid). Commonly used salient pole field of the main generator part of BSG was replaced with a prototype non-salient pole field. The main objective of the research is an investigation into the advantages and disadvantages of a cylindrical field of the main generator part of BSG over the original salient pole field. The design process of the prototype generator is presented with a focus on the electromagnetic and mechanical finite element method (FEM) analysis. The measurements of prototype and commercial BSG were conducted for the nominal speed of 8 krpm. The advantages and disadvantages of the proposed solution were established based on measurements in load and no-load conditions.

**Keywords:** autonomous systems; brushless synchronous generator; electric power generation; high speed generator

#### **1. Introduction**

Modern commercial planes are designed according to the concept of more electric aircraft (MEA) [1–7]. This concept states that the electric power system in future aircraft should replace pneumatic and hydraulic systems supplied from main turbofan engines. The electrical power system should provide the ability to control the aircraft via electromechanical and electrohydraulic actuators and also provide deicing protection and control the pressurization of the cabin. The high increase in electric power demand is the main result of such an approach. High electric power demand, in turn, requires much more powerful generators to supply it. In addition, an important role of the generator is also the ability to work as a starter to accelerate a turbofan engine [7–10]. The increase in power output of the generator should not affect the volume and weight of the generators in a significant way. One of the obvious directions is the increase in the rotational speed of the machine [11–13]. However, one of the consequences of an increase in rotational speed is a higher centrifugal force acting on the rotor parts and the limitation in the rotor and shaft diameters.

To meet the requirements of the MEA concept considers various types of machines with high power densities defined as the power to weight ratio [7,11–14]. In addition to standard, wound field brushless synchronous generator (BSG), it is also proposed to use a permanent magnet (PM) [15,16], a switched reluctance (SR) [17–19] or an induction (IM) [8,20] machine as a starter-generator for aircraft. Despite the diversity of machine types proposed in research and prototype solutions, the BSG with a wound rotor is still the most often used in commercial aircraft. BSG has been used for decades, so its design is well known, tested and relatively reliable. Among its main advantages are the ability to easily control the voltage at the terminals by changing the excitation field, a simple adjustment system, the ability to easily reconfigure the system. Unfortunately, standard BSG solutions are characterized by

low efficiency (<80%), the maximum speed limit is less than 28 krpm and, consequently, also low power density (about 2.5 kW/kg) [11]. Therefore, research and design works are underway on new, alternative constructions for aircraft generators, including PM, SR, and IM generators [11–14]. These machines are characterized by a simpler rotor design, which allows operation at higher rotational velocity with less maintenance required. Permanent magnet generators achieve the highest power densities in relation to weight (from 3.3 to 8 kW/kg), and high efficiency (up to 95%) [10,11,15,16]. The disadvantage of PM generators is the complex regulation of the excitation field and voltage at the terminals and also the problem with excitation during a short circuit fault. The main advantage of SR generators is a very simple construction of the rotor, and therefore high reliability. SR generators have lower efficiency and lower power density than PM machines; however, it is still better than BSG [17–19]. On the other hand, SR requires complex control systems and algorithms. IM machines are also a prudent alternative, especially in the multiphase, fault-tolerance design [8,20]. Simple and cheap design, proven control systems and algorithms guarantee reliable work required in the aviation industry.

Only proven and reliable solutions can be used in aviation. The use of new solutions and devices for a commercial aircraft requires a time-consuming and expensive design, verification and certification process [21]. Therefore, research is still being carried out to improve the BSG design both in the context of achieving a better power-to-weight ratio and working at higher speeds and frequencies. These requirements can be met when materials, cooling, and mechanical structure would be improved [11].

High requirements regarding the parameters of the aircraft generator can be met by optimizing their design. Therefore, both new construction solutions, as well as new methods of their design, project management and production, are being developed. A wide range of issues related to the design and optimization of electrical machines is of interest to both manufactures and academic centers [14,21–25]. The requirements for high power density and high speed mean that it is important to optimize the generator system as a whole, and not its individual components. The design of the electric machine ceases to be an independent problem, and the optimal solution from the component point of view does not necessarily mean an optimal solution for the entire drive system. The design of the generator system is, therefore, a multi-level, multi-task, multi-disciplinary problem and at the same time nonlinear with many independent variables [26–29].

Computer modelling and simulation methods are widely used in the design process of electrical machines. Optimal use of computer modelling methods and simulations allows one to accelerate design work, verify the importance of potential solutions and, as a result, achieve a product with better parameters. Product development uses both simple analytical methods as well as more complex finite element method (FEM). Obtaining a low mass of the generator can be achieved by increasing the speed or reducing the mass of active materials. High-speed operation requires special attention when designing the rotor. The basis for design should be the selection of appropriate materials and FEM simulations to achieve the expected electromagnetic and mechanical parameters [30]. The reduction of the mass of conductive materials can be achieved by increasing the density of currents. Conducting computational fluid dynamics (CFD) simulation based on a complete conjugate heat transfer (CHT) allows the design of an appropriate cooling system [31].

Even if the final optimization of the machine requires the use of advanced simulation tools (FEM, CFD, etc.), the basis of the project is an analytical method. One of the methods widely used in the pre-design of electric machines is the sizing equation method [32–34].

The brushless synchronous generator is a complex power generation device (Figure 1). Three machines are installed on one shaft: the main generator, the exciter, and the subexciter. Both the main generator and the exciter are electromagnetically excited synchronous generators and the subexciter is a permanent magnet synchronous generator. The main generator field is supplied from the armature winding of the exciter through a 6-phase diode rectifier. The exciter armature is on the rotor along with the field of the main generator and the exciter field winding is on the stator. The exciter field winding is supplied from the subexciter via a controlled rectifier. The generator

control unit (GCU) controls the voltage and the protection systems for the generator. In the variable high-speed operation of the generator, the excitation system is designed to ensure the voltage control (RMS value) requirements of the power grid.

**Figure 1.** The brushless synchronous generator—based on a three-stage electrical machine topology: the subexciter—permanent magnet generator (PMG); the brushless exciter—synchronous machine with stationary field winding, rotating armature winding and rotating diode rectifier; the main generator—synchronous machine with rotating field winding; GCU—generator control unit.

In modern aircraft such as the Boeing 787 or the Airbus A350XWB the variable frequency generation system is used, where generators are directly driven from the turbofan engines. The system voltage is regulated at 230 V and the bus frequency varies from 350 to 800 Hz. Changing the way in MEA electric power generation requires adaptation of the distribution system and the use of more power electronic converters supplying loads [4]. This affects the power quality, which is also an important issue, as it has a direct effect on the reliability and efficiency of the power system [35]. The harmonic content of produced voltage and current of the main generator can also have a significant influence on the system power quality. In the salient pole machines, the effect of unsymmetrical pole shoe saturation during loaded conditions can increase the harmonic content in induced voltage [36].

#### **2. Objectives and Scope**

The main objective of the research is an investigation into the advantages and disadvantages of a cylindrical field of the main generator part of BSG over the original salient pole field. The main motivation for the research was the analysis of the possibility of introducing a cylindrical main filed construction without any detriment to the generated power quality in reference to the original construction. The proposed construction should allow for:


The development of a modern aircraft power system is based on MEA concept. However, only repetitively proven and reliable solutions can be used in commercial avionic applications. The MOET FP6 [37] project investigated the component level power system development. The practical result of the investigation was an evolutional change in the power generation system components and

particularly the BSG. The focus of those changes was the improvement of the generator structure for higher efficiency without any detriment to its relatability. That meant the introduction of a variable frequency system with a higher rotational velocity of the generator and improved power quality. The proposal in this paper solution is aiming to investigate the possible advantages in this regard of a cylindrical structure of the main exciter of the BSG. The literature on the subject of the construction of BSG main exciter is very limited and only investigates the possibility of using claw pole design [38] for the main exciter of the machine or the possibility of using a cylindrical structure for operation with higher rotational velocity [39].

The comparative study is conducted on a commercial and redesigned generator construction. Both machines are three stage machines consisting of a permanent magnet generator (PMG) subexciter, an inverted synchronous generator exciter and the main synchronous generator. The only difference between the two machines is the construction of the main generator field (Figure 2). The commercially available machine has a salient pole main generator field construction and the proposed prototype has a cylindrical (non-salient pole) construction of the main generator field. The proposed prototype was designed as a cylindrical rotor machine in order to obtain sinusoidal distributed excitation winding MMF with the possibility of achieving adequate mechanical strength when working at a higher rotational speed. The reverse engineering approach was used on the commercial salient pole generator to determine the electromagnetic parameters of the machine. The generator performance in the steady state was measured and analyzed. The main factor for the analysis was the minimization of the high-order harmonic content of currents and voltages waveforms during the no-load and load operation.

The design of the prototype generator is based on the salient pole brushless synchronous generator type GT40PCz8 (*Sn* = 40 kVA, *Un* = 208 V, *p*. *f* . = 0.8, *nn* = 8000 rpm) used in the Russian MI-28 helicopter and is a typical construction for avionic application. The commercial generator has a salient pole field of the main generator part and the proposed and constructed prototype uses a non-salient pole field (Figure 2).

**Figure 2.** Two types of exciter of the main generator: (**a**) salient pole (original commercial structure, exciter, and subexciter on one shaft); (**b**) non-salient pole (prototype structure, exciter, and subexciter temporally removed from the shaft).

The commercial GT40PCz8 generator was the subject of the reverse engineering process using analytical calculations and FEM simulations. This approach was used to determine the value of electromagnetic parameters such as current densities in the armature and field windings and magnetic flux density in various parts of the machine. Based on those values, the new main generator field was designed and developed.

The main contribution of the research is the design of the prototype of the main generator part cylindrical exciter that provides better power quality in no-load and under load conditions. This will allow for future developments for variable frequency power system higher maximum frequency.

#### **3. Solution**

The development and the implementation of a new design of the electrical machine requires a comprehensive approach to the physical phenomena taking part in the electromechanical energy conversion process in the drive system.

The first step of the design and analysis process is an electric machine design based on analytical calculations (Figure 3). The dimensions of a machine are calculated based on the sizing equations approach [22,34]. The sizing equation describes the relationship between the output power of the machine (*P*) and its main dimensions design features, material parameters, and rotational speed:

$$P = \frac{\pi}{2} K\_I K\_P K\_E \frac{f\_s}{p} (A\_s B\_m) (D\_s^2 l\_s) \tag{1}$$

where: *fs*—the stator voltage frequency, *p*—the number of pole pairs; the machine main dimensions: *Ds*—the inside stator diameter, *ls*—the length of stator core; the machine design features: *KI*—the current waveform factor, *KP*—the electrical power waveform factor, *KE*—the EMF (electromotive force or induced voltage) factor which incorporates the winding distribution factor and the air gap magnetic flux distribution, and the material parameters: *As*—the stator electric loading, and *Bm*—the air gap flux magnitude. The values of the design feature factors depend on the type of the machine, the type of the power supply (shape of the current waveform), the air gap flux distribution and the field and armature winding distributions [22].

**Figure 3.** Design process and analysis of the electrical machine—analytical calculation, numerical field model implemented in FEM software.

The main generator stage designed process was conducted based on the sizing Equation (1) and the parameters presented in Table 1. The parameters for the design process were derived from the reverse engineering of the original generator. The design process assumed that the new exciter should allow for the higher rotational speed of the BSG. In addition, the excitation winding should generate a sinusoidal distribution of magnetomotive force (MMF) to generate a sinusoidal distribution of the air gap flux density. Part of the solution for the proposed assumption was to use a high number of rotor slots. The outer diameter of the original and new field structure of the main generator is *Dr* = 173.8 mm, the airgap length is *δ* = 0.6 mm and the core length is *l* = 73 mm. The most important difference is that in a salient pole solution the air gap length varies, and, in a non-salient solution, it is constant. An additional advantage of the uniform air gap length is the decrease in non-uniform pole shoe saturation of the field during the under load operation of the BSG.


**Table 1.** Assumptions for the main generator stage design process.

The prototype field dimensions are partly based on the original generator construction and partly on the selected manufacturing technology, in particular, the technology of field winding construction and manufacturing. In the proposed manufacturing process, a square profile wire was selected for the field winding to maximize the fill factor. The minimum tooth thickness was calculated and then verified using FEM models for the machine to sustain the centrifugal forces and also due to the maximum value of the flux density in the field core. Figure 4 shows the dimensions of the main generator stator and both the commercial and the prototype generator rotors. The original salient pole excitation winding has 42 turns per pole and the design cylindrical prototype has 40 turns per pole. The cross section of the original and prototype excitation winding wire is the same and is equal to 4.5 mm2.

**Figure 4.** Dimensions of the: (**a**) main generator stator; (**b**) original salient pole rotor; (**c**) prototype non-salient pole rotor.

The proposed design prototype apart from allowing higher rotational velocity of the machine should also decrease the high-order harmonic content of the induced EMF. Both the air gap flux density spatial distribution and the field and armature winding spatial distribution influence the harmonic

content of the induced EMF. To analyze this influence, the air gap length distribution functions and field and armature winding distribution functions have to be defined. The air gap length distribution function is defined as:

$$\delta\left(\phi\_{\delta}-\theta\_{r}\right) = \frac{1}{\mathfrak{a}\_{\delta}} \frac{1}{\mathfrak{a}\_{0} - \sum\_{k=1}^{k\_{\delta}} \mathfrak{a}\_{2k-1} \cos\left(2\left(2k-1\right)\left(\phi\_{\delta}-\theta\_{r}\right)\right)}\tag{2}$$

where *αδ* defines the average length of the air gap, *α*<sup>0</sup> and *α*2*k*−<sup>1</sup> are the relative components of the air gap length distribution function (*α*<sup>0</sup> = 1), *φ<sup>s</sup>* is the angular position along the stator and *θ<sup>r</sup>* is the displacement angle of the rotor axis in reference to the stator axis. In case of the prototype machine with cylindrical stator and rotor, this function is defined as:

$$
\delta = \frac{1}{a\_{\delta}} \tag{3}
$$

The armature winding distribution is defined based on the winding MMF distribution as:

$$N\_{\rm xs} \left(\phi\_{\rm s}\right) = \frac{1}{i\_{\rm xs}} \frac{\partial \, MMF\_{\rm xs} \left(i\_{\rm xs}, \phi\_{\rm s}\right)}{\partial \, \phi\_{\rm s}} \tag{4}$$

where *x* refers to phase *a*, *b* or *c* in three phase machine, *ixs* is the *x*'th phase armature current and the *MMFxs* (*ixs*, *φs*) is the *x*'th phase magnetomotive force distribution. This force is defined as:

$$MMF\_{\rm x5} \left( i\_{\rm x5}, \phi\_{\rm 5} \right) = \frac{N\_{\rm s}}{2} i\_{\rm x5} \sum\_{k=1}^{k\_{\rm MMF\_k}} A\_{\rm s,2k-1} \cos \left( \left( 2k - 1 \right) \left( \phi\_{\rm 5} + \theta\_{\rm x5} \right) \right) \tag{5}$$

where *Ns* represents the number of turns of equivalent sinusoidally distributed armature winding, *As*,2*k*−<sup>1</sup> are the relative amplitudes of armature winding MMF distribution (*As*,1 = 1), *<sup>θ</sup>xs* is the angular displacement of *x*'th phase in reference to phase *a* (*θas* = 0) and *kMMFs* is the number of odd harmonics used to approximate the armature MMF distribution. The same winding distribution and magnetomotive force distribution functions can be defined for field winding:

$$N\_{fd} \left(\phi\_r\right) = \frac{1}{i\_{fd}} \frac{\partial \, MMF\_{fd}\left(i\_{fd}\phi\_r\right)}{\partial \, \phi\_r} \tag{6}$$

$$\text{MMF}\_{fd}\left(\mathbf{i}\_{fd}, \boldsymbol{\phi}\_{\mathcal{V}}\right) = \frac{N\_{fd}}{2}\mathbf{i}\_{fd} \sum\_{k=1}^{k\_{MMF}} A\_{fd,2k-1} \sin\left(\left(2k-1\right)\left(\boldsymbol{\phi}\_{\mathcal{V}}\right)\right) \tag{7}$$

where *if d* is the field current, *φ<sup>r</sup>* is the angular position along the rotor, *Nf d* represents the number of turns of the equivalent sinusoidally distributed field winding, *Af d*,2*k*−<sup>1</sup> are the relative amplitudes of field winding MMF distribution (*Af d*,1 = 1) and *kMMFf d* is the number of odd harmonics used to approximate the field MMF distribution. Based on the winding distribution the MMF distribution and the air gap length distribution, one can calculate the mutual inductance between field winding and phase winding:

$$\mathcal{L}\_{\text{xsfd}}(\theta\_r) = \frac{1}{i\_{fd}} \int\_{\pi}^{2\pi} \left( N\_{\text{xs}} \left( \phi\_s \right) \int\_{\phi\_s}^{\phi\_s + \pi} \frac{MMF\_{\text{xs}} \left( i\_{\text{xs}}, \phi\_s \right)}{\delta(\mathcal{J} - \theta\_r)} d\mathcal{J}\_s \right) d\phi\_s \tag{8}$$

If Equations (3), (4) and (7) are applied to relation (8), the function for mutual inductance is defined as:

$$L\_{xsfd}(\theta\_{\boldsymbol{\tau}}) = L\_{sfd} \sum\_{k=1}^{k\_{\rm \perp\_{sfd}}} \left( (2k-1) \, a\_0 A\_{s,2k-1} A\_{fd,2k-1} \sin \left( (2k-1) \left( \theta\_{\boldsymbol{\tau}} + \theta\_{\boldsymbol{\text{xs}}} \right) \right) \right) \tag{9}$$

where:

$$L\_{sfd} = \frac{N\_s N\_{fd}}{4} \mu\_0 l \tau r \alpha\_d \tag{10}$$

where *μ*<sup>0</sup> is the magnetic permeability of vacuum, *l* is the machine length and *r* is the distance from the axis of machine to the middle of the air gap. As can be noticed, the mutual inductance harmonic components depend on the air gap length and the armature and field winding distribution high-order harmonics. In the designed prototype, the armature winding was not changed. However, both the air gap length distribution and the field winding distribution were changed. The proposed solution for more sinusoidally distributed excitation winding MMF was to use a high number of rotor slots. Because the prototype has a cylindrical construction and more sinusoidally distributed winding, then the commercial generator the resulting EMF is much more sinusoidal. Figure 5 shows the MMF distribution for both the commercial and the prototype main generator field windings.

**Figure 5.** MMF distribution for the commercial and the prototype machine main generator field windings.

Simplified mechanical calculations have been performed to approximate forces and pressures acting on the prototype main generator field core. The force acting on the weakest part of the rotor was calculated as a sum of two forces:


The centrifugal force acting on the rotor tooth was calculated as:

$$F\_{\text{toot}h} = \omega^2 \int\_0^\varepsilon \rho\_{\text{Fe}} \cdot l \left( a + \left( \mathbf{x} \frac{b - a}{\varepsilon} \right) \right) (r\_1 + \mathbf{x}) \, d\mathbf{x} \tag{11}$$

where *ω* is the mechanical rotational velocity of the rotor, *ρFe* is the density of rotor core, *l* is the length of the main generator part of the machine, *a* is tooth thickness at the bottom of the tooth, *b* is the tooth thickness at the top of the tooth, *c* is the height of the tooth and *r*<sup>1</sup> is the bottom radius of the rotor tooth. The force acting on the copper in the slot was calculated as:

$$F\_{\mathbb{C}u} = \omega^2 \int\_0^\varepsilon \rho\_{\mathbb{C}u} \cdot l \cdot d\left(r\_1 + \mathbf{x}\right) d\mathbf{x} \tag{12}$$

where *d* is the width of the copper profile wire in the rotor slot. The resulting force acting on the rotor tooth from both the tooth itself and the copper in the slot is about 23 kN. The resulting pressure in the thinnest part of the rotor tooth is about 126 MPa.

Because the resulting machine dimensions and the resistance of excitation winding did not change for the prototype and original commercial generator, the thermal analysis of the prototype was not conducted. However, if this solution was to be implemented in a commercial power system, such analysis would have to be conducted.

Three-dimensional (3D) geometric models of the machine are then developed as the parametric model using CAD software. This approach enables easy modification of the geometric models and studying of the influence of the geometrical parameters and material properties on the magnetic field distribution, power losses, integral parameters (inductances, torque), mechanical stresses, etc. (Figure 3).

#### **4. Results**

#### *4.1. FEM Simulation*

The FEM simulations (using Cedrat/FLUX2D software (version 9.3, Cedrat, Grenoble, France)) have been carried out to verify the designed BSG. Due to a lack of information about the materials used in the commercial generator, the FEM simulations were verified using measurement data. As a result of reverse engineering processes for the given dimensions and measured behavior of the commercial machine, the material parameters of the stator and rotor magnetic core were established. Those parameters were calculated using the iteration process by minimizing the error between measurements of no-load state and FEM simulation results. The initial relative permeability *μ<sup>r</sup>* = 8000 and the magnetic polarization at saturation of 1.6 T parameters were set for the commercial machine stator and rotor core. For the prototype field core of the main generator, M530-50A steel sheets have been used. Unfortunately, the manufacturing process and relatively poor stacking factor caused the necessity for additional varnish coating insulation between steel sheets. During the laser cutting process of the prototype core steel sheets, a slag caused by laser cutter needed to be removed. This process also damaged the original varnish insulation coating on the surface of the steel. Because of this decrease in the stacking factor and the damage done to the cut edge of the electrical steel, the material magnetic characteristic has changed. The field core parameters for the simulation were set to *μ<sup>r</sup>* = 500 (initial relative permeability) and the magnetic polarization at saturation of 1.55 T.

The model was implemented in the Cedrat/FLUX2D software using a dedicated Python script. For the discretization of the FEM model, the air gap was divided into three sections: the stator static mesh section, the rotor moving mesh region and the air gap re-meshing region for the mesh generation during the rotor movement. The resulting mesh and geometry of both the commercial and prototype generator are shown in Figure 6.

**Figure 6.** FEM model mesh of the: (**a**) commercial salient pole field BSG; (**b**) prototype cylindrical field BSG.

The resulting mesh generation process led to a total number of about 36 thousand nodes for the original salient pole main generator of BSG and about 53 thousand for the prototype. The simulation process was conducted with set maximum variation between the iterations of integral calculation of less than 0.1%. Table 2 shows the maximum values of magnetic field densities in the crucial parts of the commercial and designed machines.

**Table 2.** Maximum values of flux density of the main generator part of BSG.


FEM simulation results at no-load conditions for the commercial salient-pole and prototype cylindrical generators are shown on Figure 7.

**Figure 7.** FEM simulation results in no-load conditions: (**a**) commercial salient pole field BSG; (**b**) prototype cylindrical field BSG (marked direction of current in field windings).

The magnetic equipotential lines are shown for both the commercial and the prototype BSG main generator part. In this simulation, the main generator field winding has been supplied with nominal no-load current and FEM computations in magnetostatic conditions have been conducted. The amplitude of the normal component of the air gap flux density of both types of the main generator field has the same maximum value of 0.8 T (Figure 8). The prototype cylindrical generator has got more sinusoidal distribution of air gap flux density (Figure 8).

**Figure 8.** FEM simulation results at no-load conditions—air gap flux density normal component distribution (0.8T—amplitude of the fundamental component): (**a**) commercial salient pole and prototype air gap flux distribution BSG; (**b**) harmonic components of air gap flux distribution.

In addition to electromagnetic calculations, the mechanical calculations concerning the centrifugal forces have been conducted. The simulation of this force acting on one of the rotor teeth was analyzed. The stress analysis module of Autodesk Inventor software (version 2014, Autodesk, San Rafael, CA, USA) was used to perform the mechanical FEM simulation. The rotation speed of 16 krpm has been applied to the rotor. The materials used for stator (steel sheets), windings (copper), and slot wedges (bras) were defined. The stress analysis results are shown in Figure 9. In the simulation, only the copper inside of the slots is considered as the copper in the end windings will not be held by the rotor core itself.

**Figure 9.** FEM simulation stress analysis results of centrifugal force acting on the rotor.

The electromagnetic and mechanical FEM simulations verified the assumptions and results of the analytical calculations. The proposed solution delivers more sinusoidal distribution of the excitation winding MMF and the maximum tensile stress due to centrifugal force does not exceed maximum tensile stress allowed for M530-50A steel.

#### *4.2. Measurements*

Measurements have been carried out for the no-load (open circuit) and load conditions of the prototype and commercial BSG. Measurements of the no-load back EMF in the function of field current for constant speed are shown in Figure 10. This shows that the prototype's main generator field requires more than two times the no-load nominal excitation current. The main factor in this is the smaterial and manufacturing technology used for the core of the main generator field. The difference in number of turns between the prototype and the original machine is only 5% and the air gap length is the same. This also means that the armature reaction of the prototype generator is diminished.

**Figure 10.** No-load EMF measurements (line to line value) in the function of the exciter field current for commercial generator (back\_EMF\_comm) and prototype generator (back\_EMF\_prot).

The no-load back EMF waveform of the prototype has less high-order harmonic components than the commercial one (Figure 11). In both machines, the stator winding is the same which means that the back EMF waveform is dependent on the air gap flux density spatial distribution which in case of a prototype generator is more sinusoidal.

For the load conditions measurements, a load impedance consisting of 0.8 Ω active part and 0.15 Ω reactive inductive part (@400 Hz) has been used. The measurement results of the voltage and current harmonic components in load conditions are shown in Figure 12. As can be observed, the designed prototype with non-salient pole field of the main generator produces less high-order harmonics in the back EMF waveform.

**Figure 11.** No-load voltage for prototype and commercial BSG: (**a**) back EMF waveform, (**b**) harmonic content of back EMF (vas\_prot—measured for prototype BSG, vas\_comm—measured for commercial BSG).

**Figure 12.** Measurement of armature voltage and current of the salient and non-salient BSG in load conditions (balanced star connected load *R* = 0.8 Ω, *X* = 0.15 Ω (@400 Hz), *p*. *f* . = 0.98), (**a**) harmonic content of voltage waveforms, (**b**) harmonic content of current waveforms (vas\_comm, ias\_comm—results for salient pole generator, vas\_prot, ias\_prot—results for non-salient pole generator).

In addition to no-load and under load steady state test, a transient state of the symmetric short circuit was performed (Figure 13). This test was performed by shorting the armature terminals when the machine was operating at no-load condition with the nominal voltage at the terminals. As can be observed, the excitation current is higher in the prototype generator than in the commercial one and equal in steady state to the nominal no-load excitation current.

**Figure 13.** Measurement of armature and exciter currents during transient short circuit test: (**a**) armature currents; (**b**) exciter currents (ias\_comm, ifd\_ex\_comm—results for salient pole generator, ias\_prot, ifd\_ex\_prot—results for non-salient pole generator).

#### **5. Conclusions**

The comparative study of two types of BSG main generator exciter has been conducted. The salient pole and cylindrical exciter construction have been compered. The exciter and subexciter of BSG were not changed. The design process was composed of analytical calculation, electromagnetic FEM calculations, and mechanical FEM calculations. The mechanical calculations were conducted to analyze the centrifugal force impact on the designed field of the main generator at higher than nominal (16 krpm while nominal is 8 krpm) rotational velocity, whereas the electromagnetic FEM simulations have been conducted in order to verify the analytically calculated dimensions and parameters of the main generator exciter.

The BSG with the prototype non-salient field has a higher volume and mass of the magnetic core (cylindrical core—9.3 kg, salient pole core—8.1 kg). The volume and mass of the copper used for the field coils are slightly smaller in design and build a prototype (cylindrical field windings copper—2.7 kg, salient pole field windings copper—2.9 kg).

The measurements were conducted at a nominal rotational speed of 8 krpm due to the same exciter and subexciter in both generators. The results of the measurements show significant advantages of the designed solution. Because of more sinusoidal MMF distribution of the field winding, the no-load phase voltage contains less high-order harmonics than the voltage of the silent pole commercial BSG (GT40PCz8). The possibility of operation at higher rotational speeds and better quality of produced electrical energy are a good indicator for the development of the BSG for the MEA variable frequency power system.

The designed generator also has its flaws. It is slightly heavier than the original, and the cooling of the field winding will be different (lack of axial vents in the designed rotor). The selected technique of the manufacturing leaves much to be desired and will have to be modified to archive better magnetic permeability of the field core.

**Author Contributions:** Conceptualization, F.K., M.M., and G.K.; methodology, F.K. and M.M.; software, F.K. and G.K.; validation, F.K., M.M., and G.K.; formal analysis, M.M.; investigation, G.K.; resources, F.K. and M.M.; writing—original draft preparation, F.K.; writing—review and editing, M.M. and G.K.; visualization, F.K.; supervision, M.M. and G.K.; project administration, F.K.; funding acquisition, F.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Polish Ministry of Science and Higher Education Grant No. N N510 328937.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**


c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Smart-Sensors to Estimate Insulation Health in Induction Motors via Analysis of Stray Flux**

#### **Israel Zamudio-Ramirez 1, Roque Alfredo Osornio-Rios 1, Miguel Trejo-Hernandez 1, Rene de Jesus Romero-Troncoso <sup>1</sup> and Jose Alfonso Antonino-Daviu 2,\***


Received: 6 April 2019; Accepted: 25 April 2019; Published: 1 May 2019

**Abstract:** Induction motors (IMs) are essential components in industrial applications. These motors have to perform numerous tasks under a wide variety of conditions, which affects performance and reliability and gradually brings faults and efficiency losses over time. Nowadays, the industrial sector demands the necessary integration of smart-sensors to effectively diagnose faults in these kinds of motors before faults can occur. One of the most frequent causes of failure in IMs is the degradation of turn insulation in windings. If this anomaly is present, an electric motor can keep working with apparent normality, but factors such as the efficiency of energy consumption and mechanical reliability may be reduced considerably. Furthermore, if not detected at an early stage, this degradation could lead to the breakdown of the insulation system, which could in turn cause catastrophic and irreversible failure to the electrical machine. This paper proposes a novel methodology and its application in a smart-sensor to detect and estimate the healthiness of the winding insulation in IMs. This methodology relies on the analysis of the external magnetic field captured by a coil sensor by applying suitable time-frequency decomposition (TFD) tools. The discrete wavelet transform (DWT) is used to decompose the signal into different approximation and detail coefficients as a pre-processing stage to isolate the studied fault. Then, due to the importance of diagnosing stator winding insulation faults during motor operation at an early stage, this proposal introduces an indicator based on wavelet entropy (WE), a single parameter capable of performing an efficient diagnosis. A smart-sensor is able to estimate winding insulation degradation in IMs using two inexpensive, reliable, and noninvasive primary sensors: a coil sensor and an E-type thermocouple sensor. The utility of these sensors is demonstrated through the results obtained from analyzing six similar IMs with differently induced severity faults.

**Keywords:** induction motor; smart-sensor; stray flux; time-frequency transforms; wavelet entropy

#### **1. Introduction**

In the companies, electric motors have gained great importance, and have been widely used as electromechanical devices for the conversion of energy, consuming more than 60% of all the energy of any industrial nation [1]. Current quality requirements consider the use of monitoring systems and the development of incipient failure detection techniques increasingly necessary in order to enhance the reliability of these industrial systems so that production is not interrupted. Machines operating under faulty conditions consume and spend more energy, causing additional economical losses. Furthermore, some failures can remain unnoticed in motors that work continuously with apparent normality; nevertheless, if not detected in time, incipient faults can result in catastrophic and irreversible damage to the machine, and if the fault progresses it can cause collateral damages to others systems coupled to the induction motor (IM). Therefore, it is of paramount importance to study the main faults in induction motors, and there is a clear necessity to develop emergent techniques that can detect faults in the early stages, and to integrate new technologies. In this regard, some authors have adopted the concept of a smart-sensor, in which one or more primary sensors are combined with a processing unit in order to gather certain functionalities like processing, communication and integration. Smart-sensors have found an application in different research fields, including the monitoring and diagnosis of faults in distinct industrial applications [2,3], real-time high-resolution frequency measurement [4], identification of broken bars and unbalance in induction motors [5,6], among others.

Surveys on motor reliability have determined that the distribution of failures in IMs can essentially be classified into four classes: bearing faults, stator related faults, rotor related faults, and other faults (cooling, connection, terminal boxes) [7]. Some investigations have shown that most failures of electric motors can be attributed to bearings and windings [8]. Depending on the type and size of the machine, problems related to stator windings correspond to a range between 16%–36% of total reported faults [9,10], which is the second largest type of fault for IMs, just after bearings defects. Stator winding insulation failures have recently received special attention. This is mainly because the worst stator faults start from undetectable insulation degradation problems between drastically adjacent turns [11–14] that lead to the appearance of an inter-turn fault, where two or more turns become short circuited. If undetected at early stages or after its appearance, this type of fault can develop into more severe problems very quickly. Many techniques found in the literature have been proposed to detect winding faults, and focus on two main approaches: offline methods and online methods [15]. Common offline methods that are typically used in industry include insulation resistance measurement, polarization index/dielectric absorption measurement, offline partial discharge tests, and evaluation of the dissipation factor [16–18]. A disadvantage of offline tests is the necessity to remove the machine from service, a drawback that can lead to false indications caused by unrealistic operations [19]. On the other hand, online monitoring methods are desirable due to their capability to diagnose faults when a motor is in service. To this end, several techniques have been proposed to perform online diagnosis, and many physical magnitudes have been highlighted as potential sources of information, with each one having its own advantages and disadvantages, as discussed below. Vibration [20], thermographic [21,22], and partial discharge [23] are some of the online methods used to detect insulation inter-turn faults; however, most of these techniques are not yet proven to detect faults during early stages, before reaching a severe phase, and in the case of the thermographic approach, diagnosis is difficult to perform under real working environmental conditions, since optimal conditions must be met to get confident results. Other classical approaches are focused on the use of current and voltage signals: spectral analysis of the steady-state current using the Fourier transform (motor current signature analysis, MCSA) [24], analysis of the zero and negative sequence currents [25], and analysis of the zero-sequence voltage [26]. The main disadvantage of analyzing the zero-sequence voltage is that the final diagnosis can be affected by the influence of other parameters, such as voltage unbalances, measurement errors, and inherent asymmetries during the manufacturing process, which can cause false diagnoses. Although analysis of the negative sequence current overcomes these problems [27], it is required to measure three-line currents, a condition that is not always available.

Due to the need for a system that is able to automatically diagnose in an online mode and monitor the health of winding insulation in induction motors (before an irreversible fault occurs), this work introduces a smart-sensor composed of two primary sensors (a coil sensor and an E-type thermocouple sensor) and a hardware signal processing unit (HSP unit) in order to accomplish this task. The wavelet entropy (WE) of the coil sensor was used as an auxiliary parameter in the final diagnosis, since it is able to characterize the dynamism and the order/disorder of a signal using a single value [28–31]. The coil sensor was used to capture the stray flux signal, and the E-type thermocouple sensor acquired the temperature of the motor, since both of them have a non-invasive nature. In order to constantly monitor and diagnose winding insulation degradation, the smart-sensor applied a signal processing stage composed of the computation of the discrete wavelet transform (DWT) followed by the calculation of the WE. Furthermore, with the purpose of automating the full process, a trained artificial neural network (ANN) performed a regression estimation by using the wavelet entropy and the induction motor temperature as input signals. All computations were performed by a field-programmable gate array (FPGA) HSP digital unit by developing proprietary hardware cores focused on the above-mentioned tools, as described below.

#### **2. Materials and Methods**

In this section, we detail the mathematical tools and methodologies that constitute the main core of the smart-sensor. The DWT is used to obtain a representation of the frequency content for the different bands that make up the input electromotive force (*emf*) signal. Wavelet entropy is used as the main parameter that will serve to subtract relevant information about the healthiness of the winding insulation, since it is a tool capable of describing the dynamic behavior of a signal, in addition to indicating the amount of order/disorder of that signal. Furthermore, wavelet entropy shows a clear relation to the healthiness of the winding insulation, as will be shown below. Finally, a final diagnosis through an artificial neural network, whose inputs are indeed wavelet entropy and IM temperature, will indicate the healthiness of the winding insulation using an automatic process.

#### *2.1. Discrete Wavelet Transform (DWT)*

As is well known, DWT is a time-frequency analysis transform that provides significant features for the analysis of a time-variant signal, since this technique is very suitable for decomposing a signal into well-defined "wavelet signals" that cover specific frequency ranges that are known to be directly dependent on the sampling frequency used to capture the analyzed signals [32]. The DWT of a signal can be defined as

$$\mathcal{W}(i,k) = \sum \mathbf{x}(k)\boldsymbol{\psi}\_{i,k}(t) \tag{1}$$

where *i* is the decomposition level, *k* is the number of the sample, and ψ*i*,*k*(*t*) is the discrete wavelet mother function.

To compute the DWT of a signal, a Mallat's algorithm facilitates its application and improves its performance, processing time, and the computational burden that its application entails. The DWT of a signal *x*[*n*] of length *N* is calculated by applying a mathematical convolution defined by Equation (2) with a bank of high-pass filters (HPF) with impulse response *g[n]* to analyze the high frequencies, and simultaneously with a bank of low-pass filters (LPF) with impulse response *h[n]* to analyze the low frequencies.

$$y[n] = (\mathbf{x} \* h)[n] = \sum\_{k=0}^{N} \mathbf{x}[k]h[n-k] \tag{2}$$

The DWT decomposes the time-domain signal in several levels, which are limited by the sample size *N.* The frequency content of every decomposition level for both *aCi* and *dCi* is estimated by

$$a\mathbb{C}\_i \to \left[0, \frac{f\_s}{2^{i+1}}\right], d\mathbb{C}\_i \to \left[\frac{f\_s}{2^{i+1}}, \frac{f\_s}{2^1}\right] \tag{3}$$

where *fs* is the sampling frequency and *i* is the desired decomposition level.

The coefficients of the HPF and LPF are determined by the selection of a mother wavelet according to the application.

In this regard, some investigations have been developed to evaluate the performance of DWT in extracting features from the current signals of induction motors. This serves the purpose of detecting eccentricities [33], rotor-asymmetries [34], broken rotor bars [35], and other factors using the Daubechies (db), Symlet (sym), Morlet, and Meyer wavelet families, and varying the orders. These works have shown that a Daubechies family of higher order is well suited to extracting the information required for the detection of motor failures. Furthermore, studies have shown that higher order filters behave as more-ideal filters, allowing less overlap between adjacent frequency bands.

#### *2.2. Wavelet Entropy*

Due to the inherent constraints of some time-frequency transforms, there can be problems when a specific window is applied to a series of data. Such is the case of the uncertainty problem given in the DWT—if the window is too narrow, the resolution of the frequency will be poor, whereas if the window is too wide, the location during the time of the signal will be less precise. This limitation is of great importance when it comes to the analysis of signals with transient components located in time, which is the case for the great majority of signals with real physical magnitudes.

To minimize the effects of this limitation, a parameter based on the entropy of a signal has been defined from a time-frequency representation of the signal provided by the wavelet transform [36]. In this regard, the entropy based on the wavelet transform (wavelet entropy) reflects the degree of order/disorder in the signal, so it can provide additional information about the underlying dynamic processes associated with the signal [29]. This is achieved by combining the information of all the wavelet bands, since data between adjacent wavelet signals is taken and combined into one index in order to avoid focusing on just one wavelet band having its own time-frequency resolution.

The total wavelet entropy (*SWT*) is defined according to [37].

$$S\_{WT} \equiv S\_{WT}(p) = -\sum\_{i=m}^{n} p\_i \ln p\_i \tag{4}$$

where *m* and *n* are the first- and last-considered decomposition levels for analysis, respectively, and *pi* represent the relative wavelet energy normalized values, which can be computed as

$$p\_i = \frac{E\_i}{E\_{tot}}\tag{5}$$

where *Ei* (Equation (6)) and *Etot* (Equation (7)) represent the energy of wavelet level decomposition *i* and the total energy of all wavelet level decompositions, respectively.

$$E\_i = \sum\_k \left| \mathbb{C}\_i(k) \right|^2 \tag{6}$$

$$E\_{\text{tot}} = \sum\_{i} E\_i \tag{7}$$

#### *2.3. Stray Flux Analysis*

Effective analysis of the stray flux by applying suitable signal processing techniques to detect several failures in induction motors, such as broken rotor bars, static and dynamic rotor eccentricity, bearing faults, and shorted turns in stator winding, have been proven and validated in a number of works [38–41].

The external magnetic field can be analyzed by its axial and radial components [42]. The axial radial field is generated by currents in the stator end windings or rotor cage end ring. The radial field is related to air gap flux density, which is attenuated by the stator magnetic circuit and by the external machine frame.

The coil sensor can be installed in the vicinity of the motor frame in convenient positions in order to measure the electromotive forces that are indicative of the axial and radial flux components, depending on its placement. Thus, Figure 1 shows the positions A, B, and C in which the sensor can be installed to measure both fields.

**Figure 1.** Coil sensor positions in the tested induction motor.

In position A, the sensor's placement enables measurement of the axial flux; on the other hand, if the sensor is placed in position B, the result of the measurements will correspond to components of the axial flux and part of the radial flux simultaneously. In position C, the coil sensor is parallel to the longitudinal cross-section of the machine, which makes the axial field null and the radial flux essentially present.

#### *2.4. Artificial Neural Network*

Artificial neural networks (ANNs) are computational models that simulate the neurological structure of the human brain and its capability to learn and solve problems through pattern recognition [43]. As is well known, this method has exceptional characteristics enabling it to process and extract relevant information from large amounts of data. Among the most popular ANN architectures, feed-forward neural networks (FFNNs) are widely used, since they are simple, practical, and very good at approximating real-valued functions and at classifying data. Furthermore, the operation of this kind of ANN demands a very low computational burden, which makes it appropriate for implementation in digital systems. FFNNs are composed of a layered architecture possessing essentially one input layer, one or more hidden layers, and one output layer, as shown in Figure 2a. Each layer has one or more elementary processing units called neurons (see Figure 2b), whose processing capability is stored in the connections of synaptic weights, and whose adaptation depends on learning [44]. The mathematical model describing the functionality of each neuron is given by

$$y = f\left(\sum\_{i=1}^{n} w\_i \mathbf{x}\_i + b\right) \tag{8}$$

where *y*, *wi*, *xi*, *b*, *f*(·), and *n* are the output, synaptic weights, inputs, bias, activation function, and the total number of inputs, respectively. To define the network weights, a training process is carried out where pairs of input–output data are presented, then a training rule is defined for adjusting these weights.

**Figure 2.** Artificial neural network (ANN): (**a**) feed-forward neural network (FFNN) architecture, (**b**) functional structure of a neuron.

#### *2.5. Smart-Sensor*

The smart-sensor proposed here is based on a low-cost system on a chip (SoC) field-programmable gate array (FPGA) to estimate the insulation health of an induction motor. Figure 3 shows the scheme of the structure of the proposed smart-sensor. The system uses a coil sensor and a thermocouple sensor as primary sensors that can be installed on the frame of the analyzed IM to capture the stray flux and temperature of the IM, respectively. The information coming from the primary sensors is acquired in the data acquisition system (DAS) module, then the signal processing is performed in the FPGA-based HSP unit by applying suitable time-frequency decomposition (TFD) tools and by extracting an efficient indicator based on the wavelet entropy. Finally, the estimated health of the insulation is supplied to the final user using an liquid crystal display (LCD).

**Figure 3.** Block diagram of the proposed smart-sensor: primary sensor, data acquisition system (DAS), field-programmable gate array (FPGA) and liquid crystal display (LCD).

#### 2.5.1. Primary Sensors

Two sensors are used as primary sensors—one flux sensor, and one thermocouple sensor. The flux sensor is generated by 1000 turns of a coil. Its dimensions are specified in Figure 4b. The main purpose of the coil is to detect the largest amount of stray flux possible through the induced electromotive force (*emf*) in that coil. The coil is protected with a material that is able to isolate the greater amount of electromagnetic noise coming from the outside using a special meshed cable for the transmission of the induced voltage towards the DAS. On the other hand, to capture the temperature of the analyzed IM, an E-type thermocouple sensor (Figure 4a) is used, since it is a non-magnetic sensor and has a wide temperature ranging from −50 ◦C to 740 ◦C.

**Figure 4.** Primary sensors: (**a**) thermocouple sensor, (**b**) coil sensor dimensions.

#### 2.5.2. DAS and LCD

As secondary elements, the DAS and LCD modules enable the interaction between the final user and the primary sensors. On the one hand, the DAS module is constituted by a signal-conditioning submodule and an ADS7841 analog-to-digital converter. The signal-conditioning module is composed of one operational amplifier with two processing stages, the first of which is configured to sum a constant voltage to the input flux signal, while the second amplifies it by a factor of 10 in order to standardize the input-voltage range to the analog-to-digital converter used. On the other hand, the LCD is used to display the estimated winding insulation degradation to the final user.

#### *2.6. HSP Unit*

The FPGA-based HSP unit is created by processing the DWT, feature-extraction wavelet entropy, and regression FFNN, and by mapping the min–max function to normalize inputs and outputs of the FFNN, as shown in Figure 5. First, the input *emf* signal (φ) is decomposed by the DWT in multiple "wavelet signals" in order to obtain the time-frequency representation of the input signal in well-known frequency bands. Then, the feature extraction is performed by applying Equation (3). Note that the *SWT* value of a signal is a normalized parameter ranging from zero to 1, where a minimum value indicates a light disorder in the analyzed signal (that is, the signal is mainly represented by one wavelet signal having the highest amplitude). On the other hand, if the value is near 1, the analyzed signal is considered to have a high disorder, since it penetrates several wavelet signals, each one having high relative amplitudes. Next, the min–max function map normalizes the IM temperature signal (*T*) and the *SWT*, in order to perform the mathematical operations inside the FFNN in a defined closed range. Finally, the FFNN unit performs the regression diagnosis by using the normalized values of the extracted *SWT* and the temperature of the induction motor as inputs.

**Figure 5.** FPGA-based Hardware Signal Processing (HSP) unit: Analogical-Digital Converter (ADC), Discrete Wavelet Transform (DWT), Inverse Discrete Wavelet Transform (IDWT), Wavelet Entropy and Healthiness estimation module.

Figure 6 shows the *emf* and induction motor temperature signals processing the flow to the FFNN evaluation for the purpose of obtaining a health estimation of the winding insulation. The processing flow starts with the *emf* signal acquisition; then, it is computed by the DWT to a level defined by the final user. The next step is to obtain the *SWT* parameter from the detail decomposition signals after performing the DWT by applying Equation (4). Finally, wavelet entropy and induction machine temperature are used as inputs to the FFNN, so that the information of both parameters can be combined to offer an automated estimation of the health of the winding insulation. Note that the FFNN is composed of one input layer with two input neurons (the wavelet entropy and the induction machine temperature); two hidden layers with four and two neurons, respectively; and one output neuron (the estimated health of the winding insulation, a parameter shown to the final user via the LCD). To specify the health of the winding insulation, results are shown in a continuous scale ranging from 10% (indicating a severe degradation) up to 95% (implying a healthy winding insulation system).

**Figure 6.** Proposed methodology flow.

#### 2.6.1. DWT–IDWT Digital Structure

In Figure 7, the counters *n* and *rwdir* indicate the index of the sample *x*[*n*] to read and the read/write direction of the approximation and detail coefficients, respectively. Note that it is necessary to fill the *RAMemf*[*n*] with the number of samples from the coil sensor specified by the user prior to testing. When asserting *strDWT*, the module starts to compute the convolution operation defined by Equation (2), which is essentially composed of a multiply–accumulate (MAC) process. The MAC operation requires the coefficients of a filter obtained from four Read-only Memory (ROM)previously filled with the corresponding coefficients (*ROM Lo-D* and *ROM Hi-D* for low-pass and high-pass decomposition filters, respectively; and *ROM Lo-R* and *ROM Hi-R* for low-pass and high-pass reconstruction filters, respectively), as well as the *emf*[*n-k*] signal, which is obtained by passing *emf*[*n*] through a *k*-level pipeline register. Finally, when the computation process is finished, signal *rdyDWT* is set to high. The approximation and detail coefficients will be given by the output signals *aCik* and *dCik*, respectively.

**Figure 7.** Block diagram of the DWT–IDWT unit: including Random Access Memory (RAM), Multiplier Accumulator (MAC), Read-only Memory (ROM), Alternating current AC and Direct current (DC) modules.

#### 2.6.2. Wavelet Entropy Digital Structure

Figure 8 shows a block diagram of the proposed digital architecture used to obtain the wavelet entropy of a signal by applying Equation (3) where the MAC process is the main operation. Note that it is essential to fill the RAM and *RAM-pi* with the corresponding relative wavelet energies *pi* of the corresponding "wavelet signals" prior to testing. To start the *SWT* computation, the signal StrWE is asserted. After that, the counter *Rddir* selects the signal *pi* to be processed, and the log2(*pi* 2) is computed by applying the algorithm proposed in [45], since it offers an easy implementation in hardware. Next, to obtain the required log*e*(*pi* <sup>2</sup>) value, a simple multiplication factor defined by Equation (9) is applied to log2(*pi* 2).

$$\log\_e(\mathbf{x}) = \frac{\log\_2(\mathbf{x})}{\log\_2(e)} \approx 0,693147 \ast \log\_2(\mathbf{x})\tag{9}$$

**Figure 8.** Block diagram of the wavelet entropy unit.

#### 2.6.3. FFNN Digital Structure

Figure 9 shows the block diagram of the proposed digital architecture to compute a regression FFNN. When the signal *strANN* is asserted, the input signals *SWT* and *T* are stored in in the first two memory elements of the *Ni* submodule. Note that the submodule *Ni* works as a storage memory for the total number of neurons that constitute the FFNN architecture used here (that is, neurons on the input layer, neurons on the hidden layers and neurons on the output layer). Memory ROMs *ROMidx-rd*, *ROMwi*, *ROM-layer*, and *ROMb[i]* store the indices of each neuron to read/write, depending on the actual layer (*layer*), synaptic weights, and biases. This design is based on a MAC operation in order to save element resources and use only one multiplier. The inputs for the MAC operation are the synaptic weights (*wi*) and the corresponding neuron outputs (*xi*). Finally, when the MAC process is finished, its output is summed by the corresponding bias (*bi*) in order to compute the activation function tansig(x) using the piece-wise linear function defined in [46].

**Figure 9.** Block diagram of the FFNN digital unit.

To train the FFNN, the information of the wavelet entropy and motor frame temperature are used. In addition, to establish a frame of reference between a motor with healthy winding insulation and a motor with a degraded insulation system, the index of dielectric absorption of the motor is extracted by means of a megger device during the degradation tests. These data (wavelet entropy, motor frame temperature, and dielectric absorption index) are obtained during the process of induced degradation to the winding insulation system on a three-phase induction motor whose characteristics are specified in the results section. The wavelet entropy and the temperature of the motor frame are used as the input data set for the training of the neural network and as the desired output, then the interpolation between the dielectric absorption rates obtained for each test are carried out and limits are established for the purposes of this work (that is, 95% for a healthy motor, and 10% for motor with serious winding insulation degradation).

#### **3. Results**

The diagnostic procedure and functionality of the smart-sensor proposed in this paper was validated in the laboratory on six IMs with the following same characteristics: 1.1 kW, 400 V, Y-connected, 50 Hz, 4 poles, where several experiments were carried out for healthy IMs and IMs with induced winding insulation degradation, as explained below.

#### *Experimental Set-Up*

An experimental test bench was designed to develop and implement the diagnostic technique proposed here and simulate a load using a three-phase squirrel-cage induction motor connected to a Direct Current (DC) generator, as shown in Figure 10b,c, respectively. The coil sensor and the thermocouple sensor (see Figure 10a), which was connected to the encased proprietary FPGA-based HSP unit, were attached to the frame of the motor. The laboratory room where the experiments were carried out was a closed space where the ambient temperature was maintained at an approximate value of 26 ◦C. Other elements that could potentially interfere with the experiments were removed (inverters, other test benches, etc.) to ensure that no other factor might influence the results.

**Figure 10.** Laboratory test bench: (**a**) coil sensor and thermocouple sensor, (**b**) induction motor, (**c**) DC generator acting as the load, (**d**) proprietary FPGA-based HSP unit.

Two different experiments were carried out in order to probe the functionality and effectiveness of the proposed smart-sensor. In the first experiment, five IMs with the same characteristics but different health statuses where diagnosed. The smart-sensor was placed on the frame of a healthy motor, then on an IM with one or two broken bars (one of the most common failures in this type of motor), but with a healthy winding insulation. Finally, smart-sensors were added to one IM with light winding insulation degradation, and one with severe winding insulation degradation. All experiments in this first stage were developed maintaining the same operating temperature in order to keep the winding insulation temperature in a controlled range of approximately 26 ◦C.

In the second experiment, a winding insulation degradation was progressively induced to one IM, in order to fully diagnose several levels of deterioration. With the purpose of establishing a reference between a healthy motor and a faulty motor in the insulation system, the first IM used was in a healthy condition at the beginning of the tests. Afterwards, an overheating of the winding insulation was artificially created by connecting and disconnecting one of the motor supply phases in successive cycles. In that way, while one phase was disconnected, the other two were overloaded, leading to abrupt thermal increments that produced higher temperatures than those defined by the thermal class of the insulation (class F). The connection–disconnection cycles of one supply phase were repeated a large number of times, a fact that led to an accelerated degradation of the insulation due to thermal effects. It is worth noting that this experimental setup tries, for the first time, to study the thermal degradation that the insulation system of an induction motor suffers when it is in service (that is, before a short circuit occurs between turns). A total of 100 tests were carried out on the same induction motor, thus generating a premature and irreversible wear to the insulation of the winding, since the temperatures reached in the machine frame exceeded 150 ◦C. This level implies that much higher temperatures were present inside the motor that clearly exceeded the limit for class insulation (155 ◦C at the hottest point).

#### **4. Discussion**

In this section, the results obtained from testing the smart-sensor by installing it on six similar induction motors with different induced failures are shown.

Firstly, to probe the effectiveness of the smart-sensor and diagnose the winding insulation degradation over different faults, the smart-sensor was installed on five IMs at an ambient temperature (26 ◦C), all with the same constructive characteristics, but with each one possessing a special failure case, namely: minimal insulation degradation, light insulation degradation, severe insulation degradation, and an IM with one and two broken bars, but with a healthy winding insulation.

It can be clearly seen in Figure 11 how the *SWT* parameter amplitudes were highly related to the winding insulation degradation, and it is also evident that the combined failures negligibly affected negligibly the results (that is, other failures like broken bars—one of the most frequent failures in IMs—did not affect the proposed methodology).

**Figure 11.** Wavelet entropy for different Induction Motors (IM) winding insulation degradations.

Figure 12 shows the results obtained when using the smart-sensor to compute the wavelet entropy for different prematurely induced degradation stages on the winding insulation of an IM. At the top of the same figure, the temperatures reached on the IM frame for different test points are shown, in order to contrast the results of when the temperature changes drastically. Similarly, five results displayed by the LCD of the smart-sensor proposed here are included, and shown at the top of Figure 12. These results correspond to tests labeled as A, B, C, and D for different winding insulation health states. For this purpose, over 100 tests were run. In each test, the winding insulation was degraded continuously. Evidently, the more severe the winding insulation degradation, the higher the amplitude of the wavelet entropy. Furthermore, note how the *SWT* parameter is also dependent on the temperature of the motor, since numerous tests showed an increase of *SWT* amplitudes with higher temperatures (especially in frame temperatures above 130 ◦C). Additionally at the top of Figure 12, the diagnosis offered by the smart-sensor proposed here is shown. The final results ranged from 10% to 95%, indicating the healthiness of the winding insulation (where 10% indicates a severe degradation, and 95% indicates very low or null degradation).

Considering a wavelet entropy value over 0.18 at ambient temperature (26 ◦C), a threshold value of 35% or below could be set as the criterion for discriminating between healthy and severe winding insulation system (requiring immediate maintenance) conditions (see Figure 12).

**Figure 12.** Wavelet entropy for different IM winding insulation degradation.

#### **5. Conclusions**

This work has introduced a new approach to performing an online estimation of the status of winding insulation degradation in IMs (a very common failure in this type of motors). The methodology is implemented in an FPGA in order to generate a smart-sensor, which is achieved by developing the digital cores needed to compute the DWT, the *SWT* index, and the regression FFNN. These tools provide the smart-sensor with the capability to automatically diagnose the health of the winding insulation, specifically before incipient faults progress into irreversible damage to the motor, making the smart-sensor an excellent device for the online diagnosis of winding insulation degradation.

What makes the smart-sensor proposed here even more attractive, is that the signal processing tools rely on a stray-flux analysis combined with the temperature of the analyzed IM, where the signals are obtained from a coil sensor and an E-type thermocouple sensor, both of which are installed externally on the motor frame. The coil sensor complies with several characteristics that make it an excellent alternative as a source of information, including its simple design, small size, low cost, installation flexibility, and non-invasive nature.

Furthermore, the proposed methodology relies on the study of the WE parameter obtained from the stray flux captured by a coil sensor. The wavelet entropy provides a very useful and practical index to gather information related to the analyzed signal, since it can characterize and combine the dynamism and order/disorder of this signal in a single value.

In addition, it can be deduced from the results obtained here that the *SWT* was sensitive to temperature variations in the analyzed IM, so it is very important to take this fact into account when diagnosing the severity degradation of the winding insulation, which was situation-controlled by the FFNN in this proposal.

For future work, the authors propose to further research in this open field, and suggest testing different IMs in order to deepen knowledge of the relationship discovered between the IM temperature and winding insulation degradation.

**Author Contributions:** J.A.A.-D. and R.A.O.-R. conceived and designed the experiments; I.Z.-R. performed the experiments and processed data; J.A.A.-D., R.A.O.-R., and R.d.J.R.-T. analyzed the data; I.Z.-R. and M.T.-H. wrote the paper.

**Funding:** We would like to thank Consejo Nacional de Ciencia y Tecnología (CONACYT) for providing economic support in this work (scholarship). Finally, thanks to the next projects: SEP-CONACYT 222453-2013, and FOFIUAQ-FIN201812. This wok was also funded by Spanish 'Ministerio de Ciencia Innovación y Universidades' and FEDER program in the framework of the 'Proyectos de I+D de Generación de Conocimiento

del Programa Estatal de Generación de Conocimiento y Fortalecimiento Científico y Tecnológico del Sistema de I+D+i, Subprograma Estatal de Generación de Conocimiento' (ref: PGC2018-095747-B-I00).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
