A General Method to Compute the Electric Flux Lines between Two Magnet Wires in Close Contact and Its Application for the Evaluation of Partial Discharge Risks in the Slots of Electric Machines Embedded in Future Transportation Systems

The sizing of the Electrical Insulation System (EIS) is an important 
challenge in electric motors of higher specific power 
driven by faster inverters. That keeps increasing the electric stress which the winding is submitted in the stator 
slot. Consequently, Partial Discharges (PD) are more likely to occur. Nowadays, the Paschen’s criterion is widely used to 
evaluate the risk of partial discharge. It requires the knowledge of electric field lines. This paper presents a method 
to precisely compute the electric field lines in a two-dimensional (2D) 
electrostatic problem. The field of study is composed of two magnet wires in 
close contact. Such configuration is representative of the turn-to-turn interaction in an electric motor slot. The problem 
is solved using the scalar potential formulation only. The notion of flux tubes 
is used for the post process of the electric field lines in a developed numerical code on Matlab. The 
developed method is compared to a ballistic method already included on Matlab. 
The work presented here is included in an automatic tool to suppress or reduce 
the partial discharge risk in a stator slot of high power density motor destined for future transportation systems.


Introduction
The role of transportation in the world sustainable development was firstly pointed out during the 1992 United Nation's Earth Summit at Rio de Janeiro [1]. Since then, the strong role of transportation in climate change, raw material depletion, human health and ecosystems equilibrium is as important as ever. The transportation sector is responsible for 24% of world CO 2 emissions [2] per year, just for being on the move [3]. Both passengers and freight road vehicles are responsible for 74% of the total transportation CO 2 emissions, while both aviation and shipping reach 22% and rail 1.3%. The transportation sector remains the largest consumer of oil: 57% of the global demand [4] where road represents 78.1%, air 11.4%, sea 7.3% and rail 3.2%. It is also responsible for 12% -70% of the total tropospheric air pollution mix [5]. Outdoor air pollution kills more than 8 million people across the world every year. Transportations are responsible for direct and indirect damages, or major changes, on ecosystems (air, marine and earth), which are often unpredictable [6].
In order to reduce the impact of transportation on global warming, human health and environmental issues, different efforts must be undertaken or pursued in all kinds of transport: increasing even more the efficiency of existing, Internal Combustion Engine (ICE) powertrains: lower fuel consumption, use of bio-or low carbon-fuels, better decontamination of exhaust gas; increasing the electrification of ICE powertrains: hybridization or full-electrification and increasing the capacity, the efficiency, the lifetime and the hybridization of embedded power sources. Such measures are essentially dependent on R & D efforts in the field of both Electrical and Electrochemical (applied to energy) Engineering.
Aircraft manufacturers such as Airbus, Boeing and Bombardier are engaged in the competition to develop more-and full-electric aircrafts. That incoming revolution takes place in a context where more and more people and countries are expecting a much greener air transportation. The proportion of on board electric power has continuously increased in aircrafts. The electric power tends to replace more and more systems which were powered by either pneumatic or hydraulic power. Figure 1 illustrates the increasing of inboard electrical equipment demand in commercial aviation [7]. The former on board electric bus first provided both constant voltage and frequency to the aircrafts devices. An Integrated Drive Generator (IDG) was used to change the variable speed of the jet engine to constant speed [8]. Between 1936 and 1946, the voltage supply has increased from 14.25 VDC to 28 VDC [9]. In recent aircrafts such as Airbus A380, Airbus A350 and Boeing 787, there is no more IDG. A gearbox is used to directly couple the engine generator to the jet engine. An alternative voltage of 115/200 VAC is produced with a frequency range from 350 to 800 Hz [8].
Since the introduction of power electronic power supplies, that provides easy control of the machine rotational speed, the electrical insulation of inverter fed motors faces new hazards. Fast changing supply voltage, with high dV/dt, may cause the apparition of Partial Discharges (PD), that results in accelerated insulation aging [10] and often leads to premature failure of the motors. In low voltage rotating machines, the stator insulation system is multi-levels. Its first component (primary insulation) is the polymer enamel on the magnet wire, among the others: inter-phase insulation, slot insulation and impregnation varnish. Depending on the desired thermal properties, there are several types of polymers being used nowadays in enameled wires: polyamide (PA), polyamide-imide (PAI), polyester-imide (PEI) and polyimide (PI). In random-wound stators powered by power inverters, in comparison with sinusoidal power supply, the magnet wire insulation is far more endangered. Hence the context of this paper is the primary insulation. Once the voltage exceeds the Partial Discharge Inception Voltage (PDIV), electronic avalanche will take place in the EIS. That leads to an ion bombardment of the insulator surfaces and an increase in temperature in the area submitted to PD. That chemically degrades the insulators.
The Paschen's criterion is widely used to evaluate the PD risk. It is essential to get the electric field lines precise computation. This paper presents a general method to get precise electric field lines for an electrostatic problem made of two magnet wires in close contact. The advantages of that method are that it, first, only use the electric scalar formulation, second, uses the same mesh defined for solving the problem with finite elements. First, the scalar potential formulation is presented. Then, the classic ballistic method already included on Matlab [11] is described and illustrated on a simple example. It follows the introduction of the flux tubes theory and its use in the proposed method. Finally, the developed method is compared to the ballistic method on a 2D electrostatic problem made of two magnet wires in close contact.

Scalar Potential Formulation
The basis of electromagnetism is the Maxwell's equations. The scalar potential formulation used to solve the problem with the finite elements method is derived from the following Maxwell's equations: (1) Advances in Aerospace Science and Technology div ρ = D (2) With E the electric field intensity, B the magnetic flux density, D the electric flux density and ρ the volume charge density of the dielectric medium.
The fields produced by a power frequency alternating voltage are not electrostatic. However, they are considered quasi-stationary: the fields time variation is neglected. Besides, in this work both air and polymer insulation materials are considered free of charge carrier. The simplified Maxwell's equations which come into play in the considered electrostatic problem are given below: Because E is curl free, it can be expressed as a gradient of a scalar potential.
That is the electric potential V: The electric flux density D is derived from the electric field intensity E and the medium properties: With ε 0 and ε r the permittivity of the air and the dielectric constant of the medium respectively. By combining (3), (4), (5) and (6) it comes the so-called Poisson's equation: Most of commercial finite elements software solve (7). That formulation is easy to handle and gives a unique solution. The main disadvantage of that formulation is that it does not give directly the electric field lines. Additional steps are required and are presented in the following paragraphs. By analogy with the magnetic vector potential, there also exists a vector potential formulation. Such formulation is complex to execute in an electrostatic problem with multiple conductors. It requires to put in place a network of branches and cuts all across the field of study [12]. However, [13] [14] have successfully put in application that formulation to solve an electrostatic problem with multiple conductors.

Electric Field Lines Computation Derived from a Scalar Potential Formulation
In this paragraph, two methods are presented. The first one is the ballistic method. It is used by Matlab stream [11] functions. The second method is based on the definition of electric flux tubes. This concept was used to developed our method.:

Ballistic Method
This method is one way of computing field lines from a scalar potential formulation. It is presented in [15]. If l is a vector tangential to a field line, then, in a direct orthonormal system, the cross-product of l and E is null. The point Advances in Aerospace Science and Technology on a field line can be calculated step by step with the following equations: With l x and l y respectively the projections of l along x and y axis in a 2D-system. From (8) it comes: E x and E y are the field components over the mesh. These are known. Starting at any point, the field line l can be computed by incrementing (9) given an arbitrary displacement ∆ x along x axis. Thus, the vertical displacement ∆ y along y axis is expressed as follow: Let us express Equation (10) between two consecutive points M i (x i , y i ) and M i+1 (x i + 1, y i+1 ): Let us define a constant step increase ∆ step at each iteration, so that: It is thus possible to express the M i+1 (x i + 1, y i+1 ) coordinates from M i (x i , y i ) data (coordinates and field components) and ∆ step: When chosen arbitrary, the starting point may not be the start of the field line. It is then necessary to integrate (12) backward to complete the line. This method is called a ballistic method. The mesh is swept and field lines are started in elements which do not already contain a field line. The field lines are computed by integration. The integration process is stopped if one of the following condition is checked:  The field line enters a forbidden region (for instance the limit of the domain);  The field line reaches a null field;  The field line loops back onto itself;  The field line has too many segments. This is a safety measure in case the previous conditions do not work properly. By doing so, there is the same density of field line all over the model. It does not allow to represent the electric field intensity.
As an illustration, let us consider an electrostatic problem made of two infinitely long cylinder oppositely. Due to the symmetries, only a quarter of the field of study is considered as displayed on Figure 2. The mesh is represented with black crosses. It is made of rectangular elements (dotted red line). The blue line delimits the conductor contour. Let us apply the Whittaker ballistic method on that system. The algorithm noted above is implanted in a Matlab script [11]. The orientation vector of the electric field is computed on all mesh nodes, except for the ones located inside the conductor. The starting points are chosen on the contour of the domain. Figure 2 displays the computed orientation vectors and the starting points (circles). The starting points are arbitrary chosen on the limit of the field of study. Let us designate the starting point Q(x 0 , y 0 ) of a field line. Then, the next point on the line M 1 (x 1 , y 1 ) is computed by iteration from point Q. The iteration is backward because the starting points are on the end of the field lines: With E x,0 and E y,0 the electric field orientation vector components along x and y axis respectively on starting point Q. The step increase ∆ step at each iteration is taken equal to the grid spacing along x axis d x . Figure  At each iteration, the field lines extend. The ith iteration between two consecutive points M i (x i , y i ) and M i+1 (x i + 1, y i+1 ) correspond to: The computed field lines after 11 iterations are displays on Figure 3.

Flux Tubes Based Method
The objective is to represent the strength of the electric field by the density of the electric flux lines [16]. The core of the method is the choice of the field lines starting points. A predetermined amount of flux δϕ has to be present between two adjacent lines a and b: From (17) to (21) E refers to the normal electric field component.
The field lines are then built from their starting point by integrating as in Whittaker's method [16]. As the same flux quantity δϕ is present between two adjacent lines, the electric field strength is represented by the concentration of field lines. The same termination criteria as in Whittaker's method are used. However, an electric machine slot filled with conductors is a more complex case.
It mainly happens that all the field lines do not pass through a single V contour. Let us consider two conductors at potential of V 1 and V 2 respectively. A field line starting from V 1 contour is identified by  Step 3: repeat step 2 until the whole V 2 contour is swept;  Step 4: integrate backward from the first intersection point to find the remaining starting points. Figure 5 [16] illustrates the field lines computed between two conductors at potential V 1 and V 2 using the presented algorithm.

Developed Method
The method we developed is an upgrade of the method presented by Horowitz [16]. A flux function is built from the electric field components on all the nodes. The field lines correspond to iso-values of this scalar function. The proposed method does not require a uniform mesh. The calculation is done on the same mesh on which the scalar potential V is computed.

Electric Field and Electric Flux
The 2D-finite element model on Ansys Mechanical APDL [17] uses eight nodes elements. On each element, the coordinates and the scalar potential solution are expressed as a combination of each nodes data with a defined shape function.
Shape functions are expressed in the local coordinate system (u, v) of the element. It is a coordinate system attached to the element which defines the location of each node.  With x e and y e the coordinates in the global coordinate system (x,y) of the nodes of an element e. For eight nodes elements, the shape functions N are given in [18]. The electric field on the nodes of an element is derived from the voltage values on the node using (5). It requires the expression of partial derivatives in the global system. First, the partial derivatives of a function in the local system (u,v) can be expressed from its partial derivatives in the global system (x,y): With: With: The electric flux components are computed from electric field components using (6).

Equipotential Lines
The next step consists of forming equipotential lines. These are contours on which the voltage value is constant. Equipotential lines are obtained by regrouping nodes with the same voltage value. Let us see the steps to compute an equipotential line at voltage V. A test is done on each edge of each element to check whether or not it is crossed by the equipotential line V. As can be seen on Figure 6, each edge is composed of three nodes. There are three possible outcomes:  The voltage V is lower or bigger than any of the three nodes voltage. The edge is not intersected by the equipotential line V;  The voltage V is equal to one of the three nodes voltage. The equipotential line intersects the edge at the corresponding node of voltage V;  The voltage V is between any of the three nodes voltage. The edge is intersected by the equipotential line at a point which is interpolated. Figure 6 illustrates an equipotential intersecting two edges of an element. On one edge the intersection point is an already existing node (node 6). The second edge intersection point has to be interpolated (red cross). The local numbering of the nodes is rearranged compared to Ansys eight nodes reference element to facilitate the programming. A second order polynomial is used for the interpolation. The following equation gives the parametric expression of the polynomial: In the considered example, (V 0 , V 0.5 , V 1 ) are respectively the finite elements voltage solution on nodes (1, 2, 3) in Figure 6. The three polynomial coefficients (a 1,V , a 2,V , a 3,V ) are thus determined.
With V the voltage of the considered equipotential line V. The obtained parameter t 0 is injected back into (29) to determine the interpolated node coordinates, electric field components and electric flux components. As in (30), the associated polynomial coefficients are deduced from the known nodes data.

Flux Function
At this point, equipotential lines made of nodes at the same scalar potential V are determined. The coordinates and fields components are determined on all these nodes. The number of equipotential lines depends on the accuracy from which the electric scalar potential problem is solved. The finer the mesh used to solve the problem, the higher the number of equipotential lines that can be accurately determined. For each equipotential line, a geometrical reference is defined on its barycenter. The points on the line are located by using polar coordinates in this reference. A starting point Q is chosen for instance by means of the angular coordinates. Figure 7 The trace of a flux tube on the equipotential line is delimited by two points P i and P i+1 which are given by the predetermined flux per meter δϕ:

Field of Study
The problem consists of two enameled magnet wires in close contact in air. Such configuration is representative of a contact between two adjacent wires in a sta- The dielectric constant of air is 1. The dielectric constant of the enamel is greater than 1 and depends on the polymer used. In this paper, the dielectric constant of the enamel has been taken as 3.5. The wires parameters are recapped on Table 1. Advances in Aerospace Science and Technology  The model is realized on Ansys Mechanical APDL [17]. The mesh is made with plane 121 element. This is an eight nodes element. The mesh is composed of 4144 elements. Figure 9 illustrates the mesh. The finite elements solution is obtained in 12.6 seconds.

Developed Method
A voltage amplitude of 1000 V peak is applied as boundary conditions on the copper contour of the left wire ( Figure 8). A null potential is imposed on the copper contour of the right wire. The voltage drop in the domain is subdivided into one hundred equipotential lines. The starting point of each equipotential line is taken on the symmetry axis. That because the symmetry axis is an electric field line. The flux tubes contours are determined following the process presented on paragraph 4. Figure 10 shows that the electric flux crossing the equipotential lines is constant. It has been arbitrary chosen to plot twenty field lines. Each flux tubes limit is thus one nineteenth of the electric flux. That because the first field line is the symmetry axis. Figure 11 displays the one hundred equipotential lines (colored dotted lines) and the twenty electric field lines (colored solid lines). The black contours are the enamel external layer. The field lines are computed from the finite elements solution in 9.11 seconds.

Matlab Ballistic Method
A ballistic method is also used to compute the electric field lines. The voltage drop is also taken equals to 1000 V peak between the copper cores. This method is Advances in Aerospace Science and Technology   implanted on Matlab function stream2 [11]. That function requires a uniform mesh. However, the finite elements mesh on Figure 9 is not uniform.  X a and Y a results in an N uniformly spaced values.
In this paper N is taken equals to 600. The mesh generated by meshgrid is referred as Matalb mesh. The Matlab mesh is different from the one generated with the finite elements software. Multiple points of the finite elements mesh are deleted. Besides, the Matlab mesh adds some points. The voltage solution on these points is linearly interpolated from the existing ones using griddata function [11]. Then, the voltage gradient over the Matlab mesh is computed using gradient function [11]. We now have all the inputs necessary for using stream2 functions. Figure 12 displays twenty electric field lines in air (black solid lines). Only the contours of the enamel external layer are represented. The field lines are computed in 1.4 seconds with the Matlab ballistic method.

Comparison for PDIV Evaluation
In PD evaluation, only the part of the electric field line in the air gap is considered. Besides, the electric field along the field lines has to be uniform. Figure 13 displays the electric field amplitude. The field lines from Figure 12 are superposed. It can be seen that the electric field in the air gap is quite constant along the obtained field lines. So the Paschen's criterion can be applied on the part of the field lines in the air gap. To take into account the impact of the enamel layer on PDIV, the Paschen's criterion is corrected as in [19]. The secondary electron emission coefficient γ is equal to 9 × 10 −4 . The voltage on the enamel overcoat contours at the interface with the air is picked up. Figure 14 displays the computed voltage drops along each part of field lines in air for the two approaches at   evaluated PDIV (1260 V peak ). The points intersecting with the Paschen's curve indicate field lines on which PD activity is likely to occur. Both methods give close results for field lines in air longer than 10 µm. PD is evaluated to take place on field line of 37 µm length with a voltage drop in air of 561 V.

Conclusion
Partial Discharges (PD) phenomenon represents a great deal in the design of future rotating machines fed by inverter. On one hand, new faster components made out of SiC and GaN technologies will considerably improve the performances of the inverters. On the other hand, they will generate harder voltage fronts at the motor terminals which lead to higher transient voltage overshoots. PD will be more likely to appear and the insulation lifespan will be reduced due to both voltage overshoots and high switching frequency. That is the reason why it is absolutely necessary to take into account such a phenomenon when designing the motor to avoid any PD appearance, rather than searching for solutions in a motor already produced to remove them. The Paschen's criterion is used to evaluate PD activity. This criterion is corrected for taking into account the impact of the enamel layer. The Paschen's criterion requires the complete knowledge of electric field lines geometry. A numerical code has then been developed on Matlab to drawn electric field lines from a 2D-Finite Elements (2D-FEM) solution of an electrostatic problem. The originality of the proposed method is that only the scalar potential solution is required and the density of displayed lines is analogous to the field intensity. It has been compared to a ballistic approach already implanted on Matlab. It results that the Matlab function is much faster. Both methods give similar results for the considered 2D-electrostatic problem. The development of the proposed method provided better understanding on scalar potential formulation and the notions of electric flux tubes. These methods have been used to evaluate the PDIV of several magnet wires. The ultimate goal is the computation of insulation design graphs for suppressing PD risk between turns in a stator slot.